← Back to branch summary

~ubuntu-branches/ubuntu/utopic/paml/utopic

« back to all changes in this revision

Viewing changes to src/tools.c

  • Committer: Package Import Robot
  • Author(s): Andreas Tille, Pjotr Prins, Andreas Tille
  • Date: 2012-05-15 11:10:59 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20120515111059-34qte02qzkcnzq5n
Tags: 4.5-1
http://bugs.debian.org/661519
[ Pjotr Prins ]
* Improved long description.

[ Andreas Tille ]
* New upstream version
* Fixed watch file
* debian/control: Standards-Version: 3.9.3 (no changes needed)
* debian/rules: remove useless dh-make template
* debian/paml.dirs: removed because unneeded
* debian/{install,links}: install all binaries to
  /usr/lib/debian-med/bin and symlink to /usr/bin with the
  exception of evolver which has a name space conflict.  This
  is renamed to /usr/bin/paml-evolver
  Closes: #661519
* debian/{paml.docs,paml-doc.docs}:
   - Adapted to new names
   - do not duplicate the same files into both packages (paml and
     paml-doc)
* debian/profile.d/paml: Set PATH to find paml executables under
  their expected names
* debian/README.Debian: Document the name space conflict solution
* debian/upstream: Citation information
* debhelper 9 (control+compat)
* debian/patches/hardening.patch: Enable propagation of hardening
  flags

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools.c
678
678
      { if(rates) FOR(h,ls) rates[h]=1; }
679
679
   else {
680
680
      if (K>1) {
681
 
         DiscreteGamma (freqK, rK, alpha, alpha, K, DGammaMean);
 
681
         DiscreteGamma(freqK, rK, alpha, alpha, K, DGammaUseMedian);
682
682
 
683
683
         MultiNomialAliasSetTable(K, freqK, Falias, Lalias, space+5*K);
684
684
         MultiNomialAlias(ls, K, Falias, Lalias, counts);
1346
1346
   double x0=x, nroundsVirtual, small=1e-50;
1347
1347
 
1348
1348
   if(b-a<small) {
1349
 
      printf("\nimproper range x0=%.6g (%.6g, %.6g)\n", x0,a,b);
 
1349
      printf("\nimproper range x0=%.6g (%.6g, %.6g)\n", x0, a, b);
1350
1350
      exit(-1);
1351
1351
   }
1352
1352
   nroundsVirtual = floor((x-a)/(2*(b-a)));
1516
1516
   if(a0 < 1) 
1517
1517
      v *= pow(rndu(), 1/a0);
1518
1518
   if(v==0) 
1519
 
      printf("\nrndgamma returning 0.\n");
 
1519
      printf("\a\nrndgamma returning 0.\n");
1520
1520
   return v;
1521
1521
}
1522
1522
 
1803
1803
}
1804
1804
 
1805
1805
 
 
1806
int StirlingS2(int n, int k)
 
1807
{
 
1808
/* Stirling number of the second type, calculated using the recursion (loop)
 
1809
   S(n, k) = S(n - 1, k - 1) + k*S(n - 1, k).
 
1810
   This works for small numbers of n<=15.
 
1811
*/
 
1812
   int S[16]={0}, i, j;
 
1813
 
 
1814
   if((n==0 && k==0) || k==1 || k==n)
 
1815
      return 1;
 
1816
   if(k==0 || k>n)
 
1817
      return 0;
 
1818
   if(k==2)
 
1819
      return (int) ldexp(1,n-1) - 1;
 
1820
   if(k==n-1)
 
1821
      return n*(n-1)/2;
 
1822
   if(n>15)
 
1823
      error2("n>15 too large in StirlingS2()");
 
1824
 
 
1825
   S[1] = S[2] = 1;  /* start with n = 2 */
 
1826
   for(i=3; i<=n; i++) {
 
1827
      for(j=min2(k,i); j>=2; j--)
 
1828
         S[j] = S[j-1] + j*S[j];
 
1829
   }
 
1830
   return S[k];
 
1831
}
 
1832
 
 
1833
double lnStirlingS2(int n, int k)
 
