1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate and effectice population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6 Likelihood ratio test R O U T I N E S
7
8 moved out from world.c
9 Peter Beerli 2000, Seattle
10 beerli@fsu.edu
11
12 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
13 Copyright 2003-2004 Peter Beerli, Tallahassee FL
14
15 Permission is hereby granted, free of charge, to any person obtaining
16 a copy of this software and associated documentation files (the
17 "Software"), to deal in the Software without restriction, including
18 without limitation the rights to use, copy, modify, merge, publish,
19 distribute, sublicense, and/or sell copies of the Software, and to
20 permit persons to whom the Software is furnished to do so, subject
21 to the following conditions:
22
23 The above copyright notice and this permission notice shall be
24 included in all copies or substantial portions of the Software.
25
26 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
27 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
28 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
29 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
30 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
31 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
32 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
33
34 $Id: lrt.c 2158 2013-04-29 01:56:20Z beerli $
35
36 -------------------------------------------------------*/
37 /* \file lrt.c
38 Calculates likelihood ratio tests
39 */
40 #include "migration.h"
41 #include "sighandler.h"
42 #include "tools.h"
43 #include "broyden.h"
44 #include "combroyden.h"
45 #include "options.h"
46 #include "aic.h"
47 #include "migrate_mpi.h"
48 #include "pretty.h"
49
50 #ifndef LAGUERRE
51 //#include "derivatives2.h"
52 #endif
53 #include "sort.h"
54
55 //#ifdef UEP
56 //#include "uep.h"
57 //#endif
58 #ifdef DMALLOC_FUNC_CHECK
59 #include <dmalloc.h>
60 #endif
61
62 void print_lratio_test (world_fmt * world, long *Gmax);
63 void test_loci_like (nr_fmt * nr, MYREAL *param0,
64 MYREAL *param1, long df, long zeros,
65 long loci, world_fmt * world,
66 long *maxwhich, long maxnum,
67 boolean withhead, char *this_string1);
68 long set_test_param (MYREAL *param, lr_data_fmt *lrtdata, world_fmt * world,
69 long lrline, long locus, long *maxwhich,
70 long *maxnum, long *zeros);
71 void print_lrt_box(world_fmt *world, MYREAL *param0, MYREAL *param1, long zeros, long elem,
72 char * this_string1, MYREAL like0, MYREAL like1, long df);
73
74 long parse_h0(char ***box, long *size, char *thisString, MYREAL *param0, MYREAL *param1 , long elem, world_fmt *world);
75 void parse_h0part(char ***box, long *size, long *tempsize,
76 MYREAL *param, long elem, char head[], world_fmt *world);
77 void parse_h0string(char ***box, long *size, long *tempsize,
78 char *thisString);
79
80 void print_box(world_fmt* world,char **box,long size,
81 MYREAL like0, MYREAL like1,MYREAL lrt, MYREAL chiprob, MYREAL chiprob2, MYREAL aic, long df, long aicparamnum);
82
83 long simplify_lrtvalues(char * in, char** out);
84
85 boolean
mywhitespace(char ch)86 mywhitespace (char ch)
87 {
88 if (!(ch == ',' || ch == '\0' || ch == '\t' || ch == '\r' || ch == '\n' || ch == ';'))
89 {
90 return TRUE;
91 }
92 return FALSE;
93 }
94
95
96 void
print_lratio_test(world_fmt * world,long * Gmax)97 print_lratio_test (world_fmt * world, long *Gmax)
98 {
99 long c;
100 long r, locus;
101 long df;
102 long zeros;
103 int header;
104 nr_fmt *nr;
105 MYREAL *param0;
106 MYREAL *param1;
107 long *maxwhich;
108 long nparam;
109 long maxnum = 0;
110 long rep = !world->options->replicate ? 0 :
111 (world->options->replicatenum == 0 ?
112 world->options->lchains : world->options->replicatenum);
113 long repstop = world->repstop;
114 param0 = (MYREAL *) mycalloc (1, sizeof (MYREAL) * (world->numpop2 + 1));
115 param1 = (MYREAL *) mycalloc (1, sizeof (MYREAL) * (world->numpop2 + 1));
116 maxwhich = (long *) mycalloc (1, sizeof (long) * (world->numpop2 + 1));
117
118
119 if (world->options->progress)
120 FPRINTF (stdout, " Printing likelihood ratio tests\n");
121
122 nr = (nr_fmt *) mycalloc (1, sizeof (nr_fmt));
123 create_nr (nr, world, *Gmax, 0L, world->loci, world->repkind, world->rep);
124 for (locus = 0; locus < world->loci; locus++)
125 {
126 if (world->options->replicate)
127 for (r = 0; r < repstop; r++)
128 create_multiapg0 (nr->apg0[r][locus], nr, r, locus);
129 else
130 create_apg0 (nr->apg0[0][locus], nr, &world->atl[0][locus], locus);
131 }
132 if (world->loci > 1)
133 {
134 locus = world->loci;
135 rep = 0;
136 }
137 else
138 {
139 world->locus=0;
140 locus = 0;
141 }
142 PAGEFEEDWORLD;
143 nparam = world->options->gamma ? world->numpop2 + 1 : world->numpop2;
144 for (c = 0; c < world->options->lratio->counter; c++)
145 {
146 header = (c == 0) ? HEADER : NOHEADER;
147 if (world->options->lratio->data[c].type == MLE)
148 {
149 memcpy (param1, world->atl[rep][locus].param,
150 sizeof (MYREAL) * nparam);
151 df = set_test_param (param0,
152 &world->options->lratio->data[c],
153 world, 0L, -1L, maxwhich, &maxnum, &zeros);
154 test_loci_like (nr, param0, param1,
155 df, zeros, world->loci, world, maxwhich,
156 maxnum, header,
157 world->options->lratio->data[c].value1);
158 }
159 }
160 fflush (world->outfile);
161 myfree(param0);
162 myfree(param1);
163 myfree(maxwhich);
164 destroy_nr (nr, world);
165 }
166
167
168 //remember: param0 is the parameterset to test and
169 // NOT the parameterset from migrate.
170 #define BOXSIZE 50 /*defines the print width of the left box containing the H0*/
171 #define BOXSIZE2 40 /*defines the print width of the left box with the legend*/
172 void
test_loci_like(nr_fmt * nr,MYREAL * param0,MYREAL * param1,long df,long zeros,long loci,world_fmt * world,long * maxwhich,long maxnum,boolean withhead,char * this_string1)173 test_loci_like (nr_fmt * nr, MYREAL *param0, MYREAL *param1, long df,
174 long zeros, long loci, world_fmt * world, long *maxwhich,
175 long maxnum, boolean withhead, char *this_string1)
176 {
177
178 // char *teststat, temp[LRATIO_STRINGS];
179 MYREAL like1, like0;//, testval, chiprob, chiprob2;
180 // int length;
181 //
182 long numparam;
183 long i, j, g = 0;
184 long elem;
185 //long zi;
186 long z = 0, w = 0;
187 // long pop;
188 //MYREAL normd = 0.0;
189 long *which;
190 MYREAL **hess;
191 MYREAL /**values,*/ *saveparam0;
192 //MYREAL tparam;
193 // long spaces = 0;
194 char temp[100];
195 helper_fmt helper;
196 MYREAL *lparam0;
197 MYREAL *lparam1;
198 MYREAL aicfull;
199 long aicfullparamnum;
200 // char *message;
201 numparam = world->numpop2 + 1;
202 doublevec2d(&hess,numparam, numparam);
203 lparam0 = (MYREAL *) mycalloc (numparam, sizeof (MYREAL));
204 lparam1 = (MYREAL *) mycalloc (numparam, sizeof (MYREAL));
205 which = (long *) mycalloc (1, sizeof (long) * numparam);
206 // message = (char *) mycalloc (20, sizeof (char));
207 // values = (MYREAL *) mycalloc (1, sizeof (MYREAL) * (world->numpop2 + 1));
208 saveparam0 = (MYREAL *) mymalloc (sizeof (MYREAL) * numparam);
209 memcpy (saveparam0, param0, sizeof (MYREAL) * numparam);
210 for (i = 0; i < loci; i++)
211 {
212 for (j = 0; j < world->repstop; j++)
213 {
214 if (g < world->atl[j][i].T)
215 g = world->atl[j][i].T;
216 }
217 }
218 elem = world->options->gamma ? numparam : nr->numpop2;
219 nr->skiploci = world->data->skiploci;
220 helper.multilocus = world->loci == 1 ? FALSE : TRUE;
221 if (maxnum > 0)
222 {
223 for (i = 0; i < elem; i++)
224 {
225 if (i != maxwhich[z])
226 {
227 which[w] = i;
228 nr->values[w++] = param0[i];
229 }
230 else
231 {
232 if (maxnum > z + 1)
233 z++;
234 }
235 }
236 nr->profilenum = w;
237 maximize (¶m0, world, nr, hess, PROFILE, world->repkind);
238 like0 = nr->llike;
239 //xcode normd = nr->normd;
240 memcpy (param0, world->param0, sizeof (MYREAL) * nr->partsize);
241 }
242 else
243 {
244 set_logparam (lparam0, param0, elem);
245 fill_helper (&helper, param0, lparam0, world, nr);
246 like0 = CALCLIKE (&helper, param0, lparam0);
247 }
248 set_logparam (lparam1, param1, elem);
249 fill_helper (&helper, param1, lparam1, world, nr);
250 like1 = CALCLIKE (&helper, param1, lparam1);
251 // sprintf (message, " %f ", chiprob2);
252
253 if (withhead)
254 {
255 FPRINTF (world->outfile,
256 "==============================================================================\n");
257 FPRINTF (world->outfile, "Likelihood ratio tests\n");
258 FPRINTF (world->outfile,
259 "==============================================================================\n");
260 FPRINTF (world->outfile, "Over all loci\n");
261 FPRINTF (world->outfile, "Legend for the LRT tables\n");
262 print_line(world->outfile,'-',79L,CONT);
263 sprintf (temp,"Null-Hypothesis: your test model");
264 FPRINTF(world->outfile,"%-*.*s | Log(likelihood) of test model\n",BOXSIZE2,BOXSIZE2,temp);
265 sprintf (temp,"=same=");
266 FPRINTF(world->outfile,"%-*.*s | Log(likelihood) of full model\n",BOXSIZE2,BOXSIZE2, temp);
267 sprintf (temp,"full model (the model under which the");
268 FPRINTF(world->outfile,"%-*.*s | Likelihood ratio test value\n",BOXSIZE2,BOXSIZE2, temp);
269 sprintf (temp,"genealogies were sampled)");
270 FPRINTF(world->outfile,"%-*.*s | Degrees of freedom of test\n",BOXSIZE2,BOXSIZE2,temp);
271 sprintf (temp,"[Theta values are on the diagonal of the ");
272 FPRINTF(world->outfile,"%-*.*s | Probability*\n",BOXSIZE2,BOXSIZE2,temp);
273 sprintf (temp,"Migration matrix, migration rates are ");
274 FPRINTF(world->outfile,"%-*.*s | Probability**\n",BOXSIZE2,BOXSIZE2,temp);
275 sprintf (temp,"specified as %s]", world->options->usem ? "M" : "Theta * M");
276 FPRINTF(world->outfile,"%-*.*s | Akaike's Information Criterion***\n",BOXSIZE2,BOXSIZE2,temp);
277 sprintf (temp," ");
278 FPRINTF(world->outfile,"%-*.*s | Number of parameters used\n",BOXSIZE2,BOXSIZE2,temp);
279 print_line(world->outfile,'-',79L,CONT);
280 FPRINTF(world->outfile," *) Probability under the assumption that parameters have range -Inf to Inf\n");
281 FPRINTF(world->outfile," **) Probability under the assumption that parameters have range 0 to Inf\n");
282 FPRINTF(world->outfile,"***) AIC: the smaller the value the better the model\n");
283 aicfullparamnum = find_paramnum(world,NULL);
284 aicfull = -2. * like1 + 2. * aicfullparamnum;
285 FPRINTF(world->outfile," [the full model has AIC=%f, num(param)=%li]\n\n",aicfull,aicfullparamnum);
286
287 pdf_print_lrt_header(aicfull,aicfullparamnum, &world->page_width, &world->page_height);
288 }
289 print_lrt_box(world,param0, param1, zeros, elem, this_string1, like0, like1, df);
290 myfree(lparam0);
291 myfree(lparam1);
292 myfree(saveparam0);
293 myfree(which);
294 myfree(hess[0]);
295 myfree(hess);
296 }
297
print_lrt_box(world_fmt * world,MYREAL * param0,MYREAL * param1,long zeros,long elem,char * this_string1,MYREAL like0,MYREAL like1,long df)298 void print_lrt_box(world_fmt *world, MYREAL *param0, MYREAL *param1, long zeros, long elem,
299 char * this_string1, MYREAL like0, MYREAL like1, long df)
300 {
301 // Alternative plot
302 //
303 //---------------------------------------------------------------------
304 // H0: (p1,p2,p3,..., | LRT = -2(val1 - val2) = val
305 // pi,....,pn) == (v1,v2, | df = x
306 // v3,.... vi,...vn) | Prob = x.xx
307 // | Probc= x.xx
308 // | AIC = x.xxx
309 //---------------------------------------------------------------------
310 //
311 // LRT
312 MYREAL testval = -2. * (like0 - like1);
313 // standard probability assuming indpendence and range of -inf .. +inf
314 MYREAL chiprob = probchi (df, testval);
315 // probability assuming independence and range of 0 .. inf
316 MYREAL chiprob2 = probchiboundary (testval, zeros, df);
317 // AIC value penalizing for number of paramters
318 long aicparamnum = find_paramnum(world,this_string1);
319 MYREAL aic = -2. * like0 + 2. * aicparamnum;
320
321 char **box;
322 long size = HUNDRED; // we need at least 7 lines to print the right side of the table
323 /// \todo reallocation of likelihood ratio boxes needs reevaluation
324 long newsize;
325 long i;
326
327 if(world->options->progress)
328 print_line(stdout,'-',79L,CONT);
329 print_line(world->outfile,'-',79L,CONT);
330 box = (char **) mycalloc (size, sizeof (char *));
331 for(i=0; i<size; i++)
332 box[i] = (char *) mycalloc(LRATIO_STRINGS,sizeof(char));
333 newsize = parse_h0(&box,&size, this_string1,param0,param1,elem,world);
334 print_box(world,box,newsize,like0,like1,testval, chiprob, chiprob2, aic, df, aicparamnum);
335 pdf_print_LRT_box(world,box,newsize,like0,like1,testval, chiprob, chiprob2, aic, df, aicparamnum, &world->page_height, & world->page_width);
336 if(world->options->progress)
337 print_line(stdout,'-',79L,CONT);
338 print_line(world->outfile,'-',79L,CONT);
339 for(i=0; i<size; i++)
340 myfree(box[i]);
341 myfree(box);
342 }
343
parse_h0(char *** box,long * size,char * thisString,MYREAL * param0,MYREAL * param1,long elem,world_fmt * world)344 long parse_h0(char ***box, long *size, char *thisString, MYREAL *param0, MYREAL *param1 , long elem, world_fmt *world)
345 {
346 long tempsize=0;
347 parse_h0part(box, size, &tempsize, param0, elem, "H0:", world);
348 tempsize++;
349 parse_h0part(box, size, &tempsize, param1, elem, " = ", world);
350 tempsize++;
351 parse_h0string(box,size,&tempsize,thisString);
352 return tempsize;
353 }
354
parse_h0part(char *** box,long * size,long * tempsize,MYREAL * param,long elem,char head[],world_fmt * world)355 void parse_h0part(char ***box, long *size, long *tempsize,
356 MYREAL *param, long elem, char head[], world_fmt *world)
357 {
358 int length;
359 long i;
360 long zi;
361 long pop;
362 MYREAL tparam;
363 char tmp[LRATIO_STRINGS];
364 long counter = 0;
365 //xcode length = (int) MAX (0, 5 - sprintf (tmp, "%i", ((int) param[0])));
366 counter = sprintf ((*box)[*tempsize],"%3.3s", head);
367 for (i = 0; i < elem; i++)
368 {
369 zi = mml2m (i, world->numpop);
370 if (zi >= world->numpop)
371 {
372 pop = (zi - world->numpop) / (world->numpop - 1);
373 if(world->options->usem)
374 tparam = param[zi];
375 else
376 tparam = param[zi] * param[pop];
377 }
378 else
379 tparam = param[zi];
380 length = MAX (0, 5 - sprintf (tmp, "%i", ((int) tparam)));
381 counter += sprintf ((*box)[*tempsize] + counter, " %.*f", length, tparam);
382 if (counter + length > BOXSIZE)
383 {
384 (*tempsize)++;
385 if(*tempsize>= *size-1)
386 {
387 *box = (char **) myrealloc(*box,sizeof(char*) * (*tempsize+1));
388 (*box)[*tempsize] = (char *) mycalloc(LRATIO_STRINGS,sizeof(char));
389 *size = *tempsize+1;
390 }
391 counter = sprintf((*box)[*tempsize],"%3.3s"," ");
392 }
393 }
394 }
395
parse_h0string(char *** box,long * size,long * tempsize,char * thisString)396 void parse_h0string(char ***box, long *size, long *tempsize,
397 char *thisString)
398 {
399 long i;
400 long counter = 0;
401 long strsize = (long) strlen(thisString);
402 counter = sprintf ((*box)[*tempsize],"[");
403 for (i = 0; i < strsize; i++)
404 {
405 counter += sprintf ((*box)[*tempsize] + counter, "%c",thisString[i]);
406 if (counter >= BOXSIZE)
407 {
408 (*tempsize)++;
409 if(*tempsize>= *size-1)
410 {
411 *box = (char **) myrealloc(*box,sizeof(char*) * (*tempsize+1));
412 (*box)[*tempsize] = (char *) mycalloc(LRATIO_STRINGS,sizeof(char));
413 *size = *tempsize+1;
414 }
415 counter = sprintf((*box)[*tempsize]," ");
416 }
417 }
418 sprintf ((*box)[*tempsize] + counter, "]");
419 }
420
print_box(world_fmt * world,char ** box,long size,MYREAL like0,MYREAL like1,MYREAL lrt,MYREAL chiprob,MYREAL chiprob2,MYREAL aic,long df,long aicparamnum)421 void print_box(world_fmt* world,char **box,long size,
422 MYREAL like0, MYREAL like1,MYREAL lrt, MYREAL chiprob, MYREAL chiprob2, MYREAL aic, long df, long aicparamnum)
423 {
424 long i;
425 FPRINTF(world->outfile,"%-*.*s | LnL(test) = %f\n",BOXSIZE,BOXSIZE, box[0],like0);
426 FPRINTF(world->outfile,"%-*.*s | LnL(full) = %f\n",BOXSIZE,BOXSIZE, box[1],like1);
427 FPRINTF(world->outfile,"%-*.*s | LRT = %f\n",BOXSIZE,BOXSIZE, box[2],lrt);
428 FPRINTF(world->outfile,"%-*.*s | df = %li\n",BOXSIZE,BOXSIZE,box[3],(long) df);
429 FPRINTF(world->outfile,"%-*.*s | Prob = %f\n",BOXSIZE,BOXSIZE,box[4],chiprob);
430 FPRINTF(world->outfile,"%-*.*s | Probc = %f\n",BOXSIZE,BOXSIZE,box[5],chiprob2);
431 FPRINTF(world->outfile,"%-*.*s | AIC = %f\n",BOXSIZE,BOXSIZE,box[6],aic);
432 FPRINTF(world->outfile,"%-*.*s | num(param)= %li\n",BOXSIZE,BOXSIZE,box[7],aicparamnum);
433 for(i=8; i< size; i++)
434 FPRINTF(world->outfile,"%-*.*s |\n",BOXSIZE,BOXSIZE,box[i]);
435 if(world->options->progress)
436 {
437 FPRINTF(stdout,"%-*.*s | LnL(test) = %f\n",BOXSIZE,BOXSIZE, box[0],like0);
438 FPRINTF(stdout,"%-*.*s | LnL(full) = %f\n",BOXSIZE,BOXSIZE, box[1],like1);
439 FPRINTF(stdout,"%-*.*s | LRT = %f\n",BOXSIZE,BOXSIZE, box[2],lrt);
440 FPRINTF(stdout,"%-*.*s | df = %li\n",BOXSIZE,BOXSIZE,box[3],(long) df);
441 FPRINTF(stdout,"%-*.*s | Prob = %f\n",BOXSIZE,BOXSIZE,box[4],chiprob);
442 FPRINTF(stdout,"%-*.*s | Probc = %f\n",BOXSIZE,BOXSIZE,box[5],chiprob2);
443 FPRINTF(stdout,"%-*.*s | AIC = %f\n",BOXSIZE,BOXSIZE,box[6],aic);
444 FPRINTF(stdout,"%-*.*s | num(param)= %li\n",BOXSIZE,BOXSIZE,box[7],aicparamnum);
445 for(i=8; i< size; i++)
446 FPRINTF(stdout,"%-*.*s |\n",BOXSIZE,BOXSIZE,box[i]);
447 }
448 }
449
450
451 long
set_test_param(MYREAL * param,lr_data_fmt * lrtdata,world_fmt * world,long lrline,long locus,long * maxwhich,long * maxnum,long * zeros)452 set_test_param (MYREAL *param, lr_data_fmt *lrtdata, world_fmt * world, long lrline,
453 long locus, long *maxwhich, long *maxnum, long *zeros)
454 {
455 long elem = world->options->gamma ? world->numpop2 + 1 : world->numpop2;
456 long elements;
457 long repstop = !world->options->replicate ? 0 :
458 (world->options->replicatenum == 0 ?
459 world->options->lchains : world->options->replicatenum);
460 char *paramtype;
461 char *ss;
462 char **custm;
463 MYREAL *meanparam;
464 //long zzz;
465 long zi;
466 long el;
467 long df = 0;
468 long count=0;
469 MYREAL mean=0.;
470 long offset;
471 long limit;
472 long zz;
473 long numpop = world->numpop;
474 long pop1, pop2;
475
476 paramtype = (char *) mycalloc (1, sizeof (char) * elem);
477 ss = (char *) mycalloc (LONGLINESIZE, sizeof (char));
478 charvec2d(&custm,world->numpop2 + 1, LONGLINESIZE);
479 strcpy (ss, lrtdata->value1);
480
481 elements = simplify_lrtvalues(ss,custm);
482
483 *zeros = 0;
484
485 if(elem != elements)
486 warning("Not enough elements in the l-ratio specification\n");
487
488 //xcode zzz=0;
489 if (world->loci - world->skipped > 1)
490 meanparam = world->atl[0][world->loci].param;
491 else
492 meanparam = world->atl[repstop][0].param;
493
494 for(el=0; el < elements; el++)
495 {
496 zi = mml2m (el, numpop);
497 m2mm (zi, numpop, &pop1, &pop2);
498
499 switch (custm[el][0])
500 {
501 case 'x':
502 paramtype[zi] = '-';
503 lrtdata->connect[zi]='x';
504 param[zi] = meanparam[zi];
505 maxwhich[(*maxnum)++] = zi;
506 df++;
507 break;
508 case '*':
509 paramtype[zi] = '-';
510 lrtdata->connect[zi]='*';
511 param[zi] = meanparam[zi];
512 break;
513 case 't':
514 case 'm':
515 paramtype[zi] = '-';
516 lrtdata->connect[zi]='m';
517 if(custm[el][0] != world->options->custm2[zi])
518 df++;
519 mean = 0.0;
520 count = 0;
521 offset = (zi >= world->numpop) ? world->numpop : 0;
522 limit = (zi >= world->numpop) ? world->numpop2 : world->numpop;
523 for (zz = offset; zz < limit; zz++)
524 {
525 if(custm[m2mml(zz,numpop)][0] == 'm')
526 {
527 mean += meanparam[zz];
528 count++;
529 }
530 }
531 mean /= count; //limit - offset;
532 param[zi] = mean;
533 break;
534 case 'M':
535 if(custm[el][0] != world->options->custm2[zi])
536 df++;
537 mean = 0.0;
538 count = 0;
539 offset = (zi >= world->numpop) ? world->numpop : 0;
540 limit =
541 (zi >= world->numpop) ? world->numpop2 : world->numpop;
542 if (offset < numpop)
543 {
544 paramtype[zi] = '-';
545 lrtdata->connect[zi]='m';
546 for (zz = offset; zz < limit; zz++)
547 {
548 if(custm[m2mml(zz,numpop)][0] == 'm')
549 {
550 mean += meanparam[zz];
551 count++;
552 }
553 }
554 mean /= count; //limit - offset;
555 param[zi] = mean;
556 }
557 else
558 {
559 paramtype[zi] = '+';
560 lrtdata->connect[zi]='M';
561 for (zz = offset; zz < limit; zz++)
562 {
563 if(custm[m2mml(zz,numpop)][0] == 'm')
564 {
565 m2mm (zz, numpop, &pop1, &pop2);
566 mean += meanparam[zz] * meanparam[pop2];
567 count++;
568 }
569 }
570 mean /= count;
571 param[zi] = mean;
572 }
573 break;
574 case 's':
575 if(custm[el][0] != world->options->custm2[zi])
576 df++;
577 if (zi < world->numpop)
578 {
579 paramtype[zi] = '-';
580 lrtdata->connect[zi]='*';
581 param[zi] = meanparam[zi];
582 }
583 else
584 {
585 paramtype[zi] = '-';
586 lrtdata->connect[zi]='s';
587 param[zi] =
588 (meanparam[zi] +
589 meanparam[mm2m (pop2, pop1, world->numpop)]) / 2.;
590 }
591 break;
592 case 'S':
593 if(custm[el][0] != world->options->custm2[zi])
594 df++;
595 if (zi < world->numpop)
596 {
597 paramtype[zi] = '-';
598 lrtdata->connect[zi]='*';
599 param[zi] = meanparam[zi];
600 }
601 else
602 {
603 paramtype[zi] = '+';
604 lrtdata->connect[zi]='S';
605 param[zi] =
606 (meanparam[zi] * meanparam[pop2] +
607 meanparam[mm2m (pop2, pop1, world->numpop)] *
608 meanparam[pop1]) / 2.;
609 }
610 break;
611 default:
612 paramtype[zi] = '+';
613 if(custm[el][0] != world->options->custm2[zi])
614 df++;
615 if (custm[el][0] == '0' && (long) strlen(custm[el])==1)
616 {
617 if(custm[el][0] != world->options->custm2[zi])
618 (*zeros)++;
619 lrtdata->connect[zi]='0';
620 }
621 else
622 lrtdata->connect[zi]='c';
623 param[zi] = MAX (atof (custm[el]), SMALLEST_THETA);
624 break;
625 }
626
627 }
628 for (zi = world->numpop; zi < world->numpop2; zi++)
629 {
630 if (paramtype[zi] == '+')
631 {
632 zz = (zi - world->numpop) / (world->numpop - 1);
633 if(!world->options->usem)
634 param[zi] /= param[zz];
635 }
636 }
637 myfree(paramtype);
638 myfree(ss);
639 free_charvec2d(custm);
640 return df;
641 }
642
643
simplify_lrtvalues(char * in,char ** out)644 long simplify_lrtvalues(char *in, char **out)
645 {
646 char *ptr;
647 long elements=0;
648 while(in != NULL)
649 {
650
651 ptr = strsep(&in,"\n\r\t; {,}");
652 if(strlen(ptr)>0)
653 {
654 if(ptr[0]=='0' && ((atof(ptr) - SMALLEPSILON) <= 0.0))
655 sprintf(out[elements++], "0");
656 else
657 sprintf(out[elements++],"%s", ptr);
658 }
659 }
660 return elements;
661 }
662
663 // in is the string from the l-ratio option, custom is a database array that
664 // gets allocated while its filled
665 // maximal event numbers is LONGLINESIZE [search for use of simplify_lrtvalues()]
simplify_lrtvalues_old(char * in,char *** out)666 long simplify_lrtvalues_old(char * in, char*** out)
667 {
668 long elements=0;
669 long count = 0;
670 boolean white=FALSE;
671 boolean element_ready = TRUE;
672 // (*out)[0] = (char *) mycalloc(LONGLINESIZE,sizeof(char));
673 while(*in!='\0')
674 {
675 switch(*in)
676 {
677 case '\r':
678 case '\n':
679 case '\t':
680 case ' ':
681 case ',':
682 case ';':
683 white=TRUE;
684 break; // jump over "whitespace",
685 // if elements are from the list {*, m, M, s, S} read one character and then restart a new element
686 case 't':
687 case 'm':
688 case 'M':
689 case 's':
690 case 'S':
691 case '*':
692 case 'x':
693 white = FALSE;
694 // fill element and reset counter within element
695 count=0;
696 (*out)[elements][0] = *in;
697 // make new element ready
698 elements++;
699 // (*out)[elements] = (char *) mycalloc(LONGLINESIZE,sizeof(char));
700 element_ready = TRUE;
701 break;
702 default:
703 if(!element_ready)
704 {
705 (*out)[elements][count++] = *in;
706 }
707 else
708 {
709 if(white)
710 {
711 count=0;
712 elements++;
713 // (*out)[elements] = (char *) mycalloc(LONGLINESIZE,sizeof(char));
714 element_ready = FALSE;
715 (*out)[elements][count++] = *in;
716 }
717 else
718 {
719 (*out)[elements][count++] = *in;
720 }
721 }
722 white = FALSE;
723 break;
724 }
725 in++;
726 }
727 //*out = (char**) myrealloc(*out, elements * sizeof(char*));
728 return elements;
729 }
730
731 #if 0
732 long
733 set_test_param (MYREAL *param, lr_data_fmt *lrtdata, world_fmt * world, long lrline,
734 long locus, long *maxwhich, long *maxnum, long *zeros)
735 {
736 long i = 0, z = 0, zi = 0, zz = 0, zzz = 0, df = 0;
737 long count;
738 long offset = 0, limit = 0, pop, pop1, pop2;
739 long numpop = world->numpop;
740 char *tmp, *tmp2, *ss, *custm, *paramtype;
741 MYREAL *meanparam, mean;
742 long elem = world->options->gamma ? world->numpop2 + 1 : world->numpop2;
743 long repstop = !world->options->replicate ? 0 :
744 (world->options->replicatenum == 0 ?
745 world->options->lchains : world->options->replicatenum);
746 *zeros = 0;
747 ss = (char *) mycalloc (1, sizeof (char) * LOGLINESIZE);
748 custm = (char *) mycalloc (1, sizeof (char) * LONGLINESIZE);
749 tmp = (char *) mycalloc (1, sizeof (char) * LONGLINESIZE);
750 paramtype = (char *) mycalloc (1, sizeof (char) * elem);
751
752 *maxnum = 0;
753 strcpy (ss, lrtdata->value1);
754 strcpy (custm, lrtdata->value1);
755 tmp2 = custm;
756 zzz = 0;
757 while(*tmp2!='\0')
758 {
759 if(strchr("tmMsS*xc",*tmp2))
760 custm[zzz++] = *tmp2;
761 if(strchr("0",*tmp2))
762 {
763 if(!strchr(".0123456789",*(tmp2+1)))
764 custm[zzz++] = *tmp2;
765 }
766 tmp2++;
767 }
768 custm[zzz]="\0";
769 zzz=0;
770 if (world->loci - world->skipped > 1)
771 meanparam = world->atl[0][world->loci].param;
772 else
773 meanparam = world->atl[repstop][0].param;
774 while (ss[zzz] != '\0')
775 {
776 tmp[i] = ss[zzz++];
777 if (mywhitespace (tmp[i]))
778 {
779 if (tmp[i] != ' ')
780 i++;
781 }
782 else
783 {
784 tmp[i] = '\0';
785 i = 0;
786 zi = mml2m (z, numpop);
787 m2mm (zi, numpop, &pop1, &pop2);
788 switch (tmp[0])
789 {
790 case 'x':
791 paramtype[zi] = '-';
792 lrtdata->connect[zi]='x';
793 param[zi] = meanparam[zi];
794 maxwhich[(*maxnum)++] = zi;
795 z++;
796 df++;
797 break;
798 case '*':
799 paramtype[zi] = '-';
800 lrtdata->connect[zi]='*';
801 param[zi] = meanparam[zi];
802 z++;
803 break;
804 case 't':
805 case 'm':
806 paramtype[zi] = '-';
807 lrtdata->connect[zi]='m';
808 zz = atol (tmp) - 1;
809 tmp2 = tmp;
810 if(tmp[0] != world->options->custm2[zi])
811 df++;
812 if (zz < 0)
813 {
814 mean = 0.0;
815 count = 0;
816 offset = (zi >= world->numpop) ? world->numpop : 0;
817 limit = (zi >= world->numpop) ? world->numpop2 : world->numpop;
818 for (zz = offset; zz < limit; zz++)
819 {
820 if(custm[m2mml(zz,world->numpop)] == 'm')
821 {
822 mean += meanparam[zz];
823 count++;
824 }
825 }
826 mean /= count; //limit - offset;
827 param[zi] = mean;
828 }
829 else
830 {
831 param[zi] = meanparam[zz];
832 }
833 z++;
834 break;
835 case 'M':
836 zz = atol (tmp) - 1;
837 df++;
838 if (zz < 0)
839 {
840 mean = 0.0;
841 offset = (zi >= world->numpop) ? world->numpop : 0;
842 limit =
843 (zi >= world->numpop) ? world->numpop2 : world->numpop;
844 if (offset < numpop)
845 {
846 paramtype[zi] = '-';
847 lrtdata->connect[zi]='m';
848 for (zz = offset; zz < limit; zz++)
849 mean += meanparam[zz];
850 mean /= limit - offset;
851 param[zi] = mean;
852 }
853 else
854 {
855 paramtype[zi] = '+';
856 lrtdata->connect[zi]='M';
857 for (zz = offset; zz < limit; zz++)
858 {
859 m2mm (zz, numpop, &pop1, &pop2);
860 mean += meanparam[zz] * meanparam[pop2];
861 }
862 mean /= limit - offset;
863 param[zi] = mean;
864 }
865 }
866 else
867 {
868 paramtype[zi] = '-';
869 lrtdata->connect[zi]='M';
870 param[zi] = meanparam[zz];
871 }
872 z++;
873 break;
874 case 's':
875 df++;
876 if (zi < world->numpop)
877 {
878 paramtype[zi] = '-';
879 lrtdata->connect[zi]='*';
880 param[zi] = meanparam[zi];
881 z++;
882 }
883 else
884 {
885 paramtype[zi] = '-';
886 lrtdata->connect[zi]='s';
887 param[zi] =
888 (meanparam[zi] +
889 meanparam[mm2m (pop2, pop1, world->numpop)]) / 2.;
890 z++;
891 }
892 break;
893 case 'S':
894 df++;
895 if (zi < world->numpop)
896 {
897 paramtype[zi] = '-';
898 lrtdata->connect[zi]='*';
899 param[zi] = meanparam[zi];
900 z++;
901 }
902 else
903 {
904 paramtype[zi] = '+';
905 lrtdata->connect[zi]='S';
906 param[zi] =
907 (meanparam[zi] * meanparam[pop2] +
908 meanparam[mm2m (pop2, pop1, world->numpop)] *
909 meanparam[pop1]) / 2.;
910 z++;
911 }
912 break;
913 default:
914 paramtype[zi] = '+';
915 df++;
916 if (tmp[0] == '0')
917 {
918 (*zeros)++;
919 lrtdata->connect[zi]='0';
920 }
921 else
922 lrtdata->connect[zi]='c';
923 param[zi] = MAX (atof (tmp), SMALLEST_THETA);
924 z++;
925 break;
926 }
927 }
928 }
929 for (zi = world->numpop; zi < world->numpop2; zi++)
930 {
931 if (paramtype[zi] == '+')
932 {
933 pop = (zi - world->numpop) / (world->numpop - 1);
934 if(!world->options->usem)
935 param[zi] /= param[pop];
936 }
937 }
938 myfree(paramtype);
939 myfree(ss);
940 myfree(tmp);
941 myfree(custm);
942 return df;
943 }
944
945 #endif
946