sre_math.c

hmmer源程序

DAdd(double *vec1, double *vec2, int n)
{
  int x;
  for (x = 0; x < n; x++)
    vec1[x] += vec2[x];
}
void
FAdd(float *vec1, float *vec2, int n)
{
  int x;
  for (x = 0; x < n; x++)
    vec1[x] += vec2[x];
}

void
DCopy(double *vec1, double *vec2, int n)
{
  int x;
  for (x = 0; x < n; x++)
    vec1[x] = vec2[x];
}
void
FCopy(float *vec1, float *vec2, int n)
{
  int x;
  for (x = 0; x < n; x++)
    vec1[x] = vec2[x];
}

double
DDot(double *vec1, double *vec2, int n)
{
  double result = 0.;
  int x;

  for (x = 0; x < n; x++)
    result += vec1[x] * vec2[x];
  return result;
}
float
FDot(float *vec1, float *vec2, int n)
{
  float result = 0.;
  int x;

  for (x = 0; x < n; x++)
    result += vec1[x] * vec2[x];
  return result;
}

/* Functions: DMax(), FMax()
 * Date:      SRE, Fri Aug 29 11:14:08 1997 (Denver CO)
 * 
 * Purpose:   return index of maximum element in vec.
 */
int
DMax(double *vec, int n)
{
  int i;
  int best = 0;

  for (i = 1; i < n; i++)
    if (vec[i] > vec[best]) best = i;
  return best;
}
int
FMax(float *vec, int n)
{
  int i;
  int best = 0;

  for (i = 1; i < n; i++)
    if (vec[i] > vec[best]) best = i;
  return best;
}


/* 2D matrix operations
 */
float **
FMX2Alloc(int rows, int cols)
{
  float **mx;
  int     r;
  
  mx    = (float **) MallocOrDie(sizeof(float *) * rows);
  mx[0] = (float *)  MallocOrDie(sizeof(float) * rows * cols);
  for (r = 1; r < rows; r++)
    mx[r] = mx[0] + r*cols;
  return mx;
}
void
FMX2Free(float **mx)
{
  free(mx[0]);
  free(mx);
}
double **
DMX2Alloc(int rows, int cols)
{
  double **mx;
  int      r;
  
  mx    = (double **) MallocOrDie(sizeof(double *) * rows);
  mx[0] = (double *)  MallocOrDie(sizeof(double) * rows * cols);
  for (r = 1; r < rows; r++)
    mx[r] = mx[0] + r*cols;
  return mx;
}
void
DMX2Free(double **mx)
{
  free(mx[0]);
  free(mx);
}
/* Function: FMX2Multiply()
 * 
 * Purpose:  Matrix multiplication.
 *           Multiply an m x p matrix A by a p x n matrix B,
 *           giving an m x n matrix C.
 *           Matrix C must be a preallocated matrix of the right
 *           size.
 */
void
FMX2Multiply(float **A, float **B, float **C, int m, int p, int n)
{
  int i, j, k;

  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++)
      {
	C[i][j] = 0.;
	for (k = 0; k < p; k++)
	  C[i][j] += A[i][p] * B[p][j];
      }
}

/* Function: sre_random()
 * 
 * Purpose:  Return a uniform deviate from 0.0 to 1.0.
 *           sre_randseed is a static variable, set
 *           by sre_srandom(). sre_reseed is a static flag
 *           raised by sre_srandom(), saying that we need
 *           to re-initialize.
 *           [0.0 <= x < 1.0]
 *           
 *           Uses a simple linear congruential generator with
 *           period 2^28. Based on discussion in Robert Sedgewick's 
 *           _Algorithms in C_, Addison-Wesley, 1990.
 *
 *           Requires that long int's have at least 32 bits.
 *
 * Reliable and portable, but slow. Benchmarks on wol,
 * using IRIX cc and IRIX C library rand() and random():
 *     sre_random():    0.8 usec/call
 *     random():        0.3 usec/call
 *     rand():          0.3 usec/call
 */