1834
{
 
1835
/* This calculates the logarithm of the Stirling number of the second type.
 
1836
   Temme NM. 1993. Asymptotic estimates of Stirling numbers. Stud Appl Math 89:233-243.
 
1837
*/
 
1838
   int i;
 
1839
   double lnS=0, t0, x0, x, A, nk, y;
 
1840
 
 
1841
   if(k>n) error2("k<n in lnStirlingS2");
 
1842
 
 
1843
   if(n==0 && k==0)
 
1844
      return 0;
 
1845
   if(k==0)
 
1846
      return -1e300;
 
1847
   if(k==1 || k==n) 
 
1848
      return (0);
 
1849
   if(k==2)
 
1850
      return (n<50 ? log(ldexp(1,n-1) - 1) : (n-1)*0.693147);
 
1851
   if(k==n-1)
 
1852
      return log(n*(n-1)/2.0);
 
1853
   if(n<8)
 
1854
      return log((double)StirlingS2(n, k));
 
1855
   
 
1856
   nk = (double)n/k;
 
1857
   for(i=0,x0=x=1; i<10000; i++) {
 
1858
      x = (x0 + nk - nk*exp(-x0))/2;
 
1859
      if(fabs(x-x0)/(1+x) < 1e-10) break;
 
1860
      x0 = x;
 
1861
   }
 
1862
   t0 = n/(double)k - 1;
 
1863
   if(x<100)
 
1864
      A = -n * log(x) + k*log(exp(x) - 1);
 
1865
   else 
 
1866
      A = -n * log(x) + k*x;
 
1867
 
 
1868
   A += -k*t0 + (n-k)*log(t0);
 
1869
   lnS = A + (n-k)*log((double)k) + 0.5*log(t0/((1 + t0)*(x - t0)));
 
1870
   lnS += log(Binomial(n, k, &y));
 
1871
   lnS += y;
 
1872
 
 
1873
   return(lnS);
 
1874
}
 
1875
 
 
1876
 
