source: tags/arb-6.0/GDE/PHYLIP/hermite.c

Last change on this file was 5458, checked in by baderk, 17 years ago

Removed .cvsignore files from repository. Hopefully this time all svn:ignore flags were set right.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 13.9 KB
Line 
1From beerli@genetics.washington.edu Tue Mar 13 15:47:38 2001
2Received: from darwin.genetics.washington.edu (darwin.genetics.washington.edu [128.95.144.42])
3        by evolution.genetics.washington.edu (8.9.3/8.9.3) with ESMTP id PAA17484
4        for <joe@evolution.genetics.washington.edu>; Tue, 13 Mar 2001 15:47:38 -0800 (PST)
5Received: from localhost (beerli@localhost)
6        by darwin.genetics.washington.edu (8.9.3/8.9.3) with ESMTP id PAA10268
7        for <joe@evolution.genetics.washington.edu>; Tue, 13 Mar 2001 15:47:38 -0800
8X-Authentication-Warning: darwin.genetics.washington.edu: beerli owned process doing -bs
9Date: Tue, 13 Mar 2001 15:47:38 -0800 (PST)
10From: Peter Beerli <beerli@genetics.washington.edu>
11To: Felsenstein Joe <joe@evolution.genetics.washington.edu>
12Subject: hermite a fini laguerre
13Message-ID: <Pine.LNX.4.30.0103131538570.8726-100000@darwin.genetics.washington.edu>
14MIME-Version: 1.0
15Content-Type: TEXT/PLAIN; charset=US-ASCII
16Status: RO
17
18I have two files included below,
19and you certainly need to adapt them for phylip.
20
21the hermit stuff contains the word hermite
22and its main function is
23
24inthermitcat()
25
26
27I ripped the initgammacat into a simple driver.
28laguerre.c contains a simple test program that
29can be compiled using
30gcc -DLAGUERRE_TEST laguerre.c -o laguerretest -lm
31
32Peter
33-- 
34
35/* laguerre.h */
36#ifndef _LAGUERRE_H_
37#define _LAGUERRE_H_
38/*------------------------------------------------------
39 Maximum likelihood estimation
40 of migration rate  and effectice population size
41 using a Metropolis-Hastings Monte Carlo algorithm
42 -------------------------------------------------------
43 Laguerre integration  R O U T I N E S
44
45 Peter Beerli 2000, Seattle
46 beerli@genetics.washington.edu
47 $Id: hermite.c 5458 2008-07-16 15:24:20Z westram $
48-------------------------------------------------------*/
49
50extern void integrate_laguerre(long categs, double *rate,
51                                 double *probcat,
52                                 double (*func)(double, helper_fmt *),
53                                 helper_fmt *helper, double *result,
54                                 double *rmax);
55extern void initgammacat (long categs, double alpha, double theta1,
56              double *rate, double *probcat);
57
58#endif
59
60
61/* laguerre.c */
62/*
63   (mutation) rates following a Gamma distribution
64   using orthogonal polynomials for finding rates and
65   LOG(probabilities)
66   [based on initgammacat of Joe Felsenstein]
67
68   - Generalized Laguerre (routine by Joe Felsenstein 2000)
69     defining points for a Gamma distribution with
70     shape parameter alpha and location parameter beta=1/alpha
71     [mean=1, std = 1/alpha^2]
72   - Hermite (approximates a normal and is activated when
73     the shape parameter alpha is > 100.)
74
75   Part of Migrate (Lamarc package)
76   http://evolution.genetics.washington.edu/lamarc.html
77
78   Peter Beerli, Seattle 2001
79   $Id: hermite.c 5458 2008-07-16 15:24:20Z westram $
80*/
81#include "migration.h"
82#include "laguerre.h"
83#include "tools.h"
84
85#define SQRTPI 1.7724538509055160273
86#define SQRT2  1.4142135623730950488
87
88
89/*this triggers the test  main()
90and is called with
91gcc -DLAGUERRE_TEST -g laguerre.c -o laguerre -lm*/
92
93#ifdef LAGUERRE_TEST
94/* at the end is a test main to help test if the root/weights finding
95   is OK*/
96
97/* if machine has lgamma() use it otherwise use lgamma from
98   tools.h*/
99#undef LGAMMA
100#define LGAMMA lgamma
101
102
103/* for migrate this is defined in tools.h */
104double
105logfac (long n)
106{
107  /* log(n!) values were calculated with Mathematica
108     with a precision of 30 digits */
109  switch (n)
110    {
111    case 0:
112      return 0.;
113    case 1:
114      return 0.;
115    case 2:
116      return 0.693147180559945309417232121458;
117    case 3:
118      return 1.791759469228055000812477358381;
119    case 4:
120      return 3.1780538303479456196469416013;
121    case 5:
122      return 4.78749174278204599424770093452;
123    case 6:
124      return 6.5792512120101009950601782929;
125    case 7:
126      return 8.52516136106541430016553103635;
127    case 8:
128      return 10.60460290274525022841722740072;
129    case 9:
130      return 12.80182748008146961120771787457;
131    case 10:
132      return 15.10441257307551529522570932925;
133    case 11:
134      return 17.50230784587388583928765290722;
135    case 12:
136      return 19.98721449566188614951736238706;
137    default:
138      return LGAMMA(n + 1.);
139    }
140}
141#endif
142
143/* prototypes */
144double hermite(long n, double x);
145void root_hermite(long n, double *hroot);
146double  halfroot(double (*func)(long m, double x),
147                 long n, double startx, double delta);
148void hermite_weight(long n, double * hroot, double * weights);
149void inithermitcat(long categs, double alpha, double theta1,
150                   double *rate, double *probcat);
151
152double glaguerre (long m, double b, double x);
153void initlaguerrecat(long categs, double alpha, double theta1,
154                   double *rate, double *probcat);
155void roots_laguerre(long m, double b, double **lgroot);
156
157void initgammacat (long categs, double alpha, double theta1,
158                   double *rate, double *probcat);
159
160void integrate_laguerre(long categs, double *rate,
161                        double *probcat,
162                        double (*func)(double theta, helper_fmt * b),
163                        helper_fmt *helper, double *result, double *rmax);
164
165
166/*------------------------------------------------------
167  Generalized Laguerre polynomial computed recursively.
168  For use by initgammacat
169*/
170double
171glaguerre (long m, double b, double x)
172{
173  long i;
174  double gln, glnm1, glnp1;        /* L_n, L_(n-1), L_(n+1) */
175
176  if (m == 0)
177    return 1.0;
178  else
179    {
180      if (m == 1)
181        return 1.0 + b - x;
182      else
183        {
184          gln = 1.0 + b - x;
185          glnm1 = 1.0;
186          for (i = 2; i <= m; i++)
187            {
188              glnp1 =
189                ((2 * (i - 1) + b + 1.0 - x) * gln - (i - 1 + b) * glnm1) / i;
190              glnm1 = gln;
191              gln = glnp1;
192            }
193          return gln;
194        }
195    }
196}                                /* glaguerre */
197
198
199/* calculates hermite polynomial with degree n and parameter x */
200/* seems to be unprecise for n>13 -> root finder does not converge*/
201double hermite(long n, double x)
202{
203 double   h1 = 1.;
204 double   h2 =  2. * x;
205 double   xx =  2. * x;
206 long i;
207 for(i = 1; i < n; i++)
208   {
209     xx = 2. * x * h2 - 2. * (i) * h1;
210     h1 = h2;
211     h2 = xx;
212   }
213 return xx;
214}
215
216void root_hermite(long n, double *hroot)
217{
218  long z=0;
219  long ii;
220  long start;
221  if(n % 2 == 0)
222    {
223      start = n/2;
224      z = 1;
225    }
226  else
227    {
228      start = n/2 + 1;
229      z=2;
230      hroot[start-1] = 0.0;
231    }
232  for(ii=start; ii<n; ii++)
233    {
234      /* search only upwards*/
235      hroot[ii] = halfroot(hermite,n,hroot[ii-1]+EPSILON, 1./n);
236      hroot[start - z] = -hroot[ii];
237      z++;
238    }
239}
240
241/*searches from the bound (startx) only in one direction
242  (by positive or negative delta, which results in
243  other-bound=startx+delta)
244  delta should be small.
245  (*func) is a function with two arguments
246*/
247double
248halfroot(double (*func)(long m, double x),
249         long n, double startx,
250         double delta)
251{
252  double xl;
253  double xu;
254  double xm;
255  double fu;
256  double fl;
257  double fm = 100000.;
258  double gradient;
259  boolean down;
260  /* decide if we search above or below startx and escapes to trace back
261     to the starting point that most often will be
262     the root from the previous calculation*/
263  if(delta < 0)
264    {
265      xu = startx;
266      xl = xu + delta;
267    }
268  else
269    {
270      xl = startx;
271      xu = xl + delta;
272    }
273  delta = fabs(delta);
274  fu = (*func)(n, xu);
275  fl = (*func)(n, xl);
276  gradient = (fl-fu)/(xl-xu);
277
278  while(fabs(fm) > EPSILON)
279    {
280      /* is root outside of our bracket?*/
281      if((fu<0.0 && fl<0.0) || (fu>0.0 && fl > 0.0))
282        {
283          xu += delta;
284          fu = (*func)(n, xu);
285          fl = (*func)(n, xl);
286          gradient = (fl-fu)/(xl-xu);
287          down = gradient < 0 ? TRUE : FALSE;
288        }
289      else
290        {
291          xm = xl - fl / gradient;
292          fm = (*func)(n, xm);
293          if(down)
294            {
295              if(fm > 0.)
296                {
297                  xl = xm;
298                  fl = fm;
299                }
300              else
301                {
302                  xu = xm;
303                  fu = fm;
304                }
305            }
306          else
307            {
308              if(fm > 0.)
309                {
310                  xu = xm;
311                  fu = fm;
312                }
313              else
314                {
315                  xl = xm;
316                  fl = fm;
317                }
318            }
319          gradient = (fl-fu)/(xl-xu);
320        }
321    }
322  return xm;
323}
324
325
326// calculate the weights for the hermite polynomial
327// at the roots
328// using formula Abramowitz and Stegun chapter 25.4.46 p.890
329void hermite_weight(long n, double * hroot, double * weights)
330{
331  long i;
332  double hr2;
333  double nominator = exp(LOG2 * ( n-1.) + logfac(n)) * SQRTPI / (n*n);
334  for(i=0;i<n; i++)
335    {
336      hr2 = hermite(n-1, hroot[i]);
337      weights[i] = nominator / (hr2*hr2);
338    }
339}
340
341/* calculates rates and LOG(probabilities) */
342void
343inithermitcat(long categs, double alpha, double theta1,
344              double *rate, double *probcat)
345{
346  long i;
347  double *hroot;
348  double std = SQRT2 * theta1/sqrt(alpha);
349  hroot = (double *) calloc(categs+1,sizeof(double));
350  root_hermite(categs, hroot);  // calculate roots
351  hermite_weight(categs, hroot, probcat);  // set weights
352  for(i=0;i<categs;i++)  // set rates
353    {
354      rate[i] = theta1 + std * hroot[i];
355      probcat[i] = log(probcat[i]);
356    }
357  free(hroot);
358}
359
360void
361roots_laguerre(long m, double b, double **lgroot)
362{
363  /* For use by initgammacat.
364     Get roots of m-th Generalized Laguerre polynomial, given roots
365     of (m-1)-th, these are to be stored in lgroot[m][] */
366  long i;
367  double upper, lower, x, y;
368  boolean dwn=FALSE;
369  double tmp;
370  /* is function declining in this interval? */
371  if (m == 1)
372    {
373      lgroot[1][1] = 1.0 + b;
374    }
375  else
376    {
377      dwn = TRUE;
378      for (i = 1; i <= m; i++)
379        {
380          if (i < m)
381            {
382              if (i == 1)
383                lower = 0.0;
384              else
385                lower = lgroot[m - 1][i - 1];
386              upper = lgroot[m - 1][i];
387            }
388          else
389            {                        /* i == m, must search above */
390              lower = lgroot[m - 1][i - 1];
391              x = lgroot[m - 1][m - 1];
392              do
393                {
394                  x = 2.0 * x;
395                  y = glaguerre (m, b, x);
396                }
397              while ((dwn && (y > 0.0)) || ((!dwn) && (y < 0.0)));
398              upper = x;
399            }
400          while (upper - lower > 0.000000001)
401            {
402              x = (upper + lower) / 2.0;
403              if ((tmp=glaguerre (m, b, x)) > 0.0)
404                {
405                  if (dwn)
406                    lower = x;
407                  else
408                    upper = x;
409                }
410              else
411                {
412                  if (dwn)
413                    upper = x;
414                  else
415                    lower = x;
416                }
417            }
418          lgroot[m][i] = (lower + upper) / 2.0;
419          dwn = !dwn;                /* switch for next one */
420        }
421    }
422}                                /* root_laguerre */
423
424
425void
426initgammacat (long categs, double alpha, double theta1,
427              double *rate, double *probcat)
428{
429/* calculate rates and probabilities to approximate Gamma distribution
430   of rates with "categs" categories and shape parameter "alpha" using
431   rates and weights from Generalized Laguerre quadrature */
432
433  if(alpha>=100.)
434    {
435      inithermitcat(categs, alpha, theta1, rate, probcat);
436    }
437  else
438    {
439      initlaguerrecat(categs, alpha, theta1, rate, probcat);
440    }
441}
442
443void
444initlaguerrecat(long categs, double alpha, double theta1, double *rate,
445                double *probcat)
446{
447  long i;
448  double  **lgroot;                /* roots of GLaguerre polynomials */
449  double f, x, xi, y;
450
451  lgroot = (double **) calloc(categs+1,sizeof(double*));
452  lgroot[0] = (double *) calloc((categs+1)*(categs+1),sizeof(double));
453  for(i=1;i<categs+1;i++)
454    {
455      lgroot[i] = lgroot[0] + i * (categs+1);
456    }
457  lgroot[1][1] = 1.0 + alpha;
458  for (i = 2; i <= categs; i++)
459    roots_laguerre (i, alpha, lgroot);        /* get roots for L^(a)_n */
460  /* here get weights */
461  /* Gamma weights are
462     (1+a)(1+a/2) ... (1+a/n)*x_i/((n+1)^2 [L_{n+1}^a(x_i)]^2)  */
463  f = 1;
464  for (i = 1; i <= categs; i++)
465    f *= (1.0 + alpha / i);
466  for (i = 1; i <= categs; i++)
467    {
468      xi = lgroot[categs][i];
469          y = glaguerre (categs + 1, alpha, xi);
470          x = f * xi / ((categs + 1) * (categs + 1) * y * y);
471          rate[i - 1] = xi / (1.0 + alpha);
472          probcat[i - 1] = x;
473    }
474  for(i=0;i<categs;i++)
475    {
476      probcat[i] = log(probcat[i]);
477      rate[i] *= theta1;
478    }
479  free(lgroot[0]);
480  free(lgroot);
481}                                /* initgammacat */
482
483
484void integrate_laguerre(long categs, double *rate,
485                        double *probcat,
486                        double (*func)(double theta, helper_fmt * b),
487                        helper_fmt *helper, double *result, double *rmax)
488{
489  double sum=0.;
490  long i;
491  double *temp;
492  int *stemp;
493  *rmax = -DBL_MAX;
494
495  temp = (double *) calloc(categs, sizeof(double));
496  stemp = (int *) calloc(categs, sizeof(int));
497  for(i=0;i<categs;i++)
498    {
499      temp[i] = (*func)(rate[i], (void *)helper);
500      stemp[i] = helper->sign;
501      if(temp[i] > *rmax)
502        *rmax = temp[i];
503    }
504  for(i=0;i<categs;i++)
505    {
506      sum += ((double) stemp[i]) * probcat[i] * exp(temp[i] - *rmax);
507    }
508  free(temp);
509  free(stemp);
510  *result =  sum;
511}
512
513
514/* initgammacat test function*/
515#ifdef LAGUERRE_TEST
516int main()
517{
518  long categs=10;
519  double alpha;
520  double theta1;
521  double *rate;
522  double *probcat;
523  long i;
524  printf("Enter alpha, theta1,  and categs\n");
525  fscanf(stdin,"%lf%lf%li",&alpha,&theta1,&categs);
526  rate = (double *) calloc(categs+1,sizeof(double));
527  probcat = (double *) calloc(categs+1,sizeof(double));
528  initgammacat(categs, alpha, theta1, rate,probcat);
529  for(i=0;i<categs;i++)
530    {
531      printf("%20.20f %20.20f\n",rate[i],probcat[i]);
532    }
533  return 0;
534}
535#endif
536
537
538
539
540
541
542
543
544
545
546
547
Note: See TracBrowser for help on using the repository browser.