********************************************************************/
/*
- * RCS $Id: hhalignment-C.h 236 2011-04-14 11:30:04Z fabian $
+ * RCS $Id: hhalignment-C.h 316 2016-12-16 16:14:39Z fabian $
*/
#include "hhhmm.h"
#endif
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+//#include "new_new.h" /* memory tracking */
+
enum {KEEP_NOT = 0, KEEP_CONDITIONALLY, KEEP_ALWAYS};
//printf(">>>>>>>>%s:%s:%d: maxseq=%d, maxres=%d\n", __FUNCTION__, __FILE__, __LINE__, maxseq, maxres); /* (FS) */
longname = new(char[DESCLEN]);
- sname = new(char*[maxseq+2]); /* MR1 */
- seq = new(char*[maxseq+2]); /* MR1 */
- l = new(int[maxres]);
- X = new(char*[maxseq+2]); /* MR1 */
- I = new(short unsigned int*[maxseq+2]); /* MR1 */
- keep = new(char[maxseq+2]); /* MR1 */
- display = new(char[maxseq+2]); /* MR1 */
- wg = new(float[maxseq+2]); /* MR1 */
- nseqs = new(int[maxres+2]); /* MR1 */
+ sname = new char*[maxseq+2]; /* MR1 */
+ seq = new char*[maxseq+2]; /* MR1 */
+ l = new int[maxres];
+ X = new char*[maxseq+2]; /* MR1 */
+ I = new short unsigned int*[maxseq+2]; /* MR1 */
+ keep = new char[maxseq+2]; /* MR1 */
+ display = new char[maxseq+2]; /* MR1 */
+ wg = new float[maxseq+2]; /* MR1 */
+ nseqs = new int[maxres+2]; /* MR1 */
N_in=L=0;
nres=NULL; // number of residues per sequence k
first=NULL; // first residue in sequence k
int k; // Index of sequence being read currently (first=0)
char line[LINELEN]=""; // input line
//char cur_seq[MAXCOL]; // Sequence currently read in
- char *cur_seq=new(char[par.maxColCnt]);
+ char *cur_seq=new char[par.maxColCnt];
char* cur_name; // Sequence currently read in
int linenr=0; // current line number in input file
char skip_sequence=0;
}
// Create space for residues and paste new sequence in
- seq[k]=new(char[strlen(cur_seq)+2]);
+ seq[k]=new char[strlen(cur_seq)+2];
if (!seq[k]) MemoryError("array for input sequences");
- X[k]=new(char[strlen(cur_seq)+2]);
+ X[k]=new char[strlen(cur_seq)+2];
if (!X[k]) MemoryError("array for input sequences");
- I[k]=new(short unsigned int[strlen(cur_seq)+2]);
+ I[k]=new short unsigned int[strlen(cur_seq)+2];
if (!I[k]) MemoryError("array for input sequences");
strcpy(seq[k],cur_seq);
}
// store sequence name
if (v>=4) printf("Reading seq %-16.16s k=%3i n_displ=%3i display[k]=%i keep[k]=%i\n",cur_name,k,n_display,display[k],keep[k]);
- sname[k] = new(char[strlen(cur_name)+1]);
+ sname[k] = new char[strlen(cur_name)+1];
if (!sname[k]) {MemoryError("array for sequence names");}
strcpy(sname[k],cur_name);
} // end if(line contains sequence name)
if (k>=0) //if at least one sequence was read in
{
- seq[k]=new(char[strlen(cur_seq)+2]);
+ seq[k]=new char[strlen(cur_seq)+2];
if (!seq[k]) MemoryError("array for input sequences");
- X[k]=new(char[strlen(cur_seq)+2]);
+ X[k]=new char[strlen(cur_seq)+2];
if (!X[k]) MemoryError("array for input sequences");
- I[k]=new(short unsigned int[strlen(cur_seq)+2]);
+ I[k]=new short unsigned int[strlen(cur_seq)+2];
if (!I[k]) MemoryError("array for input sequences");
strcpy(seq[k],cur_seq);
}
/*static short unsigned int h[MAXSEQ];*/
/*short*/ unsigned int *h = NULL; /* short may lead to overflow for long alignments, FS, r235 -> r236 */
- h = new(/*short*/ unsigned int[N_in+2]); /* short -> overflow, FS, r235 -> r236 */
+ h = new /*short*/ unsigned int[N_in+2]; /* short -> overflow, FS, r235 -> r236 */
float *percent_gaps = NULL; /* FS, 2010-Nov */
char *match_state = NULL; /* FS, 2010-Nov */
i--;
if (L!=i && L!=/*MAXRES*/par.maxResLen-2 && !unequal_lengths) unequal_lengths=k; //sequences have different lengths
L=imin(L,i);
- }
+ } // end for (k)
if (unequal_lengths) break;
//Replace GAP with ENDGAP for all end gaps /* MR1 */
had to move declaration of float *percent_gaps out of switch()
*/
//float percent_gaps[MAXCOL]; //percentage of gaps in column k (with weighted sequences)
- percent_gaps = new(float[par.maxColCnt]);
+ percent_gaps = new float[par.maxColCnt];
//determine number of columns L in alignment
L=strlen(seq[kfirst])-1;
// Conversion to integer representation, checking for unequal lengths and initialization
- if (nres==NULL) nres=new(int[N_in]);
+ if (nres==NULL) nres=new int[N_in];
+#if 0
+//#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static), private(l)
+#endif
for (k=0; k<N_in; k++)
{
if (!keep[k]) continue;
for (k=0; k<N_in; k++) if (keep[k]) nl[ (int)X[k][l]]++;
for (a=0; a<20; a++) if(nl[a]) naa++;
if (!naa) naa=1; //naa=0 when column consists of only gaps and Xs (=ANY)
+#if 0
+//#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static)
+#endif
for (k=0; k<N_in; k++)
if (keep[k] && (X[k][l]<20) )
{
}
// Add up percentage of gaps
+#if 0
+//#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(k)
+#endif
for (l=1; l<=L; l++)
{
float res=0;
}
i++;
this->l[i]=l;
+#if 0
+//#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static)
+#endif
for (k=0; k<N_in; k++)
{
if (keep[k])
had to move declaration of float *percent_gaps out of switch()
*/
//char match_state[MAXCOL]; //1: column assigned to match state 0: insert state
- match_state = new(char[par.maxColCnt]);
+ match_state = new char[par.maxColCnt];
// Determine number of columns L in alignment
L=strlen(seq[0]+1);
if (par.mark) return n_display;
- char *dummy = new(char[N_in+1]);
+ char *dummy = new char[N_in+1];
int vtmp=v, seqid;
v=0;
n_display=0;
{
// In the beginnning, keep[k] is 1 for all regular amino acid sequences and 0 for all others (ss_conf, ss_pred,...)
// In the end, keep[k] will be 1 for all regular representative sequences kept in the alignment, 0 for all others
- char* in=new(char[N_in+1]); // in[k]=1: seq k has been accepted; in[k]=0: seq k has not yet been accepted at current seqid
- char* inkk=new(char[N_in+1]); // inkk[k]=1 iff in[ksort[k]]=1 else 0;
- int* Nmax=new(int[L+2]); // position-dependent maximum-sequence-identity threshold for filtering /* MR1, used to be called idmax*/
- int* idmaxwin=new(int[L+2]); // minimum value of idmax[i-WFIL,i+WFIL]
- int* seqid_prev=new(int[N_in+1]); // maximum-sequence-identity threshold used in previous round of filtering (with lower seqid)
- int* N=new(int[L+2]); // N[i] number of already accepted sequences at position i
+ char* in=new char[N_in+1]; // in[k]=1: seq k has been accepted; in[k]=0: seq k has not yet been accepted at current seqid
+ char* inkk=new char[N_in+1]; // inkk[k]=1 iff in[ksort[k]]=1 else 0;
+ int* Nmax=new int[L+2]; // position-dependent maximum-sequence-identity threshold for filtering /* MR1, used to be called idmax*/
+ int* idmaxwin=new int[L+2]; // minimum value of idmax[i-WFIL,i+WFIL]
+ int* seqid_prev=new int[N_in+1]; // maximum-sequence-identity threshold used in previous round of filtering (with lower seqid)
+ int* N=new int[L+2]; // N[i] number of already accepted sequences at position i
const int WFIL=25; // see previous line
int diffNmax=Ndiff; // current maximum difference of Nmax[i] and Ndiff /* MR1 */
int i; // counts residues
int n; // number of sequences accepted so far
-
// Initialize in[k]
for (n=k=0; k<N_in; k++) if (keep[k]==KEEP_ALWAYS) {in[k]=2/*KEEP_ALWAYS??*/; n++;} else in[k]=0;
// Determine first[k], last[k]?
if (first==NULL)
{
- first=new(int[N_in]);// first non-gap position in sequence k
- last =new(int[N_in]);// last non-gap position in sequence k
+ first=new int[N_in];// first non-gap position in sequence k
+ last =new int[N_in];// last non-gap position in sequence k
for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])
{
for (i=1; i<=L; i++) if (X[k][i]<NAA) break;
}
// Determine number of residues nres[k]?
- if ( (nres==NULL) || (sizeof(nres)<N_in*sizeof(int)) )
+ // KB: memory leak as sizeof(nres) just gives the size of an int pointer
+ //if ( (nres==NULL) || (sizeof(nres)<N_in*sizeof(int)) )
+ if (nres==NULL)
{
- nres=new(int[N_in]);
+ nres=new int[N_in];
for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])
{
int nr=0;
// Sort sequences according to length; afterwards, nres[ksort[kk]] is sorted by size
if (ksort==NULL)
{
- ksort=new(int[N_in]); // never reuse alignment object for new alignment with more sequences
+ ksort=new int[N_in]; // never reuse alignment object for new alignment with more sequences
for (k=0; k<N_in; k++) ksort[k]=k;
QSortInt(nres,ksort,kfirst+1,N_in-1,-1); //Sort sequences after kfirst (query) in descending order
}
if (diff>=qdiff_max) {keep[k]=KEEP_NOT; continue;} // too different from query? => reject once and for all
}
// printf(" qsc=%6.2f qid=%6.2f \n",qsc_sum/nres[k],100.0*(1.0-(float)(diff)/nres[k]));
- }
+ } // end for (k)
if (seqid1>seqid2)
{
// printf("\n");
// Loop over all candidate sequences kk (-> k)
+#if 0
+//#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(k, i, diff_min_frac, jj, j, first_kj, last_kj, cov_kj, diff_suff, diff)
+#endif
for (kk=0; kk<N_in; kk++)
{
if (inkk[kk]) continue; // seq k already accepted
if (keep[k]==KEEP_ALWAYS) {inkk[kk]=2; continue;} // accept all marked sequences (no n++, since this has been done already)
// Calculate max-seq-id threshold seqidk for sequence k (as maximum over idmaxwin[i])
- if (seqid>=100) {in[k]=inkk[kk]=1; n++; continue;}
+// if (seqid>=100) {in[k]=inkk[kk]=1; n++; continue;}
+ if (seqid>=100) {
+ in[k]=inkk[kk]=1;
+#if 0
+//#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ n++;
+ continue;
+ }
float seqidk=seqid1;
for (i=first[k]; i<=last[k]; i++)
if (idmaxwin[i]>seqidk) seqidk=idmaxwin[i];
if (jj>=kk) // did loop reach end? => accept k. Otherwise reject k (the shorter of the two)
{
in[k]=inkk[kk]=1;
+#if 0
+//#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
n++;
- for (i=first[k]; i<=last[k]; i++) N[i]++; // update number of sequences at position i
+ for (i=first[k]; i<=last[k]; i++) {
+#if 0
+//#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ N[i]++; // update number of sequences at position i
+ }
// printf("%i %20.20s accepted\n",k,sname[k]);
}
// else
const int Ndiff_core=0;
int n;
HMM qcore;
- char* coreseq=new(char[N_in]); // coreseq[k]=1 if sequence belongs to core of alignment (i.e. it is very similar to query)
+ char* coreseq=new char[N_in]; // coreseq[k]=1 if sequence belongs to core of alignment (i.e. it is very similar to query)
for (int k=0; k<N_in; k++) coreseq[k]=keep[k]; // Copy keep[] into coreseq[]
// Remove sequences with seq. identity larger than seqid percent (remove the shorter of two)
int i; // column in query alignment
int a; // amino acid (0..19)
int n=1; // number of sequences that passed filter
- float** logodds=new(float*[L+1]); // log-odds ratios for HMM qcore
+ float** logodds=new float*[L+1]; // log-odds ratios for HMM qcore
char gap; // 1: previous state in seq k was a gap 0: previous state in seq k was an amino acid
float score; // score of sequence k aligned with qcore
// Determine first[k], last[k]?
if (first==NULL)
{
- first=new(int[N_in]);// first non-gap position in sequence k
- last =new(int[N_in]);// last non-gap position in sequence k
+ first=new int[N_in];// first non-gap position in sequence k
+ last =new int[N_in];// last non-gap position in sequence k
for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])
{
for (i=1; i<=L; i++) if (X[k][i]<NAA) break;
// Determine number of residues nres[k]?
if (nres==NULL)
{
- nres=new(int[N_in]);
+ nres=new int[N_in];
for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])
{
int nr=0;
for (k=0; k<N_in; k++) X[k][0]=ENDGAP; // make sure that sequences ENTER subalignment j for j=1
for (k=0; k<N_in; k++) X[k][L+1]=ENDGAP; // does it have an influence?
-#ifdef HAVE_OPENMP
+#if 0
+//#ifdef HAVE_OPENMP
if(par.wg != 1)
{
#pragma omp parallel sections
q.sname[q.ncons]=new(char[10]);
if (!q.sname[q.ncons]) {MemoryError("array of names for displayed sequences");}
strcpy(q.sname[q.ncons],"Consensus");
- q.seq[q.ncons]=new(char[L+2]);
+ q.seq[q.ncons]=new char[L+2];
if (!q.seq[q.ncons]) {MemoryError("array of names for displayed sequences");}
}
if (par.cons)
{
// Reserve space for consensus sequence as first sequence in alignment
q.nfirst=n++; kfirst=-1;
- q.sname[q.nfirst]=new(char[strlen(name)+11]);
+ q.sname[q.nfirst]=new char[strlen(name)+11];
if (!q.sname[q.nfirst]) {MemoryError("array of names for displayed sequences");}
strcpy(q.sname[q.nfirst],name);
strcat(q.sname[q.nfirst],"_consensus");
- q.seq[q.nfirst]=new(char[L+2]);
+ q.seq[q.nfirst]=new char[L+2];
if (!q.seq[q.nfirst]) {MemoryError("array of names for displayed sequences");}
}
// Calculate consensus amino acids using similarity matrix
else if (k==kfirst) nn=q.nfirst=n++;
else nn=n++;
// strcut(sname[k]," "); // delete rest of name line beginning with two spaces " " // Why this?? Problem for pdb seqs without chain
- q.sname[nn]=new(char[strlen(sname[k])+1]);
+ q.sname[nn]=new char[strlen(sname[k])+1];
if (!q.sname[nn]) {MemoryError("array of names for displayed sequences");}
strcpy(q.sname[nn],sname[k]);
- q.seq[nn]=new(char[strlen(seq[k])+1]);
+ q.seq[nn]=new char[strlen(seq[k])+1];
if (!q.seq[nn]) {MemoryError("array of names for displayed sequences");}
strcpy(q.seq[nn],seq[k]);
}
//float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i
float *wi=NULL; // weight of sequence k in column i, calculated from subalignment i
//float Neff[MAXRES]; // diversity of subalignment i
- float *Neff = new(float[par.maxResLen]); // diversity of subalignment i
+ float *Neff = new float[par.maxResLen]; // diversity of subalignment i
int nseqi=0; // number of sequences in subalignment i
int ncol=0; // number of columns j that contribute to Neff[i]
char change; // has the set of sequences in subalignment changed? 0:no 1:yes
float fj[NAA+3]; // to calculate entropy
float sum;
- wi = new(float[N_in+2]);
+ wi = new float[N_in+2];
// Global weights?
if (par.wg==1)
- for (k=0; k<N_in; k++) wi[k]=wg[k];
+ for (k=0; k<N_in; k++)
+ wi[k]=wg[k];
// Initialization
q.Neff_HMM=0.0f;
Neff[0]=0.0; // if the first column has no residues (i.e. change==0), Neff[i]=Neff[i-1]=Neff[0]
- n = new(int*[L+2]);
- for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);
- for (j=1; j<=L; j++)
- for (a=0; a<NAA+3; a++) n[j][a]=0;
+ n = new int*[L+2];
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(a)
+#endif
+ for (j=1; j<=L; j++) {
+ n[j]=new(int[NAA+3]);
+ for (a=0; a<NAA+3; a++)
+ n[j][a]=0;
+ }
//////////////////////////////////////////////////////////////////////////////////////////////
// Main loop through alignment columns
for (i=1; i<=L; i++) // Calculate wi[k] at position i as well as Neff[i]
- {
-
+ {
if (par.wg==0)
- {
-
- change=0;
- // Check all sequences k and update n[j][a] and ri[j] if necessary
- for (k=0; k<N_in; k++)
- {
- if (!in[k]) continue;
- if (X[k][i-1]>=ANY && X[k][i]<ANY)
- { // ... if sequence k was NOT included in i-1 and has to be included for column i
- change=1;
- nseqi++;
- for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]++;
- }
- else if (X[k][i-1]<ANY && X[k][i]>=ANY)
- { // ... if sequence k WAS included in i-1 and has to be thrown out for column i
- change=1;
- nseqi--;
- for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]--;
- }
- } //end for (k)
- nseqs[i]=nseqi;
-
- // If subalignment changed: update weights wi[k] and Neff[i]
- if (change)
- {
- // Initialize weights and numbers of residues for subalignment i
- ncol=0;
- for (k=0; k<N_in; k++) wi[k]=1E-8; // for pathological alignments all wi[k] can get 0; /* MR1 */
-
- // sum wi[k] over all columns j and sequences k of subalignment
- for (j=1; j<=L; j++)
- {
- // do at least a fraction MAXENDGAPFRAC of sequences in subalignment contain an end gap in j?
- if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
- naa=0; for (a=0; a<20; a++) if(n[j][a]) naa++;
- if (naa==0) continue;
- ncol++;
- for (k=0; k<N_in; k++)
- {
- if (in[k] && X[k][i]<ANY && X[k][j]<ANY)
- {
- // if (!n[j][ (int)X[k][j]]) {fprintf(stderr,"Error: Mi=%i: n[%i][X[%i]]=0! (X[%i]=%i)\n",i,j,k,k,X[k][j]);}
- wi[k]+=1.0/float(n[j][ (int)X[k][j] ]*naa);
- }
- }
- }
-
- // Check whether number of columns in subalignment is sufficient
- if (ncol<NCOLMIN)
- // Take global weights
- for (k=0; k<N_in; k++)
- if(in[k] && X[k][i]<ANY) wi[k]=wg[k]; else wi[k]=0.0;
-
- // Calculate Neff[i]
- Neff[i]=0.0;
- for (j=1; j<=L; j++)
- {
- // do at least a fraction MAXENDGAPFRA of sequences in subalignment contain an end gap in j?
- if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
- for (a=0; a<20; a++) fj[a]=0;
- for (k=0; k<N_in; k++)
- if (in[k] && X[k][i]<ANY && X[k][j]<ANY)
- fj[ (int)X[k][j] ]+=wi[k];
- NormalizeTo1(fj,NAA);
- for (a=0; a<20; a++)
- if (fj[a]>1E-10) Neff[i]-=fj[a]*fast_log2(fj[a]);
- }
- if (ncol>0) Neff[i]=pow(2.0,Neff[i]/ncol); else Neff[i]=1.0;
-
- }
-
- else //no update was necessary; copy values for i-1
- {
- Neff[i]=Neff[i-1];
- }
- }
+ {
+ change=0;
+ // Check all sequences k and update n[j][a] and ri[j] if necessary
+ for (k=0; k<N_in; k++)
+ {
+ if (!in[k]) continue;
+ if (X[k][i-1]>=ANY && X[k][i]<ANY)
+ { // ... if sequence k was NOT included in i-1 and has to be included for column i
+ change=1;
+ nseqi++;
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static)
+#endif
+ for (int j=1; j<=L; j++)
+ n[j][ (int)X[k][j]]++;
+ }
+ else if (X[k][i-1]<ANY && X[k][i]>=ANY)
+ { // ... if sequence k WAS included in i-1 and has to be thrown out for column i
+ change=1;
+ nseqi--;
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static)
+#endif
+ for (int j=1; j<=L; j++)
+ n[j][ (int)X[k][j]]--;
+ }
+ } //end for (k)
+ nseqs[i]=nseqi;
+
+ // If subalignment changed: update weights wi[k] and Neff[i]
+ if (change)
+ {
+ // Initialize weights and numbers of residues for subalignment i
+ ncol=0;
+ for (k=0; k<N_in; k++)
+ wi[k]=1E-8; // for pathological alignments all wi[k] can get 0; /* MR1 */
+
+ // sum wi[k] over all columns j and sequences k of subalignment
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(naa, a, k)
+#endif
+ for (j=1; j<=L; j++)
+ {
+ // do at least a fraction MAXENDGAPFRAC of sequences in subalignment contain an end gap in j?
+ if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
+ naa=0;
+ for (a=0; a<20; a++)
+ if(n[j][a])
+ naa++;
+ if (naa==0) continue;
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ ncol++;
+ for (k=0; k<N_in; k++)
+ {
+ if (in[k] && X[k][i]<ANY && X[k][j]<ANY)
+ {
+ // if (!n[j][ (int)X[k][j]]) {fprintf(stderr,"Error: Mi=%i: n[%i][X[%i]]=0! (X[%i]=%i)\n",i,j,k,k,X[k][j]);}
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ wi[k]+=1.0/float(n[j][ (int)X[k][j] ]*naa);
+ }
+ }
+ } // end for (j)
+
+ // Check whether number of columns in subalignment is sufficient
+ if (ncol<NCOLMIN)
+ // Take global weights
+ for (k=0; k<N_in; k++)
+ if(in[k] && X[k][i]<ANY)
+ wi[k]=wg[k];
+ else wi[k]=0.0;
+
+ // Calculate Neff[i]
+ Neff[i]=0.0;
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(a, k, fj)
+#endif
+ for (j=1; j<=L; j++)
+ {
+ // do at least a fraction MAXENDGAPFRA of sequences in subalignment contain an end gap in j?
+ if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
+ for (a=0; a<20; a++) fj[a]=0;
+ for (k=0; k<N_in; k++)
+ if (in[k] && X[k][i]<ANY && X[k][j]<ANY)
+ fj[ (int)X[k][j] ]+=wi[k];
+ NormalizeTo1(fj,NAA);
+ for (a=0; a<20; a++)
+ if (fj[a]>1E-10)
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ Neff[i]-=fj[a]*fast_log2(fj[a]);
+ }
+ if (ncol>0)
+ Neff[i]=pow(2.0,Neff[i]/ncol);
+ else
+ Neff[i]=1.0;
+ } // end change
+
+ else //no update was necessary; copy values for i-1
+ {
+ Neff[i]=Neff[i-1];
+ }
+ } // end par.wg==0
// Calculate amino acid frequencies q.f[i][a] from weights wi[k]
- for (a=0; a<20; a++) q.f[i][a]=0;
- for (k=0; k<N_in; k++) if (in[k]) q.f[i][ (int)X[k][i] ]+=wi[k];
+ for (a=0; a<20; a++)
+ q.f[i][a]=0;
+ for (k=0; k<N_in; k++)
+ if (in[k])
+ q.f[i][ (int)X[k][i] ]+=wi[k];
NormalizeTo1(q.f[i],NAA,pb);
// Calculate transition probabilities from M state
q.tr[i][M2M]=q.tr[i][M2D]=q.tr[i][M2I]=0.0;
for (k=0; k<N_in; k++) //for all sequences
- {
- if (!in[k]) continue;
- //if input alignment is local ignore transitions from and to end gaps
- if (X[k][i]<ANY) //current state is M
- {
- if (I[k][i]) //next state is I
- q.tr[i][M2I]+=wi[k];
- else if (X[k][i+1]<=ANY) //next state is M
- q.tr[i][M2M]+=wi[k];
- else if (X[k][i+1]==GAP) //next state is D
- q.tr[i][M2D]+=wi[k];
- }
- } // end for(k)
+ {
+ if (!in[k]) continue;
+ //if input alignment is local ignore transitions from and to end gaps
+ if (X[k][i]<ANY) //current state is M
+ {
+ if (I[k][i]) //next state is I
+ q.tr[i][M2I]+=wi[k];
+ else if (X[k][i+1]<=ANY) //next state is M
+ q.tr[i][M2M]+=wi[k];
+ else if (X[k][i+1]==GAP) //next state is D
+ q.tr[i][M2D]+=wi[k];
+ }
+ } // end for(k)
// Normalize and take log
sum = q.tr[i][M2M]+q.tr[i][M2I]+q.tr[i][M2D]+FLT_MIN;
- q.tr[i][M2M]=log2(q.tr[i][M2M]/sum);
- q.tr[i][M2I]=log2(q.tr[i][M2I]/sum);
- q.tr[i][M2D]=log2(q.tr[i][M2D]/sum);
+ q.tr[i][M2M]=Log2(q.tr[i][M2M]/sum);
+ q.tr[i][M2I]=Log2(q.tr[i][M2I]/sum);
+ q.tr[i][M2D]=Log2(q.tr[i][M2D]/sum);
// for (k=0; k<N_in; k++) if (in[k]) w[k][i]=wi[k];
- }
- // DD TODO:fill in all the missing Neff values
+ } // end loop through alignment columns i
+ // DD TODO:fill in all the missing Neff values
- // end loop through alignment columns i
//////////////////////////////////////////////////////////////////////////////////////////////
delete[](wi); wi=NULL;
// delete n[][]
for (j=1; j<=L; j++){
- delete[](n[j]); (n[j]) = NULL;
+ delete[](n[j]);
+ (n[j]) = NULL;
}
delete[](n); (n) = NULL;
q.Neff_M[0]=99.999; // Neff_av[0] is used for calculation of transition pseudocounts for the start state
// Set emission probabilities of zero'th (begin) state and L+1st (end) state to background probabilities
- for (a=0; a<20; a++) q.f[0][a]=q.f[L+1][a]=pb[a];
+ for (a=0; a<20; a++)
+ q.f[0][a]=q.f[L+1][a]=pb[a];
// Assign Neff_M[i] and calculate average over alignment, Neff_M[0]
if (par.wg==1)
- {
+ {
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(a)
+#endif
for (i=1; i<=L; i++)
- {
- float sum=0.0f;
- for (a=0; a<20; a++)
- if (q.f[i][a]>1E-10) sum -= q.f[i][a]*fast_log2(q.f[i][a]);
- q.Neff_HMM+=pow(2.0,sum);
- }
+ {
+ float sum=0.0f;
+ for (a=0; a<20; a++)
+ if (q.f[i][a]>1E-10)
+ sum -= q.f[i][a]*fast_log2(q.f[i][a]);
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ q.Neff_HMM+=pow(2.0,sum);
+ }
q.Neff_HMM/=L;
float Nlim=fmax(10.0,q.Neff_HMM+1.0); // limiting Neff
- float scale=log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with inserts at specific pos
+ float scale=Log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with inserts at specific pos
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(k)
+#endif
for (i=1; i<=L; i++)
- {
- float w_M=-1.0/N_filtered;
- for (k=0; k<N_in; k++)
- if (in[k] && X[k][i]<=ANY) w_M+=wg[k];
- if (w_M<0) q.Neff_M[i]=1.0;
- else q.Neff_M[i] = Nlim - (Nlim-1.0)*fpow2(scale*w_M);
- // fprintf(stderr,"M i=%3i ncol=--- Neff_M=%5.2f Nlim=%5.2f w_M=%5.3f Neff_M=%5.2f\n",i,q.Neff_HMM,Nlim,w_M,q.Neff_M[i]);
- }
- }
+ {
+ float w_M=-1.0/N_filtered;
+ for (k=0; k<N_in; k++)
+ if (in[k] && X[k][i]<=ANY)
+ w_M+=wg[k];
+ if (w_M<0)
+ q.Neff_M[i]=1.0;
+ else
+ q.Neff_M[i] = Nlim - (Nlim-1.0)*fpow2(scale*w_M);
+ // fprintf(stderr,"M i=%3i ncol=--- Neff_M=%5.2f Nlim=%5.2f w_M=%5.3f Neff_M=%5.2f\n",i,q.Neff_HMM,Nlim,w_M,q.Neff_M[i]);
+ }
+ }
else
- {
+ {
for (i=1; i<=L; i++)
- {
- q.Neff_HMM+=Neff[i];
- q.Neff_M[i]=Neff[i];
- if (q.Neff_M[i] == 0) { q.Neff_M[i] = 1; } /* MR1 */
- }
+ {
+ q.Neff_HMM+=Neff[i];
+ q.Neff_M[i]=Neff[i];
+ if (q.Neff_M[i] == 0) {
+ q.Neff_M[i] = 1;
+ } /* MR1 */
+ }
q.Neff_HMM/=L;
- }
+ } // end par.wg==1
delete[] Neff; Neff = NULL;
//float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i
float *wi = NULL; // weight of sequence k in column i, calculated from subalignment i
//float Neff[MAXRES]; // diversity of subalignment i
- float *Neff = new(float[par.maxResLen]); // diversity of subalignment i
+ float *Neff = new float[par.maxResLen]; // diversity of subalignment i
int nseqi; // number of sequences in subalignment i
int ncol; // number of columns j that contribute to Neff[i]
float fj[NAA+3]; // to calculate entropy
float Nlim=0.0; // only for global weights
float scale=0.0; // only for global weights
- wi = new(float[N_in+2]);
+ wi = new float[N_in+2];
// Global weights?
if (par.wg==1)
{
for (k=0; k<N_in; k++) wi[k]=wg[k];
Nlim=fmax(10.0,q.Neff_HMM+1.0); // limiting Neff
- scale=log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with inserts at specific pos
+ scale=Log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with inserts at specific pos
}
// Initialization
- n = new(int*[L+2]);
+ n = new int*[L+2];
for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);
//////////////////////////////////////////////////////////////////////////////////////////////
// Normalize and take log
sum = q.tr[i][I2M]+q.tr[i][I2I];
- q.tr[i][I2M]=log2(q.tr[i][I2M]/sum);
- q.tr[i][I2I]=log2(q.tr[i][I2I]/sum);
+ q.tr[i][I2M]=Log2(q.tr[i][I2M]/sum);
+ q.tr[i][I2I]=Log2(q.tr[i][I2I]/sum);
}
// end loop through alignment columns i
//float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i
float *wi=NULL; // weight of sequence k in column i, calculated from subalignment i
//float Neff[MAXRES]; // diversity of subalignment i
- float *Neff = new(float[par.maxResLen]); // diversity of subalignment i
+ float *Neff = new float[par.maxResLen]; // diversity of subalignment i
int nseqi=0; // number of sequences in subalignment i (for DEBUGGING)
int ncol=0; // number of columns j that contribute to Neff[i]
char change; // has the set of sequences in subalignment changed? 0:no 1:yes
float Nlim=0.0; // only for global weights
float scale=0.0; // only for global weights
- wi = new(float[N_in+2]); /* FIXME: FS */
+ wi = new float[N_in+2]; /* FIXME: FS */
// Global weights?
if (par.wg==1)
{
- for (k=0; k<N_in; k++) wi[k]=wg[k];
+ for (k=0; k<N_in; k++)
+ wi[k]=wg[k];
Nlim=fmax(10.0,q.Neff_HMM+1.0); // limiting Neff
- scale=log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with dels at specific pos
+ scale=Log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with dels at specific pos
}
// Initialization
- n = new(int*[L+2]);
- for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);
- for (j=1; j<=L; j++)
- for (a=0; a<NAA+3; a++) n[j][a]=0;
-
-
+ n = new int*[L+2];
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(a)
+#endif
+ for (j=1; j<=L; j++) {
+ n[j]=new(int[NAA+3]);
+ for (a=0; a<NAA+3; a++)
+ n[j][a]=0;
+ }
//////////////////////////////////////////////////////////////////////////////////////////////
// Main loop through alignment columns
for (i=1; i<=L; i++) // Calculate wi[k] at position i as well as Neff[i]
+ {
+ if (par.wg==0) // if local weights
{
- if (par.wg==0) // if local weights
+ change=0;
+ // Check all sequences k and update n[j][a] and ri[j] if necessary
+ for (k=0; k<N_in; k++)
+ {
+ if (!in[k]) continue;
+ if (X[k][i-1]!=GAP && X[k][i]==GAP)
+ { // ... if sequence k was NOT included in i-1 and has to be included for column i
+ change=1;
+ nseqi++;
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static)
+#endif
+ for (int j=1; j<=L; j++)
+ n[j][ (int)X[k][j]]++;
+ }
+ else if (X[k][i-1]==GAP && X[k][i]!=GAP)
+ { // ... if sequence k WAS included in i-1 and has to be thrown out for column i
+ change=1;
+ nseqi--;
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static)
+#endif
+ for (int j=1; j<=L; j++)
+ n[j][ (int)X[k][j]]--;
+ }
+ } //end for (k)
+ nseqs[i]=nseqi;
+
+ // If there is no sequence in subalignment j ...
+ if (nseqi==0)
+ {
+ ncol=0;
+ Neff[i]=0.0; // effective number of sequences = 0!
+ q.tr[i][D2M]=-100000;
+ q.tr[i][D2D]=-100000;
+ continue;
+ }
+
+ // If subalignment changed: update weights wi[k] and Neff[i]
+ if (change)
+ {
+ // Initialize weights and numbers of residues for subalignment i
+ ncol=0;
+ for (k=0; k<N_in; k++)
+ wi[k]=0.0;
+
+ // sum wg[k][i] over all columns j and sequences k of subalignment
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(naa, a, k)
+#endif
+ for (j=1; j<=L; j++)
{
- change=0;
- // Check all sequences k and update n[j][a] and ri[j] if necessary
+ if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
+ naa=0;
+ for (a=0; a<20; a++)
+ if(n[j][a])
+ naa++;
+ if (naa==0) continue;
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ ncol++;
for (k=0; k<N_in; k++)
+ {
+ if (in[k] && X[k][i]==GAP && X[k][j]<ANY)
{
- if (!in[k]) continue;
- if (X[k][i-1]!=GAP && X[k][i]==GAP)
- { // ... if sequence k was NOT included in i-1 and has to be included for column i
- change=1;
- nseqi++;
- for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]++;
- }
- else if (X[k][i-1]==GAP && X[k][i]!=GAP)
- { // ... if sequence k WAS included in i-1 and has to be thrown out for column i
- change=1;
- nseqi--;
- for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]--;
- }
- } //end for (k)
- nseqs[i]=nseqi;
-
- // If there is no sequence in subalignment j ...
- if (nseqi==0)
- {
- ncol=0;
- Neff[i]=0.0; // effective number of sequences = 0!
- q.tr[i][D2M]=-100000;
- q.tr[i][D2D]=-100000;
- continue;
+ if (!n[j][ (int)X[k][j]]) {fprintf(stderr,"Error: Di=%i: n[%i][X[%i]]=0! (X[%i]=%i)\n",i,j,k,k,X[k][j]);}
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ wi[k]+=1.0/float(n[j][ (int)X[k][j] ]*naa);
}
+ }
+ } // end for (j)
+ // Check whether number of columns in subalignment is sufficient
+ if (ncol<NCOLMIN)
+ // Take global weights
+ for (k=0; k<N_in; k++)
+ if(in[k] && X[k][i]==GAP)
+ wi[k]=wg[k];
+ else
+ wi[k]=0.0;
+
+ // Calculate Neff[i]
+ Neff[i]=0.0;
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp parallel for schedule(static) private(a, k, fj)
+#endif
+ for (j=1; j<=L; j++)
+ {
+ if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
+ for (a=0; a<20; a++) fj[a]=0;
+ for (k=0; k<N_in; k++)
+ if (in[k] && X[k][i]==GAP && X[k][j]<ANY)
+ fj[ (int)X[k][j] ]+=wi[k];
+ NormalizeTo1(fj,NAA);
+ for (a=0; a<20; a++)
+ if (fj[a]>1E-10)
+//#if 0
+#ifdef HAVE_OPENMP
+#pragma omp atomic
+#endif
+ Neff[i]-=fj[a]*fast_log2(fj[a]);
+ }
+ if (ncol>0) Neff[i]=pow(2.0,Neff[i]/ncol); else Neff[i]=1.0;
- // If subalignment changed: update weights wi[k] and Neff[i]
- if (change)
- {
- // Initialize weights and numbers of residues for subalignment i
- ncol=0;
- for (k=0; k<N_in; k++) wi[k]=0.0;
-
- // sum wg[k][i] over all columns j and sequences k of subalignment
- for (j=1; j<=L; j++)
- {
- if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
- naa=0; for (a=0; a<20; a++) if(n[j][a]) naa++;
- if (naa==0) continue;
- ncol++;
- for (k=0; k<N_in; k++)
- {
- if (in[k] && X[k][i]==GAP && X[k][j]<ANY)
- {
- if (!n[j][ (int)X[k][j]]) {fprintf(stderr,"Error: Di=%i: n[%i][X[%i]]=0! (X[%i]=%i)\n",i,j,k,k,X[k][j]);}
- wi[k]+=1.0/float(n[j][ (int)X[k][j] ]*naa);
- }
- }
- }
-
- // Check whether number of columns in subalignment is sufficient
- if (ncol<NCOLMIN)
- // Take global weights
- for (k=0; k<N_in; k++)
- if(in[k] && X[k][i]==GAP) wi[k]=wg[k]; else wi[k]=0.0;
-
- // Calculate Neff[i]
- Neff[i]=0.0;
- for (j=1; j<=L; j++)
- {
- if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
- for (a=0; a<20; a++) fj[a]=0;
- for (k=0; k<N_in; k++)
- if (in[k] && X[k][i]==GAP && X[k][j]<ANY)
- fj[ (int)X[k][j] ]+=wi[k];
- NormalizeTo1(fj,NAA);
- for (a=0; a<20; a++)
- if (fj[a]>1E-10) Neff[i]-=fj[a]*fast_log2(fj[a]);
- }
- if (ncol>0) Neff[i]=pow(2.0,Neff[i]/ncol); else Neff[i]=1.0;
-
- }
+ }
- else //no update was necessary; copy values for i-1
- {
- Neff[i]=Neff[i-1];
- }
+ else //no update was necessary; copy values for i-1
+ {
+ Neff[i]=Neff[i-1];
+ }
- // Calculate transition probabilities from D state
- q.tr[i][D2M]=q.tr[i][D2D]=0.0;
- for (k=0; k<N_in; k++) //for all sequences
- {
- if (in[k] && X[k][i]==GAP) //current state is D
- {
- if (X[k][i+1]==GAP) //next state is D
- q.tr[i][D2D]+=wi[k];
- else if (X[k][i+1]<=ANY) //next state is M
- q.tr[i][D2M]+=wi[k];
- }
- } // end for(k)
+ // Calculate transition probabilities from D state
+ q.tr[i][D2M]=q.tr[i][D2D]=0.0;
+ for (k=0; k<N_in; k++) //for all sequences
+ {
+ if (in[k] && X[k][i]==GAP) //current state is D
+ {
+ if (X[k][i+1]==GAP) //next state is D
+ q.tr[i][D2D]+=wi[k];
+ else if (X[k][i+1]<=ANY) //next state is M
+ q.tr[i][D2M]+=wi[k];
}
+ } // end for(k)
+ }
- else // fast global weights?
+ else // fast global weights?
+ {
+ float w_D=-1.0/N_filtered;
+ ncol=0;
+ q.tr[i][D2M]=q.tr[i][D2D]=0.0;
+ // Calculate amino acid frequencies fj[a] from weights wg[k]
+ for (k=0; k<N_in; k++) //for all sequences
+ if (in[k] && X[k][i]==GAP) //current state is D
{
- float w_D=-1.0/N_filtered;
- ncol=0;
- q.tr[i][D2M]=q.tr[i][D2D]=0.0;
- // Calculate amino acid frequencies fj[a] from weights wg[k]
- for (k=0; k<N_in; k++) //for all sequences
- if (in[k] && X[k][i]==GAP) //current state is D
- {
- ncol++;
- w_D+=wg[k];
- if (X[k][i+1]==GAP) //next state is D
- q.tr[i][D2D]+=wi[k];
- else if (X[k][i+1]<=ANY) //next state is M
- q.tr[i][D2M]+=wi[k];
- }
- if (ncol>0)
- {
- if (w_D<0) Neff[i]=1.0;
- else Neff[i] = Nlim - (Nlim-1.0)*fpow2(scale*w_D);
- // fprintf(stderr,"D i=%3i ncol=%3i Neff_M=%5.2f Nlim=%5.2f w_D=%5.3f Neff_D=%5.2f\n",i,ncol,q.Neff_HMM,Nlim,w_D,Neff[i]);
- }
- else
- {
- Neff[i]=0.0; // effective number of sequences = 0!
- q.tr[i][D2M]=-100000;
- q.tr[i][D2D]=-100000;
- continue;
- }
+ ncol++;
+ w_D+=wg[k];
+ if (X[k][i+1]==GAP) //next state is D
+ q.tr[i][D2D]+=wi[k];
+ else if (X[k][i+1]<=ANY) //next state is M
+ q.tr[i][D2M]+=wi[k];
}
+ if (ncol>0)
+ {
+ if (w_D<0) Neff[i]=1.0;
+ else Neff[i] = Nlim - (Nlim-1.0)*fpow2(scale*w_D);
+ // fprintf(stderr,"D i=%3i ncol=%3i Neff_M=%5.2f Nlim=%5.2f w_D=%5.3f Neff_D=%5.2f\n",i,ncol,q.Neff_HMM,Nlim,w_D,Neff[i]);
+ }
+ else
+ {
+ Neff[i]=0.0; // effective number of sequences = 0!
+ q.tr[i][D2M]=-100000;
+ q.tr[i][D2D]=-100000;
+ continue;
+ }
+ }
- // Normalize and take log
- sum = q.tr[i][D2M]+q.tr[i][D2D];
- q.tr[i][D2M]=log2(q.tr[i][D2M]/sum);
- q.tr[i][D2D]=log2(q.tr[i][D2D]/sum);
+ // Normalize and take log
+ sum = q.tr[i][D2M]+q.tr[i][D2D];
+ q.tr[i][D2M]=Log2(q.tr[i][D2M]/sum);
+ q.tr[i][D2D]=Log2(q.tr[i][D2D]/sum);
- }
+ }
// end loop through alignment columns i
//////////////////////////////////////////////////////////////////////////////////////////////
q.Neff_D[0]=99.999;
// Assign Neff_D[i]
- for (i=1; i<=L; i++)
+ for (i=1; i<=L; i++) {
q.Neff_D[i]=Neff[i];
-
- delete[](wi); wi = NULL;/* FIXME: FS */
- // delete n[][]
- for (j=1; j<=L; j++){
- delete[](n[j]); (n[j]) = NULL;
+ delete[](n[i]);
+ (n[i]) = NULL;
}
+
delete[](n); (n) = NULL;
+ delete[](wi); wi = NULL;/* FIXME: FS */
delete[] Neff; Neff = NULL;
return;
N_filtered = Tali.Filter(par.max_seqid,par.coverage,par.qid,par.qsc,par.Ndiff);
// Record imatch[j]
- int* imatch=new(int[hit.j2+1]);
+ int* imatch=new int[hit.j2+1];
int step = hit.nsteps;
for (j=hit.j1; j<=hit.j2; j++)
{
iprev=i; lprev=l;
if (h>=maxcol-1000) // too few columns? Reserve double space
{
- char* new_seq=new(char[2*maxcol]);
+ char* new_seq=new char[2*maxcol];
strncpy(new_seq,cur_seq,maxcol); //////// check: maxcol-1 ????
delete[](cur_seq); (cur_seq) = NULL;
cur_seq=new_seq;
cur_seq[h++]='\0';
keep[N_in] = display[N_in] = KEEP_CONDITIONALLY;
- seq[N_in]=new(char[h]);
+ seq[N_in]=new char[h];
if (!seq[N_in]) MemoryError("array for input sequences");
strcpy(seq[N_in],cur_seq);
- X[N_in]=new(char[h]);
+ X[N_in]=new char[h];
if (!X[N_in]) MemoryError("array for input sequences");
- I[N_in]=new(short unsigned int[h]);
+ I[N_in]=new short unsigned int[h];
if (!I[N_in]) MemoryError("array for input sequences");
- sname[N_in]=new(char[strlen(Tali.sname[k])+1]);
+ sname[N_in]=new char[strlen(Tali.sname[k])+1];
if (!sname[N_in]) MemoryError("array for input sequences");
strcpy(sname[N_in],Tali.sname[k]);
N_in++;
{
int i; // position in query and target
if (L<=0) InternalError("L is not set in AddSequence()");
- X[N_in]=new(char[L+2]);
+ X[N_in]=new char[L+2];
for (i=0; i<=L+1; i++) X[N_in][i]=Xk[i];
if (Ik==NULL)
for (i=0; i<=L+1; i++) I[N_in][i]=0;
{
// Initialization
- n = new(int*[L+2]);
+ n = new int*[L+2];
for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);
for (j=1; j<=L; j++)
for (a=0; a<NAA+3; a++) n[j][a]=0;
/* @<allocate memory for sequences etc@> */
for (k = 0; k < iCnt; k++){
#define GOOD_MEASURE 1000 /* Temporary -- can be removed once rest in place */
- I[k] = new(short unsigned int[iLen+2+GOOD_MEASURE]);
- X[k] = new(char[iLen+2+GOOD_MEASURE]);
- seq[k] = new(char[iLen+2+GOOD_MEASURE]);
+ I[k] = new short unsigned int[iLen+2+GOOD_MEASURE];
+ X[k] = new char[iLen+2+GOOD_MEASURE];
+ seq[k] = new char[iLen+2+GOOD_MEASURE];
seq[k][0] = ' ';
seq[k][1] = '\0';
if (NULL == ppcProf[k]){
git://source.jalview.org / jabaws.git / blobdiff