1 /** \file combroyden.c */
2 /*------------------------------------------------------
3 Maximum likelihood estimation
4 of migration rate and effectice population size
5 using a Metropolis-Hastings Monte Carlo algorithm
6 -------------------------------------------------------
7 M A X I M I Z E R R O U T I N E S
8
9 - calculates single locus ML
10 - multilocus ML with constant mutation rate
11 - multilocus ML with gamma deviated mutation rate
12
13 using Broyden-Fletcher-Goldfarb-Shanno method
14
15
16 Peter Beerli 1996-1998, Seattle
17 beerli@fsu.edu
18
19
20 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
21 Copyright 2003-2004 Peter Beerli, Tallahassee FL
22
23 Permission is hereby granted, free of charge, to any person obtaining a copy
24 of this software and associated documentation files (the "Software"), to deal
25 in the Software without restriction, including without limitation the rights
26 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
27 of the Software, and to permit persons to whom the Software is furnished to do
28 so, subject to the following conditions:
29
30 The above copyright notice and this permission notice shall be included in all copies
31 or substantial portions of the Software.
32
33 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
34 INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
35 PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
36 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
37 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
38 OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
39
40 $Id: combroyden2.c 2169 2013-08-24 19:02:04Z beerli $
41
42 -------------------------------------------------------*/
43
44 #include "migration.h"
45 #include "sighandler.h"
46 #include "broyden.h"
47 #include "joint-chains.h"
48 #include "world.h"
49 #include "integrate.h"
50 #include "migrate_mpi.h"
51 #include "reporter.h"
52 #include "tools.h"
53 #include "laguerre.h"
54 #include "gammalike.h"
55 #ifdef DMALLOC_FUNC_CHECK
56 #include <dmalloc.h>
57 #endif
58
59
60 //if using slownet then this should be activated
61 // otherwise I will use a define in defintions.h that that
62 // replaces all occurrences with CALCLIKE with calc_loci_like
63 // for speed reasons we don't want any further indirection
64 // for which we have no use.
65 #ifdef SLOWNET
66 MYREAL (*calc_like) (helper_fmt *, MYREAL *, MYREAL *);
67 void (*calc_gradient) (nr_fmt *, helper_fmt *, MYREAL *);
68 void (*setup_param0) (world_fmt *, nr_fmt *, long, long, long, long, long,
69 boolean);
70 #endif
71
72 // external definitions -------------------------------
73 extern void debug_plot (helper_fmt * helper);
74 extern int myID;
75
76 /* prototypes ------------------------------------------- */
77 /* public function used in main.c */
78 long estimateParameter (long rep, long G, world_fmt * world,
79 option_fmt * options, MYREAL **dd, long chain,
80 char type, long kind, long repkind);
81
82 /* calculates likelihood */
83 MYREAL calc_loci_like (helper_fmt * helper, MYREAL *param, MYREAL *lparam);
84
85 /* first derivatives for ML with constant mutation rate */
86 void simple_loci_derivatives (MYREAL *d, nr_fmt * nr,
87 timearchive_fmt ** tyme, long locus);
88 /* progress-reporter */
89 void print_menu_finalestimate (boolean progress, long locus, long loci);
90
91 /* start parameter for multilocus estimator, creates simple average */
92 void calc_meanparam (world_fmt * world, long n, long repstart, long repstop);
93
94 void expected_param(MYREAL *param, world_fmt *world);
95
96 /* multilocus first derivatives with or without gamma distrib. mut. rate */
97 void combine_gradient (nr_fmt * nr, helper_fmt * helper, MYREAL *gxv);
98 #ifdef SLOWNET
99 void slow_net_combine_gradient (nr_fmt * nr, helper_fmt * helper,
100 MYREAL *gxv);
101 #endif
102 /* this helps to use the simple single locus material */
103 void copy_and_clear_d (nr_fmt * nr);
104 void add_back_d (nr_fmt * nr);
105 /* rescales parameter with theta1 for gamma calculation */
106 void set_gamma_param (MYREAL *paramn, MYREAL *paramo,
107 MYREAL *lparamn, MYREAL *lparamo,
108 MYREAL theta, nr_fmt * nr);
109 /* calculates norm, this creates a stopping criteria */
110 MYREAL norm (MYREAL *d, long size);
111 /* used in multilocus line search */
112 MYREAL psilo (MYREAL lamda, helper_fmt * helper, MYREAL *param,
113 MYREAL *lparam);
114 /* used in integral for gamma distr. mutation rate */
115 MYREAL phi (MYREAL t, helper_fmt * b);
116
117 /* calculates new multilocus parameter values */
118 void calc_loci_param (nr_fmt * nr, MYREAL *lparam, MYREAL *olparam,
119 MYREAL *dv, MYREAL lamda, long nnn);
120 /* used to reset the BFSG stuff */
121 void reset_work (MYREAL **work, long nnn);
122
123 /* linesearch of BFGS
124 and helper functions for this linesearch */
125 MYREAL line_search (MYREAL *dv, MYREAL *gxv,
126 MYREAL (*fpsi) (MYREAL lamda, helper_fmt * helper,
127 MYREAL *param, MYREAL *lparam),
128 helper_fmt * thishelper, MYREAL oldL, MYREAL *tempparam,
129 MYREAL *templparam);
130
131 void replace_with (int mode, MYREAL *a, MYREAL *b, MYREAL *c, MYREAL m,
132 MYREAL *la, MYREAL *lb, MYREAL *lc, MYREAL ll);
133 MYREAL psil (MYREAL (*fpsi) (MYREAL lamda, helper_fmt * helper), MYREAL a,
134 MYREAL b, helper_fmt * helper, MYREAL *gxv, MYREAL *dv,
135 MYREAL oldL, MYREAL *val1, MYREAL *val2);
136 MYREAL quadratic_lamda (MYREAL lamda, MYREAL a, MYREAL b, MYREAL c, MYREAL la,
137 MYREAL lb, MYREAL lc);
138
139
140 void print_reset_big_param (FILE * file, MYREAL average, MYREAL maxtheta,
141 long pop);
142
143 long do_profiles (world_fmt * world, nr_fmt * nr, MYREAL *likes,
144 MYREAL *normd, long kind, long rep, long repkind);
145
146 void correct_values (world_fmt * world, nr_fmt * nr);
147
148 void which_calc_like (long repkind);
149
150 void set_replicates (world_fmt * world, long repkind, long rep,
151 long *repstart, long *repstop);
152
153 void prepare_broyden (long kind, world_fmt * world, boolean * multilocus);
154
155 long multilocus_gamma_adjust (boolean multilocus, boolean boolgamma,
156 world_fmt * world, long repstart, long repstop);
157 void profile_setup (long profilenum, long *profiles, MYREAL *param,
158 MYREAL *value, long *nnn);
159
160 void setup_parameter0_standard (world_fmt * world, nr_fmt * nr, long repkind,
161 long repstart, long repstop, long loci,
162 long kind, boolean multilocus);
163
164 void report_broyden (MYREAL newL, MYREAL normd, long trials,
165 boolean boolgamma, MYREAL theta1, MYREAL alpha,
166 FILE * logfile);
167 boolean guard_broyden (MYREAL newL, MYREAL oldL, long trials, MYREAL normd,
168 MYREAL *normd20);
169
170 void fill_helper (helper_fmt * helper, MYREAL *param, MYREAL *lparam,
171 world_fmt * world, nr_fmt * nr);
172
173 /* maximizer routine */
174 long maximize (MYREAL **thisparam, world_fmt * world, nr_fmt * nr,
175 MYREAL **hess,
176 long analystype,
177 long repkind);
178
179 void set_both_delta (nr_fmt * nr,
180 MYREAL *delta, MYREAL *den, MYREAL *deo,
181 MYREAL *gama, MYREAL *gxv, MYREAL *gxv0, long size);
182 void set_expxv (MYREAL *expa, MYREAL *a, long size);
183 void set_expparam (MYREAL *expa, MYREAL *a, long size);
184 void set_logparam (MYREAL *loga, MYREAL *a, long size);
185 void finish_broyden (MYREAL newL, MYREAL normd, long trials,
186 world_fmt * world, nr_fmt * nr, MYREAL *expxv, long nn,
187 long analystype, long repkind);
188 void fill_profile_index (nr_fmt * nr);
189
190
191 /* Functions ++++++++++++++++++++++++++++++++++++++++++++++++*/
192 /* public functions---------------------------------------------- */
193 /* driver for maximizing single locus, multilocus functions */
194 long
estimateParameter(long rep,long G,world_fmt * world,option_fmt * options,MYREAL ** dd,long chain,char type,long kind,long repkind)195 estimateParameter (long rep, long G, world_fmt * world, option_fmt * options,
196 MYREAL **dd, long chain, char type, long kind,
197 long repkind)
198 //timearchive_fmt * tyme
199 {
200 long trials = 0;
201 MYREAL normd = 0.;
202 MYREAL likes = 0.;
203 MYREAL *param;
204 #ifndef MPI
205 long repstart, repstop;
206 #endif
207
208 boolean multilocus;
209 nr_fmt *nr;
210 nr = (nr_fmt *) mycalloc (1, sizeof (nr_fmt));
211 if (kind == MULTILOCUS)
212 {
213 world->locus = world->loci;
214 print_menu_finalestimate (world->options->progress, world->loci,
215 world->loci);
216 }
217
218 create_nr (nr, world, G, 0L, world->locus, repkind, rep);
219 set_replicates (world, world->repkind, world->rep, &nr->repstart,
220 &nr->repstop);
221 // fprintf(stdout,"%i> DEBUG sumfile: world-rep=%li world->repkind=%li nr->repstart=%li nr->repstop=%li\n",myID, world->rep, world->repkind, nr->repstart, nr->repstop);
222 prepare_broyden (kind, world, &multilocus);
223 #ifndef INTEGRATEDLIKE
224 // if(world->options->integrated_like)
225 SETUPPARAM0 (world, nr, world->repkind, nr->repstart, nr->repstop,
226 world->loci, kind, multilocus);
227 #else
228 SETUPPARAM0 (world, nr, world->repkind, nr->repstart, nr->repstop,
229 world->loci, kind, multilocus);
230 #endif
231 doublevec1d (¶m, nr->partsize + 1);
232
233 #ifdef MPI
234 if (myID == MASTER)
235 memcpy (param, world->param0, sizeof (MYREAL) * nr->numpop2);
236 #else
237 if (kind == MULTILOCUS || repkind != SINGLECHAIN)
238 {
239 set_replicates (world, world->repkind, world->rep, &repstart, &repstop);
240 calc_meanparam (world, world->numpop2, repstart, repstop);
241 memcpy (param, world->param0, sizeof (MYREAL) * nr->numpop2);
242 }
243 // {
244 // expected_param(param,world);
245 // }
246 #endif
247 memcpy (param, world->param0, sizeof (MYREAL) * nr->numpop2);
248
249 #ifdef LONGSUM
250
251 memcpy (param+nr->numpop2+
252 ((world->options->gamma && kind==MULTILOCUS) ? 1 : 0) , world->flucrates, sizeof (MYREAL) * nr->numpop*3);
253 #endif
254
255 if (strpbrk (world->options->custm2, "0c"))
256 {
257 if (kind == SINGLELOCUS
258 && world->options->custm2[world->numpop2] == 'c')
259 trials = maximize (¶m, world, nr, dd, kind, repkind);
260 else
261 trials = do_profiles (world, nr, &likes, &normd, kind, rep, repkind);
262 }
263 else
264 {
265 trials = maximize (¶m, world, nr, dd, kind, repkind);
266 }
267 if (kind == SINGLELOCUS && world->options->progress)
268 {
269 print_menu_chain (type, chain, world->atl[rep][world->locus].T,
270 world, options, rep);
271 FPRINTF (stdout, " Maximization cycles needed: %li\n Norm of first derivatives: %f\n",
272 trials, nr->normd);
273 if (world->options->logfile)
274 {
275 FPRINTF (world->options->logfile,
276 " Maximization cycles needed: %li\n Norm of first derivatives: %f\n", trials, nr->normd);
277 }
278 }
279 #ifdef MPI
280 if (myID==MASTER && world->options->plotnow)
281 #else
282 if (kind == MULTILOCUS && world->options->plotnow && world->loci > 1)
283 #endif
284 {
285 if (((world->options->replicate && repkind != SINGLECHAIN) ||
286 (!world->options->replicate && repkind == SINGLECHAIN)) && (world->numpop > 1))
287 #ifdef MPI
288 {
289 if(world->loci==1)
290 create_plot (world, world->plane[0], nr, world->loci, MULTILOCUSPLOT);
291 else
292 create_plot (world, world->plane[world->loci], nr, world->loci, MULTILOCUSPLOT);
293 }
294 #else
295 create_plot (world, world->plane[world->loci], nr, world->loci, MULTILOCUSPLOT);
296 #endif
297 }
298 #ifdef LONGSUM
299 memcpy (world->flucrates, param+nr->numpop2+
300 ((world->options->gamma && kind==MULTILOCUS) ? 1 : 0), sizeof (MYREAL) * nr->numpop*3);
301 #endif
302
303 world->trials = trials;
304 world->normd = nr->normd;
305 myfree(param);
306 destroy_nr (nr, world);
307 return trials;
308 }
309
310
311 /* calculates the likelihood over all loci for the new parameter set */
312 MYREAL
calc_loci_like(helper_fmt * helper,MYREAL * param,MYREAL * lparam)313 calc_loci_like (helper_fmt * helper, MYREAL *param, MYREAL *lparam)
314 {
315 nr_fmt *nr = helper->nr;
316 world_fmt *world = helper->nr->world;
317 MYREAL logres = 0;
318 #ifndef MPI
319 long locus;
320 MYREAL *mu_rates = world->options->mu_rates;
321 MYREAL * locusweight = world->data->locusweight;
322 #endif
323 if (helper->multilocus)
324 {
325 #ifdef MPI
326 /* must be MASTER node!
327 mutation rate is constant log(L(loci)) = Sum[log(L(locus))) */
328 // printf("%i> before mpi_likelihood_master\n",myID);
329 logres = mpi_likelihood_master (param, lparam, world, nr,
330 helper, world->who);
331 #else
332
333 if (!helper->boolgamma)
334 {
335 for (locus = 0; locus < world->loci; locus++)
336 {
337 if (nr->skiploci[locus])
338 continue;
339 nr->locilikes[locus] =
340 calc_locus_like (nr, param, lparam, locus) + mu_rates[locus];
341 logres += locusweight[locus] * nr->locilikes[locus]; //invariant loci treatment
342 }
343 }
344 else
345 {
346 logres = gamma_loci_like (helper, param, lparam, helper->weight);
347 }
348 #endif /* NOT MPI */
349
350 }
351 else /* ->singlelocus */
352 logres = calc_locus_like (nr, param, lparam, nr->world->locus);
353 nr->llike = logres;
354 return nr->llike;
355 }
356
357 #ifdef SLOWNET
358 /* calculates the likelihood over all loci for the new parameter set */
359 MYREAL
slow_net_calc_loci_like(helper_fmt * helper,MYREAL * param,MYREAL * lparam)360 slow_net_calc_loci_like (helper_fmt * helper, MYREAL *param, MYREAL *lparam)
361 {
362 long locus;
363 nr_fmt *nr = helper->nr;
364 world_fmt *world = helper->nr->world;
365 MYREAL *mu_rates = world->options->mu_rates;
366 MYREAL logres = 0;
367 MYREAL *locusweight = helper->nr->world->data->locusweight; //invariant loci treatment
368 if (helper->multilocus)
369 {
370 if (!helper->boolgamma)
371 {
372 for (locus = 0; locus < world->loci; locus++)
373 {
374 if (nr->skiploci[locus])
375 continue;
376 nr->locilikes[locus] =
377 calc_locus_like (nr, param, lparam, locus) + mu_rates[locus];
378 logres += locusweight[locus] * nr->locilikes[locus];
379 }
380 }
381 else
382 {
383 logres = gamma_loci_like (helper, param, lparam, helper->weight);
384 }
385 }
386 else /* ->singlelocus */
387 logres = calc_locus_like (nr, param, lparam, nr->world->locus);
388 nr->llike = logres;
389 return nr->llike;
390 }
391 #endif /*SLOWNET*/
392 /*
393 MYREAL
394 phi (MYREAL t, helper_fmt *b)
395 {
396 MYREAL ll, alpha, beta;
397 long locus;
398 MYREAL weight;
399 helper_fmt *helper;
400 nr_fmt *nr;
401 timearchive_fmt **atl;
402 helper = (helper_fmt *) b;
403 locus = (long) helper->locus;
404 weight = (MYREAL) helper->weight;
405 atl = (timearchive_fmt **) helper->atl;
406 nr = (nr_fmt *) helper->nr;
407
408 alpha = nr->oparam[nr->numpop2];
409 beta = nr->oparam[0] / alpha;
410 set_gamma_param (nr->param, nr->oparam, nr->lparam, nr->olparam, t, nr);
411 ll = calc_locus_like (nr, nr->param, nr->lparam, locus);
412 //#ifdef LAGUERRE
413 ll = -t / beta + (alpha - 1.) * LOG (t) + ll + weight;
414 //#else
415 // ll = EXP (-t / beta + (alpha - 1.) * LOG (t) + ll + weight);
416 //#endif
417 return ll;
418 }
419
420 */
421 /* private functions---------------------------------------------- */
422 /* derivatives */
423 void
combine_gradient(nr_fmt * nr,helper_fmt * helper,MYREAL * gxv)424 combine_gradient (nr_fmt * nr, helper_fmt * helper, MYREAL *gxv)
425 {
426 #ifndef MPI
427 long locus;
428 #endif
429
430 timearchive_fmt **tyme = NULL;
431 // MYREAL save_a, denom, fla, fha, la, ha;
432 long nnn = nr->partsize - nr->profilenum;
433 memset (gxv, 0, sizeof (MYREAL) * nnn);
434 if (helper->analystype == SINGLELOCUS)
435 {
436 simple_loci_derivatives (gxv, nr, tyme, nr->world->locus);
437 }
438 else
439 {
440 #ifdef MPI
441 if (myID == MASTER)
442 {
443 mpi_gradient_master (nr, nr->world, nr->world->who);
444 memcpy (gxv, nr->d, sizeof (MYREAL) * nr->partsize);
445 }
446 else //SLAVE
447 {
448 error ("not allowed to execute mpi_gradient_worker()\n");
449 }
450 #else
451 if (!helper->boolgamma)
452 {
453 for (locus = 0; locus < nr->world->loci; locus++)
454 {
455 if (nr->skiploci[locus])
456 continue;
457 memset (nr->d, 0, sizeof (MYREAL) * nnn);
458 simple_loci_derivatives (nr->d, nr, tyme, locus);
459 add_vector (gxv, nr->d, nnn);
460 }
461 }
462 else
463 {
464 //#ifdef LAGUERRE
465 //if(nr->normd<5.)
466 // gamma_loci_difference(helper);
467 //else
468 gamma_loci_derivative (helper);
469 memcpy (gxv, nr->d, sizeof (MYREAL) * nr->partsize);
470 /* save_a = helper->xv[nr->partsize-1];
471 la= save_a - save_a/100.;
472 helper->xv[nr->partsize-1]=la;
473 helper->expxv[nr->partsize-1]=EXP(la);
474 denom = -LGAMMA(helper->expxv[nr->partsize-1]) - (helper->xv[0]-helper->xv[nr->partsize-1]) * helper->xv[nr->partsize-1];
475 fla = gamma_loci_like (helper, helper->expxv, helper->xv, denom);
476 ha= save_a + save_a/100.;
477 helper->xv[nr->partsize-1]=ha;
478 helper->expxv[nr->partsize-1]=EXP(ha);
479 denom = -LGAMMA(helper->expxv[nr->partsize-1]) - (helper->xv[0]-helper->xv[nr->partsize-1]) * helper->xv[nr->partsize-1];
480 fha = gamma_loci_like (helper, helper->expxv, helper->xv, denom);
481 nr->d[nr->partsize-1] = (fla - fha)/(la-ha);
482 memcpy (gxv, nr->d, sizeof (MYREAL) * nr->partsize);*/
483 //#else
484 // dt (gxv, helper);
485 //endif ****LAGUERRE****
486 }
487 #endif /*MPI*/
488
489 }
490 force_symmetric_d (gxv, (long) nr->world->options->migration_model, nr, nnn);
491 grad2loggrad (helper->expxv, nr->indeks, gxv, nnn, nr->profilenum);
492 }
493
494 #ifdef SLOWNET
495 void
slow_net_combine_gradient(nr_fmt * nr,helper_fmt * helper,MYREAL * gxv)496 slow_net_combine_gradient (nr_fmt * nr, helper_fmt * helper, MYREAL *gxv)
497 {
498 long locus;
499 timearchive_fmt **tyme = NULL;
500
501 long nnn = nr->partsize - nr->profilenum;
502 memset (gxv, 0, sizeof (MYREAL) * nnn);
503 if (helper->analystype == SINGLELOCUS)
504 {
505 simple_loci_derivatives (gxv, nr, tyme, nr->world->locus);
506 force_symmetric_d (gxv, nr->world->options->migration_model, nr, nnn);
507 grad2loggrad (helper->expxv, nr->indeks, gxv, nnn, nr->profilenum);
508 }
509 else
510 {
511 if (!helper->boolgamma)
512 {
513 for (locus = 0; locus < nr->world->loci; locus++)
514 {
515 if (nr->skiploci[locus])
516 continue;
517 memset (nr->d, 0, sizeof (MYREAL) * nnn);
518 simple_loci_derivatives (nr->d, nr, tyme, locus);
519 add_vector (gxv, nr->d, nnn);
520 }
521 force_symmetric_d (gxv, nr->world->options->migration_model, nr,
522 nnn);
523 grad2loggrad (helper->expxv, nr->indeks, gxv, nnn, nr->profilenum);
524 }
525 else
526 {
527 //#ifdef LAGUERRE
528 gamma_loci_derivative (helper);
529 memcpy (gxv, nr->d, sizeof (MYREAL) * nr->partsize);
530 //#else
531 //
532 //dt (gxv, helper);
533 //#endif
534
535 force_symmetric_d (gxv, nr->world->options->migration_model,
536 nr, nnn);
537 grad2loggrad (helper->expxv, nr->indeks, gxv, nnn, nr->profilenum);
538 }
539 }
540 }
541 #endif
542
543 void
simple_loci_derivatives(MYREAL * d,nr_fmt * nr,timearchive_fmt ** tyme,long locus)544 simple_loci_derivatives (MYREAL *d, nr_fmt * nr, timearchive_fmt ** tyme,
545 long locus)
546 {
547 long g;
548 #ifndef LONGSUM
549
550 gradient(d, nr, locus);
551 #else /*LONGSUM*/
552
553 gradient_longsum (d, nr, locus);
554 #endif /*LONGSUM*/
555
556 for (g = 0; g < nr->partsize - nr->profilenum; g++)
557 d[g] /= nr->PGC[locus];
558 }
559
560 void
set_gamma_param(MYREAL * param,MYREAL * oparam,MYREAL * lparam,MYREAL * olparam,MYREAL theta,nr_fmt * nr)561 set_gamma_param (MYREAL *param, MYREAL *oparam,
562 MYREAL *lparam, MYREAL *olparam, MYREAL theta, nr_fmt * nr)
563 {
564 long i;
565 MYREAL logtheta;
566 param[0] = theta;
567 lparam[0] = logtheta = LOG (theta);
568 for (i = 1; i < nr->numpop; i++)
569 {
570 lparam[i] = olparam[i] + logtheta - olparam[0];
571 param[i] = EXP (lparam[i]);
572 }
573 for (i = nr->numpop; i < nr->numpop2; i++)
574 {
575 lparam[i] = olparam[i] + olparam[0] - logtheta;
576 param[i] = EXP (lparam[i]);
577 }
578 lparam[i] = olparam[i];
579 param[i] = EXP (olparam[i]);
580 }
581
582
583 void
copy_and_clear_d(nr_fmt * nr)584 copy_and_clear_d (nr_fmt * nr)
585 {
586 memcpy (nr->od, nr->d, sizeof (MYREAL) * nr->partsize);
587 memset (nr->d, 0, sizeof (MYREAL) * nr->partsize);
588 }
589
590 void
add_back_d(nr_fmt * nr)591 add_back_d (nr_fmt * nr)
592 {
593 long pop;
594 for (pop = 0; pop < nr->partsize; pop++)
595 {
596 nr->d[pop] += nr->od[pop];
597 }
598 }
599
600
601 void
print_menu_finalestimate(boolean progress,long locus,long loci)602 print_menu_finalestimate (boolean progress, long locus, long loci)
603 {
604 static char nowstr[STRSIZE];
605 if (progress)
606 {
607 get_time (nowstr, "%H:%M:%S");
608 if (locus < loci)
609 fprintf (stdout, "%s Parameter estimation for locus %li\n", nowstr,
610 locus);
611 else
612 fprintf (stdout, "%s Parameter estimation over all loci\n", nowstr);
613 }
614 }
615
616 void
print_reset_big_param(FILE * file,MYREAL average,MYREAL maxtheta,long pop)617 print_reset_big_param (FILE * file, MYREAL average, MYREAL maxtheta, long pop)
618 {
619 static boolean checkparam[MAXPOP];
620 static boolean done = FALSE;
621 static char nowstr[STRSIZE];
622 if (!done)
623 {
624 done = TRUE;
625 memset (&checkparam, 0, sizeof (boolean) * MAXPOP);
626 }
627 if (!checkparam[pop])
628 {
629 get_time (nowstr, "%H:%M:%S");
630 FPRINTF (file,
631 "%s Theta for population %li was reset to the average\n",
632 nowstr, pop);
633 FPRINTF (file, " of %f, it exceeded the maximum of %f \n",
634 average, maxtheta);
635 }
636 }
637
638 void
calc_meanparam(world_fmt * world,long n,long repstart,long repstop)639 calc_meanparam (world_fmt * world, long n, long repstart, long repstop)
640 {
641 long r, pop, locus, loci = world->loci;
642
643
644 memset (world->param0, 0, sizeof (MYREAL) * n);
645 for (r = repstart; r < repstop; r++)
646 {
647 for (locus = 0; locus < loci; locus++)
648 {
649 for (pop = 0; pop < n; pop++)
650 {
651 if (!world->data->skiploci[locus])
652 world->param0[pop] += world->atl[r][locus].param[pop];
653 }
654 }
655 }
656 for (pop = 0; pop < n; pop++)
657 {
658 world->param0[pop] /= (loci - world->skipped) * (repstop - repstart);
659 }
660 }
661
662
663 MYREAL
norm(MYREAL * d,long size)664 norm (MYREAL *d, long size)
665 {
666 int i;
667 MYREAL summ = 0.;
668 for (i = 0; i < (int) size; i++)
669 {
670 summ += d[i] * d[i];
671 }
672 return sqrt (summ);
673 }
674
675
676 MYREAL
psilo(MYREAL lamda,helper_fmt * helper,MYREAL * param,MYREAL * lparam)677 psilo (MYREAL lamda, helper_fmt * helper, MYREAL *param, MYREAL *lparam)
678 {
679 MYREAL like;
680 nr_fmt *nr = helper->nr;
681 calc_loci_param (nr, lparam, helper->xv, helper->dv, lamda, nr->partsize);
682 set_expparam (param, lparam, nr->partsize);
683 fill_helper (helper, param, lparam, nr->world, nr);
684 like = CALCLIKE (helper, param, lparam);
685 //printf("%f ",like);
686 return -like;
687 }
688
689 void
calc_loci_param(nr_fmt * nr,MYREAL * lparam,MYREAL * olparam,MYREAL * dv,MYREAL lamda,long nnn)690 calc_loci_param (nr_fmt * nr, MYREAL *lparam, MYREAL *olparam, MYREAL *dv,
691 MYREAL lamda, long nnn)
692 {
693 long i, ii, z = 0;
694 if (nr->profilenum == 0)
695 {
696 for (i = 0; i < nnn; i++)
697 lparam[i] = (MAX (-30., MIN (olparam[i] - lamda * dv[i], 30.)));
698 }
699 else
700 {
701 for (i = 0; i < nr->partsize - nr->profilenum; i++)
702 {
703 ii = nr->indeks[z++];
704 lparam[ii] = (MAX (-30., MIN (olparam[ii] - lamda * dv[i], 30.)));
705 }
706 }
707 param_all_adjust (lparam, nr);
708 //#ifdef LAGUERRE
709
710 if (nr->world->locus >= nr->world->loci && nr->world->options->gamma)
711 {
712 if (lparam[nr->numpop2] > 9.903487553)
713 {
714 lparam[nr->numpop2] = 9.903487553;
715 }
716 initgammacat (nr->categs, EXP (lparam[nr->numpop2]),1.0 /*EXP (lparam[0])*/,
717 nr->rate, nr->probcat);
718 //}
719 }
720 //#endif
721
722 }
723
724
725 void
reset_work(MYREAL ** work,long nnn)726 reset_work (MYREAL **work, long nnn)
727 {
728 long i, j;
729 for (i = nnn; i < nnn + 5; i++)
730 {
731 for (j = 0; j < nnn + 1; j++)
732 work[i][j] = 0.0;
733 }
734 }
735
736
737 MYREAL
calc_d_g(MYREAL * d,MYREAL * g,long n)738 calc_d_g (MYREAL *d, MYREAL *g, long n)
739 {
740 long i;
741 MYREAL sum = 0.;
742 for (i = 0; i < n; i++)
743 sum += d[i] * g[i];
744 return sum;
745 }
746
747
748 #define STARTLAMBDA 1
749 #define M_ 0
750
751 /*
752 returns LAMBDA (how far to jump in the search direction)
753 needed in the the multilocus broyden-fletcher-shanno maximizer
754
755 DV search direction
756 GXV first derivatives of function to minimize
757 FUNC the function to minimize
758 NR all additional variables needed to calculate FUNC
759 OLDL -log(likelihood) of FUNC with LAMBDA=0 (previously calculated)
760
761 needs functions: psil, replace_with, quadratic_lamda
762 */
763 MYREAL
line_search(MYREAL * dv,MYREAL * gxv,MYREAL (* fpsi)(MYREAL lamda,helper_fmt * helper,MYREAL * param,MYREAL * lparam),helper_fmt * thishelper,MYREAL oldL,MYREAL * tempparam,MYREAL * templparam)764 line_search (MYREAL *dv, MYREAL *gxv,
765 MYREAL (*fpsi) (MYREAL lamda, helper_fmt * helper, MYREAL *param,
766 MYREAL *lparam), helper_fmt * thishelper,
767 MYREAL oldL, MYREAL *tempparam, MYREAL *templparam)
768 {
769 //#ifdef MYREAL == float
770 //const MYREAL eps = FLT_EPSILON ;
771 //#else
772 const MYREAL eps = DBL_EPSILON ;
773 //#endif
774
775 long trials = 0;
776 MYREAL locallamda = STARTLAMBDA, ql;
777 MYREAL a, b, c;
778 MYREAL psiabc;
779 MYREAL ll, la, lb, lc;
780 MYREAL m = 1.0;
781 MYREAL denom;
782 MYREAL nom;
783 helper_fmt helper;
784 helper = *thishelper;
785 // partcopy_helper(thishelper, &helper);
786 a = 0.0;
787 b = 1.0;
788
789 la = -oldL;
790 lb = -oldL; //(*fpsi) (b, &helper, tempparam, templparam);
791 // ll = (*fpsi) (m, &helper, tempparam, templparam);
792 //if(ll > la)
793 // c = -1.;
794 //else
795 // c = -2.;
796 c = calc_d_g (dv, gxv, helper.nr->partsize - helper.nr->profilenum);
797 lc = (*fpsi) (c, &helper, tempparam, templparam);
798 while (trials++ < NTRIALS)
799 {
800 // printf ("*%3li> %f %f %f / %f %f %f\n", trials, a, b, c, la, lb, lc);
801 denom = (2. * (b * la - c * la - a * lb + c * lb + a * lc - b * lc));
802 if((fabs(denom)-eps)>eps)
803 {
804 nom = (c * c * (lb - la) + b * b * (la - lc) +
805 a * a * (lc - lb));
806 }
807 else
808 return 1.;
809 locallamda = nom / denom;
810 psiabc = ((la - lb) / (-a + b) + (la - lc) / (a - c)) / (-b + c);
811 if ((psiabc <= 0.0) || (locallamda >= m))
812 {
813 if ((a == m || b == m || c == m))
814 return m;
815 else
816 {
817 ll = (*fpsi) (m, &helper, tempparam, templparam);
818 // replace_with(M_,&a,&b,&c,m,&la,&lb,&lc,ll);
819 replace_with ((int) m, &a, &b, &c, m, &la, &lb, &lc, ll);
820 continue;
821 }
822 }
823 ll = (*fpsi) (locallamda, &helper, tempparam, templparam);
824 ql = quadratic_lamda (locallamda, a, b, c, la, lb, lc);
825 if ((fabs (ll - MYMIN3 (la, lb, lc)) <= BIGEPSILON)
826 || (fabs (ll - ql) <= BIGEPSILON))
827 return locallamda;
828 else
829 {
830 /* if(((a<b<c) || (c<b<a)) && (lb < MIN(la,lc)))
831 return locallamda;
832 if(((b<a<c) || (c<a<b)) && (la < MIN(lb,lc)))
833 return locallamda;
834 if(((a<c<b) || (b<c<a)) && (lc < MIN(la,lb)))
835 return locallamda; */
836 replace_with ((int) STARTLAMBDA, &a, &b, &c, locallamda, &la, &lb, &lc,
837 ll);
838 m = MYMAX3 (a, b, c);
839 }
840 }
841
842 return locallamda;
843 }
844
845 void
replace_with(int mode,MYREAL * a,MYREAL * b,MYREAL * c,MYREAL m,MYREAL * la,MYREAL * lb,MYREAL * lc,MYREAL ll)846 replace_with (int mode, MYREAL *a, MYREAL *b, MYREAL *c, MYREAL m, MYREAL *la,
847 MYREAL *lb, MYREAL *lc, MYREAL ll)
848 {
849 MYREAL ma, mb, mc;
850 if (mode == STARTLAMBDA)
851 {
852 ma = *la;
853 mb = *lb;
854 mc = *lc;
855 }
856 else
857 {
858 ma = ll - *la;
859 mb = ll - *lb;
860 mc = ll - *lc;
861 }
862 if (ma > mb)
863 {
864 if (ma > mc)
865 {
866 *a = m;
867 *la = ll;
868 }
869 else
870 {
871 *c = m;
872 *lc = ll;
873 }
874 }
875 else
876 {
877 if (mb > mc)
878 {
879 *b = m;
880 *lb = ll;
881 }
882 else
883 {
884 *c = m;
885 *lc = ll;
886 }
887 }
888 }
889
890 /* not used
891 MYREAL
892 psil (MYREAL (*fpsi) (MYREAL lamda, helper_fmt * helper), MYREAL a, MYREAL b,
893 helper_fmt * helper, MYREAL *gxv, MYREAL *dv, MYREAL oldL, MYREAL *val1,
894 MYREAL *val2)
895 {
896 long i;
897 long nn = helper->nr->partsize;
898 MYREAL value = 0.0;
899
900 if (a == 0.0 && b == 0)
901 {
902 for (i = 0; i < nn; i++)
903 {
904 value += dv[i]; // gxv[i] ;
905 }
906 *val1 = -oldL;
907 *val2 = -oldL;
908 return value;
909 }
910 else
911 {
912 if (a == 0)
913 value = (((*val2) = (*fpsi) (b, helper)) - ((*val1) = -oldL)) / b;
914 else
915 {
916 if (b == 0)
917 value =
918 (((*val2) = -oldL) - ((*val1) = (*fpsi) (a, helper))) / (-a);
919 else
920 value =
921 (((*val2) = (*fpsi) (b, helper)) - ((*val1) =
922 (*fpsi) (a,
923 helper))) / (b -
924 a);
925 }
926 return value;
927 }
928 }
929 */
930
931 MYREAL
quadratic_lamda(MYREAL lamda,MYREAL a,MYREAL b,MYREAL c,MYREAL la,MYREAL lb,MYREAL lc)932 quadratic_lamda (MYREAL lamda, MYREAL a, MYREAL b, MYREAL c, MYREAL la,
933 MYREAL lb, MYREAL lc)
934 {
935 MYREAL alpha, beta, gama;
936 MYREAL aa, bb, cc;
937 aa = a * a;
938 bb = b * b;
939 cc = c * c;
940 alpha =
941 ((c - b) * la + (a - c) * lb +
942 (b - a) * lc) / ((b - a) * (c - a) * (c - b));
943 beta =
944 ((cc - bb) * la + (aa - cc) * lb +
945 (bb - aa) * lc) / ((b - a) * (b - c) * (c - a));
946 gama =
947 (b * cc * la - bb * c * la + aa * c * lb - a * cc * lb - aa * b * lc +
948 a * bb * lc) / ((b - a) * (c - a) * (c - b));
949
950
951 return alpha * lamda * lamda + beta * lamda + gama;
952
953 }
954
955 #define TAU1 2.6180339887
956 #define TAU 1.6180339887
957 #define BSTART -2.360679775
golden_section(MYREAL (* fpsi)(MYREAL lamda,helper_fmt * helper,MYREAL * param,MYREAL * lparam),helper_fmt * thishelper,MYREAL * tempparam,MYREAL * templparam)958 MYREAL golden_section(MYREAL (*fpsi) (MYREAL lamda, helper_fmt * helper, MYREAL *param,
959 MYREAL *lparam), helper_fmt * thishelper,
960 MYREAL *tempparam, MYREAL *templparam)
961 {
962 MYREAL a = -10.;
963 MYREAL b = BSTART;
964 MYREAL c = 10.;
965 MYREAL bl;
966 MYREAL d=1.;
967 MYREAL dl = -MYREAL_MAX;
968 MYREAL ll = 0.;
969 helper_fmt helper;
970 helper = *thishelper;
971 bl = (*fpsi) (b, &helper, tempparam, templparam);
972 while(fabs(ll - dl)>EPSILON)
973 {
974 ll = dl;
975 if((c-b) > (b-a))
976 {
977 d = b + (c - b)/ TAU1;
978 dl = (*fpsi) (d, &helper, tempparam, templparam);
979 if(dl > bl)
980 c = d;
981 else
982 {
983 a = b;
984 b = d;
985 bl = dl;
986 }
987 }
988 else
989 {
990 d = a + (b - a)/ TAU1;
991 dl = (*fpsi) (d, &helper, tempparam, templparam);
992 if(dl > bl)
993 a = d;
994 else
995 {
996 c = b;
997 b = d;
998 bl = dl;
999 }
1000 }
1001 }
1002 return d;
1003 }
1004
1005 long
do_profiles(world_fmt * world,nr_fmt * nr,MYREAL * likes,MYREAL * normd,long kind,long rep,long repkind)1006 do_profiles (world_fmt * world, nr_fmt * nr, MYREAL *likes,
1007 MYREAL *normd, long kind, long rep, long repkind)
1008 {
1009 MYREAL **hess;
1010 char *p;
1011 long z = 0;
1012 long trials;
1013 long i = 0;
1014 p = world->options->custm2;
1015 while (*p != '\0')
1016 {
1017 if (*p == '0' || *p == 'c')
1018 {
1019 nr->profiles[z] = i;
1020 nr->values[z] = world->param0[i];
1021 z++;
1022 }
1023 p++;
1024 i++;
1025 }
1026 nr->profilenum = z;
1027 if(kind!=SINGLELOCUS)
1028 world->locus = world->loci;
1029 // use a private copy of the hessian so not to clobber the
1030 // version we use to print out the -fisher information
1031 z = world->numpop2 +
1032 (world->locus == world->loci ?
1033 world->options->gamma : 0);
1034 doublevec2d(&hess,z,z);
1035 trials = maximize (&world->param0, world, nr, hess, kind, repkind);
1036 myfree(hess[0]);
1037 myfree(hess);
1038 return trials;
1039 }
1040
1041 void
correct_values(world_fmt * world,nr_fmt * nr)1042 correct_values (world_fmt * world, nr_fmt * nr)
1043 {
1044 long pop, i, ii, z;
1045 MYREAL ssum;
1046 long nsum;
1047 long elem = (world->options->gamma ? nr->numpop2 : nr->partsize);
1048 for (pop = 0; pop < world->numpop; pop++)
1049 {
1050 if (world->param0[pop] >= BIGGEST_THETA)
1051 {
1052 ssum = 0;
1053 nsum = 0;
1054 for (i = 0; i < world->numpop; i++)
1055 {
1056 if (world->param0[i] <= BIGGEST_THETA)
1057 {
1058 ssum += world->param0[i];
1059 nsum++;
1060 }
1061 }
1062 if (nsum > 0 && ssum > 0.0)
1063 world->param0[pop] = ssum / (MYREAL) nsum;
1064 else
1065 {
1066 if (!(ssum > 0.0))
1067 world->param0[pop] = SMALLEST_THETA;
1068 else
1069 world->param0[pop] = BIGGEST_THETA;
1070 }
1071 if (world->options->writelog)
1072 print_reset_big_param (world->options->logfile,
1073 world->param0[pop], BIGGEST_THETA, pop);
1074 // print_reset_big_param (stdout, world->param0[pop], BIGGEST_THETA,
1075 // pop);
1076 }
1077 }
1078 z = 0;
1079 for (i = 0; i < elem - nr->profilenum; i++)
1080 {
1081 ii = (nr->profilenum > 0) ? nr->indeks[z++] : i;
1082 if (world->param0[ii] >= BIGGEST_MIGRATION)
1083 world->param0[ii] = BIGGEST_MIGRATION;
1084 if (world->param0[ii] < SMALLEST_MIGRATION)
1085 world->param0[ii] = SMALLEST_MIGRATION;
1086 }
1087 }
1088
1089 /// decide which conditional likelihood calculator to use
1090 void
which_calc_like(long kind)1091 which_calc_like (long kind)
1092 {
1093 switch (kind)
1094 {
1095 case PROFILE:
1096 #ifdef SLOWNET
1097
1098 calc_like =
1099 (MYREAL (*)(helper_fmt *, MYREAL *, MYREAL *))
1100 slow_net_calc_loci_like;
1101 calc_gradient =
1102 (void (*)(nr_fmt *, helper_fmt *, MYREAL *))
1103 slow_net_combine_gradient;
1104 setup_param0 =
1105 (void (*)
1106 (world_fmt *, nr_fmt *, long, long, long, long, long,
1107 boolean)) setup_parameter0_slowmpi;
1108 #endif
1109
1110 break;
1111 default:
1112 #ifdef SLOWNET
1113
1114 calc_like =
1115 (MYREAL (*)(helper_fmt *, MYREAL *, MYREAL *)) calc_loci_like;
1116 calc_gradient =
1117 (void (*)(nr_fmt *, helper_fmt *, MYREAL *)) combine_gradient;
1118 setup_param0 =
1119 (void (*)
1120 (world_fmt *, nr_fmt *, long, long, long, long, long,
1121 boolean)) setup_parameter0_mpi;
1122 #endif
1123
1124 break;
1125 }
1126 }
1127
1128 /// calculates the start and end of a replication step, sets repstop = repstart+1
1129 void
set_replicates(world_fmt * world,long repkind,long rep,long * repstart,long * repstop)1130 set_replicates (world_fmt * world, long repkind, long rep, long *repstart,
1131 long *repstop)
1132 {
1133 if (world->options->replicate)
1134 {
1135 *repstart = (repkind == SINGLECHAIN) ? rep : 0;
1136 *repstop = (repkind == SINGLECHAIN) ? rep + 1 : world->repstop;
1137 }
1138 else
1139 {
1140 *repstart = 0;
1141 *repstop = 1;
1142 }
1143 }
1144
1145
1146
1147 void
prepare_broyden(long kind,world_fmt * world,boolean * multilocus)1148 prepare_broyden (long kind, world_fmt * world, boolean * multilocus)
1149 {
1150 if (kind == SINGLELOCUS || (kind == PROFILE && world->loci == 1))
1151 {
1152 if (world->loci == 1)
1153 world->locus = 0;
1154 *multilocus = FALSE;
1155 }
1156 else
1157 {
1158 *multilocus = TRUE;
1159 }
1160 }
1161
1162
1163 long
multilocus_gamma_adjust(boolean multilocus,boolean boolgamma,world_fmt * world,long repstart,long repstop)1164 multilocus_gamma_adjust (boolean multilocus, boolean boolgamma,
1165 world_fmt * world, long repstart, long repstop)
1166 {
1167 long nn;
1168 if (multilocus)
1169 {
1170 nn = boolgamma ? world->numpop2 + 1 : world->numpop2;
1171 if (boolgamma)
1172 {
1173 world->param0 =
1174 (MYREAL *) myrealloc (world->param0,
1175 sizeof (MYREAL) * (nn > 0 ? nn : 1));
1176 world->param0[world->numpop2] = world->options->alphavalue;
1177
1178 }
1179 }
1180 else
1181 nn = world->numpop2;
1182 return nn;
1183 }
1184
1185 /* profiler setup--------------------------
1186 needs the passed in values generated in profiles setup (profile.c)
1187 */
1188 void
profile_setup(long profilenum,long * profiles,MYREAL * param,MYREAL * value,long * nnn)1189 profile_setup (long profilenum, long *profiles, MYREAL *param, MYREAL *value,
1190 long *nnn)
1191 {
1192 long i;
1193 if (profilenum > 0)
1194 {
1195 *nnn -= profilenum;
1196 for (i = 0; i < profilenum; i++)
1197 param[profiles[i]] = LOG (value[i]);
1198 }
1199 if (*nnn == 0)
1200 *nnn = 1;
1201 }
1202
1203
1204 void
setup_parameter0_standard(world_fmt * world,nr_fmt * nr,long repkind,long repstart,long repstop,long loci,long kind,boolean multilocus)1205 setup_parameter0_standard (world_fmt * world, nr_fmt * nr, long repkind,
1206 long repstart, long repstop, long loci, long kind,
1207 boolean multilocus)
1208 {
1209 long locus, r;
1210 if (multilocus)
1211 {
1212 for (locus = 0; locus < loci; locus++)
1213 {
1214 if (repkind == SINGLECHAIN)
1215 {
1216 for (r = repstart; r < repstop; r++)
1217 create_apg0 (nr->apg0[r][locus], nr,
1218 &world->atl[r][locus], locus);
1219 }
1220 else
1221 {
1222 if (kind != PROFILE)
1223 {
1224 for (r = repstart; r < repstop; r++)
1225 create_apg0 (nr->apg0[r][locus], nr,
1226 &world->atl[r][locus], locus);
1227 interpolate_like (nr, locus);
1228 }
1229 // else
1230 // {
1231 for (r = repstart; r < repstop; r++)
1232 create_multiapg0 (nr->apg0[r][locus], nr, r, locus);
1233 // }
1234 }
1235 }
1236 }
1237 else //single locus
1238 {
1239 if (repkind == SINGLECHAIN)
1240 {
1241 for (r = repstart; r < repstop; r++)
1242 create_apg0 (nr->apg0[r][world->locus], nr,
1243 &world->atl[r][world->locus], world->locus);
1244 }
1245 else
1246 {
1247 if (kind != PROFILE)
1248 {
1249 for (r = repstart; r < repstop; r++)
1250 create_apg0 (nr->apg0[r][world->locus], nr,
1251 &world->atl[r][world->locus], world->locus);
1252 interpolate_like (nr, world->locus);
1253 }
1254 for (r = repstart; r < repstop; r++)
1255 create_multiapg0 (nr->apg0[r][world->locus], nr, r, world->locus);
1256 }
1257 }
1258 }
1259
1260 void
report_broyden(MYREAL newL,MYREAL normd,long trials,boolean boolgamma,MYREAL theta1,MYREAL alpha,FILE * logfile)1261 report_broyden (MYREAL newL, MYREAL normd, long trials,
1262 boolean boolgamma, MYREAL theta1, MYREAL alpha,
1263 FILE * logfile)
1264 {
1265 if (boolgamma)
1266 FPRINTF (stdout, "%i:%li> Log(L)=%f Norm=%f Alpha=%f Theta_1=%f\n",
1267 myID, trials, newL, normd, alpha, theta1);
1268 else
1269 FPRINTF (stdout, "%i:%li> Log(L)=%f Norm=%f Theta_1=%f\n",
1270 myID, trials, newL, normd, theta1);
1271 if (logfile != NULL)
1272 {
1273 if (boolgamma)
1274 FPRINTF (logfile, "%i:%li> Log(L)=%f Norm=%f Alpha=%f Theta_1=%f\n",
1275 myID, trials, newL, normd, alpha, theta1);
1276 else
1277 FPRINTF (logfile, "%i:%li> Log(L)=%f Norm=%f Theta_1=%f\n",
1278 myID, trials, newL, normd, theta1);
1279 }
1280 }
1281
1282 boolean
guard_broyden(MYREAL newL,MYREAL oldL,long trials,MYREAL normd,MYREAL * normd20)1283 guard_broyden (MYREAL newL, MYREAL oldL, long trials, MYREAL normd,
1284 MYREAL *normd20)
1285 {
1286 if ((trials + 1) % 20 == 0)
1287 {
1288 if (fabs (normd - *normd20) < EPSILON)
1289 return TRUE;
1290 *normd20 = normd;
1291 }
1292 return FALSE;
1293 }
1294
1295 /* for profiling: filling of index array */
1296 void
fill_profile_index(nr_fmt * nr)1297 fill_profile_index (nr_fmt * nr)
1298 {
1299 long i, ii;
1300 long profilenum = nr->profilenum;
1301 if ((profilenum > 0) && profilenum < nr->partsize)
1302 {
1303 for (i = 0, ii = 0; i < nr->partsize; i++)
1304 {
1305 if (!find (i, nr->profiles, profilenum))
1306 nr->indeks[ii++] = i;
1307 }
1308 }
1309 }
1310
1311 // maximizer
1312 // - MPI
1313
1314 // likelihood
1315 // - MPI
1316 // - Gamma
1317 // - custom migration matrix aware
1318 // - profile aware
1319 // - lrt aware [special case of cumstom migration]
1320
1321 // maximizer sceleton
1322 // analystype = (SINGLELOCUS, MULTILOCUS, PROFILE)
1323 long
maximize(MYREAL ** thisparam,world_fmt * world,nr_fmt * nr,MYREAL ** hess,long analystype,long repkind)1324 maximize (MYREAL **thisparam, world_fmt * world, nr_fmt * nr, MYREAL **hess, long analystype,
1325 long repkind)
1326 {
1327 //#ifdef MYREAL == float
1328 //const MYREAL eps = FLT_EPSILON ;
1329 //const MYREAL numbermax = FLT_MAX ;
1330 //#else
1331 const MYREAL eps = DBL_EPSILON ;
1332 const MYREAL numbermax = MYREAL_MAX ;
1333 //#endif
1334 boolean stop;
1335 long repstart, repstop;
1336 // long Gmax;
1337 MYREAL two = 2.;
1338 MYREAL newL;
1339 MYREAL oldL;
1340 MYREAL lamda;
1341 long trials;
1342 MYREAL normd = 1000000;
1343 MYREAL normd20 = MYREAL_MAX;
1344 worldoption_fmt *wopt = world->options;
1345 long nnn = nr->partsize; // keeps track of profile or standard
1346 long nn; // keeps track of the full parameter array
1347
1348 MYREAL *param;
1349 // MYREAL **hess;
1350 MYREAL *delta;
1351 MYREAL *gama;
1352 MYREAL *dv = NULL;
1353 MYREAL *gxv;
1354 MYREAL *gxv0;
1355 MYREAL *xv0;
1356 MYREAL *xv;
1357 MYREAL *expxv;
1358 MYREAL *temp;
1359 //MYREAL *tempparam; //just pointer into temp
1360 //MYREAL *templparam; //just pointer into temp
1361 helper_fmt helper;
1362
1363 setdoublevec1d (¶m, *thisparam, nnn);
1364 // long reset_z = 0;
1365 helper.boolgamma = world->locus == world->loci ? wopt->gamma : FALSE;
1366 if (helper.boolgamma)
1367 param[nnn - 1] = world->options->alphavalue;
1368 //doublevec2d (&hess, nnn, nnn);
1369 reset_hess (hess, nnn);
1370 doublevec1d (&delta, nnn);
1371 doublevec1d (&gama, nnn);
1372 doublevec1d (&xv, nnn + 1);
1373 doublevec1d (&temp, 2 * (nnn + 1));
1374 //xcode tempparam = temp;
1375 //xcode templparam = temp + nnn + 1;
1376 helper.xv = xv;
1377 // expxv = xv + nnn;
1378 //memcpy(expxv,param,sizeof(MYREAL)*nnn);
1379 doublevec1d (&gxv, nnn);
1380 doublevec1d (&gxv0, nnn);
1381 setdoublevec1d (&expxv, param, nnn);
1382 helper.expxv = expxv;
1383 set_logparam (xv, expxv, nnn);
1384 setdoublevec1d (&xv0, xv, nnn);
1385 helper.analystype = analystype;
1386 set_replicates (world, world->repkind, world->rep, &repstart, &repstop);
1387 nr->repstart = repstart;
1388 nr->repstop = repstop;
1389 nr->normd = normd;
1390 prepare_broyden (analystype, world, &helper.multilocus);
1391 nnn = multilocus_gamma_adjust (helper.multilocus, helper.boolgamma,
1392 world, repstart, repstop);
1393 #ifdef LONGSUM
1394
1395 nnn=nr->partsize;
1396 #endif
1397
1398 nn = nnn;
1399 //if(helper.boolgamma)
1400 //{
1401 // save_a = nr->world->options->custm[nr->world->numpop2];
1402 // world->options->custm[nr->world->numpop2]='c';
1403 // world->options->custm2[nr->world->numpop2]='c';
1404 //}
1405 profile_setup (nr->profilenum, nr->profiles, xv, nr->values, &nnn);
1406 fill_profile_index (nr);
1407 helper.lamda = 0.;
1408 helper.dv = dv;
1409 calc_loci_param (nr, xv, xv0, gxv, 0., nnn);
1410 set_expparam (expxv, xv, nn);
1411 //SETUPPARAM0 (world, nr, world->repkind, repstart, repstop,
1412 // world->loci, analystype, helper.multilocus);
1413 // FPRINTF(stdout,"analystype is %li\n",analystype);
1414
1415 fill_helper (&helper, expxv, xv, world, nr);
1416 oldL = CALCLIKE (&helper, expxv, xv);
1417 // oldL= -MYREAL_MAX;
1418 fill_helper (&helper, expxv, xv, world, nr);
1419 CALCGRADIENT (nr, &helper, gxv);
1420 memcpy (gxv0, gxv, sizeof (MYREAL) * nnn);
1421 memcpy (xv0, xv, sizeof (MYREAL) * nn);
1422 setdoublevec1d (&dv, gxv, nnn);
1423 helper.dv = dv;
1424 /* Maximization cycle ---------------------------------------------*/
1425 for (trials = 0; trials < NTRIALS; trials++)
1426 {
1427 #ifdef __MWERKS__
1428 eventloop ();
1429 #endif
1430
1431 newL = -MYREAL_MAX;
1432 // if (helper.boolgamma)
1433 // lamda = golden_section (psilo, &helper, tempparam, templparam);
1434 //line_search (dv, gxv, psilo, &helper, oldL, tempparam, templparam);
1435 // else
1436 // lamda = golden_section (psilo, &helper, tempparam, templparam);
1437 //line_search (dv, gxv, psilo, &helper, oldL, tempparam, templparam);
1438 helper.lamda = lamda = 1.;
1439 while ((newL < oldL || MYISNAN (newL) || newL <= -numbermax) &&
1440 fabs (lamda) > 10. * eps)
1441 {
1442 calc_loci_param (nr, xv, xv0, dv, lamda, nn);
1443 set_expparam (expxv, xv, nn);
1444 fill_helper (&helper, expxv, xv, world, nr);
1445 newL = CALCLIKE (&helper, expxv, xv);
1446 lamda /= two;
1447 }
1448 nr->normd = normd = norm (gxv, nnn);
1449 //#ifdef LAGUERRE
1450
1451 if (nr->world->locus >= nr->world->loci && nr->world->options->gamma)
1452 {
1453 // initgammacat (nr->categs, expxv[nr->numpop2], expxv[0],
1454 // nr->rate, nr->probcat);
1455 if (fabs (oldL - newL) < 1e-9)
1456 lamda = 0.;
1457 }
1458 //#endif
1459 if (world->options->verbose && newL > oldL)
1460 {
1461 if(helper.boolgamma)
1462 report_broyden (newL, normd, trials, helper.boolgamma,
1463 expxv[0], expxv[nr->numpop2],
1464 world->options->logfile);
1465 else
1466 report_broyden (newL, normd, trials, helper.boolgamma,
1467 expxv[0], -1.0,
1468 world->options->logfile);
1469
1470 }
1471 stop = guard_broyden (newL, oldL, trials, normd, &normd20);
1472 if (normd < EPSILON || stop)
1473 break; // stopping criteria
1474 else if (oldL >= newL)
1475 lamda = 0;
1476 else
1477 oldL = newL;
1478 /* reset sec deri. matrix if lamda goes to 0 and retry */
1479 // if(reset_z > 10)
1480 //{
1481 // lamda=0;
1482 // reset_z = 0;
1483 //}
1484 if (fabs (lamda) <= eps * 10.)
1485 {
1486 reset_hess (hess, nnn);
1487 two = -two;
1488 memcpy (dv, gxv, sizeof (MYREAL) * nnn);
1489 continue;
1490 }
1491 //swap(gxv0,gxv);
1492 memcpy (gxv0, gxv, sizeof (MYREAL) * nnn);
1493 //calc_loci_param (nr, xv, xv0, dv, lamda * two, nn);
1494 //set_expparam(expxv,xv,nn);
1495 fill_helper (&helper, expxv, xv, world, nr);
1496 CALCGRADIENT (nr, &helper, gxv);
1497 set_both_delta (nr, delta, xv, xv0, gama, gxv, gxv0, nnn);
1498 calc_hessian (hess, nnn, delta, gama);
1499 calc_dv (dv, hess, gxv, nnn);
1500 memcpy (xv0, xv, sizeof (MYREAL) * nn);
1501 }
1502 /* end maximizer cycle ------------------------------------*/
1503 /*if(helper.boolgamma)
1504 {
1505 world->options->custm[nr->world->numpop2]=save_a;
1506 world->options->custm2[nr->world->numpop2]=save_a;
1507 }*/
1508
1509 finish_broyden (newL, normd, trials, world, nr, expxv, nn,
1510 analystype, repkind);
1511 memcpy ((*thisparam), expxv, sizeof (MYREAL) * nn);
1512 myfree(param);
1513 //myfree(hess[0]);
1514 //myfree(hess);
1515 myfree(delta);
1516 myfree(gama);
1517 myfree(dv);
1518 myfree(gxv);
1519 myfree(gxv0);
1520 myfree(xv0);
1521 myfree(xv);
1522 myfree(expxv);
1523 myfree(temp);
1524 fflush (stderr);
1525 return trials;
1526 }
1527
1528 void
finish_broyden(MYREAL newL,MYREAL normd,long trials,world_fmt * world,nr_fmt * nr,MYREAL * expxv,long nn,long analystype,long repkind)1529 finish_broyden (MYREAL newL, MYREAL normd, long trials,
1530 world_fmt * world, nr_fmt * nr, MYREAL *expxv, long nn,
1531 long analystype, long repkind)
1532 {
1533 long loci = world->loci;
1534
1535 timearchive_fmt **tyme = world->atl;
1536 memcpy (world->param0, expxv, sizeof (MYREAL) * nn);
1537 //correct_values (world, nr); // resets huge thetas to average
1538
1539 //if (world->options->verbose)
1540 // FPRINTF (stdout, "\n");
1541
1542 //#ifdef INTEGRATEDLIKE
1543 //if(analystype==MULTILOCUS && world->loci==1)
1544 // {
1545 // analystype=SINGLELOCUS;
1546 // world->locus=0;
1547 // }
1548 //#endif
1549
1550 switch (analystype)
1551 {
1552 case MULTILOCUS:
1553 // tyme[0][loci].param = (MYREAL *) myrealloc (tyme[0][loci].param,
1554 // sizeof (MYREAL) * nn);
1555 tyme[0][loci].param_like = newL;
1556 tyme[0][loci].normd = normd;
1557 tyme[0][loci].trials = trials;
1558 memmove (tyme[0][loci].param, expxv, sizeof (MYREAL) * nn);
1559 nr->normd = normd;
1560 //convergence_check (world, world->options->verbose);
1561 break;
1562 case SINGLELOCUS:
1563 world->param_like = newL;
1564 if (repkind != SINGLECHAIN)
1565 {
1566 // printf("%i> single locus combine replicates into rep=%li\n",myID,world->repstop);
1567 tyme[world->repstop][world->locus].param_like = newL;
1568 tyme[world->repstop][world->locus].normd = normd;
1569 tyme[world->repstop][world->locus].trials = trials;
1570 memmove (tyme[world->repstop][world->locus].param, expxv,
1571 sizeof (MYREAL) * nn);
1572 }
1573 else
1574 {
1575 //printf("%i> more loci put replicate into rep=%li\n",myID,world->rep);
1576 tyme[world->rep][world->locus].param_like = newL;
1577 tyme[world->rep][world->locus].normd = normd;
1578 tyme[world->rep][world->locus].trials = trials;
1579 memmove (tyme[world->rep][world->locus].param, expxv,
1580 sizeof (MYREAL) * nn);
1581 }
1582 if(world->options->gelman)
1583 convergence_check (world, world->options->verbose);
1584 if (world->loci == 1)
1585 {
1586 if ((!world->options->gelman &&
1587 world->param_like < world->options->lcepsilon &&
1588 world->options->plotnow && !world->options->simulation) ||
1589 (world->options->plotnow && world->options->gelman &&
1590 world->convergence->gelmanmeanmaxR[1] < GELMAN_MYSTIC_VALUE))
1591 {
1592 if ((world->options->replicate && repkind != SINGLECHAIN) ||
1593 (!world->options->replicate && repkind == SINGLECHAIN))
1594 create_plot (world, world->plane[world->locus], nr,
1595 nr->atl[world->rep][0].T, SINGLELOCUSPLOT);
1596 }
1597 }
1598 nr->normd = normd;
1599 break;
1600 case PROFILE:
1601 // memmove (nr->profparam, world->param0, sizeof (MYREAL) * nr->partsize);
1602 nr->llike = newL;
1603 nr->normd = normd;
1604 }
1605
1606 if (world->options->verbose && world->repkind == SINGLECHAIN)
1607 {
1608 if (analystype == SINGLELOCUS)
1609 print_contribution (nr, nr->atl, nr->atl[world->rep][world->locus].T);
1610 }
1611 }
1612
1613
1614 void
set_both_delta(nr_fmt * nr,MYREAL * delta,MYREAL * den,MYREAL * deo,MYREAL * gama,MYREAL * gxv,MYREAL * gxv0,long size)1615 set_both_delta (nr_fmt * nr,
1616 MYREAL *delta, MYREAL *den, MYREAL *deo,
1617 MYREAL *gama, MYREAL *gxv, MYREAL *gxv0, long size)
1618 {
1619 long i, ii;
1620 if ((nr->profilenum > 0) && nr->profilenum < nr->partsize)
1621 {
1622 for (i = 0; i < size; i++)
1623 {
1624 ii = nr->indeks[i];
1625 delta[i] = den[ii] - deo[ii];
1626 gama[i] = gxv[i] - gxv0[i];
1627 }
1628 }
1629 else
1630 {
1631 for (i = 0; i < size; i++)
1632 {
1633 delta[i] = den[i] - deo[i];
1634 gama[i] = gxv[i] - gxv0[i];
1635 }
1636 }
1637 }
1638
1639 /// calculates the EXP of an array \todo move this to tools.c
1640 void
set_expparam(MYREAL * expa,MYREAL * a,long size)1641 set_expparam (MYREAL *expa, MYREAL *a, long size)
1642 {
1643 long i;
1644 for (i = 0; i < size; i++)
1645 expa[i] = EXP (a[i]);
1646 }
1647
1648 /// calculates the log of an array \todo move this to the tools.c
1649 void
set_logparam(MYREAL * loga,MYREAL * a,long size)1650 set_logparam (MYREAL *loga, MYREAL *a, long size)
1651 {
1652 long i;
1653 for (i = 0; i < size; i++)
1654 loga[i] = LOG (a[i]);
1655 }
1656
1657 /// fills the helper structure with parameters and pointers to world and nr
1658 void
fill_helper(helper_fmt * helper,MYREAL * param,MYREAL * lparam,world_fmt * world,nr_fmt * nr)1659 fill_helper (helper_fmt * helper, MYREAL *param, MYREAL *lparam,
1660 world_fmt * world, nr_fmt * nr)
1661 {
1662 MYREAL alpha = 0;
1663 MYREAL theta1 = param[0];
1664 MYREAL denom = 0;
1665 helper->boolgamma = world->locus < world->loci ?
1666 (helper->analystype==PROFILE ? world->options->gamma : FALSE) : world->options->gamma;
1667 if (helper->boolgamma)
1668 {
1669 alpha = param[nr->numpop2];
1670 if (alpha <= 0.0)
1671 error ("no no! alpha=0\n");
1672 denom = LGAMMA (alpha) + LOG (theta1 / alpha) * alpha;
1673 // denom = LGAMMA (alpha) + LOG (1. / alpha) * alpha;
1674 }
1675 if (nr->lparam != lparam)
1676 memcpy (nr->lparam, lparam, sizeof (MYREAL) * nr->partsize);
1677 if (nr->param != param)
1678 memcpy (nr->param, param, sizeof (MYREAL) * nr->partsize);
1679
1680 helper->nr = nr;
1681 helper->ll = -MYREAL_MAX;
1682 helper->weight = denom;
1683 helper->atl = world->atl;
1684 }
1685
expected_param(MYREAL * param,world_fmt * world)1686 void expected_param(MYREAL *param, world_fmt *world)
1687 {
1688 long i, tree;
1689 //long j;
1690 tarchive_fmt *thistree;
1691 long T = world->atl[world->replicate][world->locus].T;
1692 memset(param,0,sizeof(MYREAL) * world->numpop2);
1693 for(tree=0; tree < T; tree++)
1694 {
1695 thistree = &world->atl[world->replicate][world->locus].tl[tree];
1696 for(i=0;i<world->numpop; i++)
1697 {
1698 param[i] += thistree->kt[i] / 2.;
1699 //xcode j or i? for(j=world->mstart[i];i<world->mend[i]; i++)
1700 for(i=world->mstart[i];i<world->mend[i]; i++)
1701 param[i] += 1./thistree->km[i];
1702 }
1703 }
1704 for(i=0;i<world->numpop2; i++)
1705 {
1706 param[i] /= T;
1707 }
1708 }
1709