--- /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