Mercurial > repos > clustalomega > clustalomega
annotate clustalomega/clustal-omega-0.2.0/src/squid/sre_math.c @ 0:ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
author | clustalomega |
---|---|
date | Tue, 07 Jun 2011 17:04:25 -0400 |
parents | |
children |
rev | line source |
---|---|
0
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
1 /***************************************************************** |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
2 * SQUID - a library of functions for biological sequence analysis |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
3 * Copyright (C) 1992-2002 Washington University School of Medicine |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
4 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
5 * This source code is freely distributed under the terms of the |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
6 * GNU General Public License. See the files COPYRIGHT and LICENSE |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
7 * for details. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
8 *****************************************************************/ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
9 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
10 /* sre_math.c |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
11 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
12 * Portability for and extensions to C math library. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
13 * RCS $Id: sre_math.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: sre_math.c,v 1.12 2002/10/09 14:26:09 eddy Exp) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
14 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
15 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
16 #include <stdio.h> |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
17 #include <stdlib.h> |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
18 #include <math.h> |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
19 #include "squid.h" |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
20 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
21 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
22 /* Function: Linefit() |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
23 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
24 * Purpose: Given points x[0..N-1] and y[0..N-1], fit to |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
25 * a straight line y = a + bx. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
26 * a, b, and the linear correlation coefficient r |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
27 * are filled in for return. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
28 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
29 * Args: x - x values of data |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
30 * y - y values of data |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
31 * N - number of data points |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
32 * ret_a - RETURN: intercept |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
33 * ret_b - RETURN: slope |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
34 * ret_r - RETURN: correlation coefficient |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
35 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
36 * Return: 1 on success, 0 on failure. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
37 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
38 int |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
39 Linefit(float *x, float *y, int N, float *ret_a, float *ret_b, float *ret_r) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
40 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
41 float xavg, yavg; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
42 float sxx, syy, sxy; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
43 int i; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
44 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
45 /* Calculate averages, xavg and yavg |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
46 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
47 xavg = yavg = 0.0; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
48 for (i = 0; i < N; i++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
49 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
50 xavg += x[i]; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
51 yavg += y[i]; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
52 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
53 xavg /= (float) N; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
54 yavg /= (float) N; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
55 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
56 sxx = syy = sxy = 0.0; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
57 for (i = 0; i < N; i++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
58 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
59 sxx += (x[i] - xavg) * (x[i] - xavg); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
60 syy += (y[i] - yavg) * (y[i] - xavg); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
61 sxy += (x[i] - xavg) * (y[i] - yavg); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
62 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
63 *ret_b = sxy / sxx; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
64 *ret_a = yavg - xavg*(*ret_b); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
65 *ret_r = sxy / (sqrt(sxx) * sqrt(syy)); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
66 return 1; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
67 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
68 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
69 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
70 /* Function: WeightedLinefit() |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
71 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
72 * Purpose: Given points x[0..N-1] and y[0..N-1] with |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
73 * variances (measurement errors) var[0..N-1], |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
74 * fit to a straight line y = mx + b. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
75 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
76 * Method: Algorithm from Numerical Recipes in C, [Press88]. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
77 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
78 * Return: (void) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
79 * ret_m contains slope; ret_b contains intercept |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
80 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
81 void |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
82 WeightedLinefit(float *x, float *y, float *var, int N, float *ret_m, float *ret_b) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
83 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
84 int i; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
85 double s; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
86 double sx, sy; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
87 double sxx, sxy; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
88 double delta; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
89 double m, b; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
90 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
91 s = sx = sy = sxx = sxy = 0.; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
92 for (i = 0; i < N; i++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
93 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
94 s += 1./var[i]; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
95 sx += x[i] / var[i]; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
96 sy += y[i] / var[i]; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
97 sxx += x[i] * x[i] / var[i]; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
98 sxy += x[i] * y[i] / var[i]; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
99 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
100 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
101 delta = s * sxx - (sx * sx); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
102 b = (sxx * sy - sx * sxy) / delta; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
103 m = (s * sxy - sx * sy) / delta; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
104 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
105 *ret_m = m; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
106 *ret_b = b; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
107 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
108 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
109 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
110 /* Function: Gammln() |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
111 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
112 * Returns the natural log of the gamma function of x. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
113 * x is > 0.0. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
114 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
115 * Adapted from a public domain implementation in the |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
116 * NCBI core math library. Thanks to John Spouge and |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
117 * the NCBI. (According to the NCBI, that's Dr. John |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
118 * "Gammas Galore" Spouge to you, pal.) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
119 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
120 double |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
121 Gammln(double x) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
122 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
123 int i; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
124 double xx, tx; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
125 double tmp, value; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
126 static double cof[11] = { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
127 4.694580336184385e+04, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
128 -1.560605207784446e+05, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
129 2.065049568014106e+05, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
130 -1.388934775095388e+05, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
131 5.031796415085709e+04, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
132 -9.601592329182778e+03, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
133 8.785855930895250e+02, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
134 -3.155153906098611e+01, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
135 2.908143421162229e-01, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
136 -2.319827630494973e-04, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
137 1.251639670050933e-10 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
138 }; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
139 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
140 /* Protect against x=0. We see this in Dirichlet code, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
141 * for terms alpha = 0. This is a severe hack but it is effective |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
142 * and (we think?) safe. (due to GJM) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
143 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
144 if (x <= 0.0) return 999999.; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
145 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
146 xx = x - 1.0; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
147 tx = tmp = xx + 11.0; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
148 value = 1.0; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
149 for (i = 10; i >= 0; i--) /* sum least significant terms first */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
150 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
151 value += cof[i] / tmp; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
152 tmp -= 1.0; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
153 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
154 value = log(value); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
155 tx += 0.5; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
156 value += 0.918938533 + (xx+0.5)*log(tx) - tx; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
157 return value; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
158 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
159 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
160 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
161 /* 2D matrix operations |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
162 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
163 float ** |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
164 FMX2Alloc(int rows, int cols) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
165 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
166 float **mx; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
167 int r; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
168 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
169 mx = (float **) MallocOrDie(sizeof(float *) * rows); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
170 mx[0] = (float *) MallocOrDie(sizeof(float) * rows * cols); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
171 for (r = 1; r < rows; r++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
172 mx[r] = mx[0] + r*cols; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
173 return mx; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
174 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
175 void |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
176 FMX2Free(float **mx) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
177 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
178 free(mx[0]); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
179 free(mx); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
180 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
181 double ** |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
182 DMX2Alloc(int rows, int cols) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
183 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
184 double **mx; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
185 int r; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
186 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
187 mx = (double **) MallocOrDie(sizeof(double *) * rows); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
188 mx[0] = (double *) MallocOrDie(sizeof(double) * rows * cols); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
189 for (r = 1; r < rows; r++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
190 mx[r] = mx[0] + r*cols; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
191 return mx; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
192 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
193 void |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
194 DMX2Free(double **mx) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
195 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
196 free(mx[0]); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
197 free(mx); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
198 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
199 /* Function: FMX2Multiply() |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
200 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
201 * Purpose: Matrix multiplication. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
202 * Multiply an m x p matrix A by a p x n matrix B, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
203 * giving an m x n matrix C. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
204 * Matrix C must be a preallocated matrix of the right |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
205 * size. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
206 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
207 void |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
208 FMX2Multiply(float **A, float **B, float **C, int m, int p, int n) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
209 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
210 int i, j, k; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
211 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
212 for (i = 0; i < m; i++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
213 for (j = 0; j < n; j++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
214 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
215 C[i][j] = 0.; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
216 for (k = 0; k < p; k++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
217 C[i][j] += A[i][p] * B[p][j]; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
218 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
219 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
220 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
221 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
222 /* Function: IncompleteGamma() |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
223 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
224 * Purpose: Returns 1 - P(a,x) where: |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
225 * P(a,x) = \frac{1}{\Gamma(a)} \int_{0}^{x} t^{a-1} e^{-t} dt |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
226 * = \frac{\gamma(a,x)}{\Gamma(a)} |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
227 * = 1 - \frac{\Gamma(a,x)}{\Gamma(a)} |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
228 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
229 * Used in a chi-squared test: for a X^2 statistic x |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
230 * with v degrees of freedom, call: |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
231 * p = IncompleteGamma(v/2., x/2.) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
232 * to get the probability p that a chi-squared value |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
233 * greater than x could be obtained by chance even for |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
234 * a correct model. (i.e. p should be large, say |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
235 * 0.95 or more). |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
236 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
237 * Method: Based on ideas from Numerical Recipes in C, Press et al., |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
238 * Cambridge University Press, 1988. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
239 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
240 * Args: a - for instance, degrees of freedom / 2 [a > 0] |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
241 * x - for instance, chi-squared statistic / 2 [x >= 0] |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
242 * |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
243 * Return: 1 - P(a,x). |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
244 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
245 double |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
246 IncompleteGamma(double a, double x) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
247 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
248 int iter; /* iteration counter */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
249 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
250 if (a <= 0.) Die("IncompleteGamma(): a must be > 0"); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
251 if (x < 0.) Die("IncompleteGamma(): x must be >= 0"); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
252 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
253 /* For x > a + 1 the following gives rapid convergence; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
254 * calculate 1 - P(a,x) = \frac{\Gamma(a,x)}{\Gamma(a)}: |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
255 * use a continued fraction development for \Gamma(a,x). |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
256 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
257 if (x > a+1) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
258 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
259 double oldp; /* previous value of p */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
260 double nu0, nu1; /* numerators for continued fraction calc */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
261 double de0, de1; /* denominators for continued fraction calc */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
262 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
263 nu0 = 0.; /* A_0 = 0 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
264 de0 = 1.; /* B_0 = 1 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
265 nu1 = 1.; /* A_1 = 1 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
266 de1 = x; /* B_1 = x */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
267 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
268 oldp = nu1; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
269 for (iter = 1; iter < 100; iter++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
270 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
271 /* Continued fraction development: |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
272 * set A_j = b_j A_j-1 + a_j A_j-2 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
273 * B_j = b_j B_j-1 + a_j B_j-2 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
274 * We start with A_2, B_2. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
275 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
276 /* j = even: a_j = iter-a, b_j = 1 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
277 /* A,B_j-2 are in nu0, de0; A,B_j-1 are in nu1,de1 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
278 nu0 = nu1 + ((double)iter - a) * nu0; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
279 de0 = de1 + ((double)iter - a) * de0; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
280 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
281 /* j = odd: a_j = iter, b_j = x */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
282 /* A,B_j-2 are in nu1, de1; A,B_j-1 in nu0,de0 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
283 nu1 = x * nu0 + (double) iter * nu1; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
284 de1 = x * de0 + (double) iter * de1; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
285 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
286 /* rescale */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
287 if (de1 != 0.) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
288 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
289 nu0 /= de1; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
290 de0 /= de1; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
291 nu1 /= de1; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
292 de1 = 1.; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
293 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
294 /* check for convergence */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
295 if (fabs((nu1-oldp)/nu1) < 1.e-7) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
296 return nu1 * exp(a * log(x) - x - Gammln(a)); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
297 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
298 oldp = nu1; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
299 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
300 Die("IncompleteGamma(): failed to converge using continued fraction approx"); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
301 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
302 else /* x <= a+1 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
303 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
304 double p; /* current sum */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
305 double val; /* current value used in sum */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
306 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
307 /* For x <= a+1 we use a convergent series instead: |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
308 * P(a,x) = \frac{\gamma(a,x)}{\Gamma(a)}, |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
309 * where |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
310 * \gamma(a,x) = e^{-x}x^a \sum_{n=0}{\infty} \frac{\Gamma{a}}{\Gamma{a+1+n}} x^n |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
311 * which looks appalling but the sum is in fact rearrangeable to |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
312 * a simple series without the \Gamma functions: |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
313 * = \frac{1}{a} + \frac{x}{a(a+1)} + \frac{x^2}{a(a+1)(a+2)} ... |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
314 * and it's obvious that this should converge nicely for x <= a+1. |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
315 */ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
316 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
317 p = val = 1. / a; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
318 for (iter = 1; iter < 10000; iter++) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
319 { |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
320 val *= x / (a+(double)iter); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
321 p += val; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
322 |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
323 if (fabs(val/p) < 1.e-7) |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
324 return 1. - p * exp(a * log(x) - x - Gammln(a)); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
325 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
326 Die("IncompleteGamma(): failed to converge using series approx"); |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
327 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
328 /*NOTREACHED*/ |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
329 return 0.; |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
330 } |
ff1768533a07
Migrated tool version 0.2 from old tool shed archive to new tool shed repository
clustalomega
parents:
diff
changeset
|
331 |