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 (&param, 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 (&param, 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 (&param, 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 (&param, *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