1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate and effectice population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6 P R O F I L E L I K E L I H O O D R O U T I N E S
7
8 Peter Beerli 1997, Seattle
9 beerli@fsu.edu
10
11 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
12 Copyright 2003-2004 Peter Beerli, Tallahassee FL
13
14 Permission is hereby granted, free of charge, to any person obtaining
15 a copy of this software and associated documentation files (the
16 "Software"), to deal in the Software without restriction, including
17 without limitation the rights to use, copy, modify, merge, publish,
18 distribute, sublicense, and/or sell copies of the Software, and to
19 permit persons to whom the Software is furnished to do so, subject
20 to the following conditions:
21
22 The above copyright notice and this permission notice shall be
23 included in all copies or substantial portions of the Software.
24
25 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
26 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
27 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
28 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
29 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
30 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
31 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
32
33 $Id: profile.c 2158 2013-04-29 01:56:20Z beerli $
34
35 -------------------------------------------------------*/
36 /*! \file profile.c
37 Calculation of profile likelihoods and printing of profile likelihood tables
38 */
39
40 #include "migration.h"
41
42 #include "world.h"
43 #include "laguerre.h"
44 #include "tools.h"
45 #include "broyden.h"
46 #include "combroyden.h"
47 #include "spline.h"
48 #include "joint-chains.h"
49 #include "sighandler.h"
50
51 #include "migrate_mpi.h"
52
53 #include "profile.h"
54 #ifdef PRETTY
55 #include "pretty.h"
56 #endif
57
58 #ifdef DMALLOC_FUNC_CHECK
59 #include <dmalloc.h>
60 #endif
61
62 extern int myID;
63
64 long find_profilelike (MYREAL testlike, MYREAL prob, long whichprob, MYREAL minparam,
65 MYREAL maxparam, MYREAL **testparam, MYREAL *llike,
66 long which, world_fmt * world, boolean QD, boolean QD2,
67 MYREAL *mlparam, MYREAL *normd, nr_fmt * nr);
68 void calc_profile_likelihood (char method, long which, MYREAL *likes,
69 MYREAL **param, world_fmt * world, long nn,
70 nr_fmt * nr);
71 void print_profile_likelihood (long which, world_fmt * world, long *gmaxptr);
72 void prepare_header (char *var, long which, world_fmt * world);
73 void print_profile_title (world_fmt * world);
74 long print_profile_table (char method, long which, char *var, MYREAL likes[],
75 MYREAL **param, world_fmt * world);
76 long print_profile_percentile (world_fmt * world);
77 void sprint_percentile_header (char *outfile, boolean first);
78 void allocate_profile_percentiles (world_fmt * world);
79 void destroy_profile_percentiles (world_fmt * world);
80 void print_menu_profile (long which, long nn);
81 MYREAL interpolate (MYREAL xl, MYREAL xc, MYREAL xh, MYREAL low,
82 MYREAL center, MYREAL high, MYREAL testlike);
83
84 long warp (long ii);
85
86 void prognose_profile_end (time_t starttime, long numpop2, long nn);
87
88
89 void profile_max_percentiles (long which, MYREAL *likes, MYREAL **param,
90 world_fmt * world, long nn, nr_fmt * nr);
91
92 void profile_max_spline (char method, long which, MYREAL *likes,
93 MYREAL **param, world_fmt * world, nr_fmt * nr,
94 long nn);
95
96 int calc_spline (MYREAL *param, MYREAL *like, long nn, long *constr,
97 MYREAL *diff, MYREAL *diff2, long *diagn, MYREAL *work,
98 long nwork);
99 void setup_spline (spline_fmt * spline);
100
101 void destroy_spline (spline_fmt * spline);
102
103 void prepare_spline_nodes (long n, long which, MYREAL **rawparam,
104 MYREAL *rawlikes, long *indeks, MYREAL *param,
105 MYREAL *likes);
106
107 void prepare_spline_first (MYREAL *param, MYREAL *likes);
108
109 void prepare_spline_last (long n, MYREAL *param, MYREAL *likes);
110
111 void calc_spline_diff (MYREAL *diff, MYREAL *param, MYREAL *likes, long n);
112
113 MYREAL recalc_profilelike (MYREAL testparam, MYREAL *param, MYREAL *like,
114 long nn, MYREAL *diff, MYREAL *diff2, MYREAL *work,
115 long nwork);
116
117
118 long set_df (long which, char *custm2, long numpop, long numpop2);
119
120
121 long set_profile_param (MYREAL *param, long which, long *allwhich,
122 MYREAL xval, MYREAL *xvals, world_fmt * world);
123 void addinvar (MYREAL *param, long *fixed, long *allwhich, MYREAL *xvals,
124 long start, long len);
125
126 long sprint_nice_param (MYREAL parameter, MYREAL bottom, MYREAL top, boolean failed,
127 char *file);
128
129 #ifdef PRIVATE
130 void master_gridder(world_fmt *world, long *gmaxptr);
131 void calc_grid(FILE *alldim, world_fmt * world, MYREAL *testparam, MYREAL *ltestparam, MYREAL *mlparam, nr_fmt *nr, long which);
132 #endif
133
134
135
136
choose_profile_parameters(world_fmt * world,long * profilelist)137 long choose_profile_parameters (world_fmt * world, long *profilelist)
138 {
139 static boolean onetheta = FALSE;
140 static boolean onemig = FALSE;
141 long i, j;
142 long z=0;
143 long which;
144 for(which=0;which<world->numpop2;which++)
145 {
146 switch (world->options->custm2[which])
147 {
148 case 'C':
149 case 'c':
150 case '0':
151 break;
152 case 'M':
153 case 'm':
154 if (which < world->numpop && onetheta)
155 continue;
156 if (which < world->numpop)
157 {
158 profilelist[z++] = which;
159 onetheta = TRUE;
160 }
161 if (which >= world->numpop && onemig)
162 continue;
163 if (which >= world->numpop)
164 {
165 profilelist[z++] = which;
166 onemig = TRUE;
167 }
168 break;
169 case 'S': // Nm is symmetrical
170 j = (which - world->numpop) % (world->numpop - 1);
171 i = (which - world->numpop) / (world->numpop - 1);
172 if(world->options->usem)
173 {
174 profilelist[z++] = which;
175 continue;
176 }
177 else
178 {
179 if (i <= j)
180 {
181 profilelist[z++] = which;
182 continue;
183 }
184 }
185 break;
186 case 's': // M is symmetric
187 j = (which - world->numpop) % (world->numpop - 1);
188 i = (which - world->numpop) / (world->numpop - 1);
189 if(!world->options->usem)
190 {
191 profilelist[z++] = which;
192 }
193 else
194 {
195 if (i <= j)
196 {
197 profilelist[z++] = which;
198 }
199 }
200 break;
201 default:
202 profilelist[z++] = which;
203 break;
204 }
205 }
206 if(world->options->gamma)
207 profilelist[z++]=world->numpop2;
208 return z;
209 }
210
211 boolean
print_profile_likelihood_driver(long which,world_fmt * world,long * gmaxptr)212 print_profile_likelihood_driver (long which, world_fmt * world, long *gmaxptr)
213 {
214 static boolean onetheta = FALSE;
215 static boolean onemig = FALSE;
216 long i, j;
217 //world->options->verbose = FALSE;
218 world->locus = 0;
219 if (which == world->numpop2 && world->options->gamma) //gamma deviated mutation rate
220 print_profile_likelihood (which, world, gmaxptr);
221 else
222 {
223 if(which>world->numpop2)
224 return 0;
225 switch (world->options->custm2[which])
226 {
227 case 'C':
228 case 'c':
229 case '0':
230 return 0;
231 case 'M':
232 case 'm':
233 if (which < world->numpop && onetheta)
234 return 0;
235 if (which < world->numpop)
236 {
237 print_profile_likelihood (which, world, gmaxptr);
238 onetheta = TRUE;
239 }
240 if (which >= world->numpop && onemig)
241 return 0;
242 if (which >= world->numpop)
243 {
244 print_profile_likelihood (which, world, gmaxptr);
245 onemig = TRUE;
246 }
247 return 1;
248 case 'S': // Nm is symmetrical
249 j = (which - world->numpop) % (world->numpop - 1);
250 i = (which - world->numpop) / (world->numpop - 1);
251 if(world->options->usem)
252 {
253 print_profile_likelihood (which, world, gmaxptr);
254 return 1;
255 }
256 else
257 {
258 if (i <= j)
259 {
260 print_profile_likelihood (which, world, gmaxptr);
261 return 1;
262 }
263 }
264 return 0;
265 case 's': // M is symmetrical
266 j = (which - world->numpop) % (world->numpop - 1);
267 i = (which - world->numpop) / (world->numpop - 1);
268 if(!world->options->usem)
269 {
270 print_profile_likelihood (which, world, gmaxptr);
271 return 1;
272 }
273 else
274 {
275 if (i <= j)
276 {
277 print_profile_likelihood (which, world, gmaxptr);
278 return 1;
279 }
280 }
281 return 0;
282 default:
283 print_profile_likelihood (which, world, gmaxptr);
284 }
285 }
286 return 1;
287 }
288
289 void
print_profile_likelihood(long which,world_fmt * world,long * gmaxptr)290 print_profile_likelihood (long which, world_fmt * world, long *gmaxptr)
291 {
292 long i, ii;
293 long bufsize = 0;
294 nr_fmt *nr;
295 char *var;
296 char method = world->options->profilemethod;
297 MYREAL **param;
298 MYREAL *likes;
299 long elem = world->options->gamma ? world->numpop * world->numpop + 1 :
300 world->numpop * world->numpop;
301 #ifdef LONGSUM
302
303 elem += world->numpop * 3;
304 #endif /*LONGSUM*/
305
306 nr = (nr_fmt *) mycalloc (1, sizeof (nr_fmt));
307 create_nr (nr, world, *gmaxptr, which, world->loci,
308 world->repkind, world->rep);
309
310 var = (char *) mycalloc(LINESIZE, sizeof(char));
311 doublevec1d(&likes, GRIDSIZE);
312 //doublevec2d(¶m,GRIDSIZE, elem);
313 param = (MYREAL **) mycalloc(GRIDSIZE,sizeof(MYREAL));
314 for(i=0;i<GRIDSIZE;i++)
315 param[i] = (MYREAL*) mycalloc(elem+1,sizeof(MYREAL));
316
317 if (world->options->progress)
318 {
319 print_menu_profile (which,
320 world->numpop2 + (long) world->options->gamma);
321 }
322 //#ifndef INTEGRATEDLIKE
323 SETUPPARAM0 (world, nr, world->repkind, nr->repstart, nr->repstop,
324 world->loci, PROFILE, world->loci>1 ? TRUE : FALSE);
325 //#endif
326 calc_profile_likelihood (method, which, likes, param, world, GRIDSIZE, nr);
327 prepare_header (var, which, world);
328 if (world->options->printprofile)
329 {
330 bufsize = print_profile_table (method, which, var, likes, param, world);
331 world->allocbufsize = bufsize;
332 #ifdef PRETTY
333 prepare_profile_table(method, which, likes, param, world, bufsize);
334 #endif
335 }
336 if (world->options->printprofsummary)
337 {
338 strcpy (world->quantiles[which].name, var);
339 for (ii = 0; ii < GRIDSIZE; ii++)
340 {
341 i = warp (ii);
342 world->quantiles[which].param[i] = param[i][which];
343 }
344 }
345 destroy_nr (nr, world);
346 myfree(var);
347 myfree(likes);
348 //myfree(param[0]);
349 myfree(param);
350 }
351
352
353 ///
354 /// print the profile likelihood table
355 /// depending on the profile option
356 /// prints a mark to the precentiles that did not succeed to converge to the
357 /// proper percentile
358 /// \return bufsize size of buffer
359 long
print_profile_table(char method,long which,char * var,MYREAL likes[],MYREAL ** param,world_fmt * world)360 print_profile_table (char method, long which, char *var, MYREAL likes[],
361 MYREAL **param, world_fmt * world)
362 {
363 const MYREAL probabilities[] = SHOWGRID;
364 long failed = 0;
365 char star;
366 char methodstring[LINESIZE];
367
368 long i, ii, j, likei = 0;
369 long bufsize = 1;
370 long elem =
371 world->options->gamma ? world->numpop * world->numpop +
372 1 : world->numpop * world->numpop;
373
374 char fp[LONGLINESIZE]; //large local string buffer
375 char **buffer = &world->buffer; // buffer for printing whole table
376 char *temp, *var2, temp2[LINESIZE], likestr[LINESIZE], *paramstr;
377 MYREAL likemax = -MYREAL_MAX;
378 #ifdef LONGSUM
379
380 elem += world->numpop * 3;
381 #endif /*LONGSUM*/
382
383 //memset (*buffer, 0, sizeof (char) * (long) strlen (*buffer));
384 (*buffer)[0] = '\0';
385 var2 = (char *) mycalloc (1, sizeof (char) * MAX (LINESIZE, (long) (elem / 10. * LINESIZE)));
386 paramstr = (char *) mycalloc (1, sizeof (char) * MAX (LINESIZE,(long) (elem / 10. * LINESIZE)));
387 temp = var2;
388 for (j = 0; j < elem; j++)
389 {
390 prepare_header (temp2, j, world);
391 if (which == j)
392 star = '*';
393 else
394 star = ' ';
395 if (((j + 1) % 5) == 0)
396 {
397 sprintf (temp, "\n %c%.7s%c ", star,
398 temp2, star);
399 temp += 28;
400 }
401 else
402 sprintf (temp, "%c%.7s%c ", star, temp2, star);
403 temp += 11;
404 }
405 switch (method)
406 {
407 case 'p':
408 strcpy (methodstring,
409 "Parameters are evaluated at percentiles\nusing bisection method (slow, but exact).");
410 break;
411 case 's':
412 strcpy (methodstring,
413 "Parameters are evaluated at percentiles\nusing cubic splines of profiled parameter\n(faster, but not so exact).");
414 break;
415 case 'd':
416 strcpy (methodstring,
417 "Parameters are evaluated at pre-defined values\n");
418 break;
419 case 'q':
420 strcpy (methodstring,
421 "Parameters are evaluated assuming complete independence\n");
422 break;
423 case 'f':
424 strcpy (methodstring,
425 "Parameters are evaluated assuming complete independence\n and then once maximized at the found profiled parameter value");
426 break;
427 }
428 sprintf (fp, "\n\nProfile likelihood for parameter %s\n", var);
429 bufsize += (long) strlen(fp);
430 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
431 strcat (*buffer, fp);
432 sprintf (fp, "%s\n", methodstring);
433 bufsize += (long) strlen(fp);
434 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
435 strcat (*buffer, fp);
436 sprint_line (fp, '-', 79, CONT);
437 bufsize += (long) strlen(fp);
438 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
439 strcat (*buffer, fp);
440 if (method == 'd')
441 sprintf (fp, " Ln(L) %7.7s %s\n", var, var2);
442 else
443 sprintf (fp, "Per. Ln(L) %7.7s %s\n", var, var2);
444 bufsize += (long) strlen(fp);
445 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
446 strcat (*buffer, fp);
447 sprint_line (fp, '-', 79, START);
448 bufsize += (long) strlen(fp);
449 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
450 strcat (*buffer, fp);
451 for (i = 0; i < GRIDSIZE; i++)
452 {
453 if (likemax < MIN (likes[i], MYREAL_MAX))
454 {
455 likemax = likes[i];
456 likei = i;
457 }
458 }
459 for (ii = 0; ii < GRIDSIZE; ii++)
460 {
461 i = warp (ii);
462 if (likes[i] >= MYREAL_MAX - EPSILON)
463 sprintf (likestr, " - ");
464 else
465 sprintf (likestr, "% 6.3f%c", likes[i], likei == i ? '*' : ' ');
466 temp = paramstr;
467 for (j = 0; j < elem; j++)
468 {
469 if (((j + 1) % 5) == 0)
470 {
471 sprintf (temp, "\n ");
472 temp += 26;
473 }
474 if (param[i][j] <= 0)
475 sprintf (temp, " - ");
476 else
477 {
478 if (param[i][j] < 0.000001)
479 sprintf (temp, " % 6.3e ", param[i][j]);
480 else
481 sprintf (temp, "% 10.6f ", param[i][j]);
482 }
483
484 temp += 11;
485 }
486 if (param[i][which] <= 0)
487 {
488 if (method == 'd')
489 sprintf (fp, " - - %s\n", paramstr);
490 else
491 sprintf (fp, "%3.3f - - %s\n",
492 probabilities[i], paramstr);
493 bufsize += (long) strlen(fp);
494 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
495 strcat (*buffer, fp);
496 }
497 else
498 {
499 if (method == 'd')
500 sprintf (fp, " %8.8s % 10.6g %s\n", likestr,
501 param[i][which], paramstr);
502 else
503 {
504 if (probabilities[i] == 0.5)
505 sprintf (fp, "MLE %8.8s % 10.6g %s\n",
506 likestr, param[i][which], paramstr);
507 else
508 {
509 if(world->percentile_failed[which][i])
510 {
511 sprintf (fp, "**** %8.8s % 10.6g %s\n",
512 likestr,
513 param[i][which], paramstr);
514 failed += 1;
515 }
516 else
517 sprintf (fp, "%4.3f %8.8s % 10.6g %s\n",
518 probabilities[i], likestr,
519 param[i][which], paramstr);
520 }
521 }
522 bufsize += (long) strlen(fp);
523 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
524 strcat (*buffer, fp);
525 }
526 }
527 sprint_line (fp, '-', 79, STOP);
528 bufsize += (long) strlen(fp);
529 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
530 strcat (*buffer, fp);
531 if(failed>0)
532 {
533 sprintf (fp, "If the percentile column is marked with **** then the convergence to the percentile value failed\nThe values are still _correct_ but not at the percentile of the profile parameter\n");
534 bufsize += (long) strlen(fp);
535 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
536 strcat (*buffer, fp);
537 }
538 myfree(paramstr);
539 myfree(var2);
540 return bufsize;
541 }
542
543
544 void
print_profile_title(world_fmt * world)545 print_profile_title (world_fmt * world)
546 {
547 print_line (world->outfile, '=', 79, CONT);
548 FPRINTF (world->outfile, "Profile likelihood tables [Summary is at the end of the file]\n");
549 print_line (world->outfile, '=', 79, CONT);
550 }
551
552
553 ///
554 /// set the profile method for ascii output
ascii_method_set(FILE * outfile,char method)555 void ascii_method_set(FILE *outfile, char method)
556 {
557 switch (method)
558 {
559 case 'p':
560 fprintf(outfile, "Parameters are evaluated at percentiles using bisection method (slow, but exact).\n");
561 break;
562 case 's':
563 fprintf(outfile,"Parameters are evaluated at percentiles\n");
564 fprintf(outfile,"using cubic splines of profiled parameter (faster, but not so exact).\n");
565 break;
566 case 'd':
567 fprintf(outfile,"Parameters are evaluated at pre-defined values\n");
568 break;
569 case 'q':
570 fprintf(outfile, "Parameters are evaluated assuming complete independence\n");
571 break;
572 case 'f':
573 fprintf(outfile, "Parameters are evaluated assuming complete independence\n");
574 fprintf(outfile,"and then once maximized at the found profiled parameter value\n");
575 break;
576 }
577 fprintf(outfile,"\n");
578 }
579
580 /*++++++++++++++++++++++++++ replacement for ascii table, not finished *************/
581 #if 0
582 ///
583 /// generic ascii table generator
584 /// \param float left_margin left edge of table
585 /// \param float * page_height page height
586 /// \param int cols number of columns in table if there are too many columns then new line
587 /// \param int rows number of rows in table, if a new page is needed then header is repeated
588 /// \param char *** elements all table elements, currently no formatting of these
589 /// \param char ** header header row
590 /// \param char *position position of all elements: L=left, R=right, C=center
591 /// \param int col_overflow when there are too many columns this is the column to restart
592 long ascii_table(int cols, int rows, char *buffer, char *position, int col_overflow, float separator)
593 {
594 boolean new_page=FALSE;
595
596 int row;
597 int col;
598 long location=0;
599
600 char pos;
601
602 float *col_widths;
603 float left_margin = 0.;
604 float page_width = 130.;
605 float right_margin = 130.;
606 float *col_starts;
607 float oldcol = HUGE;
608
609 char **header;
610 char ***elements;
611
612 col_widths = mycalloc(cols,sizeof(float));
613 col_starts = mycalloc(cols,sizeof(float));
614
615 charvec2d(&header, cols, LINESIZE);
616 elements = (char ***) mycalloc(rows, sizeof(char **));
617 for(row=0; row < rows; row++)
618 {
619 charvec2d(&(elements[row]),cols, LINESIZE);
620 }
621
622 location = translate_buffer_table(cols, rows, buffer, header, elements);
623 ascii_find_col_width(cols, rows, elements, header, col_widths);
624 define_col_start(cols, col_widths, col_overflow, position, left_margin, right_margin, separator, col_starts);
625
626 ascii_print_table_header(page_height, position, cols, col_starts,header);
627 sprint_line (fp, '-', 130, CONT);
628
629 for(row=0; row< rows; row++)
630 {
631 for(col=0;col < cols; col++)
632 {
633 if(col_starts[col] < oldcol)
634 {
635 fprintf(outfile,"\n");
636 }
637 if(position==NULL)
638 pos='L';
639 else
640 pos=position[col];
641
642 fprintf(outfilecol_starts[col],*page_height, pos, "%s",elements[row][col]);
643 oldcol = col_starts[col];
644 }
645 }
646 pdf_advance(page_height);
647 myfree(col_widths);
648 myfree(col_starts);
649 free_charvec2d(header);
650 for(row=0; row < rows; row++)
651 free_charvec2d(elements[row]);
652 myfree(elements);
653 return location;
654 }
655
656
657 ///
658 /// Printing PROFILE table for ASCII
659 /// \param method profile method
660 /// \param world master parameter structure
661 void ascii_print_profile_table(char method, world_fmt *world)
662 {
663 char *thisparam;
664 char *savedparam;
665 char *ptr;
666 char *buffer = world->buffer;
667 long position = 0;
668 FILE *outfile = world->outfile;
669
670 boolean failed = world->percentile_some_failed;
671 long thesize = WORDSIZE * world->numpop2;
672
673 thisparam = (char *) mycalloc(thesize ,sizeof(char));
674 savedparam = thisparam;
675 while(buffer[position] != '\0')
676 {
677 strncpy(thisparam, buffer, thesize-1);
678 strsep(&thisparam,"&"); // remove first element
679 strsep(&thisparam,"&"); // remove second element
680 ptr = strsep(&thisparam,"&"); // keep third
681
682 fprintf(outfile,"\n\n");
683 fprintf(outfile, "Profile likelihood table for parameter %s\n", ptr);
684 ascii_method_set(outfile, method);
685 position += pdf_table(left_margin, page_height,world->numpop2 + 3, 9, buffer + position,
686 NULL, 4, 6.);
687 pdf_table_footnote(left_margin, page_height, failed);
688 }
689 myfree(savedparam);
690 }
691 #endif /*++++++++++++++++++++++++++ replacement for ascii table, not finished *************/
692
693
694 long
warp(long ii)695 warp (long ii)
696 {
697 long i;
698 if (ii < 4)
699 i = ii;
700 else
701 {
702 if (ii == 4)
703 {
704 i = GRIDSIZE - 1;
705 }
706 else
707 {
708 i = GRIDSIZE - ii + 3;
709 }
710 }
711 return i;
712 }
713
714
715 ///
716 /// print profile tables for percentiles, returns the number
717 /// of failed convergences to a specific percentile.
718 long
print_profile_percentile(world_fmt * world)719 print_profile_percentile (world_fmt * world)
720 {
721 const MYREAL probabilities[] = SHOWGRID; //debug GRID2
722 long i, j, jj;
723 boolean first = TRUE;
724 boolean last = FALSE;
725 char fp[LONGLINESIZE];
726 char **buffer = &world->buffer;
727 long bufsize = 1;
728 long failed=0;
729
730
731 (*buffer)[0]='\0';
732
733 sprintf (fp, "\n\n");
734 bufsize += (long) strlen(fp);
735 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
736 strcat (*buffer, fp);
737
738 sprint_line (fp, '=', 79, CONT);
739 bufsize += (long) strlen(fp);
740 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
741 strcat (*buffer, fp);
742 sprintf (fp,
743 "Summary of profile likelihood percentiles of all parameters\n");
744 bufsize += (long) strlen(fp);
745 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
746 strcat (*buffer, fp);
747 sprint_line (fp, '=', 79, CONT);
748 bufsize += (long) strlen(fp);
749 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
750 strcat (*buffer, fp);
751 sprintf (fp, "\n");
752 bufsize += (long) strlen(fp);
753 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
754 strcat (*buffer, fp);
755 sprint_percentile_header (fp, first);
756 bufsize += (long) strlen(fp);
757 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
758 strcat (*buffer, fp);
759 for (i = 0; i < world->numpop2 + (long) world->options->gamma; i++)
760 {
761 if (world->quantiles[i].name[0] == '\0')
762 continue; /* this variable was not estimated, so return */
763 sprintf (fp, "%-10.10s ", world->quantiles[i].name);
764 bufsize += (long) strlen(fp);
765 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
766 strcat (*buffer, fp);
767 for (jj = 0; jj < GRIDSIZE; jj++)
768 {
769 j = warp (jj);
770 if (probabilities[j] > 0.5)
771 continue;
772 if (world->quantiles[i].param[j] < 0)
773 sprintf (fp, " - ");
774 else
775 failed += sprint_nice_param (world->quantiles[i].param[j], 0.000001,
776 999.99999, world->percentile_failed[i][j], fp);
777 bufsize += (long) strlen(fp);
778 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
779 strcat (*buffer, fp);
780 }
781 sprintf (fp, "\n");
782 bufsize += (long) strlen(fp);
783 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
784 strcat (*buffer, fp);
785 }
786 sprintf (fp, "\n\n");
787 bufsize += (long) strlen(fp);
788 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
789 strcat (*buffer, fp);
790 sprint_percentile_header (fp, last);
791 bufsize += (long) strlen(fp);
792 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
793 strcat (*buffer, fp);
794 for (i = 0; i < world->numpop2 + (long) world->options->gamma; i++)
795 {
796 if (world->quantiles[i].name[0] == '\0')
797 continue; /* this variable was not estimated, so return */
798 sprintf (fp, "%-10.10s ", world->quantiles[i].name);
799 bufsize += (long) strlen(fp);
800 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
801 strcat (*buffer, fp);
802 for (jj = 0; jj < GRIDSIZE; jj++)
803 {
804 j = warp (jj);
805 if (probabilities[j] < 0.5)
806 continue;
807 if (world->quantiles[i].param[j] < 0)
808 sprintf (fp, " - ");
809 else
810 failed += sprint_nice_param (world->quantiles[i].param[j], 0.000001,
811 999.99999, world->percentile_failed[i][j], fp);
812 bufsize += (long) strlen(fp);
813 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
814 strcat (*buffer, fp);
815 }
816 sprintf (fp, "\n");
817 bufsize += (long) strlen(fp);
818 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
819 strcat (*buffer, fp);
820 }
821
822 sprintf (fp, "\n");
823 bufsize += (long) strlen(fp);
824 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
825 strcat (*buffer, fp);
826 sprint_line (fp, '-', 80, CONT);
827 bufsize += (long) strlen(fp);
828 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
829 strcat (*buffer, fp);
830 if(failed > 0)
831 {
832 sprintf (fp, "Values with a * are NOT at the percentiles!\nThe convergence to the percentile value FAILED.\n");
833 bufsize += (long) strlen(fp);
834 *buffer = (char *) myrealloc(*buffer, sizeof(char) * bufsize);
835 strcat (*buffer, fp);
836 }
837 return bufsize;
838 }
839
840 long
sprint_nice_param(MYREAL parameter,MYREAL bottom,MYREAL top,boolean failed,char * file)841 sprint_nice_param (MYREAL parameter, MYREAL bottom, MYREAL top, boolean failed, char *file)
842 {
843 long failedcount=0;
844 if (parameter > bottom && parameter < top)
845 {
846 if(failed)
847 {
848 failedcount = 1;
849 sprintf (file, " %11.6f* ", parameter);
850 }
851 else
852 sprintf (file, " %11.6f ", parameter);
853 }
854 else
855 {
856 if(failed)
857 {
858 failedcount = 1;
859 sprintf (file, " %11.6g* ", parameter);
860 }
861 else
862 sprintf (file, " %11.6g ", parameter);
863 }
864 return failedcount;
865 }
866
867 void
sprint_percentile_header(char * outfile,boolean first)868 sprint_percentile_header (char *outfile, boolean first)
869 {
870 const MYREAL probabilities[] = GRID2;
871 long i;
872 char fp[LONGLINESIZE];
873 if (first)
874 sprintf (fp, "Parameter Lower percentiles\n");
875 else
876 sprintf (fp, "Parameter Upper percentiles\n");
877 strcat (outfile, fp);
878 sprintf (fp, " ");
879 strcat (outfile, fp);
880 sprint_line (fp, '-', 68, CONT);
881 strcat (outfile, fp);
882 sprintf (fp, " ");
883 strcat (outfile, fp);
884 for (i = 0; i < GRIDSIZE - 1; i++)
885 {
886 if (first)
887 {
888 if (probabilities[i] >= 0.5)
889 break;
890 }
891 else
892 {
893 if (probabilities[i] < 0.5)
894 continue;
895 }
896 if (probabilities[i] == 0.5)
897 sprintf (fp, " MLE ");
898 else
899 sprintf (fp, " %5.3f ", probabilities[i]);
900 strcat (outfile, fp);
901 }
902 if (probabilities[i] == 0.5)
903 sprintf (fp, " MLE\n");
904 else
905 sprintf (fp, " %5.3f\n", probabilities[i]);
906 strcat (outfile, fp);
907 sprint_line (fp, '-', 80, CONT);
908 strcat (outfile, fp);
909 }
910
911
912 void
prepare_header(char * var,long which,world_fmt * world)913 prepare_header (char *var, long which, world_fmt * world)
914 {
915 long i, j, from, to;
916 #ifdef LONGSUM
917
918 long len = world->numpop2 + world->options->gamma;
919 len += world->numpop * 3;
920 #endif /*LONGSUM*/
921
922 if (which < world->numpop)
923 {
924 sprintf (var, "Theta_%li", from = which + 1);
925 }
926 else
927 {
928 if (world->numpop > 1)
929 {
930 i = (which - world->numpop) / (long) (world->numpop - 1);
931 j = which - world->numpop - i * (world->numpop - 1);
932 to = i + 1;
933 from = (i <= j) ? j + 2 : j + 1;
934 if (world->options->profileparamtype == PLOT4NM)
935 sprintf (var, "4Nm_%li%li", from, to);
936 else
937 sprintf (var, " M_%li%li", from, to);
938 }
939 if (which == world->numpop2 && world->options->gamma)
940 sprintf (var, " Alpha ");
941 #ifdef LONGSUM
942
943 if(which >= len && world->options->fluctuate)
944 sprintf(var," Rate_%li%li ", (which - len) / world->numpop, ((which - len) % world->numpop)+1);
945 #endif /*LONGSUM*/
946
947 }
948 }
949
950 void
calc_profile_likelihood(char method,long which,MYREAL * likes,MYREAL ** param,world_fmt * world,long nn,nr_fmt * nr)951 calc_profile_likelihood (char method, long which, MYREAL *likes,
952 MYREAL **param, world_fmt * world, long nn,
953 nr_fmt * nr)
954 {
955 long i, j;
956 switch (method)
957 {
958 case 'p' /*percentiles */ :
959 case 'q': /*quick and dirty: no correlation between parametes */
960 case 'f': /*mixture of quick and precise */
961 profile_max_percentiles (which, likes, param, world, nn, nr);
962 break;
963 case 's': /* spline */
964 case 'd' /*discrete */ :
965 profile_max_spline (method, which, likes, param, world, nr, nn);
966 break;
967 }
968 if (world->options->profileparamtype == PLOT4NM)
969 {
970 for (i = 0; i < GRIDSIZE; i++)
971 {
972 for (j = world->numpop; j < world->numpop2; j++)
973 {
974 param[i][j] *=
975 param[i][(long) (j - world->numpop) / (world->numpop - 1)];
976 }
977 }
978 }
979 //else
980 // printf (" Using M (m/mu) instead of 4Nm\n");
981 }
982
983
984 void
prognose_profile_end(time_t starttime,long numpop2,long nn)985 prognose_profile_end (time_t starttime, long numpop2, long nn)
986 {
987 #ifndef NOTIME_FUNC
988 static long sumtime = 0;
989 static long tally = 0;
990 char * nowstr;
991 time_t endtime;
992 time_t nowbin;
993 struct tm *nowstruct;
994 nowstr = (char *) mycalloc(STRSIZE,sizeof(char));
995
996 if (tally % nn == 0)
997 {
998 tally++;
999 time (&endtime);
1000 sumtime = (long) (endtime - starttime);
1001 #ifdef MPI
1002
1003 nowbin = starttime + (sumtime / tally * (nn * numpop2));
1004 #else
1005
1006 nowbin = starttime + (sumtime / tally * (nn * numpop2));
1007 #endif
1008
1009 nowstruct = localtime (&nowbin);
1010
1011 strftime (nowstr, STRSIZE, "%H:%M %B %d %Y", nowstruct);
1012 FPRINTF (stdout, " Prognosed end of run: %s\n", nowstr);
1013 }
1014 else
1015 tally++;
1016 myfree(nowstr);
1017 #endif
1018 }
1019
1020
1021 void
profile_max_percentiles(long which,MYREAL * likes,MYREAL ** param,world_fmt * world,long nn,nr_fmt * nr)1022 profile_max_percentiles (long which, MYREAL *likes, MYREAL **param,
1023 world_fmt * world, long nn, nr_fmt * nr)
1024 {
1025 const MYREAL probabilities[] = GRID;
1026
1027 long i;
1028 //, trials = 0;
1029 long df;
1030 long len;
1031 boolean QD = FALSE, QD2 = FALSE;
1032 MYREAL prob, normd, maxlike, testlike, minparam = 0, maxparam = 0;
1033 MYREAL *mlparam;
1034 long rep = world->options->replicate ? (world->loci > 1 ? 0 : world->repstop) : 0;
1035 df = (world->options->df > 0) ? world->options->df : set_df (which,
1036 world->options->
1037 custm2,
1038 world->numpop,
1039 world->numpop2 +
1040 (long) world->
1041 options->gamma);
1042
1043 /* QD assumes that there is NO correlation between parameters
1044 that is certainly wrong but sometimes close and much faster
1045 */
1046 if (world->options->profilemethod == 'q')
1047 QD = TRUE;
1048 else
1049 {
1050 if (world->options->profilemethod == 'f')
1051 {
1052 QD = TRUE;
1053 QD2 = TRUE;
1054 }
1055 else
1056 QD = FALSE;
1057 }
1058 /* the minimum */
1059 minparam = which < world->numpop ? SMALLEST_THETA : SMALLEST_MIGRATION;
1060 /* the maximum */
1061 if (world->loci > 1)
1062 {
1063 mlparam = world->atl[rep][world->loci].param;
1064 maxlike = world->atl[rep][world->loci].param_like;
1065 }
1066 else
1067 {
1068 mlparam = world->atl[rep][world->locus].param;
1069 maxlike = world->atl[rep][world->locus].param_like;
1070 }
1071
1072 len = world->numpop2 + (long) world->options->gamma;
1073 #ifdef LONGSUM
1074 len += world->numpop * 3;
1075 #endif /*LONGSUM*/
1076
1077 memcpy (param[nn - 1], mlparam, sizeof (MYREAL) * len);
1078
1079 likes[nn - 1] = maxlike;
1080 maxparam = param[nn - 1][which];
1081
1082 for (i = 0; i < nn - 1; i++)
1083 {
1084 if (probabilities[i] > 0.5)
1085 prob = 1 - probabilities[i];
1086 else
1087 prob = probabilities[i];
1088 testlike = maxlike - 0.5 * find_chi (df, prob);
1089 if (i > 0)
1090 {
1091 minparam = param[i - 1][which];
1092 if (minparam < 0)
1093 minparam = which < world->numpop ?
1094 SMALLEST_THETA : SMALLEST_MIGRATION;
1095 }
1096 #ifdef __MWERKS__
1097 eventloop ();
1098 #endif
1099
1100 len = world->numpop2 + (long) world->options->gamma;
1101 #ifdef LONGSUM
1102
1103 len += world->numpop * 3;
1104 #endif /*LONGSUM*/
1105
1106 memcpy (param[i], mlparam, sizeof (MYREAL) * len);
1107
1108 //xcode trials =
1109 find_profilelike (testlike, prob, i, minparam, maxparam, ¶m[i],
1110 &likes[i], which, world, QD, QD2, mlparam, &normd,
1111 nr);
1112 world->percentile_failed[which][i] = ((fabs(likes[i] - testlike)) > PERCENTILETOLERANCE);
1113 if (world->options->progress)
1114 prognose_profile_end (world->starttime, world->numpop2 +
1115 (long) world->options->gamma, nn);
1116 }
1117 }
1118
1119 void
profile_max_spline(char method,long which,MYREAL * likes,MYREAL ** param,world_fmt * world,nr_fmt * nr,long nn)1120 profile_max_spline (char method, long which, MYREAL *likes, MYREAL **param,
1121 world_fmt * world, nr_fmt * nr, long nn)
1122 {
1123
1124 MYREAL deviate[] = DEVIATE;
1125 //xcode long errc;
1126 long indeks[] = INDEX;
1127 const MYREAL probabilities[] = GRID;
1128 long i;
1129 //, trials = 0;
1130 long panic, df;
1131 long len;
1132 MYREAL prob, normd, maxlike, testlike;
1133 //, maxparam = 0;
1134 //minparam = 0,
1135 MYREAL *mlparam;
1136 MYREAL **hess;
1137 MYREAL tmp, x, xx, xlow, xhigh, mid, low, high, value;
1138 long rep =
1139 world->options->replicate ? (world->loci > 1 ? 0 : world->repstop) : 0;
1140 spline_fmt *spline;
1141 /* set the hessian for hte maximizer*/
1142 doublevec2d(&hess, world->numpop2 + (long) world->options->gamma,
1143 world->numpop2 + (long) world->options->gamma);
1144 /* set degree of freedom */
1145 df =
1146 (world->options->df > 0) ? world->options->df : set_df (which,
1147 world->options->
1148 custm2,
1149 world->numpop,
1150 world->numpop2 +
1151 (long) world->
1152 options->gamma);
1153 /* the minimum */
1154 //xcode minparam = which < world->numpop ? SMALLEST_THETA : SMALLEST_MIGRATION;
1155 /* the maximum */
1156 if (world->loci > 1)
1157 {
1158 mlparam = world->atl[rep][world->loci].param;
1159 maxlike = world->atl[rep][world->loci].param_like;
1160 }
1161 else
1162 {
1163 mlparam = world->atl[rep][world->locus].param;
1164 maxlike = world->atl[rep][world->locus].param_like;
1165 }
1166 len = world->numpop2 + (long) world->options->gamma;
1167 #ifdef LONGSUM
1168
1169 len += world->numpop * 3;
1170 #endif /*LONGSUM*/
1171
1172 memcpy (param[nn - 1], mlparam, sizeof (MYREAL) * len);
1173 likes[nn - 1] = maxlike;
1174 //xcode maxparam = param[nn - 1][which];
1175 /* calculate likelihood at GRIDNODES */
1176 for (i = 0; i < nn - 1; i++)
1177 {
1178 if (which >= world->numpop && mlparam[which] < EPSILON)
1179 value = deviate[i];
1180 else
1181 value = deviate[i] * mlparam[which];
1182 memcpy (param[i], mlparam, sizeof (MYREAL) * len);
1183 nr->profilenum = set_profile_param (param[i], which, nr->profiles,
1184 value, nr->values, world);
1185 //do_profiles(nr->world, nr, likes, &normd, PROFILE,
1186 // nr->world->rep, nr->world->repkind);
1187 maximize (¶m[i], world, nr, hess, PROFILE, world->repkind);
1188 likes[i] = nr->llike;
1189 normd = nr->normd;
1190 memcpy (param[i], world->param0, sizeof (MYREAL) * nr->partsize);
1191 //xcode trials = 1;
1192 if (method != 's')
1193 prognose_profile_end (world->starttime,
1194 world->numpop2 + (long) world->options->gamma,
1195 nn);
1196 }
1197 if (method == 's')
1198 {
1199 spline = (spline_fmt *) mycalloc (1, sizeof (spline_fmt));
1200 setup_spline (spline);
1201 prepare_spline_nodes (nn, which, param, likes, indeks, spline->param,
1202 spline->like);
1203 prepare_spline_first (spline->param, spline->like);
1204 prepare_spline_last (nn, spline->param, spline->like);
1205 calc_spline_diff (spline->diff, spline->param, spline->like, nn);
1206 calc_spline (spline->param, spline->like, nn, spline->constr,
1207 spline->diff, spline->diff2, spline->diagn, spline->work,spline->nwork);
1208 for (i = 0; i < nn - 1; i++)
1209 {
1210 if (probabilities[i] < 0.5)
1211 {
1212 xhigh = spline->param[(nn + 2) / 2];
1213 high = spline->like[(nn + 2) / 2];
1214 xlow = spline->param[0];
1215 low = spline->like[0];
1216 prob = probabilities[i];
1217 }
1218 else
1219 {
1220 xlow = spline->param[nn + 1];
1221 low = spline->like[nn + 1];
1222 xhigh = spline->param[(nn + 2) / 2];
1223 high = spline->like[(nn + 2) / 2];
1224 prob = 1. - probabilities[i];
1225 }
1226 testlike = maxlike - 0.5 * find_chi (df, prob);
1227 x = (xhigh + xlow) / 2.;
1228 panic = 0;
1229 mid = testlike - 1000.;
1230 while (panic++ < MAX_PROFILE_TRIALS
1231 && fabs (mid - testlike) > BIGEPSILON)
1232 {
1233 mid =
1234 recalc_profilelike (x, spline->param, spline->like, nn,
1235 spline->diff, spline->diff2, spline->work,
1236 spline->nwork);
1237 if (mid < low || mid > high)
1238 break;
1239 if (testlike < mid)
1240 {
1241 high = mid;
1242 tmp = x;
1243 x = (xlow + tmp) / 2;
1244 xhigh = tmp;
1245 }
1246 else
1247 {
1248 low = mid;
1249 tmp = x;
1250 x = (tmp + xhigh) / 2;
1251 xlow = tmp;
1252 }
1253 }
1254 xx = EXP (x);
1255 nr->profilenum = set_profile_param (param[i], which, nr->profiles,
1256 xx, nr->values, world);
1257 do_profiles (nr->world, nr, likes, &normd, PROFILE,
1258 nr->world->rep, nr->world->repkind);
1259 // maximize(param[i], world, nr, PROFILE, world->repkind);
1260 likes[i] = nr->llike;
1261 normd = nr->normd;
1262 memcpy (param[i], world->param0, sizeof (MYREAL) * nr->partsize);
1263 //xcode x = LOG (xx);
1264 prognose_profile_end (world->starttime,
1265 world->numpop2 + (long) world->options->gamma,
1266 nn);
1267 }
1268 destroy_spline (spline);
1269 myfree(hess[0]);
1270 myfree(hess);
1271 }
1272 }
1273
1274 // setup spline: allocs all the necessary arrays and fill with
1275 // values when necessary.
1276 void
setup_spline(spline_fmt * spline)1277 setup_spline (spline_fmt * spline)
1278 {
1279 const long degree = 3;
1280 long i;
1281 spline->ntab = 1;
1282 spline->nwork = (GRIDSIZE + 2) * 9 + 5 + degree * (degree + 11) / 2 + 9;
1283 spline->param = (MYREAL *) mycalloc (1, sizeof (MYREAL) * (GRIDSIZE + 2));
1284 spline->like = (MYREAL *) mycalloc (1, sizeof (MYREAL) * (GRIDSIZE + 2));
1285 spline->diff = (MYREAL *) mycalloc (1, sizeof (MYREAL) * (GRIDSIZE + 2));
1286 spline->diff2 = (MYREAL *) mycalloc (1, sizeof (MYREAL) * (GRIDSIZE + 2));
1287 spline->constr = (long *) mycalloc (1, (GRIDSIZE + 2) * sizeof (long));
1288 spline->diagn = (long *) mycalloc (1, (GRIDSIZE + 2) * sizeof (long));
1289 spline->work = (MYREAL *) mycalloc (1, spline->nwork * sizeof (MYREAL));
1290
1291 for (i = 0; i < GRIDSIZE + 2; i++)
1292 spline->constr[i] = 3;
1293 }
1294
1295 void
destroy_spline(spline_fmt * spline)1296 destroy_spline (spline_fmt * spline)
1297 {
1298 myfree(spline->param);
1299 myfree(spline->like);
1300 myfree(spline->diff);
1301 myfree(spline->diff2);
1302 myfree(spline->constr);
1303 myfree(spline->diagn);
1304 myfree(spline->work);
1305 }
1306
1307 /* calc_spline() calculates the spline function derivatives,
1308 and adds two additional
1309 points, these two points bracket the real data points:
1310 so that param[0]=0 and param[n+1] = 1000 * the last real param,
1311 likes[0] is set -MYREAL_MAX and likes[n+1] is set using a linear interpolation
1312 using the last two real data points.
1313
1314 Attention the arrays param, likes, and yy need to have a size
1315 of n+2 and not n
1316 */
1317 int
calc_spline(MYREAL * param,MYREAL * like,long nn,long * constr,MYREAL * diff,MYREAL * diff2,long * diagn,MYREAL * work,long nwork)1318 calc_spline (MYREAL *param, MYREAL *like, long nn, long *constr, MYREAL *diff,
1319 MYREAL *diff2, long *diagn, MYREAL *work, long nwork)
1320 {
1321 /* call to spline.c routines */
1322 long degree = 3;
1323 long smoothness = 1;
1324 long opt = 414; //convex, monotone, derivatives constraints
1325 MYREAL d0 = 0, dnp = 0, d20 = 0, d2np = 0; //not used
1326 MYREAL eps = 0.0001;
1327 long kmax = 10; //not used
1328 long maxstp = 10; //not used
1329 long errc;
1330 dbvssc_ (param, like, &nn, °ree, &smoothness, &opt, &d0, &dnp, &d20,
1331 &d2np, constr, &eps, NULL, NULL, NULL, NULL, &kmax, &maxstp, &errc,
1332 diff, diff2, diagn, work, &nwork);
1333 /* end call to spline.c routines */
1334 return (long) errc;
1335 }
1336
1337 void
prepare_spline_nodes(long n,long which,MYREAL ** rawparam,MYREAL * rawlikes,long * indeks,MYREAL * param,MYREAL * likes)1338 prepare_spline_nodes (long n, long which, MYREAL **rawparam, MYREAL *rawlikes,
1339 long *indeks, MYREAL *param, MYREAL *likes)
1340 {
1341 long i;
1342 for (i = 0; i < n; i++)
1343 {
1344 param[i + 1] = LOG (rawparam[indeks[i]][which]);
1345 likes[i + 1] = rawlikes[indeks[i]];
1346 }
1347 }
1348
1349 void
prepare_spline_first(MYREAL * param,MYREAL * likes)1350 prepare_spline_first (MYREAL *param, MYREAL *likes)
1351 {
1352 const long i = 2;
1353 const long i0 = 1;
1354 param[0] = LOG (SMALLEST_THETA);
1355
1356 likes[0] =
1357 ((likes[i] - likes[i0]) / (param[i] - param[i0]) * param[0]) +
1358 (likes[i0] + param[i0] * (likes[i] + likes[i0]) / (param[i] - param[i0]));
1359 }
1360
1361 void
prepare_spline_last(long n,MYREAL * param,MYREAL * likes)1362 prepare_spline_last (long n, MYREAL *param, MYREAL *likes)
1363 {
1364 long n0 = n - 1;
1365 long n1 = n + 1;
1366 param[n1] = 100000. + param[n];
1367 likes[n1] =
1368 ((likes[n] - likes[n0]) / (param[n] - param[n0]) * param[n1]) +
1369 (likes[n0] + param[n0] * (likes[n] + likes[n0]) / (param[n] - param[n0]));
1370 }
1371
1372
1373 void
calc_spline_diff(MYREAL * diff,MYREAL * param,MYREAL * likes,long n)1374 calc_spline_diff (MYREAL *diff, MYREAL *param, MYREAL *likes, long n)
1375 {
1376 long i;
1377 diff[0] = ((likes[1] - likes[0]) / (param[1] - param[0]));
1378 diff[n] = ((likes[n] - likes[n + 1]) / (param[n] - param[n + 1]));
1379 for (i = 1; i < n; i++)
1380 {
1381 diff[i] =
1382 ((likes[i + 1] - likes[i - 1]) / (param[i + 1] - param[i - 1]));
1383 }
1384 diff[GRIDMIDDLE] = 0.0;
1385 }
1386
1387
1388 /* recalc_profilelike uses the spline derivatives and returns a new
1389 loglikelihood value,
1390 testparam is THE LOG(testparam) !
1391 */
1392 MYREAL
recalc_profilelike(MYREAL testparam,MYREAL * param,MYREAL * like,long nn,MYREAL * diff,MYREAL * diff2,MYREAL * work,long nwork)1393 recalc_profilelike (MYREAL testparam, MYREAL *param, MYREAL *like, long nn,
1394 MYREAL *diff, MYREAL *diff2, MYREAL *work, long nwork)
1395 {
1396 long degree = 3;
1397 long smoothness = 1;
1398 long errc;
1399 long sbopt = 2;
1400 long ntab = 1;
1401 long yy0 = 1;
1402 long yy1 = 1;
1403 long yy2 = 1;
1404 long erre;
1405 MYREAL tmp;
1406
1407 MYREAL *y0tab, *y1tab, *y2tab;
1408 //MYREAL *gaga;
1409 //long i,gagatab=10;
1410 //gaga= (MYREAL*)calloc(1,sizeof(MYREAL)*10);
1411 y0tab = (MYREAL *) mycalloc (1, sizeof (MYREAL) * 10);
1412 y1tab = (MYREAL *) mycalloc (1, sizeof (MYREAL) * 10);
1413 y2tab = (MYREAL *) mycalloc (1, sizeof (MYREAL) * 10);
1414 //gaga[0]= -5;
1415 //for(i=1;i<gagatab;i++)
1416 //gaga[i] = gaga[i-1] + 0.5;
1417
1418 dbvsse_ (param, like, &nn, °ree, &smoothness,
1419 //gaga, &gagatab,
1420 &testparam, &ntab, &sbopt, &yy0, &yy1, &yy2, &errc, diff, diff2,
1421 y0tab, y1tab, y2tab, &erre, work, &nwork);
1422 //for(i=0;i<gagatab;i++)
1423 //printf("%f %f\n",gaga[i],y0tab[i]);
1424 //exit(0);
1425
1426 tmp = y0tab[0];
1427 myfree(y0tab);
1428 myfree(y1tab);
1429 myfree(y2tab);
1430 return tmp;
1431
1432 }
1433
1434
1435 long
set_profile_param(MYREAL * param,long which,long * allwhich,MYREAL xval,MYREAL * xvals,world_fmt * world)1436 set_profile_param (MYREAL *param, long which, long *allwhich,
1437 MYREAL xval, MYREAL *xvals, world_fmt * world)
1438 {
1439 boolean has_fixed = FALSE;
1440
1441 long z = 0, i = 0, zz = 0;//zz is the counter for # fixed elements
1442 //xcode zzz = 0;
1443 long numpop = world->numpop;
1444 long numpop2 = world->numpop2;
1445 char *p;
1446 char *custm2 = world->options->custm2;
1447 long *fixed = NULL;
1448 long profnum = 0;
1449 MYREAL thetax, thetay;
1450 // find all real zero parameters and adds them to
1451 // the list of parameters not to maximize.
1452 // [inefficient, then this will be called several times
1453 // and recalculated for nothing]
1454 if (strpbrk (world->options->custm, "0c"))
1455 {
1456 has_fixed = TRUE;
1457 fixed = (long *) mycalloc (1, sizeof (long) * (numpop2 + (long) world->options->gamma));
1458 p = world->options->custm2;
1459 while (*p != '\0')
1460 {
1461 if ((*p == '0' || *p == 'c') && i != which)
1462 {
1463 fixed[zz]= i;
1464 zz++;
1465 }
1466 p++;
1467 i++;
1468
1469 }
1470 }
1471 switch (custm2[which])
1472 {
1473 case '0':
1474 profnum = 0;
1475 break;
1476 case '*':
1477 param[which] = xval;
1478 allwhich[0] = which, xvals[0] = xval;
1479 if (has_fixed)
1480 {
1481 //xcode zzz = 0;
1482 addinvar (param, fixed, allwhich, xvals, 1, zz)
1483 ;
1484 myfree(fixed)
1485 ;
1486 }
1487 profnum = 1 + zz;
1488 break;
1489 case 'S':
1490 z = 0;
1491 while (z < world->options->sym2n &&
1492 which != world->options->sym2param[z][0] &&
1493 which != world->options->sym2param[z][1])
1494 z++;
1495 thetax = param[world->options->sym2param[z][2]];
1496 thetay = param[world->options->sym2param[z][3]];
1497 if (which == world->options->sym2param[z][0])
1498 {
1499 allwhich[0] = which;
1500 allwhich[1] = world->options->sym2param[z][1];
1501 param[world->options->sym2param[z][0]] = xval;
1502 param[world->options->sym2param[z][1]] = xval * thetax / thetay;
1503 xvals[0] = xval;
1504 xvals[1] = xval * thetax / thetay;
1505 }
1506 else
1507 {
1508 allwhich[1] = world->options->sym2param[z][0];
1509 allwhich[0] = which;
1510 param[world->options->sym2param[z][0]] = xval * thetay / thetax;
1511 param[world->options->sym2param[z][1]] = xval;
1512 xvals[1] = xval * thetay / thetax;
1513 xvals[0] = xval;
1514 }
1515 if (has_fixed)
1516 {
1517 addinvar (param, fixed, allwhich, xvals, 2, zz)
1518 ;
1519 myfree(fixed)
1520 ;
1521 }
1522 profnum = 2 + zz;
1523 // profnum = 1 + zz;
1524 break;
1525 case 's':
1526 z = 0;
1527 while (z < world->options->symn &&
1528 which != world->options->symparam[z][0] &&
1529 which != world->options->symparam[z][1])
1530 z++;
1531 param[world->options->symparam[z][0]] = xval;
1532 param[world->options->symparam[z][1]] = xval;
1533 if (which == world->options->symparam[z][0])
1534 {
1535 allwhich[0] = which;
1536 allwhich[1] = world->options->symparam[z][1];
1537 }
1538 else
1539 {
1540 allwhich[0] = world->options->symparam[z][0];
1541 allwhich[1] = which;
1542 }
1543 xvals[0] = xval;
1544 xvals[1] = xval;
1545 if (has_fixed)
1546 {
1547 addinvar (param, fixed, allwhich, xvals, 2, zz)
1548 ;
1549 myfree(fixed)
1550 ;
1551 }
1552 profnum = 2 + zz;
1553 break;
1554 case 'm':
1555 if (which < numpop)
1556 {
1557 for (i = 0; i < numpop; i++)
1558 {
1559 param[i] = xval;
1560 xvals[i] = xval;
1561 allwhich[i] = i;
1562 }
1563 if (has_fixed)
1564 {
1565 addinvar (param, fixed, allwhich, xvals, numpop, zz)
1566 ;
1567 myfree(fixed)
1568 ;
1569 }
1570 profnum = numpop + zz;
1571 break;
1572 }
1573 else
1574 {
1575 for (i = numpop; i < numpop2 + (long) world->options->gamma; i++)
1576 {
1577 param[i] = xval;
1578 xvals[i - numpop] = xval;
1579 allwhich[i - numpop] = i;
1580 }
1581 profnum = numpop2 + (long) world->options->gamma - numpop;
1582 break;
1583 }
1584 }
1585 return profnum;
1586 }
1587
1588 // adds the invariants to the profiles list so that we do not
1589 // evaluate them and run unnecessary cycles in the maximizer
1590 // fixed = the index of the fixed parameters
1591 // allwhich = the variables in the profile list
1592 // start = last element filled element in allwhich
1593 // len = length of the fixed array
1594 // PB 06/30/99
1595 void
addinvar(MYREAL * param,long * fixed,long * allwhich,MYREAL * xvals,long start,long len)1596 addinvar (MYREAL *param, long *fixed, long *allwhich, MYREAL *xvals,
1597 long start, long len)
1598 {
1599 long i, z = 0;
1600 for (i = start; i < len + start; i++)
1601 {
1602 allwhich[i] = fixed[z]
1603 ;
1604 xvals[i] = param[fixed[z]]
1605 ;
1606 ++z;
1607 }
1608
1609 }
1610
1611 long
find_profilelike(MYREAL testlike,MYREAL prob,long whichprob,MYREAL minparam,MYREAL maxparam,MYREAL ** testparam,MYREAL * llike,long which,world_fmt * world,boolean QD,boolean QD2,MYREAL * mlparam,MYREAL * normd,nr_fmt * nr)1612 find_profilelike (MYREAL testlike, MYREAL prob, long whichprob, MYREAL minparam,
1613 MYREAL maxparam, MYREAL **testparam, MYREAL *llike,
1614 long which, world_fmt * world, boolean QD, boolean QD2,
1615 MYREAL *mlparam, MYREAL *normd, nr_fmt * nr)
1616 {
1617 const MYREAL probabilities[] = GRID;
1618 MYREAL tmp;
1619 MYREAL xlow, low;
1620 MYREAL x, mid = -MYREAL_MAX;
1621 MYREAL xhigh, high;
1622 long panic = 0;
1623 helper_fmt helper;
1624 MYREAL *ltestparam;
1625 MYREAL **hess;
1626 long len = world->numpop2 + (long) world->options->gamma;
1627 doublevec2d(&hess,len,len);
1628 #ifdef LONGSUM
1629
1630 len += world->numpop * 3;
1631 #endif /*LONGSUM*/
1632 world->locus = world->loci; // in case we have used gamma
1633 // used to set gamma option correctly, overly complicated, need
1634 // revamping, because the gamma test is all over the place PB022403
1635
1636 ltestparam = (MYREAL *) mycalloc (len, sizeof (MYREAL));
1637 *normd = -1;
1638 prepare_broyden (PROFILE, world, &helper.multilocus);
1639 if (probabilities[whichprob] < 0.5)
1640 {
1641 xhigh = maxparam;
1642 xlow = whichprob > 0 ? minparam : SMALLEST_THETA / 10000.;
1643 }
1644 else
1645 {
1646 if (probabilities[whichprob] > 0.5)
1647 {
1648 xhigh = maxparam;
1649 if (maxparam < EPSILON)
1650 xlow = probabilities[whichprob] >= 0.98 ? 1000000.0 : minparam;
1651 else
1652 xlow =
1653 probabilities[whichprob] >= 0.98 ? maxparam * 10000 : minparam;
1654 }
1655 else
1656 {
1657 xlow = minparam;
1658 xhigh = maxparam;
1659 }
1660 }
1661
1662
1663 // QD=TRUE -> calculate only CALCLIKE for parameter
1664 // QD2=TRUE -> assumes QD and then adds a maximizer step
1665 // if both are false maximize all parameters
1666 if (QD)
1667 {
1668 memcpy (*testparam, mlparam, sizeof (MYREAL) * len);
1669 nr->profilenum = set_profile_param (*testparam, which, nr->profiles,
1670 xlow, nr->values, world);
1671 (*testparam)[which] = xlow;
1672 set_logparam (ltestparam, *testparam, nr->partsize);
1673 //SETUPPARAM0 (world, nr, world->repkind, nr->repstart, nr->repstop,
1674 // world->loci, PROFILE, helper.multilocus);
1675 if(world->options->gamma)
1676 initgammacat (nr->categs, (*testparam)[nr->numpop2],1./* EXP (lparam[0])*/,
1677 nr->rate, nr->probcat);
1678 fill_helper (&helper, *testparam, ltestparam, nr->world, nr);
1679 low = CALCLIKE (&helper, *testparam, ltestparam);
1680 memcpy (*testparam, mlparam, sizeof (MYREAL) *
1681 (world->numpop2 + (long) world->options->gamma));
1682 nr->profilenum = set_profile_param (*testparam, which, nr->profiles,
1683 xhigh, nr->values, world);
1684 (*testparam)[which] = xhigh;
1685 set_logparam (ltestparam, *testparam, nr->partsize);
1686 if(world->options->gamma)
1687 initgammacat (nr->categs, (*testparam)[nr->numpop2],1./* EXP (lparam[0])*/,
1688 nr->rate, nr->probcat);
1689 fill_helper (&helper, *testparam, ltestparam, world, nr);
1690 high = CALCLIKE (&helper, *testparam, ltestparam);
1691 // testlike = high - 0.5 * find_chi (1, prob);
1692 }
1693 else
1694 {
1695 memcpy (*testparam, mlparam, sizeof (MYREAL) * len);
1696 nr->profilenum = set_profile_param (*testparam, which, nr->profiles,
1697 xlow, nr->values, world);
1698 maximize (testparam, world, nr, hess, PROFILE, world->repkind);
1699 low = nr->llike;
1700 memcpy (*testparam, mlparam, sizeof (MYREAL) *len);
1701 nr->profilenum = set_profile_param (*testparam, which, nr->profiles,
1702 xhigh, nr->values, world);
1703 maximize (testparam, world, nr, hess, PROFILE, world->repkind);
1704 high = nr->llike;
1705
1706 // testlike = high - 0.5 * find_chi (1, prob);
1707
1708 //if (fabs(high - world->atl[rep][loc].param_like)>EPSILON*10)
1709 // printf
1710 // ("%i> PROBLEM?: high=%f ml=%f (high=?=ml) xigh=%f mlparam[which]=%f testlike=%f\n",
1711 // myID, high, world->atl[rep][loc].param_like, xhigh,
1712 // mlparam[which], testlike);
1713 }
1714 panic = 0;
1715 x = (xhigh + xlow) / 2.;
1716 while (panic++ < MAX_PROFILE_TRIALS && fabs (mid - testlike) > BIGEPSILON)
1717 {
1718 if (QD)
1719 {
1720 memcpy (*testparam, mlparam, sizeof (MYREAL) * len);
1721 nr->profilenum = set_profile_param (*testparam, which,
1722 nr->profiles, x,
1723 nr->values, world);
1724 (*testparam)[which] = x;
1725 if(nr->world->options->gamma)
1726 initgammacat (nr->categs, (*testparam)[nr->numpop2],1.,nr->rate, nr->probcat);
1727 set_logparam (ltestparam, *testparam, nr->partsize);
1728 fill_helper (&helper, *testparam, ltestparam, world, nr);
1729 mid = CALCLIKE (&helper, *testparam, ltestparam);
1730 }
1731 else
1732 {
1733 memcpy (*testparam, mlparam, sizeof (MYREAL) * nr->partsize);
1734 nr->profilenum = set_profile_param (*testparam, which, nr->profiles,
1735 x, nr->values, world);
1736 maximize (testparam, world, nr, hess, PROFILE, world->repkind);
1737 mid = *llike = nr->llike;
1738 *normd = nr->normd;
1739 }
1740 if ((testlike < low) || (testlike > high))
1741 return MAX_PROFILE_TRIALS;
1742 if (testlike < mid)
1743 {
1744 high = mid;
1745 tmp = x;
1746 x = (xlow + tmp) / 2;
1747 xhigh = tmp;
1748 }
1749 else
1750 {
1751 low = mid;
1752 tmp = x;
1753 x = (tmp + xhigh) / 2;
1754 xlow = tmp;
1755 }
1756 }
1757 if (QD)
1758 {
1759
1760 *llike = mid;
1761 if (QD2)
1762 {
1763 x = (*testparam)[which];
1764 memcpy (*testparam, mlparam,
1765 sizeof (MYREAL) * (nr->numpop2 +
1766 (long) world->options->gamma));
1767 nr->profilenum =
1768 set_profile_param (*testparam, which, nr->profiles, x, nr->values,
1769 world);
1770 maximize (testparam, world, nr, hess, PROFILE, world->repkind);
1771 *llike = nr->llike;
1772 *normd = nr->normd;
1773 memcpy (*testparam, world->param0, sizeof (MYREAL) * nr->partsize);
1774 }
1775 }
1776 myfree(ltestparam);
1777 myfree(hess[0]);
1778 myfree(hess);
1779 return panic;
1780 }
1781
1782
1783 MYREAL
interpolate(MYREAL xl,MYREAL xc,MYREAL xh,MYREAL low,MYREAL center,MYREAL high,MYREAL testlike)1784 interpolate (MYREAL xl, MYREAL xc, MYREAL xh, MYREAL low, MYREAL center,
1785 MYREAL high, MYREAL testlike)
1786 {
1787 MYREAL xc2, xl2, xh2, a, b, c, x;
1788
1789 xc2 = xc * xc;
1790 xl2 = xl * xl;
1791 xh2 = xh * xh;
1792
1793 a =
1794 (-(high * xc) + low * xc + center * xh - low * xh - center * xl +
1795 high * xl) / ((-xc + xh) * (xh - xl) * (-xc + xl));
1796
1797 b =
1798 (high * xc2 - low * xc2 - center * xh2 + low * xh2 + center * xl2 -
1799 high * xl2) / ((-xc + xh) * (xc - xl) * (-xh + xl));
1800
1801 c =
1802 -((low * xc2 * xh - low * xc * xh2 - high * xc2 * xl + center * xh2 * xl +
1803 high * xc * xl2 - center * xh * xl2) / ((-xc + xh) * (xc - xl) * (xh -
1804 xl)));
1805
1806
1807 x = (-b - sqrt (b * b - 4 * a * (c - testlike))) / (2 * a);
1808
1809 return x;
1810 }
1811
1812
1813 void
allocate_profile_percentiles(world_fmt * world)1814 allocate_profile_percentiles (world_fmt * world)
1815 {
1816 long i;
1817 long len = world->numpop2 + (world->options->gamma ? 1 : 0);
1818 #ifdef LONGSUM
1819 len += world->numpop*3;
1820 #endif /*LONGSUM*/
1821
1822 world->quantiles =
1823 (quantile_fmt *) mycalloc (1, sizeof (quantile_fmt) *len);
1824 world->percentile_failed = (boolean **) mycalloc (len, sizeof (boolean *));
1825 for (i = 0; i < len; i++)
1826 {
1827 world->quantiles[i].name = (char *) mycalloc (1, sizeof (char) * 20);
1828 world->quantiles[i].param =
1829 (MYREAL *) mycalloc (1, sizeof (MYREAL) * GRIDSIZE);
1830 world->percentile_failed[i] =
1831 (boolean *) mycalloc (GRIDSIZE, sizeof (boolean));
1832 }
1833 }
1834
1835 ///
1836 /// free profile related quantitites: quantiles, and indicator for failed
1837 /// convergence to the percentiles
1838 void
destroy_profile_percentiles(world_fmt * world)1839 destroy_profile_percentiles (world_fmt * world)
1840 {
1841 long i;
1842 long len = world->numpop2 + (world->options->gamma ? 1 : 0);
1843 #ifdef LONGSUM
1844
1845 len += world->numpop*3;
1846 #endif /*LONGSUM*/
1847
1848 for (i = 0; i < len; i++)
1849 {
1850 myfree(world->quantiles[i].name);
1851 myfree(world->quantiles[i].param);
1852 myfree(world->percentile_failed[i]);
1853 }
1854 myfree(world->quantiles);
1855 myfree(world->percentile_failed);
1856 }
1857
1858
1859 void
print_menu_profile(long which,long nn)1860 print_menu_profile (long which, long nn)
1861 {
1862 char * nowstr;
1863 nowstr = (char*) mycalloc(LINESIZE,sizeof(char));
1864
1865 get_time (nowstr, "%H:%M:%S");
1866 FPRINTF (stdout,
1867 "%s Calculating profile likelihoods for parameter %2li out of %2li\n",
1868 nowstr, 1 + which, nn);
1869 myfree(nowstr);
1870 }
1871
1872
1873 long
set_df(long which,char * custm2,long numpop,long numpop2)1874 set_df (long which, char *custm2, long numpop, long numpop2)
1875 {
1876 char ch;
1877 long i;
1878 long zeros = 0;
1879 long ms = 0;
1880 long ss = 0;
1881
1882 for (i = 0; i < numpop2; i++)
1883 {
1884 ch = custm2[i];
1885 switch (ch)
1886 {
1887 case '0':
1888 zeros++;
1889 break;
1890 case 'm':
1891 ms++;
1892 break;
1893 case 's':
1894 ss++;
1895 break;
1896 }
1897 }
1898 if (numpop2 - zeros == ms)
1899 return 2;
1900 else
1901 {
1902 if (custm2[0] == 'm')
1903 return numpop2 - numpop + 1 - zeros - ss / 2;
1904 else
1905 return numpop2 - zeros - ss / 2;
1906 }
1907 return -1;
1908 }
1909
1910 #ifdef PRIVATE
master_gridder(world_fmt * world,long * gmaxptr)1911 void master_gridder(world_fmt *world, long *gmaxptr)
1912 {
1913 FILE *alldim;
1914 nr_fmt *nr;
1915 long rep = world->options->replicate ? (world->loci > 1 ? 0 : world->repstop) : 0;
1916 MYREAL *mlparam;
1917 MYREAL *testparam;
1918 MYREAL *ltestparam;
1919 nr = (nr_fmt *) mycalloc (1, sizeof (nr_fmt));
1920 create_nr (nr, world, *gmaxptr, 0, world->loci,
1921 world->repkind, world->rep);
1922 testparam = (MYREAL *) mycalloc ( nr->partsize, sizeof (MYREAL));
1923 ltestparam = (MYREAL *) mycalloc ( nr->partsize, sizeof (MYREAL));
1924
1925 // specify the ML parameters
1926 if (world->loci > 1)
1927 {
1928 mlparam = world->atl[rep][world->loci].param;
1929 //maxlike = world->atl[rep][world->loci].param_like;
1930 }
1931 else
1932 {
1933 mlparam = world->atl[rep][world->locus].param;
1934 //maxlike = world->atl[rep][world->locus].param_like;
1935 }
1936 memcpy (testparam, mlparam, sizeof (MYREAL) *
1937 (world->numpop2 + (long) world->options->gamma));
1938 set_logparam (ltestparam, testparam, nr->partsize);
1939 alldim = fopen("dimensions","w+");
1940 calc_grid(alldim, world, testparam, ltestparam, mlparam, nr, 0);
1941 fclose(alldim);
1942 myfree(testparam);
1943 myfree(ltestparam);
1944 destroy_nr(nr, world);
1945 }
1946
calc_grid(FILE * alldim,world_fmt * world,MYREAL * testparam,MYREAL * ltestparam,MYREAL * mlparam,nr_fmt * nr,long which)1947 void calc_grid(FILE *alldim, world_fmt * world, MYREAL *testparam, MYREAL *ltestparam, MYREAL *mlparam, nr_fmt *nr, long which)
1948 {
1949 long i;
1950 MYREAL deviate[] = DEVIATE2;
1951 //MYREAL maxlike;
1952 helper_fmt helper;
1953 MYREAL saveparam, lsaveparam;
1954 // for(j=0; j<nr->numpop2; j++)
1955 // {
1956 if(which<world->numpop2)
1957 {
1958 saveparam = testparam[which];
1959 lsaveparam = ltestparam[which];
1960 for(i=0; i<GRIDSIZE; i++)
1961 {
1962 testparam[which] = deviate[i] * mlparam[which];
1963 ltestparam[which] = log(testparam[which]);
1964 if(which==nr->numpop2-1)
1965 {
1966 if(world->options->gamma)
1967 initgammacat (nr->categs, testparam[nr->numpop2],1., nr->rate, nr->probcat);
1968 fill_helper (&helper, testparam, ltestparam, world, nr);
1969 //FPRINTF(alldim,"%li %li ",which,i);
1970 //for(z=0;z<world->numpop2;z++)
1971 // FPRINTF(alldim, "%g ", testparam[z]);
1972 FPRINTF(alldim, "%f ", CALCLIKE (&helper, testparam, ltestparam));
1973 }
1974 calc_grid(alldim, world,testparam,ltestparam, mlparam, nr,which+1);
1975
1976 }
1977 FPRINTF(alldim,"\n");
1978 testparam[which] = saveparam;
1979 ltestparam[which] =lsaveparam;
1980 }
1981 // }
1982 // FPRINTF(world->mathfile, "\n");
1983 }
1984 #endif
1985
1986
1987
1988