序列比对（22）中间字符串分支定界方法中更紧的界

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#define MAXSEQ 1000
#define max(a, b) ((a) > (b) ? (a) : (b))

struct Sequence {
    char* s;
    int* a;      // 该向量存储每个字符在字符集中的序号
    int l;
};
typedef struct Sequence* Seq;

const int nalpha = 4;     // 字符种数
const char sortedAlpha[] = {'A', 'C', 'G', 'T'};    // 排过序的字符集

char* sec2time(time_t t);
void strUpper(char *s);
int search(char c);                      /* 在排序后的字符集中寻找字符c，如果找到返回序号；找不到返回-1 */
Seq* readSeq(char* filename, const int t);          /* 从文件中读取多条序列 */
Seq create(char* s, const int n, const int i);  
void destroy(Seq seq);
int nextVertex(int* a, int d, const int l);     /* 得到下一个顶点 */
int byPass(int* a, int d, const int l);           /* 跳过某个分支，得到下一个顶点 */
int getDistance(Seq* mulSeq, const int t, int* a, const int l, const int d);        /* 计算一个l-元组和DNA的最小总距离 */
int* getStartPoints(Seq* mulSeq, const int t, int* a, const int l);    /* 获取起始位点向量 */
int getLowerBound(Seq* mulSeq, const int t, int* a, int* b, const int l);      /* 计算TotalDistance(v, DNA)的下界 */
void findMedianStr(Seq* mulSeq, const int t, const int l);           /* 寻找中间字符串 */

int main(void) {
    time_t start, end;
    char* tstr;
    int i, t, l;
    char filename[30];
    Seq* mulSeq;
    start = time(NULL);  // 开始时间
    printf("t: ");
    scanf(" %d", &t);
    printf("l: ");
    scanf(" %d", &l);
    printf("seq_file path: ");
    scanf(" %s", filename);
    mulSeq = readSeq(filename, t);
    findMedianStr(mulSeq, t, l);
    for (i = 0; i < t; i++)
        destroy(mulSeq[i]);
    free(mulSeq);
    end = time(NULL);   // 结束时间
    tstr = sec2time(end - start);
    printf("time spent: %s\n", tstr);
    free(tstr);
    return 0;
}

char* sec2time(time_t t) {
/* 将秒数转化为时间字符串 */
    int h, m, s;
    char* tstr;
    if ((tstr = malloc(sizeof(char) * 50)) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(8);
    }
    h = t / 3600;
    m = (t % 3600) / 60;
    s = t % 60;
    if (sprintf(tstr, "%02d:%02d:%02d", h, m, s) > 0)
        return tstr;
    else {
        sprintf(tstr, "more than %d hours", h);
        return tstr;
    }
}

void strUpper(char *s) {
    while (*s != '\0') {
        if (*s >= 'a' && *s <= 'z')
            *s -= 32;
        s++;
    }
}

int search(char c) {
/* 在排序后的字符集中寻找字符c，如果找到返回序号；找不到返回-1 */
    int i;
    int th = 5;
    int h, l, m;
    if (nalpha < th) {
        for (i = 0; i < nalpha; i++)
            if (sortedAlpha[i] == c)
                return i;
        return -1;
    } else {
        l = 0;
        h = nalpha - 1;
        while (l <= h) {
            m = (l + h) / 2;
            if (sortedAlpha[m] == c)
                return m;
            else if (sortedAlpha[m] < c)
                l = m + 1;
            else 
                h = m - 1;
        }
        return -1;
    }
}

Seq* readSeq(char* filename, const int t) {
/* 从文件中读取多条序列 */
    int i, n;
    char buf[MAXSEQ];
    FILE* fp;
    Seq* mulSeq;
    if ((fp = fopen(filename, "r")) == NULL) {
        fprintf(stderr, "Error: cannot open '%s'\n", filename);
        exit(1);
    }
    if ((mulSeq = (Seq*) malloc(sizeof(Seq) * t)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);        
    }
    for (i = 0; i < t; i++) {
        if (fgets(buf, MAXSEQ, fp) == NULL) {
            fprintf(stderr, "Error: reading seq from '%s'\n", filename);
            exit(2);
        }
        n = strlen(buf);
        if (buf[n - 1] != '\n' && ! feof(fp)) {
            fprintf(stderr, "Error: line %d is too long.\n", i + 1);
            exit(3);
        }
        if (! feof(fp)) {
            buf[n - 1] = '\0';
            mulSeq[i] = create(buf, n - 1, i);
        } else {
            mulSeq[i] = create(buf, n, i);
        }
    }
    fclose(fp);
    return mulSeq;
}