#define RANGE 268435456		/* 2^28        */
#define DIV   16384		/* sqrt(RANGE) */
#define MULT  72530821		/* my/Cathy's birthdays, x21, x even (Knuth)*/
float 
sre_random(void)
{
  static long  rnd;
  static int   firsttime = 1;
  long         high1, low1;
  long         high2, low2; 

  if (sre_reseed || firsttime) 
    {
      sre_reseed = firsttime = 0;
      if (sre_randseed <= 0) sre_randseed = 666; /* seeds of zero break me */
      high1 = sre_randseed / DIV;  low1  = sre_randseed % DIV;
      high2 = MULT / DIV;          low2  = MULT % DIV;
      rnd = (((high2*low1 + high1*low2) % DIV)*DIV + low1*low2) % RANGE;
    }
   high1 = rnd / DIV;  low1  = rnd % DIV;
   high2 = MULT / DIV; low2  = MULT % DIV;
   rnd = (((high2*low1 + high1*low2) % DIV)*DIV + low1*low2) % RANGE;

   return ((float) rnd / (float) RANGE);  
}
#undef RANGE
#undef DIV
#undef MULT


/* Function: sre_srandom()
 * 
 * Purpose:  Initialize with a random seed. Seed can be
 *           any integer.
 */
void
sre_srandom(int seed)
{
  if (seed < 0) seed = -1 * seed;
  sre_reseed   = 1;
  sre_randseed = seed;
}


/* Functions: DChoose(), FChoose()
 *
 * Purpose:   Make a random choice from a normalized distribution.
 *            DChoose() is for double-precision vectors;
 *            FChoose() is for single-precision float vectors.
 *            Returns the number of the choice.
 */
int
DChoose(double *p, int N)
{
  double roll;                  /* random fraction */
  double sum;                   /* integrated prob */
  int    i;                     /* counter over the probs */

  roll    = sre_random();
  sum     = 0.0;
  for (i = 0; i < N; i++)
    {
      sum += p[i];
      if (roll < sum) return i;
    }
  SQD_DASSERT2((fabs(1.0 - sum) < 1e-14)); /* a verification at level 2 */
  return (int) (sre_random() * N);         /* bulletproof */
}
int
FChoose(float *p, int N)
{
  float roll;                   /* random fraction */
  float sum;			/* integrated prob */
  int   i;                      /* counter over the probs */

  roll    = sre_random();
  sum     = 0.0;
  for (i = 0; i < N; i++)
    {
      sum += p[i];
      if (roll < sum) return i;
    }
  SQD_DASSERT2((fabs(1.0f - sum) < 1e-6f));  /* a verification at level 2 */
  return (int) (sre_random() * N);           /* bulletproof */
}

/* Functions: DLogSum(), FLogSum()
 * 
 * Calculate the sum of a log vector
 * *in normal space*, and return the log of the sum.
 */
double 
DLogSum(double *logp, int n)
{
  int x;
  double max, sum;
  
  max = logp[0];
  for (x = 1; x < n; x++)
    if (logp[x] > max) max = logp[x];
  sum = 0.0;
  for (x = 0; x < n; x++)
    if (logp[x] > max - 50.)
      sum += exp(logp[x] - max);
  sum = log(sum) + max;
  return sum;
}
float
FLogSum(float *logp, int n)
{
  int x;
  float max, sum;
  
  max = logp[0];
  for (x = 1; x < n; x++)
    if (logp[x] > max) max = logp[x];
  sum = 0.0;
  for (x = 0; x < n; x++)
    if (logp[x] > max - 50.)
      sum += exp(logp[x] - max);
  sum = log(sum) + max;
  return sum;
}


