Commit 01c9ee8d authored by Andreas Tille's avatar Andreas Tille
Browse files

New upstream version 1.5+dfsg

parent 56aba896

.gitignore

+7 −0
Original line number Diff line number Diff line 
*.a

*.o

*.help



bin/

doc/man/

include/external_libs.h

+1 −1
Original line number Diff line number Diff line
@@ -33,7 +33,7 @@ 
#else



#ifdef VECLIB

#include <vecLib/clapack.h>

#include <Accelerate/Accelerate.h>

#define LAPACK_INT __CLPK_integer

#define LAPACK_DOUBLE __CLPK_doublereal

#else

include/misc.h

+11 −0
Original line number Diff line number Diff line
@@ -87,6 +87,17 @@ int int_pow(int x, int y) { 
  return retval;

}



/** Take x,y & max as int and cast x & y as double

 *  * raise x to the power y. return min(x^y,max) as int*/



static PHAST_INLINE

int pow_bounded(int x, int y, int max) {

  double z = pow((double) x, (double) y);

  return ( z > (double) max ? max : (int)z);

}







/** \name Log based calculation functions

\{ */

src/lib/base/matrix.c

+22 −1
Original line number Diff line number Diff line
@@ -293,7 +293,8 @@ int mat_invert(Matrix *M_inv, Matrix *M) { 
   same dimension, and C is diagonal.  C is described by a vector

   representing its diagonal elements.  */

void mat_mult_diag(Matrix *A, Matrix *B, Vector *C, Matrix *D) {

  int i, j, k;



/*  int i, j, k;

  for (i = 0; i < C->size; i++) {

    for (j = 0; j < C->size; j++) {

      A->data[i][j] = 0;
@@ -302,7 +303,27 @@ void mat_mult_diag(Matrix *A, Matrix *B, Vector *C, Matrix *D) { 
    }

  }

}

*/



A->data[0][0] = B->data[0][0] * C->data[0] * D->data[0][0] + B->data[0][1] * C->data[1] * D->data[1][0] + B->data[0][2] * C->data[2] * D->data[2][0] + B->data[0][3] * C->data[3] * D->data[3][0];

A->data[0][1] = B->data[0][0] * C->data[0] * D->data[0][1] + B->data[0][1] * C->data[1] * D->data[1][1] + B->data[0][2] * C->data[2] * D->data[2][1] + B->data[0][3] * C->data[3] * D->data[3][1];

A->data[0][2] = B->data[0][0] * C->data[0] * D->data[0][2] + B->data[0][1] * C->data[1] * D->data[1][2] + B->data[0][2] * C->data[2] * D->data[2][2] + B->data[0][3] * C->data[3] * D->data[3][2];

A->data[0][3] = B->data[0][0] * C->data[0] * D->data[0][3] + B->data[0][1] * C->data[1] * D->data[1][3] + B->data[0][2] * C->data[2] * D->data[2][3] + B->data[0][3] * C->data[3] * D->data[3][3];

A->data[1][0] = B->data[1][0] * C->data[0] * D->data[0][0] + B->data[1][1] * C->data[1] * D->data[1][0] + B->data[1][2] * C->data[2] * D->data[2][0] + B->data[1][3] * C->data[3] * D->data[3][0];

A->data[1][1] = B->data[1][0] * C->data[0] * D->data[0][1] + B->data[1][1] * C->data[1] * D->data[1][1] + B->data[1][2] * C->data[2] * D->data[2][1] + B->data[1][3] * C->data[3] * D->data[3][1];

A->data[1][2] = B->data[1][0] * C->data[0] * D->data[0][2] + B->data[1][1] * C->data[1] * D->data[1][2] + B->data[1][2] * C->data[2] * D->data[2][2] + B->data[1][3] * C->data[3] * D->data[3][2];

A->data[1][3] = B->data[1][0] * C->data[0] * D->data[0][3] + B->data[1][1] * C->data[1] * D->data[1][3] + B->data[1][2] * C->data[2] * D->data[2][3] + B->data[1][3] * C->data[3] * D->data[3][3];

A->data[2][0] = B->data[2][0] * C->data[0] * D->data[0][0] + B->data[2][1] * C->data[1] * D->data[1][0] + B->data[2][2] * C->data[2] * D->data[2][0] + B->data[2][3] * C->data[3] * D->data[3][0];

A->data[2][1] = B->data[2][0] * C->data[0] * D->data[0][1] + B->data[2][1] * C->data[1] * D->data[1][1] + B->data[2][2] * C->data[2] * D->data[2][1] + B->data[2][3] * C->data[3] * D->data[3][1];

A->data[2][2] = B->data[2][0] * C->data[0] * D->data[0][2] + B->data[2][1] * C->data[1] * D->data[1][2] + B->data[2][2] * C->data[2] * D->data[2][2] + B->data[2][3] * C->data[3] * D->data[3][2];

A->data[2][3] = B->data[2][0] * C->data[0] * D->data[0][3] + B->data[2][1] * C->data[1] * D->data[1][3] + B->data[2][2] * C->data[2] * D->data[2][3] + B->data[2][3] * C->data[3] * D->data[3][3];

A->data[3][0] = B->data[3][0] * C->data[0] * D->data[0][0] + B->data[3][1] * C->data[1] * D->data[1][0] + B->data[3][2] * C->data[2] * D->data[2][0] + B->data[3][3] * C->data[3] * D->data[3][0];

A->data[3][1] = B->data[3][0] * C->data[0] * D->data[0][1] + B->data[3][1] * C->data[1] * D->data[1][1] + B->data[3][2] * C->data[2] * D->data[2][1] + B->data[3][3] * C->data[3] * D->data[3][1];

A->data[3][2] = B->data[3][0] * C->data[0] * D->data[0][2] + B->data[3][1] * C->data[1] * D->data[1][2] + B->data[3][2] * C->data[2] * D->data[2][2] + B->data[3][3] * C->data[3] * D->data[3][2];

A->data[3][3] = B->data[3][0] * C->data[0] * D->data[0][3] + B->data[3][1] * C->data[1] * D->data[1][3] + B->data[3][2] * C->data[2] * D->data[2][3] + B->data[3][3] * C->data[3] * D->data[3][3];





}

int mat_equal(Matrix *A, Matrix *B) {

  int i, j;

  if (A->nrows != B->nrows || A->ncols != B->ncols) return 0;

src/lib/msa/maf.c

+4 −8
Original line number Diff line number Diff line
@@ -139,12 +139,10 @@ MSA *maf_read_cats_subset(FILE *F, FILE *REFSEQF, int tuple_size, 
      msa->alloc_len = msa->length = refseqlen;  //this may still not be big enough because of gaps in refseq

    else 

      msa->alloc_len = 50000;

    max_tuples =  max(1000000, (int)pow(strlen(msa->alphabet)+strlen(msa->missing)+1, 

				   2 * msa->nseqs * tuple_size));

    max_tuples = pow_bounded(strlen(msa->alphabet)+strlen(msa->missing)+1, 2 * msa->nseqs * tuple_size, 1000000);

  }

  else 

    max_tuples = min(1000000,

		     (int)pow(strlen(msa->alphabet)+strlen(msa->missing)+1, 2 * msa->nseqs * tuple_size));

    max_tuples = pow_bounded(strlen(msa->alphabet)+strlen(msa->missing)+1, 2 * msa->nseqs * tuple_size, 1000000);



  if (max_tuples > 10000000 || max_tuples < 0) max_tuples = 10000000;

  if (max_tuples < 1000000) max_tuples = 1000000;
@@ -596,15 +594,13 @@ MSA *maf_read_unsorted(FILE *F, FILE *REFSEQF, int tuple_size, char *alphabet, 


  if (store_order) {

    msa->length = map != NULL ? map->msa_len : refseqlen;

    max_tuples = min(msa->length,

                     (int)pow(strlen(msa->alphabet)+strlen(msa->missing)+1, msa->nseqs * tuple_size));

    max_tuples = pow_bounded(strlen(msa->alphabet)+strlen(msa->missing)+1, msa->nseqs * tuple_size, msa->length);

    if (max_tuples < 0) max_tuples = msa->length;

    if (max_tuples > 1000000) max_tuples = 1000000; 

  }

  else {

    msa->length = 0;

    max_tuples = min(50000,

                     (int)pow(strlen(msa->alphabet)+strlen(msa->missing)+1, msa->nseqs * tuple_size));

    max_tuples = pow_bounded(strlen(msa->alphabet)+strlen(msa->missing)+1, msa->nseqs * tuple_size, 50000);

    if (max_tuples < 0) max_tuples = 50000;

  }