+++ /dev/null
-/* -*- 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 */
git://source.jalview.org / jabaws.git / blobdiff