1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate  and effectice population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6  Gammalike and derivatives  R O U T I N E S
7 
8  Peter Beerli 2001, Seattle
9  beerli@fsu.edu
10 
11  Copyright 2001-2003 Peter Beerli and Joseph Felsenstein
12  Copyright 2004-2005 Peter Beerli
13 
14 Permission is hereby granted, free of charge, to any person obtaining a copy
15 of this software and associated documentation files (the "Software"), to deal
16 in the Software without restriction, including without limitation the rights
17 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
18 of the Software, and to permit persons to whom the Software is furnished to do
19 so, subject to the following conditions:
20 
21 The above copyright notice and this permission notice shall be included in all copies
22 or substantial portions of the Software.
23 
24 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
25 INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
26 PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
27 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
28 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
29 OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
30 
31 $Id: gammalike.c 2120 2013-01-03 05:17:21Z beerli $
32 */
33 /*! \file gammalike.c
34 
35 Calculates the gamma distributed rate variation among loci
36 and derivatives
37 */
38 
39 
40 #include "migration.h"
41 #include "tools.h"
42 #include "broyden.h"
43 #include "combroyden.h"
44 #include "laguerre.h"
45 #include "sighandler.h"
46 
47 MYREAL gamma_loci_like (helper_fmt *helper, MYREAL *oparam,
48                         MYREAL *olparam, MYREAL denom);
49 MYREAL gamma_locus_like (nr_fmt * nr, MYREAL *oparam,
50                          MYREAL *olparam, MYREAL denom, long locus);
51 MYREAL loggf (MYREAL r, MYREAL alpha, MYREAL theta1);
52 void gamma_loci_derivative (helper_fmt * helper);
53 void gamma_locus_derivative (helper_fmt * helper, long locus);
54 //MYREAL gamma_theta1_derivative (MYREAL theta, MYREAL *param, nr_fmt * nr,
55 //                              long locus, int *sign);
56 //MYREAL gamma_other_derivative (long which, MYREAL theta, MYREAL *param,
57 //                             nr_fmt * nr, long locus, int *sign);
58 //MYREAL gamma_alpha_derivative (MYREAL theta, MYREAL *param, nr_fmt * nr,
59 //                             long locus, int *sign);
60 void gamma_parts_gradient (MYREAL x, MYREAL *param, nr_fmt * nr, long locus);
61 MYREAL dt_locus_gsum (long which, MYREAL x, MYREAL *param, nr_fmt * nr,
62                       timearchive_fmt ** atl, long locus);
63 
64 
65 MYREAL
gamma_loci_like(helper_fmt * helper,MYREAL * oparam,MYREAL * olparam,MYREAL denom)66 gamma_loci_like(helper_fmt *helper, MYREAL *oparam, MYREAL *olparam, MYREAL denom)
67 {
68     long locus;
69     MYREAL summ = 0.0;
70     nr_fmt *nr = helper->nr;
71     MYREAL *locusweight = helper->nr->world->data->locusweight; //invariant loci treatment
72     for (locus = 0; locus < nr->world->loci; locus++)
73     {
74         if (nr->skiploci[locus])
75             continue;
76         nr->locilikes[locus] =
77             gamma_locus_like (nr, oparam, olparam, denom, locus);
78         summ += locusweight[locus] * nr->locilikes[locus];
79     }
80     return summ;
81 }
82 
83 
84 MYREAL find_gamma_summ(MYREAL la, MYREAL theta1, nr_fmt *nr, MYREAL *oparam, MYREAL *olparam,
85                        long locus, helper_fmt *helper);
86 
87 
88 MYREAL
gamma_loci_like_iterative(helper_fmt * helper,MYREAL * oparam,MYREAL * olparam,MYREAL denom)89 gamma_loci_like_iterative(helper_fmt * helper, MYREAL *oparam, MYREAL *olparam, MYREAL denom)
90 {
91     // hack to get around the problem of not being able to get a clean derivative of alpha
92     // uses a simple maxima finder for alpha
93     long locus=helper->nr->world->locus;
94     MYREAL fa = 0.0;
95     MYREAL fla, fha;
96     MYREAL la  = -10;
97     MYREAL ha = +4;
98     MYREAL a;
99     MYREAL fm, m;
100     MYREAL theta1 = helper->nr->param[0];
101     nr_fmt *nr = helper->nr;
102     char save_a = nr->world->options->custm[nr->world->numpop2];
103     //  if(save_a == 'c')
104     // {
105     // return  gamma_loci_like_old (helper, oparam, olparam, denom);
106     // }
107     nr->world->options->custm[nr->world->numpop2]='c';
108     nr->world->options->custm2[nr->world->numpop2]='c';
109     fla = find_gamma_summ( la,  theta1, nr,  oparam, olparam,  locus, helper);
110     fha = find_gamma_summ( ha,  theta1, nr,  oparam, olparam,  locus, helper);
111     a = (ha + la)/2.;
112     fa = find_gamma_summ( a,  theta1, nr,  oparam, olparam,  locus, helper);
113     m = (la + a) / 2.;
114     fm = find_gamma_summ( m,  theta1, nr,  oparam, olparam,  locus, helper);
115     while(fabs(fla-fha)> EPSILON && fabs(ha-la) > EPSILON)
116     {
117         if(fm < fa)
118         {
119             fla = fm;
120             la = m;
121         }
122         else
123         {
124             fha = fa;
125             ha  = a;
126             a = m;
127            //xcode  fa = fm;
128         }
129         m = (ha + la)/2.;
130         fa = find_gamma_summ( m,  theta1, nr,  oparam, olparam,  locus, helper);
131 
132     }
133     nr->param[nr->numpop2] = EXP(m);
134     nr->lparam[nr->numpop2] = m;
135     nr->world->options->custm[nr->world->numpop2]=save_a;
136     nr->world->options->custm2[nr->world->numpop2]=save_a;
137     return fa;
138 }
139 
140 
141 MYREAL
gamma_loci_like_new_notworking(helper_fmt * helper,MYREAL * oparam,MYREAL * olparam,MYREAL denom)142 gamma_loci_like_new_notworking(helper_fmt * helper, MYREAL *oparam, MYREAL *olparam, MYREAL denom)
143 {
144     long locus=0;
145     long trials=0;
146     MYREAL logalpha= -MYREAL_MAX;
147     MYREAL psiabc;
148     MYREAL theta1 = helper->nr->param[0];
149     nr_fmt *nr = helper->nr;
150     char save_a = nr->world->options->custm[nr->world->numpop2];
151 
152     MYREAL a, b, c,  m=0.5;
153     MYREAL la, lb, lc, ll=0, ql;
154     //MYREAL lm;
155     a= 0.00000000001;
156     b= 1.0;
157     c=2.0;
158     nr->world->options->custm[nr->world->numpop2]='c';
159     nr->world->options->custm2[nr->world->numpop2]='c';
160     la = find_gamma_summ( a,  theta1, nr,  oparam, olparam,  locus, helper);
161     lb = find_gamma_summ( b,  theta1, nr,  oparam, olparam,  locus, helper);
162     lc = find_gamma_summ( c,  theta1, nr,  oparam, olparam,  locus, helper);
163 
164     while(trials++ < NTRIALS)
165     {
166         logalpha =
167             (c * c * (lb - la) + b * b * (la - lc) +
168              a * a * (lc - lb)) / (2. * (b * la - c * la - a * lb + c * lb +
169                                          a * lc - b * lc));
170         psiabc = ((la - lb) / (-a + b) + (la - lc) / (a - c)) / (-b + c);
171         if ((psiabc <= 0.0) || (logalpha >= m))
172         {
173             if (a == m)
174             {
175                 ll = la;
176                 logalpha = m;
177                 break;
178             }
179             if (b == m)
180             {
181                 ll = lb;
182                 logalpha = m;
183                 break;
184             }
185             if (c == m)
186             {
187                 ll = lc;
188                 logalpha = m;
189                 break;
190             }
191             ll = find_gamma_summ( m,  theta1, nr,  oparam, olparam,  locus, helper);
192             replace_with ((int) m, &a, &b, &c, m, &la, &lb, &lc, ll);
193         }
194         ll = find_gamma_summ( logalpha,  theta1, nr,  oparam, olparam,  locus, helper);
195         ql = quadratic_lamda (logalpha, a, b, c, la, lb, lc);
196         if ((fabs (ll - MYMIN3 (la, lb, lc)) <= BIGEPSILON)
197                 || (fabs (ll - ql) <= BIGEPSILON))
198         {
199             break;
200         }
201         else
202         {
203             replace_with (1, &a, &b, &c, logalpha, &la, &lb, &lc,
204                           ll);
205             m = MYMAX3 (a, b, c);
206         }
207     }
208     nr->param[nr->numpop2] = EXP(logalpha);
209     nr->lparam[nr->numpop2] = logalpha;
210     nr->world->options->custm[nr->world->numpop2]=save_a;
211     nr->world->options->custm2[nr->world->numpop2]=save_a;
212     return ll;
213 }
214 
find_gamma_summ(MYREAL lla,MYREAL theta1,nr_fmt * nr,MYREAL * oparam,MYREAL * olparam,long locus,helper_fmt * helper)215 MYREAL find_gamma_summ(MYREAL lla, MYREAL theta1, nr_fmt *nr, MYREAL *oparam, MYREAL *olparam,
216                        long locus, helper_fmt *helper)
217 {
218     MYREAL la = EXP(lla);
219     MYREAL denom =  -LGAMMA(la) - la * log(theta1/la);
220     MYREAL summ=0.0;
221     MYREAL so = oparam[nr->numpop2];
222     MYREAL slo = olparam[nr->numpop2];
223     MYREAL *locusweight = helper->nr->world->data->locusweight; //invariant loci treatment
224     oparam[nr->numpop2] = la;
225     olparam[nr->numpop2] = lla;
226     if (olparam[nr->numpop2] > 9.903487553)
227     {
228         olparam[nr->numpop2] = 9.903487553;
229     }
230     initgammacat (nr->categs, oparam[nr->numpop2],1./* EXP (lparam[0])*/,
231                   nr->rate, nr->probcat);
232     //calc_loci_param (nr, olparam, oparam, helper->dv, helper->lamda, nr->partsize);
233     for (locus = 0; locus < nr->world->loci; locus++)
234     {
235         if (nr->skiploci[locus])
236             continue;
237         nr->locilikes[locus] = gamma_locus_like (nr, oparam, olparam, denom, locus);
238         summ += locusweight[locus] * nr->locilikes[locus]; //invariant loci treatment
239     }
240     oparam[nr->numpop2] = so;
241     olparam[nr->numpop2] = slo;
242     return summ;
243 }
244 
245 MYREAL
gamma_locus_like(nr_fmt * nr,MYREAL * oparam,MYREAL * olparam,MYREAL denom,long locus)246 gamma_locus_like (nr_fmt * nr, MYREAL *oparam, MYREAL *olparam, MYREAL denom,
247                   long locus)
248 {
249     //  nr_fmt *nr = helper->nr;
250     long r;
251     MYREAL tempmax = -MYREAL_MAX;
252     MYREAL *temp;
253     MYREAL summ = 0.;
254     MYREAL alpha = oparam[nr->numpop2];
255     MYREAL theta1 = oparam[0];
256     temp = (MYREAL *) mycalloc (GAMMA_INTERVALS, sizeof (MYREAL));
257     for (r = 0; r < GAMMA_INTERVALS; r++)
258     {
259         set_gamma_param (nr->param, oparam,
260                          nr->lparam, olparam, nr->rate[r], nr);
261         temp[r] = nr->probcat[r] + loggf (nr->rate[r], alpha, theta1) - denom +
262                   calc_locus_like (nr, nr->param, nr->lparam, locus);
263         //      printf("%f ",temp[r]);
264         if (temp[r] > tempmax)
265             tempmax = temp[r];
266     }
267     //  printf("\n");
268     for (r = 0; r < GAMMA_INTERVALS; r++)
269         summ += EXP (temp[r] - tempmax);
270     myfree(temp);
271     return LOG (summ) + tempmax;
272 }
273 
274 MYREAL
loggf(MYREAL r,MYREAL alpha,MYREAL theta1)275 loggf (MYREAL r, MYREAL alpha, MYREAL theta1)
276 {
277     return -r * alpha / theta1 + LOG (r) * (alpha - 1.);
278 }
279 
280 void
gamma_loci_difference(helper_fmt * helper)281 gamma_loci_difference (helper_fmt * helper)
282 {
283     long i;
284     nr_fmt *nr = helper->nr;
285     long size = nr->partsize;
286     MYREAL mle, high, low, param1;
287     MYREAL *lparam;
288     MYREAL *param;
289     MYREAL *xv = helper->xv;
290     MYREAL *expxv = helper->expxv;
291     MYREAL epsilon = EPSILON;
292     //  MYREAL logepsilon = log(epsilon);
293     memset (nr->d, 0, size * sizeof (MYREAL));
294     setdoublevec1d (&lparam, helper->xv, size);
295     setdoublevec1d (&param, helper->expxv, size);
296     mle = gamma_loci_like (helper, param, lparam, helper->weight);
297     for (i = 0; i < size; i++)
298     {
299         memcpy (param, expxv, sizeof (MYREAL) * size);
300         memcpy (lparam, xv, sizeof (MYREAL) * size);
301         param1 = (param[i] += epsilon);
302         lparam[i] = LOG (param[i]);
303         high = gamma_loci_like (helper, param, lparam, helper->weight);
304         param[i] -= 2. * epsilon;
305         if (param[i] < 0.0)
306         {
307             low = mle;
308             param[i] = expxv[i];
309         }
310         else
311         {
312             lparam[i] = LOG (param[i]);
313             low = gamma_loci_like (helper, param, lparam, helper->weight);
314         }
315         nr->d[i] = (high - low) / (param1 - param[i]);
316     }
317     myfree(lparam);
318     myfree(param);
319 }
320 
321 void
gamma_loci_derivative(helper_fmt * helper)322 gamma_loci_derivative (helper_fmt * helper)
323 {
324     long locus;
325     memset (helper->nr->d, 0, helper->nr->partsize * sizeof (MYREAL));
326     for (locus = 0; locus < helper->nr->world->loci; locus++)
327     {
328         if (helper->nr->skiploci[locus])
329             continue;
330         //adds up in nr->d
331         gamma_locus_derivative (helper, locus);
332     }
333 }
334 
335 void
gamma_locus_derivative(helper_fmt * helper,long locus)336 gamma_locus_derivative (helper_fmt * helper, long locus)
337 {
338     long r, i, ii, z;
339     MYREAL **temp, *tempmax;
340     MYREAL ll;
341     MYREAL sumi;
342     int **stemp;
343     //MYREAL *summll;
344     //MYREAL summllsum;
345     //MYREAL summllmax;
346     MYREAL rate;
347     MYREAL lograte;
348     MYREAL tmp;
349     MYREAL denom = helper->weight;
350     nr_fmt *nr = helper->nr;
351     MYREAL alpha = helper->expxv[nr->numpop2];
352     MYREAL theta1 = helper->expxv[0];
353     long partsize = nr->partsize;
354     //summll = (MYREAL *) mycalloc(GAMMA_INTERVALS, sizeof(MYREAL));
355     tempmax = (MYREAL *) mycalloc (partsize, sizeof (MYREAL));
356     temp = (MYREAL **) mycalloc (partsize, sizeof (MYREAL *));
357     temp[0] = (MYREAL *) mycalloc (partsize * GAMMA_INTERVALS, sizeof (MYREAL));
358     stemp = (int **) mycalloc (partsize, sizeof (int *));
359     stemp[0] = (int *) mycalloc (partsize * GAMMA_INTERVALS, sizeof (int));
360     tempmax[0] = -MYREAL_MAX;
361     for (r = 1; r < partsize; r++)
362     {
363         temp[r] = temp[0] + r * GAMMA_INTERVALS;
364         stemp[r] = stemp[0] + r * GAMMA_INTERVALS;
365         tempmax[r] = -MYREAL_MAX;
366     }
367     //summllmax = -MYREAL_MAX;
368     //summllsum=0.0;
369     for (r = 0; r < GAMMA_INTERVALS; r++)
370     {
371         rate = nr->rate[r];
372         lograte = LOG (rate);
373         set_gamma_param (nr->param, helper->expxv,
374                          nr->lparam, helper->xv, rate, nr);
375         ll = nr->probcat[r] + calc_locus_like (nr, nr->param, nr->lparam, locus);
376         ll +=  -rate * alpha / theta1 + lograte * (alpha - 1.);
377         //ll +=  -rate * alpha + lograte * (alpha - 1.);
378         gamma_parts_gradient (rate, helper->expxv, nr, locus); //changes nr->parts
379         //first of the derivatives: theta1
380         //     tmp =  nr->parts[0];
381         tmp = rate * alpha / (theta1 * theta1) + nr->parts[0];
382 
383         stemp[0][r] = tmp < 0. ? -1 : 1;
384         temp[0][r] = ll + LOG (fabs (tmp));
385         if (tempmax[0] < temp[0][r])
386             tempmax[0] = temp[0][r];
387         //all others except alpha
388         for (i = 1; i < nr->numpop2; i++)
389         {
390             stemp[i][r] = (nr->parts[i] < 0.) ? -1 : 1;
391             temp[i][r] = ll + LOG (fabs (nr->parts[i]));
392             if (tempmax[i] < temp[i][r])
393                 tempmax[i] = temp[i][r];
394         }
395         // last derivative: alpha
396         tmp = -rate/theta1 + lograte;
397         //tmp = - theta1*lograte + rate;
398         stemp[i][r] = tmp < 0. ? -1 : 1;
399         temp[i][r] = ll + LOG (fabs (tmp));
400         if (tempmax[i] < temp[i][r])
401             tempmax[i] = temp[i][r];
402     }
403     z = 0;
404     for (ii = 0; ii < nr->partsize - nr->profilenum; ii++)
405     {
406         i = (nr->profilenum > 0) ? nr->indeks[z++] : ii;
407         sumi = 0;
408         //   denom = -LGAMMA(alpha) - alpha * log(1./alpha);
409         denom = -LGAMMA(alpha) - alpha * log(theta1/alpha);
410         for (r = 0; r < GAMMA_INTERVALS; r++)
411             sumi +=
412                 ((MYREAL) stemp[i][r]) * EXP (temp[i][r] -
413                                               tempmax[i]  +  denom);
414         if (i == 0)
415         {
416             sumi =
417                 //     sumi * EXP (tempmax[i] - nr->locilikes[locus]);
418                 sumi * EXP (tempmax[i] - nr->locilikes[locus]) - alpha / theta1;
419         }
420         else
421         {
422             if (i == nr->numpop2)
423             {
424                 sumi = sumi *(EXP (tempmax[i] - nr->locilikes[locus]))
425                        - polygamma (0L, alpha) - LOG (theta1 / alpha) +1.;
426             }
427             else
428                 sumi = sumi * EXP (tempmax[i] - nr->locilikes[locus]);
429         }
430         nr->d[ii] += sumi;
431     }
432     myfree(temp[0]);
433     myfree(temp);
434     myfree(stemp[0]);
435     myfree(stemp);
436     //myfree(summll);
437 }
438 
439 void
gamma_parts_gradient(MYREAL x,MYREAL * param,nr_fmt * nr,long locus)440 gamma_parts_gradient (MYREAL x, MYREAL *param, nr_fmt * nr, long locus)
441 {
442     long g;
443 
444     //  MYREAL copysum=0;
445     long pop, i;
446     MYREAL expapg;  //, summ;
447     //MYREAL summm = 0.0;
448     long msta, msto;
449     MYREAL waitmig;
450     MYREAL pointmig;
451     tarchive_fmt *tl;
452     MYREAL *part;   //, df = 0.0;
453     //  MYREAL value ;
454     MYREAL theta1 = param[0];
455     MYREAL *geo = nr->world->data->geo;
456     long rep;
457     MYREAL theta_ratio = theta1 / x;
458     part = (MYREAL *) mymalloc (sizeof (MYREAL) * nr->partsize);
459     memset (nr->parts, 0, sizeof (MYREAL) * nr->partsize);
460     for (rep = nr->repstart; rep < nr->repstop; rep++)
461     {
462         for (g = 0; g < nr->atl[rep][locus].T; g++)
463         {
464             if (nr->apg[rep][locus][g] > -40.)
465             {
466                 tl = nr->atl[rep][locus].tl;
467                 //copies = (g > 0) ? tl[g].copies : tl[g].copies - 1;
468                 expapg = /*copies  */ EXP (nr->apg[rep][locus][g]);
469                 if (expapg == 0.0)
470                     continue;
471                 //      summm += expapg;
472                 // population 1 is treated differently
473                 part[0] = 0.0;
474                 msta = nr->mstart[0];
475                 msto = nr->mend[0];
476                 for (i = msta; i < msto; i++)
477                 {
478                     part[0] += -geo[i] * param[i] * tl[g].km[0] / x +
479                                tl[g].mindex[i] / theta1;
480                     if (param[i] > 0.)
481                         part[i] = ((tl[g].mindex[i] / param[i])
482                                    - geo[i] * tl[g].km[0] * theta_ratio);
483                 }
484                 // all other populations
485                 for (pop = 1; pop < nr->numpop; pop++)
486                 {
487                     part[pop] = -tl[g].p[pop] / param[pop] +
488                                 tl[g].kt[pop] * theta1 / (x * param[pop] * param[pop]);
489 
490                     msta = nr->mstart[pop];
491                     msto = nr->mend[pop];
492                     //              z = 0;
493                     waitmig = 0;
494                     pointmig = 0;
495                     for (i = msta; i < msto; i++)
496                     {
497                         waitmig += geo[i] * param[i];
498                         pointmig += tl[g].mindex[i];
499                         if (param[i] > 0.)
500                             part[i] = ((tl[g].mindex[i] / param[i])
501                                        - geo[i] * tl[g].km[pop] * theta_ratio);
502                     }
503                     part[0] += pointmig / theta1 - waitmig * tl[g].km[pop] / x +
504                                tl[g].p[pop] / theta1 - tl[g].kt[pop] / (param[pop] * x);
505                 }
506                 for (pop = 0; pop < nr->world->numpop2; pop++)
507                     nr->parts[pop] += expapg * part[pop];
508             }
509         }
510     }
511     for (pop = 0; pop < nr->world->numpop2; pop++)
512         nr->parts[pop] /= nr->PGC[locus];
513     // printf("nr->PGC[%li]=%f == %f\n",locus,nr->PGC[locus],summm);
514     myfree(part);
515 }
516