Seq create(char* s, const int n, const int i) {
/* n: s的长度，不包括末尾的\0 */
    Seq seq;
    int j;
    if ((seq = (Seq) malloc(sizeof(struct Sequence))) == NULL || \
        (seq->s = (char*) malloc(sizeof(char) * (n + 1))) == NULL || \
        (seq->a = (int*) malloc(sizeof(int) * n)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    strcpy(seq->s, s);
    // 将序列中所有字符变为大写
    strUpper(seq->s);
    // 找到每个字符在字符集中的序号。如果是在后期做这个工作，那么就有很多重复，降低效率。
    for (j = 0; j < n; j++)
        if ((seq->a[j] = search(seq->s[j])) < 0) {
            fprintf(stderr, "Error: char %d of seq %d is %c, which is out of range.\n", j + 1, i + 1, seq->s[j]);
            exit(16);
        }
    seq->l = n;
    return seq;
}

void destroy(Seq seq) {
    free(seq->a);
    free(seq->s);
    free(seq);
}

int nextVertex(int* a, int d, const int l) {
/* 得到下一个顶点 */
    int i;
    if (d < l - 1) {
        a[d + 1] = 0;
        return d + 1;
    }
    for (i = l - 1; i >= 0; i--) {
        if (a[i] < nalpha - 1) {
            a[i]++;
            return i;
        }
    }
    return -1;
}

int byPass(int* a, int d, const int l) {
/* 跳过某个分支，得到下一个顶点 */
    int i;
    for (i = d; i >= 0; i--) {
        if (a[i] < nalpha - 1) {
            a[i]++;
            return i;
        }
    }
    return -1;
}        

int getDistance(Seq* mulSeq, const int t, int* a, const int l, const int d) {
/* 计算一个l-元组和DNA的最小总距离 */
    int i, j, k;
    int dist, minDist, totalDist;
    for (i = 0, totalDist = 0; i < t; i++) {
        for (j = 0, minDist = l; j <= mulSeq[i]->l - l; j++) {
            for (k = 0, dist = 0; k <= d; k++)
                dist += a[k] == mulSeq[i]->a[j + k] ? 0 : 1;
            if (dist < minDist)
                minDist = dist;
        }
        totalDist += minDist;
    }
    return totalDist;
}       

int* getStartPoints(Seq* mulSeq, const int t, int* a, const int l) {
/* 获取起始位点向量 */
    int* st;
    int i, j, k;
    int dist, minDist;
    if ((st = (int*) malloc(sizeof(int) * t)) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    for (i = 0; i < t; i++) {
        for (j = 0, st[i] = 0, minDist = l; j <= mulSeq[i]->l - l; j++) {
            for (k = 0, dist = 0; k < l; k++)
                dist += a[k] == mulSeq[i]->a[j + k] ? 0 : 1;
            if (dist < minDist) {
                st[i] = j;
                minDist = dist;
            }
        }
    }
    return st;
}

int getLowerBound(Seq* mulSeq, const int t, int* a, int* b, const int l) {
/* 计算TotalDistance(v, DNA)的下界 */
    int i, d;
    int dist, minDist;
    // 遍历所有可能的l-元组，寻找最小距离
    for (i = 0; i < l; i++)
        a[i] = 0;
    minDist = l * t;
    d = 0;
    while (d >= 0) {
        if (d < l - 1) {
            dist = getDistance(mulSeq, t, a, l, d);
            if (dist + b[l - 2 - d] > minDist)
                d = byPass(a, d, l);
            else 
                d = nextVertex(a, d, l);
        } else {
            dist = getDistance(mulSeq, t, a, l, l - 1);
            if (dist < minDist) {  // 注意最小值对应的l-元组可能有多种，这里只保留了其中一种。其实还可以设置如果minDist = 0时终止搜索。
                minDist = dist;
                if (minDist == 0)
                    return 0;
            }
            d = nextVertex(a, d, l);
        }
    }
    return minDist;
}

void findMedianStr(Seq* mulSeq, const int t, const int l) {
/* 寻找中间字符串 */
    int* a;     // l-元组
    int* ma;    // 最小距离对应的l-元组
    int i, d, j;
    int dist, minDist;
    char* cs;    // consensus.
    int* st;     // start point vector
    int* b;     // 各长度的l-元组的最小距离
    int md, lb, ib;
    int ml = 6;     // Minimum length of l-tuple tried.
    if ((a = (int*) malloc(sizeof(int) * l)) == NULL || \
        (ma = (int*) malloc(sizeof(int) * l)) == NULL || \
        (cs = (char*) malloc(sizeof(char) * (l + 1))) == NULL || \
        (b = (int*) malloc(sizeof(int) * (l - 1))) == NULL) {
        fputs("Error: no space!\n", stderr);
        exit(8);
    }
    // 下面一段是为了找出每个长度的L-元组的总距离的下界
    ib = -1;     // 总距离最大的L-元组的长度 - 1。如果小于0，说明所有可能的L-元组其总距离都小于0。
    if (l > 2) {
        md = l - 2 >= ml ? max(ml, (int) sqrt((double) l)) : (int) sqrt((double) l);
        for (i = 0; i < l - 1; i++) { // 获取各个长度l-元组的最小距离
            b[i] = getLowerBound(mulSeq, t, a, b, i + 1);
            if (i >= md && b[i] > 0)
                break;
        }
        if (i < l - 1) 
            ib = i;
    } else if (l == 2) {
        b[0] = getLowerBound(mulSeq, t, a, b, 1);
        if (b[0] > 0)
            ib = 0;
    }
    // 遍历所有可能的l-元组，寻找最小距离
    for (i = 0; i < l; i++)
        ma[i] = a[i] = 0;
    minDist = l * t;
    d = 0;
    while (d >= 0) {
        if (d < l - 1) {   // 如果l == 1的话这一段是不会执行的。
            dist = getDistance(mulSeq, t, a, l, d);
            if (ib >= 0) {     // 如果将word拆分为u+v，此处计算TotalDistance(v, DNA)的下界
                j = (l - 1 - d) % (ib + 1);
                lb = l - 2 - d <= ib ? b[l - 2 - d] : b[ib] * (l - 1 - d) / (ib + 1) + (j == 0 ? 0 : b[j - 1]);
            } else {
                lb = 0;
            }
            if (dist + lb > minDist)
                d = byPass(a, d, l);
            else 
                d = nextVertex(a, d, l);
        } else {
            dist = getDistance(mulSeq, t, a, l, l - 1);
            if (dist < minDist) {  // 注意最小值对应的l-元组可能有多种，这里只保留了其中一种。其实还可以设置如果minDist = 0时终止搜索。
                minDist = dist;
                for (i = 0; i < l; i++)
                    ma[i] = a[i];
            }
            d = nextVertex(a, d, l);
        }
    }
    // 获取共有序列
    for (i = 0; i < l; i++)
        cs[i] = sortedAlpha[ma[i]];
    cs[l] = '\0';
    // 获取最优起始位点向量
    st = getStartPoints(mulSeq, t, ma, l);
    // 输出结果
    printf("Min distance: %d\n", minDist);
    printf("Starting point vector: ");
    for (i = 0; i < t; i++)
        printf("%d ", st[i] + 1);
    printf("\n");
    printf("Target seq:\n");
    for (i = 0; i < t; i++) {
        printf("\t");
        for (j = st[i]; j < st[i] + l; j++)
            printf("%c", mulSeq[i]->s[j]);
        printf("\n");
    }
    printf("Consensus seq: %s\n", cs);
    // 释放内存
    free(a);
    free(ma);
    free(st);
    free(cs);
    free(b);
}