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 (¶m, 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