1806
1877
double LnGamma (double x)
1807
1878
{
1808
1879
/* returns ln(gamma(x)) for x>0, accurate to 10 decimal places.
1991
2062
}
1992
2063
 
1993
2064
 
1994
 
int DiscreteGamma (double freqK[], double rK[], double alpha, double beta, int K, int dgammamean)
1995
 
{
1996
 
/* discretization of gamma distribution with equal proportions in each 
1997
 
   category.
1998
 
*/
1999
 
   int i;
2000
 
   double t, factor=alpha/beta*K, lnga1;
2001
 
 
2002
 
   if (dgammamean==0) {   /* mean */
 
2065
int DiscreteBeta (double freq[], double x[], double p, double q, int K, int UseMedian)
 
2066
{
 
2067
/* discretization of beta(p, q), with equal proportions in each category.
 
2068
*/
 
2069
   int i;
 
2070
   double mean=p/(p+q), lnbeta, lnbeta1, t;
 
2071
 
 
2072
   lnbeta = LnBeta(p, q);
 
2073
   if(UseMedian) {   /* median */
 
2074
      for(i=0,t=0; i<K; i++)
 
2075
         t += x[i] = QuantileBeta((i+0.5)/K, p, q, lnbeta);
 
2076
      for(i=0; i<K; i++)
 
2077
         x[i] *= mean*K/t;
 
2078
 
 
2079
      /* printf("\nmedian  ");  for(i=0; i<K; i++) printf("%9.5f", x[i]); */
 
2080
   }
 
2081
   else {            /* mean */
 
2082
      for(i=0; i<K-1; i++) /* cutting points */
 
2083
         freq[i] = QuantileBeta((i+1.0)/K, p, q, lnbeta);
 
2084
      freq[K-1] = 1;
 
2085
 
 
2086
      /* printf("\npoints  ");  for(i=0; i<K; i++) printf("%9.5f", freq[i]); */
 
2087
      lnbeta1 = lnbeta - log(1 + q/p);
 
2088
      for(i=0; i<K-1; i++) /* CDF */
 
2089
         freq[i] = CDFBeta(freq[i], p+1, q, lnbeta1);
 
2090
 
 
2091
      x[0] = freq[0]*mean*K;
 
2092
      for (i=1; i<K-1; i++)  x[i] = (freq[i] - freq[i-1])*mean*K;
 
2093
      x[K-1] = (1 - freq[K-2])*mean*K;
 
2094
 
 
2095
      /* printf("\nmean    ");  for(i=0; i<K; i++) printf("%9.5f", x[i]); */
 
2096
      for(i=0,t=0; i<K; i++) t += x[i]/K;
 
2097
   }
 
2098
 
 
2099
   for (i=0; i<K; i++) freq[i] = 1.0/K;
 
2100
   return (0);
 
2101
}
 
2102
 
 
2103
int DiscreteGamma (double freqK[], double rK[], double alpha, double beta, int K, int UseMedian)
 
2104
{
 
2105
/* discretization of G(alpha, beta) with equal proportions in each category.
 
2106
*/
 
2107
   int i;
 
2108
   double t, mean=alpha/beta, lnga1;
 
2109
 
 
2110
   if(UseMedian) {   /* median */
 
2111
      for(i=0; i<K; i++) rK[i] = QuantileGamma((i*2.+1)/(2.*K), alpha, beta);
 
2112
      for(i=0,t=0; i<K; i++) t += rK[i];
 
2113
      for(i=0; i<K; i++) rK[i] *= mean*K/t;   /* rescale so that the mean is alpha/beta. */
 
2114
   }
 
2115
   else {            /* mean */
2003
2116
      lnga1 = LnGamma(alpha+1);
2004
2117
      for (i=0; i<K-1; i++) /* cutting points, Eq. 9 */
2005
2118
         freqK[i] = QuantileGamma((i+1.0)/K, alpha, beta);
2006
2119
      for (i=0; i<K-1; i++) /* Eq. 10 */
2007
2120
         freqK[i] = IncompleteGamma(freqK[i]*beta, alpha+1, lnga1);
2008
 
      rK[0] = freqK[0]*factor;
2009
 
      rK[K-1] = (1-freqK[K-2])*factor;
2010
 
      for (i=1; i<K-1; i++)  rK[i] = (freqK[i] - freqK[i-1])*factor;
2011
 
   }
2012
 
   else {                 /* median */
2013
 
      for(i=0; i<K; i++) rK[i] = QuantileGamma((i*2.+1)/(2.*K), alpha, beta);
2014
 
      for(i=0,t=0; i<K; i++) t += rK[i];
2015
 
      for(i=0; i<K; i++) rK[i] *= factor/t;
2016
 
   }
 
2121
      rK[0] = freqK[0]*mean*K;
 
2122
      for (i=1; i<K-1; i++)  rK[i] = (freqK[i] - freqK[i-1])*mean*K;
 
2123
      rK[K-1] = (1-freqK[K-2])*mean*K;
 
2124
   }
 
2125
 
2017
2126
   for (i=0; i<K; i++) freqK[i] = 1.0/K;
2018
2127
   
2019
2128
   return (0);
2055
2164
      }
2056
2165
   }
2057
2166
 
2058
 
   DiscreteGamma (freqK, rK, alpha, alpha, K, DGammaMean);
 
2167
   DiscreteGamma(freqK, rK, alpha, alpha, K, DGammaUseMedian);
2059
2168
 
2060
2169
   for (i=0,v1=*rho1=0; i<K; i++) {
2061
2170
      v1+=rK[i]*rK[i]*freqK[i];
2586
2695
      q = pin;
2587
2696
   }
2588
2697
 
2589
 
   if(lnbeta==0) lnbeta = LnGamma(p) + LnGamma(q) - LnGamma(p+q);
 
2698
   if(lnbeta==0) lnbeta = LnBeta(p, q);
2590
2699
 
