1 /*! \file aic.c */
2 /*------------------------------------------------------
3 Maximum likelihood estimation
4 of migration rate and effectice population size
5 using a Metropolis-Hastings Monte Carlo algorithm
6 -------------------------------------------------------
7 AIC model test R O U T I N E S
8
9
10 Copyright 1997-2002 Peter Beerli and Joseph Felsenstein
11 Copyright 2003-2008 Peter Beerli
12 Copyright 2009-2012 Peter Beerli and Michal Palczewski
13
14 Permission is hereby granted, free of charge, to any person obtaining a copy
15 of this software and associated documentation files (the "Software"), to deal
16 in the Software without restriction, including without limitation the rights
17 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
18 of the Software, and to permit persons to whom the Software is furnished to do
19 so, subject to the following conditions:
20
21 The above copyright notice and this permission notice shall be included in all copies
22 or substantial portions of the Software.
23
24 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
25 INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
26 PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
27 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
28 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
29 OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
30
31 $Id: aic.c 2067 2012-07-27 20:59:32Z beerli $
32 -------------------------------------------------------*/
33
34 #include "migration.h"
35 #include "tools.h"
36 #include "broyden.h"
37 #include "combroyden.h"
38 #include "options.h"
39 #include "migrate_mpi.h"
40 #include "aic.h"
41 #include "sighandler.h"
42
43 #ifdef DMALLOC_FUNC_CHECK
44 #include <dmalloc.h>
45 #endif
46
47
48 static void print_progressheader_aic(boolean progress, FILE *file, MYREAL mleaic, MYREAL numparam);
49
50 static void aic_score (aic_fmt ** aicvec, long *aicnum, nr_fmt * nr,
51 long zero, long which, char *temppattern, MYREAL *param0,
52 int migtype);
53
54 static boolean legal_pattern (char *matrix, long numpop);
55
56 void fast_aic_score (aic_fmt ** aicvec, long *aicnum, nr_fmt * nr,
57 long zero, long which, char *temppattern,
58 MYREAL *param0, int migtype);
59
60 static void add_aiclist (int migtype, long numparam, long freeparam, long remainnum, long zero, long m,
61 MYREAL *param0, char *temppattern, char *custm2,
62 long *aicnum, aic_fmt ** aicvec, nr_fmt * nr);
63
64 static void add_aicvec(MYREAL aic, aic_fmt **aicvec, long *aicnum, nr_fmt *nr, long numparam, long remainnum);
65
66 static void aic_print_driver(aic_fmt *aicvec, MYREAL mleaic, long aicnum, nr_fmt *nr, world_fmt *world);
67
68 static void aic_score_driver(aic_struct *aic, MYREAL mle, MYREAL mleaic, nr_fmt *nr, world_fmt *world, char *temppattern);
69
70 //void akaike_information (world_fmt * world, long *Gmax);
71
72 static void print_header_aic (nr_fmt * nr, MYREAL mleaic);
73
74 static void print_aicfile(aic_fmt **aicvec, long *aicnum, nr_fmt *nr);
75
76 static void print_progress_aic(boolean add
77 , aic_fmt **aicvec, long *aicnum, nr_fmt *nr, long numparam, long freeparam);
78
79 //static long find_paramnum(world_fmt *world, char *connect);
80
81 ///
82 /// changes a linear matrix of ddd mm mm mm
83 /// to a diagonal matrix dmm mdm mmd
84 /// it destroys pattern
reshuffle(char * pattern,char * origpattern,long numpop)85 static char * reshuffle (char *pattern, char *origpattern, long numpop)
86 {
87 long space = 0;
88 long i, j, z = numpop;
89
90 for (i = 0; i < numpop; i++)
91 {
92 for (j = 0; j < numpop; j++)
93 {
94 if (i == j)
95 pattern[i * numpop + j + space] = 'x';
96 else
97 pattern[i * numpop + j + space] = origpattern[z++];
98 }
99 pattern[i * numpop + j + space] = ' ';
100 space++;
101 }
102 return pattern;
103 }
104
105 /// print the progress header for AIC based model selection
print_progressheader_aic(boolean progress,FILE * file,MYREAL mleaic,MYREAL numparam)106 void print_progressheader_aic(boolean progress, FILE *file, MYREAL mleaic, MYREAL numparam)
107 {
108 if (progress)
109 {
110 FPRINTF (file, "\n\n");
111 FPRINTF (file, " Selecting the best migration model for this run,\n");
112 FPRINTF (file, " This may take a while!\n");
113 FPRINTF (file, " Checking all parameter combinations\n");
114 FPRINTF (file, " with (AIC= -2 Log(L(Param)+n_param)<%f+%f\n",
115 mleaic, numparam);
116 }
117 }
118
119 ///
120 /// print the header for the AIC based model selection
121 /// the routine reports also progress to screen and logfile
122 void
print_header_aic(nr_fmt * nr,MYREAL mleaic)123 print_header_aic (nr_fmt * nr, MYREAL mleaic)
124 {
125 world_fmt *world = nr->world;
126
127 print_progressheader_aic(world->options->progress, stdout,mleaic, world->options->aicmod*world->numpop2);
128 print_progressheader_aic(world->options->writelog, stdout,mleaic, world->options->aicmod*world->numpop2);
129
130 FPRINTF (world->outfile, "\n\n\n Akaike's Information Criterion (AIC)\n");
131 FPRINTF (world->outfile, "=========================================\n\n");
132 FPRINTF (world->outfile, "[Linearized migration matrix, x=diagonal]\n");
133 }
134
135 ///
136 /// Calculates AIC score
137 /// \callgraph
aic_score_driver(aic_struct * aic,MYREAL mle,MYREAL mleaic,nr_fmt * nr,world_fmt * world,char * temppattern)138 void aic_score_driver(aic_struct *aic, MYREAL mle, MYREAL mleaic, nr_fmt *nr, world_fmt *world, char *temppattern)
139 {
140 char * savecustm;
141 char *savecustm2;
142
143 savecustm = (char *) mycalloc (2 * world->numpop2, sizeof (char));
144 savecustm2 = savecustm + world->numpop2;
145 memcpy (savecustm, world->options->custm2, sizeof (char) * nr->numpop2);
146 memcpy (savecustm2, world->options->custm, sizeof (char) * nr->numpop2);
147
148 aic->aicnum = 1;
149 aic->aicvec = (aic_fmt *) mycalloc (aic->aicnum, sizeof (aic_fmt));
150 aic->aicvec[0].mle = mle;
151 aic->aicvec[0].aic = mleaic;
152 aic->aicvec[0].lrt = 0.0;
153 aic->aicvec[0].prob = 1.0;
154 aic->aicvec[0].probcorr = 1.0;
155 aic->aicvec[0].numparam = nr->partsize;
156 aic->aicvec[0].pattern = (char *) mycalloc (nr->partsize + 1, sizeof (char));
157 memcpy (aic->aicvec[0].pattern, world->options->custm2, sizeof (char) * nr->partsize);
158 if (world->options->fast_aic)
159 {
160 fast_aic_score (&aic->aicvec, &aic->aicnum, nr, 0L, world->numpop,
161 temppattern, aic->param0, world->options->aictype[0]);
162 memcpy (world->options->custm2, savecustm, sizeof (char) * nr->numpop2);
163 if(world->options->aictype[1]!='\0')
164 fast_aic_score (&aic->aicvec, &aic->aicnum, nr, 0L, world->numpop,
165 temppattern, aic->param0, world->options->aictype[1]);
166 }
167 else
168 {
169 //find aic scores in a branch-and-bound fashion
170 // with some parameters set to zero, this needs more
171 // investigation because of boundary problems
172 aic_score (&aic->aicvec, &aic->aicnum, nr, 0L, world->numpop,
173 temppattern, aic->param0, world->options->aictype[0]);
174 memcpy (world->options->custm2, savecustm, sizeof (char) * nr->numpop2);
175 // aic scores based on averaging M (not 4Nm)
176 if(world->options->aictype[1]!='\0')
177 {
178 aic_score (&aic->aicvec, &aic->aicnum, nr, 0L, world->numpop, temppattern,
179 aic->param0, world->options->aictype[1]);
180 memcpy (world->options->custm2, savecustm, sizeof (char) * nr->numpop2);
181 }
182 }
183 memcpy (world->options->custm2, savecustm, sizeof (char) * nr->numpop2);
184 memcpy (world->options->custm, savecustm2, sizeof (char) * nr->numpop2);
185 myfree(savecustm); // savecustom2 is also freed, because it's part of savecustm
186 }
187
188 ///
189 /// print aic results, assume they are already ordered by qsort()
190 /// printing is:
191 /// pattern {migration matrix}
192 /// aic value, number of parameters, MLE of pattern, Prob of standard LRT, Prob of weighted chisquare
aic_print_driver(aic_fmt * aicvec,MYREAL mleaic,long aicnum,nr_fmt * nr,world_fmt * world)193 void aic_print_driver(aic_fmt *aicvec, MYREAL mleaic, long aicnum, nr_fmt *nr, world_fmt *world)
194 {
195 //#ifdef MYREAL == float
196 //const MYREAL eps = FLT_EPSILON ;
197 //#else
198 const MYREAL eps = DBL_EPSILON ;
199 //#endif
200 long i;
201 char *temppattern;
202 boolean mldone=FALSE;
203 char dashes[101] = "----------------------------------------------------------------------------------------------------";
204 temppattern = (char *) mycalloc (world->numpop2 + 1 + world->numpop, sizeof (char));
205 FPRINTF (world->outfile,
206 "%-*.*s AIC #param Log(L) LRT Prob Probc\n",
207 (int) MAX (18, nr->partsize + nr->numpop),
208 (int) (nr->partsize + nr->numpop), "Pattern");
209 FPRINTF (world->outfile,
210 "%-*.*s ---------- ---- ---------- ---------- ------- -------\n",
211 (int) MAX (18, nr->partsize + nr->numpop),
212 (int) (nr->partsize + nr->numpop), dashes);
213
214 for (i = 0; i < aicnum; i++)
215 {
216 FPRINTF (world->outfile, "%-*.*s %10.5f %4li % 10.5f % 10.5f %6.5f %6.5f\n",
217 (int) MAX (18, nr->partsize + nr->numpop),
218 (int) (nr->partsize + nr->numpop),
219 reshuffle (temppattern, aicvec[i].pattern, nr->numpop),
220 aicvec[i].aic, aicvec[i].numparam, aicvec[i].mle,
221 aicvec[i].lrt, aicvec[i].prob, aicvec[i].probcorr);
222 if ((aicvec[i].aic - mleaic > eps) && !mldone)
223 {
224 mldone = TRUE;
225 FPRINTF (world->outfile, "%-*.*s %53.53s\n",
226 (int) MAX (18, nr->partsize + nr->numpop),
227 (int) (nr->partsize + nr->numpop),dashes,
228 dashes);
229 }
230 myfree(aicvec[i].pattern);
231 }
232 FPRINTF (world->outfile, "\n\n\n");
233 myfree(temppattern);
234 }
235
236 ///
237 /// this drives the aic calculation and is called from main.c
238 /// \callgraph
239 void
akaike_information(world_fmt * world,long * Gmax)240 akaike_information (world_fmt * world, long *Gmax)
241 {
242 aic_struct aic;
243 nr_fmt *nr;
244 long kind = world->loci > 1 ? MULTILOCUS : SINGLELOCUS;
245 long repstop=0;
246 long repstart=0;
247
248 MYREAL mleaic;
249 MYREAL mle;
250 boolean multilocus=FALSE;
251
252 char *temppattern;
253
254 prepare_broyden (kind, world, &multilocus);
255 world->options->migration_model = MATRIX_ARBITRARY;
256
257 temppattern = (char *) mycalloc (world->numpop2 + 1 + world->numpop, sizeof (char));
258 aic.param0 = (MYREAL *) mycalloc (world->numpop2 + 1, sizeof (MYREAL));
259
260 set_replicates (world, world->repkind, world->rep, &repstart, &repstop);
261
262 if (kind == MULTILOCUS)
263 {
264 mle = world->atl[0][world->loci].param_like;
265 mleaic = -2. * mle + 2. * find_paramnum(world,NULL);
266 }
267 else
268 {
269 mle = world->atl[repstop == 1 ? 0 : repstop][0].param_like;
270 mleaic = -2. * mle + 2. * find_paramnum(world,NULL);
271 }
272 nr = (nr_fmt *) mycalloc (1, sizeof (nr_fmt));
273
274 create_nr (nr, world, *Gmax, 0L, world->loci, world->repkind, world->rep);
275
276 SETUPPARAM0 (world, nr, world->repkind, repstart, repstop,
277 world->loci, kind, multilocus);
278
279 print_header_aic (nr, mleaic);
280
281 if (kind == MULTILOCUS)
282 memcpy (aic.param0, nr->world->atl[0][nr->world->loci].param, sizeof (MYREAL) * nr->numpop2);
283 else
284 memcpy (aic.param0, nr->world->atl[repstop ==1 ? 0 : repstop][nr->world->locus].param, sizeof (MYREAL) * nr->numpop2);
285
286 aic_score_driver(&aic, mle, mleaic, nr, world, temppattern);
287
288 qsort ((void *) aic.aicvec, (unsigned long) aic.aicnum, sizeof (aic_fmt), aiccmp);
289
290 aic_print_driver(aic.aicvec, mleaic, aic.aicnum, nr, world);
291
292 myfree(aic.aicvec);
293 (void) fflush (world->outfile);
294 myfree(aic.param0);
295 myfree(temppattern);
296 destroy_nr (nr, world);
297 //myfree(nr);
298 }
299
300 ///
301 /// Does check all enumeration on one level and then picks only the best model to continue to
302 /// the next level
303 /// [this is different to the "branch-bound" algorithm]
fast_aic_score(aic_fmt ** aicvec,long * aicnum,nr_fmt * nr,long zero,long which,char * temppattern,MYREAL * param0,int migtype)304 void fast_aic_score (aic_fmt ** aicvec, long *aicnum, nr_fmt * nr,
305 long zero, long which, char *temppattern,
306 MYREAL *param0, int migtype)
307 {
308 //#ifdef MYREAL == float
309 //MYREAL numbermax = FLT_MAX ;
310 //#else
311 const MYREAL numbermax = MYREAL_MAX ;
312 //#endif
313
314 long m, ii;
315 MYREAL likes = 0;
316 MYREAL normd = 0;
317 MYREAL aic;
318 MYREAL borderaic = (*aicvec)[0].aic + nr->world->options->aicmod * nr->numpop2;
319 char savecustm2;
320 long remainnum = 0;
321 // boolean mylegal;
322 char *custm2 = nr->world->options->custm2;
323 char *scustm2;
324 long numparam;
325 long freeparam;
326 aic_fmt *best;
327
328 scustm2 = (char *) mycalloc (nr->partsize, sizeof (char));
329 memcpy (scustm2, custm2, sizeof (char) * nr->partsize);
330 if (migtype == 'm')
331 remainnum = 1;
332 numparam = zero;
333 freeparam = (nr->numpop2 - numparam - 1 + remainnum);
334 best = (aic_fmt *) mycalloc (nr->partsize, sizeof (aic_fmt));
335 for (m = nr->numpop; m < nr->numpop2; m++)
336 {
337 best[m].aic = numbermax;
338 best[m].numparam = m;
339 if (scustm2[m] == migtype)
340 continue;
341 savecustm2 = custm2[m];
342 custm2[m] = migtype;
343 memcpy (nr->world->param0, param0, sizeof (MYREAL) * nr->numpop2);
344 resynchronize_param (nr->world);
345 if ((legal_pattern (nr->world->options->custm2, nr->numpop)))
346 {
347 (void) do_profiles (nr->world, nr, &likes, &normd, PROFILE,
348 nr->world->rep, nr->world->repkind);
349 aic = -2. * nr->llike + 2. * freeparam;
350 best[m].aic = aic;
351 best[m].numparam = m;
352 add_aiclist(migtype, numparam, freeparam, remainnum, zero, m,
353 param0, temppattern, custm2, aicnum, aicvec, nr);
354 }
355 else
356 { // illegal combination of parameters
357 if (nr->world->options->progress)
358 FPRINTF (stdout, " F %s %20s\n",
359 reshuffle (temppattern, custm2, nr->numpop), "-----");
360 (void) fflush (stdout);
361 if (nr->world->options->writelog)
362 FPRINTF (nr->world->options->logfile, " F %s %20s\n",
363 reshuffle (temppattern, custm2, nr->numpop), "-----");
364 }
365 custm2[m] = savecustm2;
366 }
367 for (ii = nr->numpop; ii < nr->partsize; ii++)
368 {
369 if (best[ii].aic < borderaic && custm2[best[ii].numparam] != migtype)
370 {
371 custm2[best[ii].numparam] = migtype;
372 fast_aic_score (aicvec, aicnum, nr, zero + 1, best[ii].numparam,
373 temppattern, param0, migtype);
374 }
375 }
376 myfree(best);
377 myfree(scustm2);
378 }
379
380 /// add model and AIC score to the list of models for the final printing
add_aicvec(MYREAL aic,aic_fmt ** aicvec,long * aicnum,nr_fmt * nr,long numparam,long remainnum)381 void add_aicvec(MYREAL aic, aic_fmt **aicvec, long *aicnum, nr_fmt *nr, long numparam, long remainnum)
382 {
383 MYREAL lrt = -2. * (nr->llike - (*aicvec)[0].mle);
384
385 *aicvec = (aic_fmt *)
386 myrealloc (*aicvec, sizeof (aic_fmt) * (*aicnum + 1));
387 (*aicvec)[*aicnum].pattern = (char *)
388 mycalloc (nr->partsize + 1, sizeof (char));
389 (*aicvec)[*aicnum].aic = aic;
390 (*aicvec)[*aicnum].mle = nr->llike;
391 (*aicvec)[*aicnum].numparam = nr->numpop2 - numparam - 1 + remainnum;
392 memcpy ((*aicvec)[*aicnum].pattern, nr->world->options->custm2, sizeof (char) * nr->partsize);
393
394 (*aicvec)[*aicnum].lrt = lrt;
395 (*aicvec)[*aicnum].prob = probchi (numparam, (*aicvec)[*aicnum].lrt);
396 (*aicvec)[*aicnum].probcorr = probchiboundary ((*aicvec)[*aicnum].lrt, numparam, numparam);
397 }
398
399 /// print AIC model selection results
400 void
print_aicfile(aic_fmt ** aicvec,long * aicnum,nr_fmt * nr)401 print_aicfile(aic_fmt **aicvec, long *aicnum, nr_fmt *nr)
402 {
403 long i;
404 if (nr->world->options->aicfile)
405 {
406 FPRINTF (nr->world->options->aicfile, "%f %f %li %f %f ",
407 (*aicvec)[*aicnum].aic, (*aicvec)[*aicnum].lrt,
408 (*aicvec)[*aicnum].numparam,
409 (*aicvec)[*aicnum].prob,
410 (*aicvec)[*aicnum].probcorr);
411
412 for (i = 0; i < nr->partsize; i++)
413 FPRINTF (nr->world->options->aicfile, "%f ", nr->world->param0[i]);
414 FPRINTF (nr->world->options->aicfile, "\n");
415 }
416 }
417
418 /// print progress report for AIC model selection
419 void
print_progress_aic(boolean add,aic_fmt ** aicvec,long * aicnum,nr_fmt * nr,long numparam,long freeparam)420 print_progress_aic(boolean add
421 , aic_fmt **aicvec, long *aicnum, nr_fmt *nr, long numparam, long freeparam)
422 {
423 MYREAL mle, aic, lrt, prob, probcorr;
424 char *temppattern;
425 char *custm2 = nr->world->options->custm2;
426 temppattern = (char *) mycalloc (nr->world->numpop2 + 1 + nr->world->numpop, sizeof (char));
427 if(add
428 )
429 {
430 if (nr->world->options->progress)
431 FPRINTF (stdout,
432 " + %s %20.5f %3li %8.4f %8.4f %6.4f %6.4f\n",
433 reshuffle (temppattern, custm2, nr->numpop), (*aicvec)[*aicnum].aic,
434 freeparam, (*aicvec)[*aicnum].mle,
435 (*aicvec)[*aicnum].lrt, (*aicvec)[*aicnum].prob,
436 (*aicvec)[*aicnum].probcorr);
437 if (nr->world->options->writelog)
438 FPRINTF (nr->world->options->logfile,
439 " + %s %20.5f %3li %8.4f %8.4f %6.4f %6.4f\n",
440 reshuffle (temppattern, custm2, nr->numpop), (*aicvec)[*aicnum].aic,
441 freeparam, (*aicvec)[*aicnum].mle,
442 (*aicvec)[*aicnum].lrt, (*aicvec)[*aicnum].prob,
443 (*aicvec)[*aicnum].probcorr);
444 }
445 else
446 {
447 mle = nr->llike;
448 aic = -2. * mle + 2. * freeparam;
449 lrt = -2. * (mle - (*aicvec)[0].mle);
450 prob = probchi (numparam, lrt);
451 probcorr = probchiboundary (lrt, numparam, numparam);
452 if (nr->world->options->progress)
453 FPRINTF (stdout,
454 " - %s %20.5f %3li %8.4f %8.4f %6.4f %6.4f\n",
455 reshuffle (temppattern, custm2, nr->numpop), aic,
456 freeparam, mle, lrt,
457 prob,
458 probcorr);
459 if (nr->world->options->writelog)
460 FPRINTF (nr->world->options->logfile,
461 " + %s %20.5f %3li %8.4f %8.4f %6.4f %6.4f\n",
462 reshuffle (temppattern, custm2, nr->numpop), aic,
463 freeparam, mle, lrt,
464 prob,
465 probcorr);
466 }
467 (void) fflush (stdout);
468 myfree(temppattern);
469 }
470
471 /// add the AIC models and scores to the result list
472 void
add_aiclist(int migtype,long numparam,long freeparam,long remainnum,long zero,long m,MYREAL * param0,char * temppattern,char * custm2,long * aicnum,aic_fmt ** aicvec,nr_fmt * nr)473 add_aiclist (int migtype, long numparam, long freeparam, long remainnum, long zero, long m,
474 MYREAL *param0, char *temppattern, char *custm2,
475 long *aicnum, aic_fmt ** aicvec, nr_fmt * nr)
476 {
477 MYREAL aic = -2. * nr->llike + 2. * freeparam;
478
479 if (aic < (*aicvec)[0].aic + nr->world->options->aicmod * freeparam)
480 {
481 if (migtype != 'm' || (migtype == 'm' && nr->world->options->mmn > 1))
482 {
483 add_aicvec(aic, aicvec, aicnum, nr, numparam, remainnum);
484 print_aicfile(aicvec,aicnum,nr);
485 print_progress_aic(TRUE, aicvec,aicnum,nr, numparam, freeparam);
486 (*aicnum)++;
487 aic_score (aicvec, aicnum, nr, zero + 1, m + 1,
488 temppattern, param0, migtype);
489 }
490 }
491 else
492 {
493 print_progress_aic(FALSE, aicvec,aicnum,nr, numparam, freeparam);
494 }
495 }
496
497 /// calculates AIC scores
498 void
aic_score(aic_fmt ** aicvec,long * aicnum,nr_fmt * nr,long zero,long which,char * temppattern,MYREAL * param0,int migtype)499 aic_score (aic_fmt ** aicvec, long *aicnum, nr_fmt * nr,
500 long zero, long which, char *temppattern, MYREAL *param0,
501 int migtype)
502 {
503 long m;
504 MYREAL likes = 0;
505 MYREAL normd = 0;
506 char savecustm2;
507 long remainnum = 0;
508 // boolean mylegal;
509 char *custm2 = nr->world->options->custm2;
510 long numparam = 0;
511 long freeparam;
512
513 switch (migtype)
514 {
515 case '0':
516 numparam = nr->world->options->zeron;
517 remainnum = 0;
518 break;
519 case 'm':
520 numparam = nr->world->options->mmn;
521 remainnum = 1;
522 if (nr->world->options->custm2[which] == 'm')
523 return;
524 break;
525 }
526 freeparam = (nr->numpop2 - numparam - 1 + remainnum);
527 for (m = which; m < nr->numpop2; m++)
528 {
529 savecustm2 = custm2[m];
530 custm2[m] = migtype;
531 memcpy (nr->world->param0, param0, sizeof (MYREAL) * nr->numpop2);
532 resynchronize_param (nr->world);
533
534 if ((legal_pattern (nr->world->options->custm2, nr->numpop)))
535 {
536 (void) do_profiles (nr->world, nr, &likes, &normd, PROFILE,
537 nr->world->rep, nr->world->repkind);
538
539 add_aiclist(migtype, numparam, freeparam, remainnum, zero, m, param0,
540 temppattern, nr->world->options->custm2, aicnum, aicvec, nr);
541 }
542 else
543 {
544 if (nr->world->options->progress)
545 {
546 FPRINTF (stdout, " F %s %20s\n", reshuffle (temppattern, custm2, nr->numpop), "-----");
547 fflush (stdout);
548 }
549 if (nr->world->options->writelog)
550 FPRINTF (nr->world->options->logfile, " F %s %20s\n", reshuffle (temppattern, custm2, nr->numpop), "-----");
551 }
552 custm2[m] = savecustm2;
553 }
554 }
555
556 /* not used
557 /// check how many parameter are set to zero or are set to m (average)
558 static boolean
559 check_numparam (long which, long migtype, worldoption_fmt * options,
560 long *numparam, long *remainnum)
561 {
562 boolean rc = FALSE;
563 switch (migtype)
564 {
565 case '0':
566 *numparam = options->zeron;
567 *remainnum = 0;
568 break;
569 case 'm':
570 *numparam = options->mmn;
571 *remainnum = 1;
572 if (options->custm2[which] == 'm')
573 rc = TRUE;
574 break;
575 }
576 return rc;
577 }
578 end not used */
579
580 ///
581 /// check whether the migration pattern in matrix is legal or not
582 /// a legal migration pattern, needs to connect all populations
583 /// eahc population needs to be connected to at least one other population
584 /// in one or both direction of the migration route
585 boolean
legal_pattern(char * matrix,long numpop)586 legal_pattern (char *matrix, long numpop)
587 {
588 //#ifdef MYREAL == float
589 //const MYREAL eps = FLT_EPSILON ;
590 //#else
591 const MYREAL eps = DBL_EPSILON ;
592 //#endif
593 long from=0, to=0, i;
594 MYREAL summ = -1;
595 long oldto;
596 for (i = 0; i < numpop; i++)
597 {
598 if (matrix[i] == '0')
599 return FALSE;
600 }
601 oldto = -1;
602 for (i = numpop; i < numpop * numpop; i++)
603 {
604 m2mm (i, numpop, &from, &to);
605 if (oldto != to)
606 {
607 if (fabs(summ)< eps )
608 return FALSE;
609 oldto = to;
610 summ = 0;
611 }
612 summ += (MYREAL) (matrix[i] != '0') + (MYREAL) (matrix[mm2m (to, from, numpop)] != '0');
613 }
614 if (summ < eps)
615 return FALSE;
616 return TRUE;
617 }
618
find_paramnum(world_fmt * world,char * connect)619 long find_paramnum(world_fmt *world, char *connect)
620 {
621 long cc;
622 long tm;
623 long mm;
624 // long mmm;
625 long ms;
626 long mss;
627 long total;
628 char *temp, *temp1, *temp2, *t;
629
630 // checkout if GAMMA mutation rate
631 total = world->numpop * world->numpop + (world->options->gamma ? 1 : 0) ;
632
633 if(connect==NULL)
634 {
635 // checkout zeroes and constants in custom migration matrix
636 // removes 1 for each zero and (constant>0)
637 total -= world->options->zeron + world->options->constn;
638 // checkout mean values
639 // removes number of parameters that are part of the average and adds
640 // 1 for average migrations and 1 for average thetas
641 total -= (world->options->tmn > 0 ? world->options->tmn - 1 : 0) + (world->options->mmn > 0 ? world->options->mmn - 1 : 0);
642 // checkout symmetric values
643 // removes 1 parameter for each migration pair
644 total -= world->options->symn + world->options->sym2n;
645 }
646 else
647 {
648 temp1 = (char *) mycalloc(strlen(connect)+1, sizeof(char));
649 //temp2 = (char *) mycalloc(strlen(connect)+1, sizeof(char));
650 temp = (char *) mycalloc(strlen(connect)+1, sizeof(char));
651 t = temp;
652 strcpy(temp1,connect);
653 temp2 = strtok(temp1," ;,:\n\t\r");
654 while(temp2!= NULL)
655 {
656 switch(temp2[0])
657 {
658 case '*': *t++ = '*'; break;
659 case 'c': *t++ = 'c'; break;
660 case 'm': *t++ = 'm'; break;
661 case 'M': *t++ = 'M'; break;
662 case 's': *t++ = 's'; break;
663 case 'S': *t++ = 'S'; break;
664 case ' ' : break;
665 default:
666 *t++ = 'c';
667 }
668 temp2 = strtok(NULL," ;,:\n\t\r");
669 }
670 cc = scan_connect (temp, 0L, world->numpop2, 'c');
671 tm = scan_connect (temp, 0L, world->numpop, 'm');
672 mm = scan_connect (temp, world->numpop, world->numpop2, 'm');
673 // mmm = scan_connect (temp, world->numpop, world->numpop2, 'M');
674 ms = scan_connect (temp, world->numpop, world->numpop2, 's')/2;
675 mss = scan_connect (temp, world->numpop, world->numpop2, 'S')/2;
676 total -= cc + (tm > 0 ? tm-1 : 0);
677 total -= (mm > 0 ? mm-1 : 0 );
678 total -= ms;
679 total -= mss;
680 myfree(temp);
681 myfree(temp1);
682 //myfree(temp2);
683 }
684 return total;
685 }
686