/* Function: IncompleteGamma()
 * 
 * Purpose:  Returns 1 - P(a,x) where:
 *           P(a,x) = \frac{1}{\Gamma(a)} \int_{0}^{x} t^{a-1} e^{-t} dt
 *                  = \frac{\gamma(a,x)}{\Gamma(a)}
 *                  = 1 - \frac{\Gamma(a,x)}{\Gamma(a)}
 *                  
 *           Used in a chi-squared test: for a X^2 statistic x
 *           with v degrees of freedom, call:
 *                  p = IncompleteGamma(v/2., x/2.) 
 *           to get the probability p that a chi-squared value
 *           greater than x could be obtained by chance even for
 *           a correct model. (i.e. p should be large, say 
 *           0.95 or more).
 *           
 * Method:   Based on ideas from Numerical Recipes in C, Press et al.,
 *           Cambridge University Press, 1988. 
 *           
 * Args:     a  - for instance, degrees of freedom / 2     [a > 0]
 *           x  - for instance, chi-squared statistic / 2  [x >= 0] 
 *           
 * Return:   1 - P(a,x).
 */          
double
IncompleteGamma(double a, double x)
{
  int iter;			/* iteration counter */

  if (a <= 0.) Die("IncompleteGamma(): a must be > 0");
  if (x <  0.) Die("IncompleteGamma(): x must be >= 0");

  /* For x > a + 1 the following gives rapid convergence;
   * calculate 1 - P(a,x) = \frac{\Gamma(a,x)}{\Gamma(a)}:
   *     use a continued fraction development for \Gamma(a,x).
   */
  if (x > a+1) 
    {
      double oldp;		/* previous value of p    */
      double nu0, nu1;		/* numerators for continued fraction calc   */
      double de0, de1;		/* denominators for continued fraction calc */

      nu0 = 0.;			/* A_0 = 0       */
      de0 = 1.;			/* B_0 = 1       */
      nu1 = 1.;			/* A_1 = 1       */
      de1 = x;			/* B_1 = x       */

      oldp = nu1;
      for (iter = 1; iter < 100; iter++)
	{
	  /* Continued fraction development:
	   * set A_j = b_j A_j-1 + a_j A_j-2
	   *     B_j = b_j B_j-1 + a_j B_j-2
           * We start with A_2, B_2.
	   */
				/* j = even: a_j = iter-a, b_j = 1 */
				/* A,B_j-2 are in nu0, de0; A,B_j-1 are in nu1,de1 */
	  nu0 = nu1 + ((double)iter - a) * nu0;
	  de0 = de1 + ((double)iter - a) * de0;

				/* j = odd: a_j = iter, b_j = x */
				/* A,B_j-2 are in nu1, de1; A,B_j-1 in nu0,de0 */
	  nu1 = x * nu0 + (double) iter * nu1;
	  de1 = x * de0 + (double) iter * de1;

				/* rescale */
	  if (de1) 
	    { 
	      nu0 /= de1; 
	      de0 /= de1;
	      nu1 /= de1;
	      de1 =  1.;
	    }
				/* check for convergence */
	  if (fabs((nu1-oldp)/nu1) < 1.e-7)
	    return nu1 * exp(a * log(x) - x - Gammln(a));

	  oldp = nu1;
	}
      Die("IncompleteGamma(): failed to converge using continued fraction approx");
    }
  else /* x <= a+1 */
    {
      double p;			/* current sum               */
      double val;		/* current value used in sum */

      /* For x <= a+1 we use a convergent series instead:
       *   P(a,x) = \frac{\gamma(a,x)}{\Gamma(a)},
       * where
       *   \gamma(a,x) = e^{-x}x^a \sum_{n=0}{\infty} \frac{\Gamma{a}}{\Gamma{a+1+n}} x^n
       * which looks appalling but the sum is in fact rearrangeable to
       * a simple series without the \Gamma functions:
       *   = \frac{1}{a} + \frac{x}{a(a+1)} + \frac{x^2}{a(a+1)(a+2)} ...
       * and it's obvious that this should converge nicely for x <= a+1.
       */
      
      p = val = 1. / a;
      for (iter = 1; iter < 10000; iter++)
	{
	  val *= x / (a+(double)iter);
	  p   += val;
	  
	  if (fabs(val/p) < 1.e-7)
	    return 1. - p * exp(a * log(x) - x - Gammln(a));
	}
      Die("IncompleteGamma(): failed to converge using series approx");
    }
  /*NOTREACHED*/
  return 0.;
}