2591
2700
   if ((p + q) * y / (p + 1) < eps) {  /* tail approximation */
2592
2701
      ans = 0;
2685
2794
   if (prob == 0 || prob == 1)
2686
2795
      return prob;
2687
2796
 
2688
 
   if(lnbeta==0) lnbeta=LnGamma(p)+LnGamma(q)-LnGamma(p+q);
 
2797
   if(lnbeta==0) lnbeta = LnBeta(p, q);
2689
2798
 
2690
2799
   /* change tail if necessary;  afterwards   0 < a <= 1/2    */
2691
2800
   if (prob <= 0.5) {
3877
3986
/* this approximates the integral Nintegrate[fun[x], {x,a,b}].
3878
3987
   npoints is 10, 20, 32 or 64 nodes for legendre.");
3879
3988
*/
3880
 
   int j;
3881
 
   double *x=NULL, *w=NULL, s=0, t[2];
 
3989
   int j, ixw;
 
3990
   double *x=NULL, *w=NULL, sign, s=0, t;
3882
3991
 
 
3992
   if(npoints%2 != 0)
 
3993
      error2("this assumes even number of points.");
3883
3994
   GaussLegendreRule(&x, &w, npoints);
3884
3995
 
3885
 
   /* this goes through the natural order from a to b */
3886
 
   for(j=npoints/2-1; j>=0; j--) {
3887
 
      t[1] = (a+b)/2 - (b-a)/2*x[j];
3888
 
      s += w[j]*fun(t[1]);
3889
 
   }
3890
 
   for(j=0; j<npoints/2; j++) {
3891
 
      t[0] = (a+b)/2 + (b-a)/2*x[j];
3892
 
      s += w[j]*fun(t[0]);
3893
 
   }
3894
 
 
3895
 
/*
3896
 
   for(j=0,s=0; j<npoints/2; j++) {
3897
 
      t[0] = (a+b)/2 + (b-a)/2*x[j];
3898
 
      t[1] = (a+b)/2 - (b-a)/2*x[j];
3899
 
      for(i=0; i<2; i++)
3900
 
         s += w[j]*fun(t[i]);
3901
 
   }
3902
 
*/
3903
 
   return s *= (b-a)/2;
 
3996
   /* x changes monotonically from a to b. */
 
3997
   for(j=0; j<npoints; j++) {
 
3998
      if(j<npoints/2) { ixw = npoints/2-1-j;  sign=-1; }
 
3999
      else            { ixw = j-npoints/2;    sign=1; }
 
4000
      t = (a+b)/2 + sign*(b-a)/2*x[ixw];
 
4001
      s += w[ixw]*fun(t);
 
4002
   }
 
4003
   return s *= (b - a)/2;
3904
4004
}
3905
4005
 
3906
4006
 
4305
4405
 
4306
4406
      P(t) = (I + Qt/m + (Qt/m)^2/2)^m, with m = 2^TimeSquare.
4307
4407
 
4308
 
   See equation (2.22) and the discussion below in Yang (2006).
 
4408
   See equation (2.22) in Yang (2006) and the discussion below it.
4309
4409
   T[it=0] is the current matrix, and T[it=1] is the squared result matrix,
4310
4410
   used to avoid copying matrices.
4311
4411
   Use an even TimeSquare to avoid one round of matrix copying.
4330
4430
      matby (T[1-it], T[1-it], T[it], n, n, n);
4331
4431
   }
4332
4432
   if(it==1) 
4333
 
      for(i=0;i<n*n;i++) Q[i]=T[1][i];
 
4433
      for(i=0;i<n*n;i++) Q[i] = T[1][i];
4334
4434
   return(0);
4335
4435
}
4336
4436
 
