Mercurial > repos > clustalomega > clustalomega
diff clustalomega/clustal-omega-1.0.2/src/clustal/ktuple_pair.c @ 1:bc707542e5de
Uploaded
author | clustalomega |
---|---|
date | Thu, 21 Jul 2011 13:35:08 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/clustalomega/clustal-omega-1.0.2/src/clustal/ktuple_pair.c Thu Jul 21 13:35:08 2011 -0400 @@ -0,0 +1,842 @@ +/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */ + +/********************************************************************* + * Clustal Omega - Multiple sequence alignment + * + * Copyright (C) 2010 University College Dublin + * + * Clustal-Omega is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * This file is part of Clustal-Omega. + * + ********************************************************************/ + +/* + * RCS $Id: ktuple_pair.c 230 2011-04-09 15:37:50Z andreas $ + * + * + * K-Tuple code for pairwise alignment (Wilbur and Lipman, 1983; PMID + * 6572363). Most code taken from showpair.c (Clustal 1.83) + * DD: some functions now have lots of parameters as static variables + * were removed to make code OpenMP-friendly + * + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + + +#include <stdio.h> +#include <string.h> +#include <ctype.h> +#include <stdlib.h> +#include <math.h> +#include <assert.h> + +#ifdef HAVE_OPENMP +#include <omp.h> +#endif + +#include "squid/squid.h" +#include "util.h" +#include "symmatrix.h" +#include "ktuple_pair.h" +#include "log.h" +#include "progress.h" + +#define END_MARK -3 /* see interface.c in 1.83 */ +#define NUMRES 32 /* max size of comparison matrix */ + +/* see notes below */ +#undef SORT_LAST_ELEMENT_AS_WELL + +/* gap_pos1 = NUMRES-2; /@ code for gaps inserted by clustalw @/ */ +static const int GAP_POS2 = NUMRES-1; /* code for gaps already in alignment */ +static bool DNAFLAG = FALSE; + +static const char *AMINO_ACID_CODES = "ABCDEFGHIKLMNPQRSTUVWXYZ-"; +static const char *NUCLEIC_ACID_CODES = "ACGTUN-"; +/* As far as I understand the gap symbol should not be necessary here, + * because we use isgap for testing later anyway. But changing this, + * will affect max_res_code and max_nuc as well. So I leave it for now + * as it is. AW + */ + +static bool percent = TRUE; + +static void make_ptrs(int *tptr, int *pl, const int naseq, const int l, const int ktup, const int max_res_code, char **seq_array); +static void put_frag(const int fs, const int v1, const int v2, const int flen, const int curr_frag, int *next, int *maxsf, int **accum); +static bool frag_rel_pos(int a1, int b1, int a2, int b2, int ktup); +static void des_quick_sort(int *array1, int *array2, const int array_size); +static void pair_align(int seq_no, int l1, int l2, int max_res_code, ktuple_param_t *aln_param, + char **seq_array, int *maxsf, int **accum, int max_aln_length, + int *zza, int *zzb, int *zzc, int *zzd); +static void encode(char *seq, char *naseq, int l, const char *res_codes); +static int res_index(const char *lookup, char c); + + +typedef struct { + int i1; + int i2; +} two_ints_t; + + + +/* default ktuple pairwise alignment parameters + * + */ +/* protein + */ +/* designated initializer */ +const ktuple_param_t default_protein_param = { + .ktup = 1, + .wind_gap = 3, + .signif = 5, + .window = 5, +}; +/* dna + */ +/* designated initializer */ +const ktuple_param_t default_dna_param = { + .ktup = 2, + .wind_gap = 5, + .signif = 4, + .window = 4, +}; + + +/** + * note: naseq should be unit-offset + */ +static void +encode(char *seq, char *naseq, int l, const char *res_codes) +{ + /* code seq as ints .. use GAP_POS2 for gap */ + register int i; + bool seq_contains_unknown_char = FALSE; + /*LOG_DEBUG("seq=%s naseq=%p l=%d", &(seq[1]), naseq, l); */ + + + for (i=1; i<=l; i++) { + char res = toupper(seq[i]); + if (isgap(res)) { + naseq[i] = GAP_POS2; /* gap in input */ + } else { + naseq[i] = res_index(res_codes, res); + } + + /*LOG_DEBUG("Character '%c' at pos %d", res, i);*/ + if (-1 == naseq[i]) { + seq_contains_unknown_char = TRUE; + /*LOG_DEBUG("Unknown character '%c' at pos %d", res, i);*/ + } + /*LOG_DEBUG("na_seq[%d]=%d", i, naseq[i]);*/ + } + + if (TRUE == seq_contains_unknown_char) + Log(&rLog, LOG_WARN, "Unknown character in seq '%s'", &(seq[1])); + + naseq[i] = END_MARK; + + return; +} +/* end of encode */ + + +/** + * + */ +static int +res_index(const char *t, char c) +{ + register int i; + for (i=0; t[i] && t[i] != c; i++) + ; + if (t[i]) { + return (i); + } else { + return -1; + } +} +/* end of res_index */ + + +/** + * + */ +static void +make_ptrs(int *tptr, int *pl, const int naseq, const int l, const int ktup, const int max_res_code, char **seq_array) +{ + /* FIXME make 10 a constant and give it a nice name */ + static int a[10]; + int i, j, code, flag; + char residue; + const int limit = (int) pow((double)(max_res_code+1),(double)ktup); + + for (i=1;i<=ktup;i++) + a[i] = (int) pow((double)(max_res_code+1),(double)(i-1)); + + for (i=1; i<=limit; ++i) + pl[i]=0; + for (i=1; i<=l; ++i) + tptr[i]=0; + + for (i=1; i<=(l-ktup+1); ++i) { + code=0; + flag=FALSE; + for (j=1; j<=ktup; ++j) { + /* Log(&rLog, LOG_FORCED_DEBUG, "naseq=%d i=%d j=%d seq_array[naseq]=%p", + * naseq, i, j, seq_array[naseq]); + */ + residue = seq_array[naseq][i+j-1]; + /* Log(&rLog, LOG_FORCED_DEBUG, "residue = %d", residue); */ + if ((residue<0) || (residue > max_res_code)){ + flag=TRUE; + break; + } + code += ((residue) * a[j]); + } + if (flag) + continue; + ++code; + if (0 != pl[code]) + tptr[i] =pl[code]; + pl[code] = i; + } + + return; +} +/* end of make_ptrs */ + + +/** + * + * FIXME Why hardcoding of 5? + */ +static void +put_frag(const int fs, const int v1, const int v2, const int flen, const int curr_frag, int *next, int *maxsf, int **accum) +{ + int end; + accum[0][curr_frag]=fs; + accum[1][curr_frag]=v1; + accum[2][curr_frag]=v2; + accum[3][curr_frag]=flen; + + if (!*maxsf) { + *maxsf=1; + accum[4][curr_frag]=0; + return; + } + + if (fs >= accum[0][*maxsf]) { + accum[4][curr_frag]=*maxsf; + *maxsf=curr_frag; + return; + } else { + *next=*maxsf; + while (TRUE) { + end=*next; + *next=accum[4][*next]; + if (fs>=accum[0][*next]) + break; + } + accum[4][curr_frag]=*next; + accum[4][end]=curr_frag; + } + + return; +} +/* end of put_frag */ + + +/** + * + */ +static bool +frag_rel_pos(int a1, int b1, int a2, int b2, int ktup) +{ + if (a1-b1 == a2-b2) { + if (a2<a1) { + return TRUE; + } + } else { + if (a2+ktup-1<a1 && b2+ktup-1<b1) { + return TRUE; + } + } + return FALSE; +} +/* end of frag_rel_pos */ + + + + +/** + * + * @note: This is together with des_quick_sort most time consuming + * routine according to gprof on r110. Tried to replace it with qsort + * and/or QSortAndTrackIndex(), which is always slower! So we keep the + * original. + * + * Original doc: Quicksort routine, adapted from chapter 4, page 115 + * of software tools by Kernighan and Plauger, (1986). Sort the + * elements of array1 and sort the elements of array2 accordingly + * + * There might be a bug here. The original function apparently never + * touches the last element and keeps it as is. Tried to fix this (see + * SORT_LAST_ELEMENT_AS_WELL) which gives slightly worse performance + * (-0.5% on BB). My fix might not be working or it's not a bug at + * all... + * + * + * + */ +static void +des_quick_sort(int *array1, int *array2, const int array_size) +{ + int temp1, temp2; + int p, pivlin; + int i, j; + int lst[50], ust[50]; /* the maximum no. of elements must be*/ + /* < log(base2) of 50 */ + +#if 0 + for (i=1; i<=array_size; i++) { + Log(&rLog, LOG_FORCED_DEBUG, "b4 sort array1[%d]=%d array2[%d]=%d", i, array1[i], i, array2[i]); + } +#endif + lst[1] = 1; + +#ifdef SORT_LAST_ELEMENT_AS_WELL + ust[1] = array_size; +#else + /* original */ + ust[1] = array_size-1; +#endif + p = 1; + + + while (p > 0) { + if (lst[p] >= ust[p]) { + p--; + } else { + i = lst[p] - 1; + j = ust[p]; + pivlin = array1[j]; + while (i < j) { + for (i=i+1; array1[i] < pivlin; i++) + ; + for (j=j-1; j > i; j--) + if (array1[j] <= pivlin) break; + if (i < j) { + temp1 = array1[i]; + array1[i] = array1[j]; + array1[j] = temp1; + + temp2 = array2[i]; + array2[i] = array2[j]; + array2[j] = temp2; + } + } + + j = ust[p]; + + temp1 = array1[i]; + array1[i] = array1[j]; + array1[j] = temp1; + + temp2 = array2[i]; + array2[i] = array2[j]; + array2[j] = temp2; + + if (i-lst[p] < ust[p] - i) { + lst[p+1] = lst[p]; + ust[p+1] = i - 1; + lst[p] = i + 1; + + } else { + lst[p+1] = i + 1; + ust[p+1] = ust[p]; + ust[p] = i - 1; + } + p = p + 1; + } + } + +#if 0 + for (i=1; i<=array_size; i++) { + Log(&rLog, LOG_FORCED_DEBUG, "after sort array1[%d]=%d array2[%d]=%d", i, array1[i], i, array2[i]); + } +#endif + + return; +} +/* end of des_quick_sort */ + + + +/** + * + * FIXME together with des_quick_sort most time consuming routine + * according to gprof on r110 + * + */ +static void +pair_align(int seq_no, int l1, int l2, int max_res_code, ktuple_param_t *aln_param, + char **seq_array, int *maxsf, int **accum, int max_aln_length, + int *zza, int *zzb, int *zzc, int *zzd) +{ + int next; /* forrmerly static */ + int pot[8],i, j, l, m, flag, limit, pos, vn1, vn2, flen, osptr, fs; + int tv1, tv2, encrypt, subt1, subt2, rmndr; + char residue; + int *diag_index; + int *displ; + char *slopes; + int curr_frag; + const int tl1 = (l1+l2)-1; + + assert(NULL!=aln_param); + + /* + Log(&rLog, LOG_FORCED_DEBUG, "DNAFLAG=%d seq_no=%d l1=%d l2=%d window=%d ktup=%d signif=%d wind_gap=%d", + DNAFLAG, seq_no, l1, l2, window, ktup, signif, + wind_gap); + */ + + slopes = (char *) CKCALLOC(tl1+1, sizeof(char)); + displ = (int *) CKCALLOC(tl1+1, sizeof(int)); + diag_index = (int *) CKMALLOC((tl1+1) * sizeof(int)); + + for (i=1; i<=tl1; ++i) { + /* unnecessary, because we calloced: slopes[i] = displ[i] = 0; */ + diag_index[i] = i; + } + + for (i=1;i<=aln_param->ktup;i++) + pot[i] = (int) pow((double)(max_res_code+1),(double)(i-1)); + limit = (int) pow((double)(max_res_code+1),(double)aln_param->ktup); + + + + /* increment diagonal score for each k_tuple match */ + + for (i=1; i<=limit; ++i) { + vn1=zzc[i]; + while (TRUE) { + if (!vn1) break; + vn2 = zzd[i]; + while (0 != vn2) { + osptr = vn1-vn2+l2; + ++displ[osptr]; + vn2 = zzb[vn2]; + } + vn1=zza[vn1]; + } + } + + + /* choose the top SIGNIF diagonals + */ + +#ifdef QSORT_REPLACEMENT + /* This was an attempt to replace des_quick_sort with qsort(), + * which turns out to be much slower, so don't use this + */ + + /* FIXME: if we use this branch, we don't need to init diag_index + * before, because that is done in QSortAndTrackIndex() + * automatically. + */ +#if 0 + for (i=1; i<=tl1; i++) { + Log(&rLog, LOG_FORCED_DEBUG, "b4 sort disp[%d]=%d diag_index[%d]=%d", i, diag_index[i], i, displ[i]); + } +#endif + + QSortAndTrackIndex(&(diag_index[1]), &(displ[1]), tl1, 'a', TRUE); + +#if 0 + for (i=1; i<=tl1; i++) { + Log(&rLog, LOG_FORCED_DEBUG, "after sort disp[%d]=%d diag_index[%d]=%d", i, diag_index[i], i, displ[i]); + } +#endif + +#else + + des_quick_sort(displ, diag_index, tl1); + +#endif + + j = tl1 - aln_param->signif + 1; + + if (j < 1) { + j = 1; + } + + /* flag all diagonals within WINDOW of a top diagonal */ + + for (i=tl1; i>=j; i--) { + if (displ[i] > 0) { + pos = diag_index[i]; + l = (1 > pos - aln_param->window) ? + 1 : pos - aln_param->window; + m = (tl1 < pos + aln_param->window) ? + tl1 : pos + aln_param->window; + for (; l <= m; l++) + slopes[l] = 1; + } + } + + for (i=1; i<=tl1; i++) { + displ[i] = 0; + } + + curr_frag=*maxsf=0; + + for (i=1; i<=(l1-aln_param->ktup+1); ++i) { + encrypt=flag=0; + for (j=1; j<=aln_param->ktup; ++j) { + residue = seq_array[seq_no][i+j-1]; + if ((residue<0) || (residue>max_res_code)) { + flag=TRUE; + break; + } + encrypt += ((residue)*pot[j]); + } + if (flag) { + continue; + } + ++encrypt; + + vn2=zzd[encrypt]; + + flag=FALSE; + while (TRUE) { + if (!vn2) { + flag=TRUE; + break; + } + osptr=i-vn2+l2; + if (1 != slopes[osptr]) { + vn2=zzb[vn2]; + continue; + } + flen=0; + fs=aln_param->ktup; + next=*maxsf; + + /* + * A-loop + */ + + while (TRUE) { + if (!next) { + ++curr_frag; + if (curr_frag >= 2*max_aln_length) { + Log(&rLog, LOG_VERBOSE, "(Partial alignment)"); + goto free_and_exit; /* Yesss! Always wanted to + * use a goto (AW) */ + } + displ[osptr]=curr_frag; + put_frag(fs, i, vn2, flen, curr_frag, &next, maxsf, accum); + + } else { + tv1=accum[1][next]; + tv2=accum[2][next]; + + if (frag_rel_pos(i, vn2, tv1, tv2, aln_param->ktup)) { + if (i-vn2 == accum[1][next]-accum[2][next]) { + if (i > accum[1][next]+(aln_param->ktup-1)) { + fs = accum[0][next]+aln_param->ktup; + } else { + rmndr = i-accum[1][next]; + fs = accum[0][next]+rmndr; + } + flen=next; + next=0; + continue; + + } else { + if (0 == displ[osptr]) { + subt1=aln_param->ktup; + } else { + if (i > accum[1][displ[osptr]]+(aln_param->ktup-1)) { + subt1=accum[0][displ[osptr]]+aln_param->ktup; + } else { + rmndr=i-accum[1][displ[osptr]]; + subt1=accum[0][displ[osptr]]+rmndr; + } + } + subt2=accum[0][next] - aln_param->wind_gap + aln_param->ktup; + if (subt2>subt1) { + flen=next; + fs=subt2; + } else { + flen=displ[osptr]; + fs=subt1; + } + next=0; + continue; + } + } else { + next=accum[4][next]; + continue; + } + } + break; + } + /* + * End of Aloop + */ + + vn2=zzb[vn2]; + } + } + +free_and_exit: + CKFREE(displ); + CKFREE(slopes); + CKFREE(diag_index); + + return; +} +/* end of pair_align */ + + + +/** + * + * Will compute ktuple scores and store in tmat + * Following values will be set: tmat[i][j], where + * istart <= i <iend + * and + * jstart <= j < jend + * i.e. zero-offset + * tmat data members have to be preallocated + * + * if ktuple_param_t *aln_param == NULL defaults will be used + */ +void +KTuplePairDist(symmatrix_t *tmat, mseq_t *mseq, + int istart, int iend, + int jstart, int jend, + ktuple_param_t *param_override, + progress_t *prProgress, + unsigned long int *ulStepNo, unsigned long int ulTotalStepNo) +{ + /* this first group of variables were previously static + and hence un-parallelisable */ + char **seq_array; + int maxsf; + int **accum; + int max_aln_length; + /* divide score with length of smallest sequence */ + int *zza, *zzb, *zzc, *zzd; + int private_step_no = 0; + + int i, j, dsr; + double calc_score; + int max_res_code = -1; + + int max_seq_len; + int *seqlen_array; + /* progress_t *prProgress; */ + /* int uStepNo, uTotalStepNo; */ + ktuple_param_t aln_param = default_protein_param; + bool bPrintCR = (rLog.iLogLevelEnabled<=LOG_VERBOSE) ? FALSE : TRUE; + + + if(prProgress == NULL) { + NewProgress(&prProgress, LogGetFP(&rLog, LOG_INFO), + "Ktuple-distance calculation progress", bPrintCR); + } + + /* conversion to old style data types follows + * + */ + + seqlen_array = (int*) CKMALLOC((mseq->nseqs+1) * sizeof(int)); + for (i=0; i<mseq->nseqs; i++) { + seqlen_array[i+1] = mseq->sqinfo[i].len; + } + + /* setup alignment parameters + */ + if (SEQTYPE_PROTEIN == mseq->seqtype) { + DNAFLAG = FALSE; + max_res_code = strlen(AMINO_ACID_CODES)-2; + aln_param = default_protein_param; + + } else if (SEQTYPE_RNA == mseq->seqtype || SEQTYPE_DNA == mseq->seqtype) { + DNAFLAG = TRUE; + max_res_code = strlen(NUCLEIC_ACID_CODES)-2; + aln_param = default_dna_param; + + } else { + Log(&rLog, LOG_FATAL, "Internal error in %s: Unknown sequence type.", __FUNCTION__); + } + + if (NULL!=param_override) { + aln_param.ktup = param_override->ktup; + aln_param.wind_gap = param_override->wind_gap; + aln_param.signif = param_override->signif; + aln_param.window = param_override->window; + } + + /*LOG_DEBUG("DNAFLAG = %d max_res_code = %d", DNAFLAG, max_res_code);*/ + + /* convert mseq to clustal's old-style int encoded sequences (unit-offset) + */ + max_aln_length = 0; + max_seq_len = 0; + seq_array = (char **) CKMALLOC((mseq->nseqs+1) * sizeof(char *)); + seq_array[0] = NULL; + /* FIXME check that non of the seqs is smaller than ktup (?). + * Otherwise segfault occurs + */ + for (i=0; i<mseq->nseqs; i++) { + seq_array[i+1] = (char *) CKMALLOC((seqlen_array[i+1]+2) * sizeof (char));; + } + for (i=0; i<mseq->nseqs; i++) { + /*LOG_DEBUG("calling encode with seq_array[%d+1] len=%d and seq=%s", + i, seqlen_array[i+1], mseq->seq[i]);*/ + if (TRUE == DNAFLAG) { + encode(&(mseq->seq[i][-1]), seq_array[i+1], + seqlen_array[i+1], NUCLEIC_ACID_CODES); + } else { + encode(&(mseq->seq[i][-1]), seq_array[i+1], + seqlen_array[i+1], AMINO_ACID_CODES); + } + + if (seqlen_array[i+1]>max_seq_len) { + max_seq_len = seqlen_array[i+1]; + } + } + max_aln_length = max_seq_len * 2; + /* see sequence.c in old source */ + + /* FIXME: short sequences can cause seg-fault + * because max_aln_length can get shorter + * than (max_res_code+1)^k + * FS, r222->r223 */ + max_aln_length = max_aln_length > pow((max_res_code+1), aln_param.ktup)+1 ? + max_aln_length : pow((max_res_code+1), aln_param.ktup)+1; + + /* + * + * conversion to old style clustal done (in no time) */ + + + accum = (int **) CKCALLOC(5, sizeof (int *)); + for (i=0;i<5;i++) { + accum[i] = (int *) CKCALLOC((2*max_aln_length+1), sizeof(int)); + } + zza = (int *) CKCALLOC( (max_aln_length+1), sizeof(int)); + zzb = (int *) CKCALLOC( (max_aln_length+1), sizeof(int)); + zzc = (int *) CKCALLOC( (max_aln_length+1), sizeof(int)); + zzd = (int *) CKCALLOC( (max_aln_length+1), sizeof(int)); + + /* estimation of total number of steps (if istart and jstart are + * both 0) (now handled in the calling routine) + */ + /* uTotalStepNo = iend*jend - iend*iend/2 + iend/2; + uStepNo = 0; */ + /*LOG_DEBUG("istart=%d iend=%d jstart=%d jend=%d", istart, iend, jstart, jend);*/ + + for (i=istart+1; i<=iend; ++i) { + /* by definition a sequence compared to itself should give + a score of 0. AW */ + SymMatrixSetValue(tmat, i-1, i-1, 0.0); + make_ptrs(zza, zzc, i, seqlen_array[i], aln_param.ktup, max_res_code, seq_array); + +#ifdef HAVE_OPENMP + #pragma omp critical(ktuple) +#endif + { + ProgressLog(prProgress, *ulStepNo, ulTotalStepNo, FALSE); + } + + for (j=MAX(i+1, jstart+1); j<=jend; ++j) { + (*ulStepNo)++; + private_step_no++; + /*LOG_DEBUG("comparing pair %d:%d", i, j);*/ + + make_ptrs(zzb, zzd, j, seqlen_array[j], aln_param.ktup, max_res_code, seq_array); + pair_align(i, seqlen_array[i], seqlen_array[j], max_res_code, &aln_param, + seq_array, &maxsf, accum, max_aln_length, zza, zzb, zzc, zzd); + + if (!maxsf) { + calc_score=0.0; + } else { + calc_score=(double)accum[0][maxsf]; + if (percent) { + dsr=(seqlen_array[i]<seqlen_array[j]) ? + seqlen_array[i] : seqlen_array[j]; + calc_score = (calc_score/(double)dsr) * 100.0; + } + } + + /* printf("%d %d %d\n", i-1, j-1, (100.0 - calc_score)/100.0); */ + SymMatrixSetValue(tmat, i-1, j-1, (100.0 - calc_score)/100.0); + + /* the function allows you not to compute the full matrix. + * here we explicitely make the resulting matrix a + * rectangle, i.e. we always set full rows. in other + * words, if we don't complete the full matrix then we + * don't have a full symmetry. so only use the defined + * symmetric part. AW + */ + /*LOG_DEBUG("setting %d : %d = %f", j, i, tmat[i][j]);*/ + /* not needed anymore since we use symmatrix_t + if (j<=iend) { + tmat[j][i] = tmat[i][j]; + } + */ +#ifdef HAVE_OPENMP + #pragma omp critical(ktuple) +#endif + { + Log(&rLog, LOG_DEBUG, "K-tuple distance for sequence pair %d:%d = %lg", + i, j, SymMatrixGetValue(tmat, i-1, j-1)); + } + } + } + /* + Log(&rLog, LOG_FORCED_DEBUG, "uTotalStepNo=%d for istart=%d iend=%d jstart=%d jend=%d", uStepNo, istart, iend, jstart, jend); + Log(&rLog, LOG_FORCED_DEBUG, "Fabian = %d", iend*jend - iend*iend/2 + iend/2); + */ + +/* printf("\n\n%d\t%d\t%d\t%d\n\n", omp_get_thread_num(), uStepNo, istart, iend); */ + + for (i=0;i<5;i++) { + CKFREE(accum[i]); + } + CKFREE(accum); + +#ifdef HAVE_OPENMP + #pragma omp critical(ktuple) +#if 0 + { + printf("steps: %d\n", private_step_no); + } +#endif +#endif + + CKFREE(zza); + CKFREE(zzb); + CKFREE(zzc); + CKFREE(zzd); + + free(seqlen_array); + + for (i=1; i<=mseq->nseqs; i++) { + CKFREE(seq_array[i]); + } + CKFREE(seq_array); +} +/* end of KTuplePairDist */