5465
5565
   for(i=0; i<n; i++) {
5466
5566
      for (j=i; j<n; j++)  {
5467
5567
         for(k=0; k<n; k++) x1[k] = x[k];
5468
 
         x1[i] += h[i];    x1[j] += h[j];
 
5568
         x1[i] += h[i]; 
 
5569
         x1[j] += h[j];
5469
5570
         fpp = (*fun)(x1,n);                  /* (+hi, +hj) */
5470
 
         x1[i] -= h[i]*2;  x1[j] -= h[j]*2;
 
5571
         x1[i] -= h[i]*2;
 
5572
         x1[j] -= h[j]*2;
5471
5573
         fmm = (*fun)(x1,n);                  /* (-hi, -hj) */
5472
5574
         if (i==j)  {
5473
 
             H[i*n+i] = (fpp+fmm-2*f0)/(4*h[i]*h[i]);
 
5575
             H[i*n+i] = (fpp + fmm - 2*f0)/(4*h[i]*h[i]);
5474
5576
             g[i] = (fpp - fmm)/(h[i]*4);
5475
5577
         }
5476
5578
         else {
5477
5579
            x1[i] += 2*h[i];                     fpm = (*fun)(x1,n);  /* (+hi, -hj) */
5478
5580
            x1[i] -= 2*h[i];   x1[j] += 2*h[j];  fmp = (*fun)(x1,n);  /* (-hi, +hj) */
5479
 
            H[i*n+j] = H[j*n+i] = (fpp+fmm-fpm-fmp)/(4*h[i]*h[j]);
 
5581
            H[i*n+j] = H[j*n+i] = (fpp + fmm - fpm - fmp)/(4*h[i]*h[j]);
5480
5582
         }
5481
5583
      }
5482
5584
   }
6007
6109
   double *x0=space, *x1=space+n, eh0=Small_Diff, eh;  /* eh0=1e-6 || 1e-7 */
6008
6110
 
6009
6111
   FOR(i,n) {
6010
 
      eh=eh0*(fabs(x[i])+1);
 
6112
      eh = eh0*(fabs(x[i])+1);
6011
6113
      if (xmark[i]==0 && (AlwaysCenter || SIZEp<1)) {   /* central */
6012
 
         FOR (j, n)  x0[j]=x1[j]=x[j];
6013
 
         eh=pow(eh,.67);  x0[i]-=eh; x1[i]+=eh;
6014
 
         g[i]=((*fun)(x1,n)-(*fun)(x0,n))/(eh*2.0);
 
6114
         for(j=0; j<n; j++)  x0[j] = x1[j] = x[j];
 
6115
         eh = pow(eh, .67);  x0[i] -= eh;  x1[i] += eh;
 
6116
         g[i] = ((*fun)(x1,n) - (*fun)(x0,n))/(eh*2.0);
6015
6117
      }
6016
6118
      else  {                         /* forward or backward */
6017
 
         FOR (j, n)  x1[j]=x[j];
6018
 
         if (xmark[i]) eh*=-xmark[i];
 
6119
         for(j=0; j<n; j++)  x1[j] = x[j];
 
6120
         if (xmark[i]) eh *= -xmark[i];
6019
6121
         x1[i] += eh;
6020
 
         g[i]=((*fun)(x1,n)-f0)/eh;
 
6122
         g[i] = ((*fun)(x1,n) - f0)/eh;
6021
6123
      }
6022
6124
   }
6023
6125
   return(0);
6191
6293
      FOR (i,nfree)  FOR (j,nfree)  
6192
6294
         H[i*nfree+j] += s[i]*s[j]/v - z[i]*z[j]/w;
6193
6295
#else /* BFGS */
6194
 
      for (i=0,w=v=0.; i<nfree; i++) {
6195
 
         for (j=0,z[i]=0.; j<nfree; j++) z[i]+=H[i*nfree+j]*y[j];
6196
 
         w+=y[i]*z[i];    v+=y[i]*s[i];
 
6296
      for(i=0,w=v=0.; i<nfree; i++) {
 
6297
         for(j=0,z[i]=0.; j<nfree; j++) z[i] += H[i*nfree+j]*y[j];
 
6298
         w += y[i]*z[i];    v += y[i]*s[i];
6197
6299
      }
6198
6300
      if (fabs(v)<small)   { identity(H,nfree); fail=1; continue; }
6199
6301
      FOR (i,nfree)  FOR (j,nfree)
6203
6305
   }    /* for (Iround,maxround)  */
6204
6306
 
6205
6307
   /* try to remove this after updating LineSearch2() */
6206
 
   *f=(*fun)(x,n);  
 
6308
   *f = (*fun)(x,n);
6207
6309
   if(noisy>2) FPN(F0);
6208
6310
 
6209
 
   if (Iround==maxround) {
 
6311
   if(Iround==maxround) {
6210
6312
      if (fout) fprintf (fout,"\ncheck convergence!\n");
6211
6313
      return(-1);
6212
6314
   }