1 /* \File options.c */
2 /*------------------------------------------------------
3 Maximum likelihood estimation
4 of migration rate and effectice population size
5 using a Metropolis-Hastings Monte Carlo algorithm
6 -------------------------------------------------------
7 O P T I O N S R O U T I N E S
8
9 creates options structures,
10 reads options from parmfile if present
11
12 prints options,
13 and finally helps to destroy itself.
14
15 Peter Beerli 1996-2006
16 beerli@fsu.edu
17
18 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
19 Copyright 2003-2013 Peter Beerli, Tallahassee FL
20
21 Permission is hereby granted, free of charge, to any person obtaining
22 a copy of this software and associated documentation files (the
23 "Software"), to deal in the Software without restriction, including
24 without limitation the rights to use, copy, modify, merge, publish,
25 distribute, sublicense, and/or sell copies of the Software, and to
26 permit persons to whom the Software is furnished to do so, subject
27 to the following conditions:
28
29 The above copyright notice and this permission notice shall be
30 included in all copies or substantial portions of the Software.
31
32 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
33 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
34 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
35 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
36 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
37 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
38 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
39
40
41 $Id: options.c 2174 2013-12-16 05:27:21Z beerli $
42 -------------------------------------------------------*/
43 //#include <stdio.h>
44 #include <time.h>
45 #include "migration.h"
46 #include "sighandler.h"
47 #include <stdarg.h>
48 #include "fst.h"
49 #include "tools.h"
50 #include "migrate_mpi.h"
51 #include "bayes.h"
52 #include "menu.h"
53 #include "random.h"
54 #ifdef DMALLOC_FUNC_CHECK
55 #include <dmalloc.h>
56 #endif
57 #include "options.h"
58
59 extern char * generator;
60
61 /* parmfile parameter specifications and keywords */
62 //#define LONGLINESIZE 10000
63 #define NUMBOOL 31
64 #define BOOLTOKENS {"menu","interleaved","print-data","mixfile",\
65 "moving-steps","freqs-from-data","useroldtree", \
66 "autocorrelation", "simulation","plot", "weights",\
67 "read-summary","write-summary","NODATA", "include-unknown", "print-fst",\
68 "distfile","geo","gelman-convergence", "randomtree",\
69 "fast-likelihood", "aic-modeltest", "use-M", "use-Nm","bayes-update", "bayes-allfile", "bayes-file", "divergence-times","tipdate-file","heated-swap","auto-tune"}
70 #define NUMNUMBER 65
71 #define NUMBERTOKENS {"ttratio","short-chains",\
72 "short-steps","short-inc","long-chains",\
73 "long-steps", "long-inc", "theta", \
74 "nmlength","random-seed","migration","mutation",\
75 "datatype", "categories","rates","prob-rates", \
76 "micro-max", "micro-threshold", "delimiter","burn-in",\
77 "infile", "outfile", "mathfile", "title", \
78 "long-chain-epsilon","print-tree","progress","l-ratio",\
79 "fst-type","profile","custom-migration","sumfile","short-sample",\
80 "long-sample", "replicate","cpu","logfile", "seqerror-rate","uep", \
81 "uep-rates","uep-bases", "mu-rates","heating","fluctuate", "resistance",\
82 "bayes-updatefreq", "bayesfile","bayes-prior", "usertree", "bayes-posteriorbins",\
83 "mig-histogram", "bayes-posteriormaxtype", "pdf-outfile",\
84 "bayes-allfileinterval", "bayes-priors","skyline","rates-gamma", "bayes-proposals", \
85 "generation-per-year", "mutationrate-per-year", "inheritance-scalars", "micro-submodel","random-subset", "population-relabel","analyze-loci"};
86
87 // myID is a definition for the executing node (master or worker)
88 extern int myID;
89
90 fpos_t thePos; // used to track pposition in the parmfile
91
92 /* prototypes ------------------------------------------- */
93 void set_usem_related (option_fmt * options);
94 //void create_options (option_fmt ** options);
95 //void init_options (option_fmt * options);
96 //void set_param (world_fmt * world, data_fmt * data, option_fmt * options,long locus);
97 //void set_profile_options (option_fmt * options);
98 //void print_menu_options (world_fmt * world, option_fmt * options,
99 // data_fmt * data);
100 //void print_options (FILE * file, world_fmt * world, option_fmt * options,
101 // data_fmt * data);
102 //void decide_plot (worldoption_fmt * options, long chain, long chains,
103 // char type);
104 //void destroy_options (option_fmt * options);
105
106 //void read_options_master (option_fmt * options);
107 //void read_options_worker (char **buffer, option_fmt * options);
108 ///* private functions */
109 boolean booleancheck (option_fmt * options, char *var, char *value);
110 //long boolcheck (char ch);
111 boolean numbercheck (option_fmt * options, char *var, char *value);
112 //void reset_oneline (option_fmt * options, long position);
113 void reset_oneline (option_fmt * options, fpos_t * position);
114 void read_theta (option_fmt * options);
115 void read_mig (option_fmt * options);
116 #ifdef MPI
117 void read_theta_worker (char **buffer, option_fmt * options);
118 void read_mig_worker (char **buffer, option_fmt * options);
119 char skip_sspace (char **buffer);
120 void read_custom_migration_worker (char **buffer, option_fmt * options,
121 char *value, long customnumpop);
122 #endif
123 char skip_space (option_fmt * options);
124 //void read_custom_migration (FILE * file, option_fmt * options, char *value,
125 // long customnumpop);
126 //void synchronize_param (world_fmt * world, option_fmt * options);
127 //void resynchronize_param (world_fmt * world);
128 void specify_migration_type (option_fmt * options);
129 //void print_options_nomcmc (FILE * file, option_fmt * options,
130 // world_fmt * world);
131 //void destroy_options (option_fmt * options);
132 //void set_partmean_mig (long **mmparam, MYREAL *param, char *custm2, long migm,
133 // long numpop2);
134 //long scan_connect (char *custm2, long start, long stop, int check);
135 //void set_plot (option_fmt * options);
136 void set_plot_values (MYREAL **values, MYREAL plotrange[], long intervals, int type);
137 //void set_grid_param (world_fmt * world, long gridpoints);
138 void print_arbitrary_migration_table (FILE * file, world_fmt * world, option_fmt *options,
139 data_fmt * data);
140 void print_distance_table (FILE * file, world_fmt * world,
141 option_fmt * options, data_fmt * data);
142
143 void fillup_custm (long len, world_fmt * world, option_fmt * options);
144 //
145
146 //long save_options_buffer (char **buffer, long *allocbufsize, option_fmt * options, data_fmt *data);
147 //void print_parm_delimiter(long *bufsize, char **buffer, long *allocbufsize);
148 //void print_parm_br(long *bufsize, char **buffer, long *allocbufsize);
149 //
150 void prior_consistency(option_fmt *options, int paramtype, char priortype);
151 //
152 ///*======================================================*/
153 /// initialize filename
init_filename(char ** filename,char initstring[])154 void init_filename (char **filename, char initstring[])
155 {
156 int len=LINESIZE;
157 // len = strlen(initstring);
158 *filename = (char *) mycalloc(len+2, sizeof(char));
159 sprintf(*filename, "%-*s", len, initstring);
160 remove_trailing_blanks(filename);
161 }
162
163 void
create_options(option_fmt ** options)164 create_options (option_fmt ** options)
165 {
166 (*options) = (option_fmt *) mycalloc (1, sizeof (option_fmt));
167 }
168
169 /// initialize options
init_options(option_fmt * options)170 void init_options (option_fmt * options)
171 {
172 unsigned long i;
173 unsigned long timeseed;
174 /* General options --------------------------------------- */
175 //
176 // name length
177 options->nmlength = DEFAULT_NMLENGTH;
178 options->popnmlength = DEFAULT_POPNMLENGTH;
179 options->allelenmlength = DEFAULT_ALLELENMLENGTH;
180 //
181 // custom migration matrix setup
182 options->custm = (char *) mycalloc (1, sizeof (char) * 1000); /// \todo needs
183 //options->custm2 is allocated later
184 //options->symn = 0;
185 //options->sym2n = 0;
186 //options->zeron = 0;
187 //options->constn = 0;
188 //options->mmn = 0;
189 //options->tmn = 0;
190
191 /* input/output options ---------------------------------- */
192 options->menu = TRUE;
193 options->progress = TRUE;
194 options->verbose = FALSE;
195 options->writelog = FALSE;
196 options->geo = FALSE;
197 options->div = FALSE;
198 #ifdef UEP
199 options->uep = FALSE;
200 options->ueprate = 1.0;
201 options->uepmu = 0.9999;
202 options->uepnu = 0.0001;
203 options->uepfreq0=0.5;
204 options->uepfreq1=0.5;
205 // options->uep_last = FALSE;
206 #endif
207 options->printdata = FALSE;
208 options->printcov = FALSE;
209 options->usertree = FALSE;
210 options->usertreewithmig = FALSE;
211 options->randomtree = TRUE; //default changed May 29 2006 [convergence statistic]
212 options->fastlike = FALSE;
213 options->treeprint = _NONE;
214 options->treeinc = 1;
215 options->printfst = FALSE;
216 options->fsttype = THETAVARIABLE;
217 options->usem = TRUE;
218 options->migvar = PLOT4NM;
219 fst_type (options->fsttype);
220 options->plot = FALSE;
221 options->plotmethod = PLOTALL; /* outfile and mathematica file */
222 options->plotvar = PLOT4NM;
223 options->plotscale = PLOTSCALELOG;
224 options->plotrange[0] = PLANESTART; /* start x axis */
225 options->plotrange[1] = PLANEEND; /*end x axis */
226 options->plotrange[2] = PLANESTART; /*start y axis */
227 options->plotrange[3] = PLANEEND; /* end y axis */
228 options->plotintervals = PLANEINTERVALS;
229 options->simulation = FALSE;
230 options->movingsteps = FALSE;
231 options->acceptfreq = 0.0;
232 options->mighist = FALSE;
233 options->mighist_all = FALSE;
234 options->mighist_counter = 0;
235 options->mighist_increment = 1;
236 options->skyline = FALSE;
237 options->eventbinsize = (float) 0.001;
238 options->mixplot = FALSE;
239 init_filename( &options->parmfilename, PARMFILE);
240 init_filename( &options->infilename, INFILE);
241 init_filename( &options->outfilename, OUTFILE);
242 #ifdef PRETTY
243 init_filename( &options->pdfoutfilename, PDFOUTFILE);
244 #endif
245 #ifdef THERMOCHECK
246 init_filename( &options->mixfilename, MIXFILE);
247 #endif
248 init_filename( &options->logfilename, LOGFILE);
249 init_filename( &options->mathfilename, MATHFILE);
250 init_filename( &options->sumfilename, SUMFILE);
251 init_filename( &options->treefilename, TREEFILE);
252 init_filename( &options->utreefilename, UTREEFILE);
253 init_filename( &options->catfilename, CATFILE);
254 init_filename( &options->weightfilename, WEIGHTFILE);
255
256 init_filename( &options->mighistfilename, MIGHISTFILE);
257 init_filename( &options->skylinefilename, SKYLINEFILE);
258
259 init_filename( &options->distfilename, DISTFILE);
260 init_filename( &options->geofilename, GEOFILE);
261 init_filename( &options->divfilename, DIVFILE);
262 init_filename( &options->bootfilename, BOOTFILE);
263 init_filename( &options->aicfilename, AICFILE);
264 init_filename( &options->seedfilename, SEEDFILE);
265 #ifdef UEP
266 init_filename( &options->uepfilename, UEPFILE);
267 #endif
268 init_filename( &options->bayesfilename, BAYESFILE);
269 #ifdef ZNZ
270 init_filename( &options->bayesmdimfilename, BAYESMDIMFILE);
271 options->use_compressed=1;
272 #else
273 init_filename( &options->bayesmdimfilename, BAYESMDIMFILE2);
274 options->use_compressed=0;
275 #endif
276 init_filename( &options->datefilename, TIPDATEFILE);
277 // allocation of prior parameter settings
278 options->bayespriortheta = (prior_fmt *) mycalloc (3, sizeof (prior_fmt));
279 options->bayespriorm = options->bayespriortheta + 1;
280 options->bayespriorrate = options->bayespriortheta + 2;
281 // prior theta
282 options->bayespriortheta->min = SMALLEST_THETA;
283 options->bayespriortheta->mean = DNA_GUESS_THETA;
284 options->bayespriortheta->max = 10. * DNA_GUESS_THETA;
285 options->bayespriortheta->bins = BAYESNUMBIN;
286 options->bayespriortheta->delta = (options->bayespriortheta->max - options->bayespriortheta->min)/10.;
287 // prior M
288 options->bayespriorm->min = SMALLEST_MIGRATION;
289 options->bayespriorm->mean = DNA_GUESS_MIG;
290 options->bayespriorm->max = 10. * DNA_GUESS_MIG;
291 options->bayespriorm->bins = BAYESNUMBIN;
292 options->bayespriorm->delta = (options->bayespriorm->max - options->bayespriorm->min)/10.;
293 // prior rate
294 options->bayespriorrate->min = 1.0;
295 options->bayespriorrate->mean = 1.0;
296 options->bayespriorrate->max = 1.0;
297 options->bayespriorrate->bins = BAYESNUMBIN;
298 options->bayespriorrate->delta = (options->bayespriorrate->max - options->bayespriorrate->min)/10.;
299 for(i=0;i<PRIOR_SIZE;i++)
300 options->slice_sampling[i] = FALSE;// option for proposal setting
301 #ifdef PRETTY
302 options->bayespretty = PRETTY_P100;
303 #endif
304 strcpy (options->title, "\0");
305 options->lratio = (lratio_fmt *) mycalloc (1, sizeof (lratio_fmt));
306 options->lratio->alloccounter = 1;
307 options->lratio->data =
308 (lr_data_fmt *) mycalloc (1, sizeof (lr_data_fmt) *options->lratio->alloccounter);
309 options->lratio->data[0].value1 =
310 (char *) mycalloc (1, sizeof (char) * LONGLINESIZE);
311 options->lratio->data[0].value2 =
312 (char *) mycalloc (1, sizeof (char) * LONGLINESIZE);
313 options->lratio->data[0].connect =
314 (char *) mycalloc (1, sizeof (char) * LONGLINESIZE);
315 options->aic = FALSE;
316 options->fast_aic = FALSE;
317 strcpy(options->aictype,"0m");
318 options->profile = ALL;
319 options->profilemethod = 'p';
320 options->df = 1;
321 options->qdprofile = FALSE;
322 options->printprofsummary = TRUE;
323 options->printprofile = TRUE;
324 if (options->usem)
325 options->profileparamtype = options->plotvar = options->migvar = PLOTM;
326 else
327 options->profileparamtype = options->plotvar = options->migvar = PLOT4NM;
328 /* data options ------------------------------------------ */
329 options->datatype = 's';
330 options->include_unknown = FALSE;
331 options->migration_model = MATRIX;
332 options->thetag = (MYREAL *) mycalloc (1, sizeof (MYREAL) * NUMPOP);
333 options->mg = (MYREAL *) mycalloc (1, sizeof (MYREAL) * NUMPOP);
334 options->gamma = FALSE;
335 options->alphavalue = START_ALPHA;
336 options->murates = FALSE;
337 options->bayesmurates = FALSE;
338 #ifdef LONGSUM
339 options->fluctuate=FALSE;
340 options->flucrates = (MYREAL *) mycalloc(NUMPOP * 3, sizeof(MYREAL));
341 #endif
342 options->mu_rates = NULL;
343 options->inheritance_scalars = NULL;
344 options->newpops = NULL;
345 options->prioralone=FALSE;
346 /* EP data */
347 options->dlm = '\0';
348 /* microsat data */
349 options->micro_threshold = MICRO_THRESHOLD;
350 options->micro_stepnum = MAX_MICROSTEPNUM;
351 options->msat_option = SINGLESTEP;
352 options->msat_tuning[0] = 0.0;
353 options->msat_tuning[1] = 0.5;
354 /*sequence data */
355 options->seqerror = 0.0;
356 options->interleaved = FALSE;
357 options->ttratio = (MYREAL *) mycalloc (1, sizeof (MYREAL) * 2);
358 options->ttratio[0] = 2.0;
359 options->freqsfrom = TRUE;
360 options->categs = ONECATEG;
361 options->rate = (MYREAL *) mycalloc (1, sizeof (MYREAL));
362 options->rate[0] = 1.0;
363 options->rcategs = 1;
364 options->rrate = (MYREAL *) mycalloc (1, sizeof (MYREAL));
365 options->probcat = (MYREAL *) mycalloc (1, sizeof (MYREAL));
366 options->autocorr = FALSE;
367 options->rrate[0] = 1.0;
368 options->probcat[0] = 1.0;
369
370 options->probsum = 0.0;
371
372 options->lambda = 1.0;
373 options->weights = FALSE;
374 /* random number options --------------------------------- */
375 options->autoseed = AUTO;
376 options->autoseed = AUTO;
377 //#ifndef MAC
378 timeseed = (unsigned long) time (NULL) / 4;
379 //#else
380 //timeseed = (unsigned long) clock () / 4;
381 //#endif
382 options->inseed = (long) timeseed + 1;
383 /* mcmc options ------------------------------------------ */
384 options->thetaguess = FST;
385 options->migrguess = FST;
386 for (i = 0; i < NUMPOP; i++)
387 {
388 options->thetag[i] = 1.0;
389 options->mg[i] = 1.0;
390 options->custm[i] = '*';
391 }
392 options->custm[i] = '\0';
393 options->custm2 = (char * ) strdup(options->custm);
394 options->numthetag = options->nummg = 0;
395 options->bayes_infer = TRUE; //DEFAULT TO BAYESIAN INFERERENCE July 2012
396 #ifdef INTEGRATEDLIKE
397 options->integrated_like = TRUE;
398 #else
399 options->integrated_like = FALSE;
400 #endif
401 options->updateratio = HALF;
402 options->bayesmdiminterval = 1;
403 options->schains = 10;
404 options->sincrement = 100;
405 options->ssteps = 500;
406 options->lchains = 1;
407 options->lincrement = 100;
408 options->lsteps = 5000;
409 options->burn_in = BURNINPERIOD;
410 options->burnin_autostop = ' ';
411 //
412 // heating options
413 options->heating = 0; /* no heating */
414 options->adaptiveheat=NOTADAPTIVE;
415 options->heating_interval = 0;
416 options->heatedswap_off=FALSE;
417 options->heat[0] = COLD;
418 options->heat[1] = WARM;
419 options->heat[2] = HOT;
420 options->heat[3] = VERYHOT;
421 options->heated_chains = 1;
422 //
423 // obscure options to accept more different trees
424 options->lcepsilon = LONGCHAINEPSILON;
425 options->gelman = FALSE; // changed default May 29 2006 changed back Feb 2013
426 options->gelmanpairs = FALSE;
427 options->pluschain = PLUSCHAIN;
428 //
429 // replication options
430 options->replicate = FALSE;
431 options->replicatenum = 0;
432 options->gridpoints = 0;
433 /* genealogy summary options----------------------------------- */
434 options->readsum = FALSE;
435 options->writesum = FALSE;
436 /*threading over loci */
437 options->cpu = 1;
438 //
439 // fatal attraction to zero resistance
440 options->minmigsumstat = MINMIGSUMSTAT;
441 //
442 options->segregs=NULL;
443 options->wattersons=NULL;
444
445 options->onlyvariable=FALSE; //analyzes only variable sequence loci:default: analyze all
446 options->has_variableandone = FALSE; // analyze only one invariant locus (not more) and reweight
447 options->firstinvariant = -1; // locus ID of invariant locus analyzed
448
449 options->mutationrate_year = (MYREAL*) mycalloc(1, sizeof(MYREAL));
450 options->mutationrate_year_numalloc = 1;
451 options->generation_year = 1.0;
452 options->randomsubset = 0;
453 options->randomsubsetseed = -1;
454 options->inheritance_scalars = (MYREAL*) mycalloc(1, sizeof(MYREAL));
455 options->inheritance_scalars[0] = 1.0;
456 options->inheritance_scalars_numalloc = 1;
457 options->newpops = (long*) mycalloc(1, sizeof(long));
458 options->newpops[0] = 1;
459 options->newpops_numalloc = 1;
460
461 options->slice_sticksizes = (MYREAL*) mycalloc(1, sizeof(MYREAL));
462 options->slice_sticksizes[0] = 1.0;
463 options->heatedswap_off=FALSE;
464 options->has_autotune = TRUE;
465 options->autotune=AUTOTUNEDEFAULT;
466 options->totalsites = 0;
467 }
468
469 void
read_options_master(option_fmt * options)470 read_options_master (option_fmt * options)
471 {
472 const unsigned long linecpmax = LINESIZE -1;
473
474 long counter = 0;
475 //long position = 0;
476 // to make windows reading unix parmfile happy
477 // fpos_t thePos;
478
479 char varvalue[LONGLINESIZE];
480 char parmvar[LONGLINESIZE];
481 char input[LONGLINESIZE];
482 char *p, *tmp;
483
484 options->parmfile = fopen(options->parmfilename, "r");
485 if (options->parmfile)
486 {
487 counter = 0;
488 //position = ftell (options->parmfile);
489 fgetpos(options->parmfile, &thePos);
490 FGETS (input, LINESIZE, options->parmfile);// read first set of ###
491 FGETS (input, LINESIZE, options->parmfile);// read "# Parmfile for Migrate"
492 if(strncmp(input,"# Parmfile for Migrate", 20L))
493 {
494 usererror("This file is not a parameter file (parmfile) for migrate\nSyntax for calling from the commandline is\nmigrate-n parmfile\n and not (!!!!)\nmigrate-n datafile\n");
495 exit(-1);
496 }
497 fprintf(stdout,"Reading parmfile \"%s\"....\n",options->parmfilename);
498 while (FGETS (input, LINESIZE, options->parmfile) != EOF)
499 {
500 counter++;
501 if ((input[0] == '#') || isspace ((int) input[0])
502 || input[0] == ';')
503 {
504 //position = ftell (options->parmfile);
505 fgetpos(options->parmfile, &thePos);
506 continue;
507 }
508 else
509 {
510 if (!(isalnum ((int) input[0]) || strchr ("{}", input[0])))
511 {
512 usererror ("The parmfile contains an error on line %li\n",
513 counter);
514 }
515 }
516 if ((p = strchr (input, '#')) != NULL)
517 *p = '\n';
518 if (!strncmp (input, "end", 3))
519 break;
520 tmp = strtok (input, "=");
521 if (tmp != NULL)
522 strncpy (parmvar, tmp, linecpmax);
523 else
524 {
525 if(input[0]!='\0')
526 fprintf (stderr,
527 "WARNING: error in parmfile on line %li with %s\n",
528 counter, input);
529 continue;
530 }
531
532 // DEBUG check that we really read stuff
533 #ifdef DEBUG
534 fprintf(stdout,"%s\n", parmvar);
535 #endif
536 if (!strncmp (parmvar, "theta", 5))
537 {
538 //reset_oneline (options, position);
539 reset_oneline (options, &thePos);
540 read_theta (options);
541 //position = ftell (options->parmfile);
542 fgetpos(options->parmfile, &thePos);
543 continue;
544 }
545 if (!strncmp (parmvar, "migration", 5))
546 {
547 //reset_oneline (options, position);
548 reset_oneline (options, &thePos);
549 read_mig (options);
550 //position = ftell (options->parmfile);
551 fgetpos(options->parmfile, &thePos);
552 continue;
553 }
554 tmp = strtok (NULL, "\n");
555 if (tmp != NULL)
556 {
557 if(strlen(tmp)<LONGLINESIZE)
558 strcpy (varvalue, tmp);
559 else
560 strncpy (varvalue, tmp,LONGLINESIZE-1);
561 }
562 if (!booleancheck (options, parmvar, varvalue))
563 {
564 if (!numbercheck (options, parmvar, varvalue))
565 {
566 warning ("Inappropiate entry in parmfile: %s ignored\n",
567 input);
568 }
569 }
570 //position = ftell (options->parmfile);
571 fgetpos(options->parmfile, &thePos);
572 }
573 }
574 }
575
576 #ifdef MPI
577 void
read_options_worker(char ** buffer,option_fmt * options)578 read_options_worker (char **buffer, option_fmt * options)
579 {
580 const unsigned long linecpmax = LINESIZE - 1;
581
582 long counter = 0;
583 char *position;
584 char varvalue[LONGLINESIZE];
585 char parmvar[LONGLINESIZE];
586 char input[LONGLINESIZE];
587 char *p, *tmp;
588
589 options->buffer = buffer;
590 options->parmfile = NULL;
591 if (strlen (*buffer) > 0)
592 {
593 counter = 0;
594 position = *buffer;
595 while (sgets (input, LINESIZE, buffer) != NULL)
596 {
597 counter++;
598 if ((input[0] == '#') || isspace ((int) input[0])
599 || input[0] == ';')
600 {
601 position = *buffer;
602 continue;
603 }
604 else
605 {
606 if (!(isalnum ((int) input[0]) || strchr ("{}", input[0])))
607 {
608 usererror ("The parmfile contains an error on line %li\n",
609 counter);
610 }
611 }
612 if ((p = strchr (input, '#')) != NULL)
613 *p = '\n';
614 if (!strncmp (input, "end", 3))
615 break;
616 tmp = strtok (input, "=");
617 if (tmp != NULL)
618 strncpy (parmvar, tmp, linecpmax);
619 else
620 {
621 if(input[0]!='\0')
622 fprintf (stderr,
623 "WARNING: error in parmfile on line %li with %s\n",
624 counter, input);
625 continue;
626 }
627 if (!strncmp (parmvar, "theta", 5))
628 {
629 *buffer = position;
630 read_theta_worker (buffer, options);
631 position = *buffer;
632 continue;
633 }
634 if (!strncmp (parmvar, "migration", 5))
635 {
636 *buffer = position;
637 read_mig_worker (buffer, options);
638 position = *buffer;
639 continue;
640 }
641 tmp = strtok (NULL, "\n");
642 if (tmp != NULL)
643 strncpy (varvalue, tmp, LINESIZE - 1);
644 if (!booleancheck (options, parmvar, varvalue))
645 {
646 if (!numbercheck (options, parmvar, varvalue))
647 {
648 warning ("Inappropiate entry in parmfile: %s ignored\n",
649 input);
650 }
651 }
652 position = *buffer;
653 }
654 }
655 if(options->bayes_infer)
656 options->schains=0;
657 #ifdef DEBUG_MPI
658 printf("%i> custm=%s custm2=%s\n",myID,options->custm, options->custm2);
659 printf("%i> generation/year=%f mutationrate/year=%f\n",myID, options->generation_year, options->mutationrate_year[0]);
660 #endif
661
662 }
663 #endif
664
665 void
print_menu_options(world_fmt * world,option_fmt * options,data_fmt * data)666 print_menu_options (world_fmt * world, option_fmt * options, data_fmt * data)
667 {
668 if (options->numpop > data->numpop)
669 usererror ("Inconsistency between your Menu/Parmfile and your datafile\n \
670 Check the number of populations!\n");
671 if (options->progress)
672 {
673 print_options (stdout, world, options, data);
674 if (options->writelog)
675 print_options (options->logfile, world, options, data);
676 }
677 }
678
679 ///
680 /// prints the value of minimum value of the prior distribution
show_priormin(char * tmp,prior_fmt * prior,int priortype)681 char * show_priormin(char *tmp, prior_fmt *prior, int priortype)
682 {
683 sprintf(tmp, "%f ", prior->min);
684 return tmp;
685 }
686
687 ///
688 /// prints the value of mean value of the prior distribution
689 /// or for the gamma proposal the alpha value, and the ,multipler proposal the multiplicator
show_priormean(char * tmp,prior_fmt * prior,int priortype)690 char * show_priormean(char *tmp, prior_fmt *prior, int priortype)
691 {
692 sprintf(tmp, "%f ", prior->mean);
693 return tmp;
694 }
695
696 /// prints the value of maximum value of the prior distribution
show_priormax(char * tmp,prior_fmt * prior,int priortype)697 char * show_priormax(char *tmp, prior_fmt *prior, int priortype)
698 {
699 sprintf(tmp, "%f ", prior->max);
700 return tmp;
701 }
702
703 /// prints the alpha value of the gamma prior distribution
show_prioralpha(char * tmp,prior_fmt * prior,int priortype)704 char * show_prioralpha(char *tmp, prior_fmt *prior, int priortype)
705 {
706 sprintf(tmp, "%f ", prior->alpha);
707 return tmp;
708 }
709
710 ///
711 /// prints the value of delta value of the prior distribution
712 /// for some distribution this has no meaning
show_priordelta(char * tmp,prior_fmt * prior,int priortype)713 char * show_priordelta(char *tmp, prior_fmt *prior, int priortype)
714 {
715 switch(priortype)
716 {
717 case EXPPRIOR:
718 sprintf(tmp, "-"); break;
719 default:
720 sprintf(tmp, "%f ", prior->delta);
721 break;
722 }
723 return tmp;
724 }
725
726 ///
727 /// prints the bins used for the prior/posterior distribution
show_priorbins(char * tmp,prior_fmt * prior,int priortype)728 char * show_priorbins(char *tmp, prior_fmt *prior, int priortype)
729 {
730 sprintf(tmp, "%li ", prior->bins);
731 return tmp;
732 }
733
734
735 void
print_options(FILE * file,world_fmt * world,option_fmt * options,data_fmt * data)736 print_options (FILE * file, world_fmt * world, option_fmt * options,
737 data_fmt * data)
738 {
739 int count;
740 long i, j, tt;
741 char mytext[LINESIZE];
742 char mytext1[LINESIZE];
743 char mytext2[LINESIZE];
744 char mytext3[LINESIZE];
745 char mytext4[LINESIZE];
746 char seedgen[LINESIZE], spacer[LINESIZE];
747 char paramtgen[LINESIZE], parammgen[LINESIZE];
748 if (options->datatype != 'g')
749 {
750 switch ((short) options->autoseed)
751 {
752 case AUTO:
753 strcpy (seedgen, "with internal timer");
754 strcpy (spacer, " ");
755 break;
756 case NOAUTOSELF:
757 strcpy (seedgen, "from parmfile");
758 strcpy (spacer, " ");
759 break;
760 case NOAUTO:
761 strcpy (seedgen, "from seedfile");
762 strcpy (spacer, " ");
763 break;
764 default:
765 strcpy (seedgen, "ERROR");
766 strcpy (spacer, " ");
767 break;
768 }
769 switch (options->thetaguess)
770 {
771 case OWN:
772 strcpy (paramtgen, "from guessed values");
773 break;
774 case FST:
775 strcpy (paramtgen, "from an FST-calculation");
776 break;
777 // case PARAMGRID:
778 // strcpy (paramtgen, "GRID values around a center");
779 // break;
780 case NRANDOMESTIMATE:
781 strcpy (paramtgen, "RANDOM start value from N(mean,std) or U(min,max)");
782 break;
783 case URANDOMESTIMATE:
784 strcpy (paramtgen, "RANDOM start value from U(min,max)");
785 break;
786 default:
787 strcpy (paramtgen, "ERROR");
788 break;
789 }
790 switch (options->migrguess)
791 {
792 case OWN:
793 strcpy (parammgen, "from guessed values");
794 break;
795 case FST:
796 strcpy (parammgen, "from the FST-calculation");
797 break;
798 // case PARAMGRID:
799 // strcpy (parammgen, "GRID values around a center");
800 // break;
801 case NRANDOMESTIMATE:
802 strcpy (parammgen, "RANDOM start value from N(mean,std)");
803 break;
804 case URANDOMESTIMATE:
805 strcpy (parammgen, "RANDOM start value from U(min,max)");
806 break;
807 default:
808 strcpy (parammgen, "ERROR");
809 break;
810 }
811 }
812 fprintf (file, "Options in use:\n");
813 fprintf (file, "---------------\n\n");
814 if(options->bayes_infer)
815 {
816 fprintf (file, "Analysis strategy is BAYESIAN INFERENCE\n\n");
817 fprintf (file, "Proposal distribution:\n");
818 fprintf (file, "Parameter group Proposal type\n");
819 fprintf (file, "----------------------- -------------------\n");
820 fprintf (file, "Population size (Theta) %20s\n",
821 is_proposaltype(options->slice_sampling[THETAPRIOR]));
822 if(world->numpop > 1)
823 {
824 fprintf (file, "Migration rate %7.7s %20s\n",
825 options->usem ? "(M)" : "(xNm)",
826 is_proposaltype(options->slice_sampling[MIGPRIOR]));
827 }
828 if(world->bayes->mu)
829 {
830 fprintf (file, "Mutation rate modifier %20s\n",
831 is_proposaltype(options->slice_sampling[RATEPRIOR]));
832 }
833 fprintf(file,"\n\n");
834 if(world->options->has_autotune)
835 fprintf (file, "Prior distribution (Proposal-delta will be tuned to acceptance frequence %f):\n",world->options->autotune);
836 else
837 {
838 fprintf (file, "Prior distribution:\n");
839 }
840 fprintf (file, "Parameter group Prior type Minimum Mean(*) Maximum Delta\n");
841 fprintf (file, "----------------------- ------------ ---------- ---------- ---------- ----------\n");
842 fprintf (file, "Population size (Theta) %12s %10.10s %10.10s %10.10s %10.10s\n",
843 is_priortype(options->bayesprior[THETAPRIOR]),
844 show_priormin(mytext1, options->bayespriortheta,options->bayesprior[THETAPRIOR]),
845 show_priormean(mytext2, options->bayespriortheta,options->bayesprior[THETAPRIOR]),
846 show_priormax(mytext3, options->bayespriortheta,options->bayesprior[THETAPRIOR]),
847 show_priordelta(mytext4, options->bayespriortheta,options->bayesprior[THETAPRIOR]));
848 if(world->numpop > 1)
849 {
850 fprintf (file, "Migration rate %7.7s %12s %10.10s %10.10s %10.10s %10.10s\n",
851 options->usem ? "(M)" : "(xNm)",
852 is_priortype(options->bayesprior[MIGPRIOR]),
853 show_priormin(mytext1,options->bayespriorm,options->bayesprior[MIGPRIOR]),
854 show_priormean(mytext2,options->bayespriorm,options->bayesprior[MIGPRIOR]),
855 show_priormax(mytext3,options->bayespriorm,options->bayesprior[MIGPRIOR]),
856 show_priordelta(mytext4,options->bayespriorm,options->bayesprior[MIGPRIOR]));
857 }
858 if(world->bayes->mu)
859 {
860 fprintf (file, "Mutation rate modifier %12s %10.10s %10.10s %10.10s %10.10s\n",
861 is_priortype(options->bayesprior[RATEPRIOR]),
862 show_priormin(mytext1,options->bayespriorrate,options->bayesprior[RATEPRIOR]),
863 show_priormean(mytext2,options->bayespriorrate,options->bayesprior[RATEPRIOR]),
864 show_priormax(mytext3,options->bayespriorrate,options->bayesprior[RATEPRIOR]),
865 show_priordelta(mytext4,options->bayespriorrate,options->bayesprior[RATEPRIOR]));
866 }
867 fprintf(file,"\n\n\n");
868 }
869 else
870 {
871 fprintf (file, "Analysis strategy is MAXIMUM LIKELIHOOD\n\n\n\n");
872 }
873
874 switch (options->datatype)
875 {
876 case 'a':
877 fprintf (file, "Datatype: Allelic data\n");
878 fprintf (file, "Missing data is %s\n",options->include_unknown ? "included" : "not included");
879 break;
880 case 'b':
881 fprintf (file, "Datatype: Microsatellite data [Brownian motion]\n");
882 fprintf (file, "Missing data is %s\n",options->include_unknown ? "included" : "not included");
883 break;
884 case 'm':
885 if(options->msat_option==SINGLESTEP)
886 fprintf (file, "Datatype: Microsatellite data [Singlestep model]\n");
887 else
888 fprintf (file, "Datatype: Microsatellite data [Multistep model (Tune=%f, P_increase=%f)]\n",options->msat_tuning[0], options->msat_tuning[1]);
889 fprintf (file, "Missing data is %s\n",options->include_unknown ? "included" : "not included");
890 break;
891 case 's':
892 fprintf (file, "Datatype: DNA sequence data\n");
893 break;
894 case 'n':
895 fprintf (file, "Datatype: Single nucleotide polymorphism data\n");
896 break;
897 case 'h':
898 fprintf (file, "Datatype: Hapmap Single nucleotide polymorphism data\n");
899 break;
900 case 'u':
901 fprintf (file,
902 "Datatype: Single nucleotide polymorphism data (PANEL)\n");
903 break;
904 case 'f':
905 fprintf (file, "Datatype: Ancestral state method\n");
906 break;
907 case 'g':
908 fprintf (file, "Datatype: Genealogy summary of an older run\n");
909 break;
910 }
911
912 count=0;
913
914 fprintf (file, "\nInheritance scalers in use for Thetas (specified scalars=%li)\n",options->inheritance_scalars_numalloc);
915
916 for(i=0;i<data->loci;i++)
917 {
918 if(i<options->inheritance_scalars_numalloc)
919 fprintf (file, "%2.2f ", options->inheritance_scalars[i]);
920 else
921 fprintf (file, "%2.2f ", options->inheritance_scalars[options->inheritance_scalars_numalloc-1]);
922 if(count++ > 7)
923 {
924 fprintf(file,"\n");
925 count=0;
926 }
927 else
928 {
929 count++;
930 }
931 }
932 fprintf (file, "\n[Each Theta uses the (true) ineritance scalar of the first locus as a reference]\n\n");
933
934 if(options->randomsubset > 0)
935 {
936 if(options->randomsubsetseed>0)
937 fprintf (file, "\nData set was subsampled: used a random sample of size: %li\nwith private random number seed %li\n\n", options->randomsubset, options->randomsubsetseed);
938 else
939 fprintf (file, "\nData set was subsampled: used a random sample of size: %li\n\n", options->randomsubset);
940 }
941
942 if (options->datatype != 'g')
943 {
944 fprintf (file, "\n%-80s\n", generator);
945 #ifndef QUASIRANDOM
946 fprintf (file, "Random number seed (%s)%s%20li\n", seedgen, " ",
947 options->saveseed);
948 #endif
949 fprintf (file, "\nStart parameters:\n");
950 fprintf(file, " First genealogy was started using a %s\n",
951 options->usertree ? "user tree" : (options->randomtree ? "random tree" :
952 (options->dist ? "tree from userdistances" :
953 "UPGMA-tree")));
954 fprintf(file, " Theta values were generated ");
955 fprintf (file, " %s\n", paramtgen);
956 if (options->thetaguess == OWN)
957 {
958 fprintf (file, " Theta = ");
959 for (i = 0; i < options->numthetag - 1; i++)
960 {
961 fprintf (file, "%.5f,", options->thetag[i]);
962 }
963 fprintf (file, "%.5f\n", options->thetag[i]);
964 }
965 fprintf (file, " M values were generated %s\n", parammgen);
966 if (options->migrguess == OWN)
967 {
968 tt = 0;
969 if (options->usem)
970 fprintf (file, " M-matrix: ");
971 else
972 fprintf (file, " 4Nm-matrix: ");
973 if (options->nummg == 1)
974 {
975 fprintf (file, "%5.2f [all are the same]\n", options->mg[0]);
976 }
977 else
978 {
979 for (i = 0; i < world->numpop; i++)
980 {
981 for (j = 0; j < world->numpop; j++)
982 {
983 if (i != j)
984 {
985 fprintf (file, "%5.2f ", options->mg[tt]);
986 if(tt<options->nummg)
987 tt++;
988 }
989 else
990 fprintf (file, "----- ");
991 if (j > 10)
992 fprintf (file, "\n ");
993 }
994 fprintf (file, "\n ");
995 }
996 fprintf (file, "\n");
997 }
998 }
999 }
1000 fprintf(file,"\n");
1001 print_arbitrary_migration_table (file, world, options, data);
1002 print_distance_table (file, world, options, data);
1003 fprintf(file,"\n");
1004 // mutation related material
1005 if (options->gamma)
1006 {
1007 fprintf (file, "Mutation rate among loci is Gamma-distributed\n");
1008 fprintf (file, "Initial scale parameter alpha = %f\n",
1009 options->alphavalue);
1010 if (options->custm[world->numpop2] == 'c')
1011 fprintf (file, "and is constant [will not be estimated]\n");
1012 }
1013 else /*Gamma*/
1014 {
1015 if (options->murates && world->loci > 1)
1016 {
1017 fprintf (file, "Mutation rate among loci is varying with\n");
1018 fprintf (file, " Rates per locus: ");
1019 for (i = 0; i < world->loci - 1; i++)
1020 {
1021 fprintf (file, "%.3f, ", options->mu_rates[i]);
1022 if ((i + 1) % 6 == 0)
1023 fprintf (file, "\n ");
1024 }
1025 fprintf (file, "%.3f\n", options->mu_rates[i]);
1026 if(options->murates_fromdata)
1027 fprintf (file, "[Estimated from the data using the Watterson estimator (ignoring migration)]");
1028 else
1029 fprintf (file, "[User defined]");
1030 }
1031 else
1032 fprintf (file, "Mutation rate is constant %s\n", world->loci > 1 ?
1033 "for all loci" : "");
1034 }
1035 #ifdef UEP
1036 if (options->uep)
1037 {
1038 fprintf (file, "0/1 polymorphism analysis, with 0/1 data in file:%s\n",
1039 options->uepfilename);
1040 fprintf (file, " with forward mutation rate %f*mu\n",
1041 options->uepmu);
1042 fprintf (file, " with back mutation rate %f*mu\n",
1043 options->uepnu);
1044 fprintf (file, " with base frequencies \"0\"=%f and \"1\"=%f\n",
1045 options->uepfreq0,options->uepfreq1);
1046 }
1047 #endif /*UEP*/
1048 if (options->datatype != 'g')
1049 {
1050 fprintf (file, "\nMarkov chain settings:\n");
1051 if(!options->bayes_infer)
1052 {
1053 fprintf (file, " Short chains (short-chains): %20li\n",
1054 options->schains);
1055 fprintf (file, " Trees sampled (short-inc*samples):%20li\n",
1056 options->sincrement * options->ssteps);
1057 fprintf (file, " Trees recorded (short-sample): %20li\n",
1058 options->ssteps);
1059 }
1060 fprintf (file, " Long chains (long-chains): %20li\n",
1061 options->lchains);
1062 if(options->bayes_infer)
1063 {
1064 if(options->replicate)
1065 {
1066 fprintf (file, " Steps sampled (inc*samples*rep): %20li\n",
1067 options->lincrement * options->lsteps * options->replicatenum);
1068 fprintf (file, " Steps recorded (sample*rep): %20li\n",
1069 options->lsteps*options->replicatenum);
1070 }
1071 else
1072 {
1073 fprintf (file, " Steps sampled (long-inc*samples): %20li\n",
1074 options->lincrement * options->lsteps);
1075 fprintf (file, " Steps recorded (long-sample): %20li\n",
1076 options->lsteps);
1077 }
1078 }
1079 else
1080 {
1081 fprintf (file, " Trees sampled (long-inc*samples): %20li\n",
1082 options->lincrement * options->lsteps);
1083 fprintf (file, " Trees recorded (long-sample): %20li\n",
1084 options->lsteps);
1085 }
1086 if (options->replicate)
1087 {
1088 if (options->replicatenum == 0)
1089 fprintf (file, " Averaging over long chains\n");
1090 else
1091 {
1092 if(!options->bayes_infer)
1093 fprintf (file, " Averaging over replicates: %20li\n",
1094 options->replicatenum);
1095 else
1096 fprintf (file, " Combining over replicates: %20li\n",
1097 options->replicatenum);
1098 }
1099 }
1100 if (options->heating > 0)
1101 {
1102 fprintf (file,
1103 " %s heating scheme\n %li chains with %s temperatures\n ",
1104 options->adaptiveheat!=NOTADAPTIVE ? (options->adaptiveheat==STANDARD ? "Adaptive_Standard" : "Bounded_adaptive") : "Static", options->heated_chains, options->adaptiveheat ? "start" : "" );
1105 for (i = 0; i < options->heated_chains - 1; i++)
1106 fprintf (file, "%5.2f,", options->heat[i]);
1107 if(!options->heatedswap_off)
1108 fprintf (file, "%5.2f\n Swapping interval is %li\n",
1109 options->heat[i], options->heating_interval);
1110 else
1111 fprintf (file, "%5.2f\n No swapping\n",
1112 options->heat[i]);
1113 }
1114 if (options->movingsteps)
1115 {
1116 fprintf (file, " Forcing at least this percentage of new genealogies:%6.2f\n",
1117 (MYREAL) options->acceptfreq);
1118 }
1119 if (options->burn_in > 0)
1120 {
1121 switch (options->burnin_autostop)
1122 {
1123 case 'a':
1124 sprintf(mytext,"(var)%10li", options->burn_in);
1125 break;
1126 case 't':
1127 sprintf(mytext,"(acc)%10li", options->burn_in);
1128 break;
1129 case 'e':
1130 sprintf(mytext,"(ESS)%10li", options->burn_in);
1131 break;
1132 case ' ':
1133 default:
1134 sprintf(mytext,"%10li", options->burn_in * options->lincrement);
1135 }
1136 if (options->bayes_infer)
1137 fprintf (file, " Burn-in per replicate (samples*inc): %20.20s\n", mytext);
1138 else
1139 fprintf (file, " Burn-in per chain: %35.35s\n", mytext);
1140 }
1141 if (options->lcepsilon < LONGCHAINEPSILON)
1142 {
1143 fprintf (file, " Parameter-likelihood epsilon: %20.5f\n",
1144 options->lcepsilon);
1145 }
1146 }
1147 fprintf (file, "\nPrint options:\n");
1148 if (options->datatype != 'g')
1149 {
1150 fprintf (file, " Data file: %46.46s\n", options->infilename);
1151 #ifdef PRETTY
1152 fprintf (file, " Output file (ASCII text): %31.31s\n", options->outfilename);
1153 fprintf (file, " Output file (PDF): %31.31s\n", options->pdfoutfilename);
1154 #else
1155 fprintf (file, " Output file: %44.44s\n", options->outfilename);
1156 #endif
1157 if(options->bayes_infer)
1158 {
1159 fprintf (file, " Posterior distribution: %33.33s\n", options->bayesfilename);
1160 if(options->has_bayesmdimfile)
1161 fprintf (file, " All values of Post.Dist:%33.33s\n", options->bayesmdimfilename);
1162 }
1163 fprintf (file, " Print data: %45.45s\n",
1164 options->printdata ? "Yes" : "No");
1165 switch (options->treeprint)
1166 {
1167 case _NONE:
1168 fprintf (file, " Print genealogies: %38.38s\n", "No");
1169 break;
1170 case ALL:
1171 fprintf (file, " Print genealogies: %38.38s\n", "Yes, all");
1172 break;
1173 case LASTCHAIN:
1174 fprintf (file, " Print genealogies: %32.32s%li\n",
1175 "Yes, only those in last chain, every ", options->treeinc);
1176 break;
1177 case BEST:
1178 fprintf (file, " Print genealogies: %38.38s\n",
1179 "Yes, only the best");
1180 break;
1181 }
1182 }
1183 if (options->plot)
1184 {
1185 switch (options->plotmethod)
1186 {
1187 case PLOTALL:
1188 sprintf (mytext, "Yes, to outfile and %s", options->mathfilename);
1189 break;
1190 default:
1191 strcpy (mytext, "Yes, to outfile");
1192 break;
1193 }
1194 fprintf (file, " Plot data: %-46.46s\n", mytext);
1195 fprintf (file,
1196 " Parameter: %s, Scale: %s, Intervals: %li\n",
1197 options->plotvar == PLOT4NM ? "{Theta, 4Nm}" : "{Theta, M}",
1198 options->plotscale == PLOTSCALELOG ? "Log10" : "Standard",
1199 options->plotintervals);
1200 fprintf (file, " Ranges: X-%5.5s: %f - %f\n",
1201 options->plotvar == PLOT4NM ? "4Nm" : "M",
1202 options->plotrange[0], options->plotrange[1]);
1203 fprintf (file, " Ranges: Y-%5.5s: %f - %f\n", "Theta",
1204 options->plotrange[2], options->plotrange[3]);
1205 }
1206 else
1207 {
1208 fprintf (file, " Plot data: %-46.46s\n", "No");
1209 }
1210
1211 if (options->mighist)
1212 {
1213 if(options->mighist_all)
1214 sprintf(mytext,"Yes: All events");
1215 else
1216 sprintf(mytext,"Yes: Migration events");
1217 fprintf (file,
1218 " Frequency histogram of events %26.26s\n", mytext);
1219 fprintf (file,
1220 " Time of events are saved in file %24.24s\n",
1221 options->mighistfilename);
1222 }
1223 if (options->mighist && options->skyline)
1224 {
1225 sprintf(mytext,"%3s","Yes");
1226 fprintf (file,
1227 " Histogram of the parameter values through time %10s\n", mytext);
1228 fprintf (file,
1229 " Parameters through time are saved in file %15.15s\n",
1230 options->skylinefilename);
1231 }
1232 if(!options->bayes_infer)
1233 {
1234 switch (options->profile)
1235 {
1236 case _NONE:
1237 strcpy (mytext, "No");
1238 break;
1239 case ALL:
1240 strcpy (mytext, "Yes, tables and summary");
1241 break;
1242 case TABLES:
1243 strcpy (mytext, "Yes, tables");
1244 break;
1245 case SUMMARY:
1246 strcpy (mytext, "Yes, summary");
1247 break;
1248 }
1249 fprintf (file, " Profile likelihood: %-36.36s\n", mytext);
1250 if (options->profile != _NONE)
1251 {
1252 switch (options->profilemethod)
1253 {
1254 case 'p':
1255 fprintf (file, " Percentile method\n");
1256 break;
1257 case 'q':
1258 fprintf (file, " Quick method\n");
1259 break;
1260 case 'f':
1261 fprintf (file, " Fast method\n");
1262 break;
1263 case 'd':
1264 fprintf (file, " Discrete method\n");
1265 break;
1266 case 's':
1267 fprintf (file, " Spline method\n");
1268 break;
1269 default:
1270 fprintf (file, " UNKOWN method????\n");
1271 break;
1272 }
1273 fprintf (file, " with df=%li and for Theta and %s\n\n\n\n",
1274 options->df, options->profileparamtype ? "M=m/mu" : "4Nm");
1275 }
1276 }
1277 }
1278
1279 /// \brief sets the theta parameters for random start parameter setting
1280 ///
1281 /// Sets the startparameter to a random number derived from a uniform distribution
1282 /// with minimum options->thetag[0] and maximum in options->thetag[1]
set_theta_urandomstart(world_fmt * world,option_fmt * options)1283 void set_theta_urandomstart(world_fmt *world, option_fmt *options)
1284 {
1285 long i;
1286 long pos=0;
1287 char temp[SUPERLINESIZE];
1288 if(world->options->progress)
1289 pos = sprintf(temp,"Random start Theta values:");
1290 for (i = 0; i < world->numpop; i++)
1291 {
1292 world->param0[i] = options->thetag[0] + RANDUM() * (options->thetag[1]-options->thetag[0]);
1293 if(world->options->progress)
1294 pos += sprintf(temp+pos," %f",world->param0[i]);
1295 }
1296 if(world->options->progress)
1297 FPRINTF(stdout,"[%3i] %s\n",myID, temp);
1298 if(world->options->writelog)
1299 FPRINTF(world->options->logfile,"[%3i] %s\n",myID, temp);
1300 }
1301
1302 /// \brief sets the theta parameters for random start parameter setting
1303 ///
1304 /// Sets the startparameter to a random number derived from a normal distribution
1305 /// with men options->thetag[0] and standard deviation options-.thetag[1]
set_theta_nrandomstart(world_fmt * world,option_fmt * options)1306 void set_theta_nrandomstart(world_fmt *world, option_fmt *options)
1307 {
1308 long i;
1309 long pos=0;
1310 char temp[SUPERLINESIZE];
1311 if(world->options->progress)
1312 pos = sprintf(temp,"Random start Theta values:");
1313 for (i = 0; i < world->numpop; i++)
1314 {
1315 world->param0[i] = rannor (options->thetag[0], options->thetag[1]);
1316 while (world->param0[i] < 0)
1317 {
1318 world->param0[i] =
1319 rannor (options->thetag[0], options->thetag[1]);
1320 }
1321 if(world->options->progress)
1322 pos += sprintf(temp+pos," %f",world->param0[i]);
1323 }
1324 if(world->options->progress)
1325 FPRINTF(stdout,"%i> %s\n",myID, temp);
1326 if(world->options->writelog)
1327 FPRINTF(world->options->logfile,"%i> %s\n",myID, temp);
1328 }
1329
1330
1331 /// \brief sets the theta parameters for OWN start parameter setting
1332 ///
1333 /// Sets the startparameter to a user-defined value
set_theta_ownstart(world_fmt * world,option_fmt * options)1334 void set_theta_ownstart(world_fmt *world, option_fmt *options)
1335 {
1336 long i, ii;
1337 for (i = 0; i < world->numpop; i++)
1338 {
1339 if (i < options->numthetag - 1)
1340 ii = i;
1341 else
1342 ii = options->numthetag - 1;
1343 if (options->thetag[ii] == 0.0)
1344 world->param0[i] = SMALLEST_THETA;
1345 else
1346 {
1347 world->param0[i] = options->thetag[ii];
1348 }
1349 }
1350 }
1351
1352
1353
1354 /// \brief sets the theta parameters to a start parameter using an FST value
1355 ///
1356 /// Sets the startparameter to a value from the FST calculation
set_theta_fststart(world_fmt * world,option_fmt * options,long locus)1357 void set_theta_fststart(world_fmt *world, option_fmt *options, long locus)
1358 {
1359 long i;
1360 for (i = 0; i < world->numpop; i++)
1361 {
1362 if (world->fstparam[locus][i] > SMALLEST_THETA)
1363 {
1364 if (world->fstparam[locus][i] > 100)
1365 world->param0[i] = 1.0;
1366 else
1367 world->param0[i] = world->fstparam[locus][i];
1368 }
1369 else
1370 {
1371 if (strchr (SEQUENCETYPES, options->datatype))
1372 {
1373 world->param0[i] = 0.01;
1374 }
1375 else
1376 {
1377 world->param0[i] = 1.0;
1378 }
1379 }
1380 }
1381 }
1382
1383
1384 /// \brief sets the migration parameters for random start parameter setting
1385 ///
1386 /// Sets the startparameter to a random number derived from a uniform distribution
1387 /// with minimum options->mg[0] and maximum in options->mg[1]
set_mig_urandomstart(world_fmt * world,option_fmt * options)1388 void set_mig_urandomstart(world_fmt *world, option_fmt *options)
1389 {
1390 long i;
1391 long pos=0;
1392 char temp[SUPERLINESIZE];
1393 if(world->options->progress)
1394 pos = sprintf(temp,"Random start M values:");
1395
1396 for (i = world->numpop; i < world->numpop2; i++)
1397 {
1398 world->param0[i] = options->mg[0] + RANDUM() * (options->mg[1]-options->mg[0]);
1399 if(world->options->progress)
1400 pos += sprintf(temp+pos," %f",world->param0[i]);
1401 }
1402 if(world->options->progress)
1403 FPRINTF(stdout,"%i> %s\n",myID, temp);
1404 if(world->options->writelog)
1405 FPRINTF(world->options->logfile,"%i> %s\n",myID, temp);
1406 }
1407
1408 /// \brief sets the migration parameters for random start parameter setting
1409 ///
1410 /// Sets the startparameter to a random number derived from a normal distribution
1411 /// with men options->mg[0] and standard deviation options->mg[1]
set_mig_nrandomstart(world_fmt * world,option_fmt * options)1412 void set_mig_nrandomstart(world_fmt *world, option_fmt *options)
1413 {
1414 long i;
1415 long pos=0;
1416 char temp[SUPERLINESIZE];
1417 if(world->options->progress)
1418 pos = sprintf(temp,"Random start M values:");
1419
1420 for (i = world->numpop; i < world->numpop2; i++)
1421 {
1422 world->param0[i] = rannor (options->mg[0], options->mg[1]);
1423 while (world->param0[i] <= 0.0)
1424 {
1425 world->param0[i] = rannor (options->mg[0], options->mg[1]);
1426 }
1427 if(world->options->progress)
1428 pos += sprintf(temp+pos," %f",world->param0[i]);
1429 }
1430 if(world->options->progress)
1431 FPRINTF(stdout,"%i> %s\n",myID, temp);
1432 if(world->options->writelog)
1433 FPRINTF(world->options->logfile,"%i> %s\n",myID, temp);
1434 }
1435
1436 /// \brief sets the migration parameters for OWN start parameter setting
1437 ///
1438 /// Sets the startparameter to a user-defined value
set_mig_ownstart(world_fmt * world,option_fmt * options,data_fmt * data)1439 void set_mig_ownstart(world_fmt *world, option_fmt *options, data_fmt *data)
1440 {
1441 long i;
1442 long j;
1443 long ii;
1444 long iitest;
1445 long numpop = world->numpop;
1446 MYREAL tmp;
1447 for (i = 0; i < numpop; i++)
1448 {
1449 for (j = 0; j < numpop; j++)
1450 {
1451 if (i == j)
1452 continue;
1453 if ((iitest = mm2m (j, i, numpop) - numpop) < options->nummg)
1454 ii = iitest;
1455 else
1456 {
1457 ii = options->nummg - 1;
1458 }
1459 //@@if (options->usem)
1460 tmp = options->mg[ii];
1461 //@@else
1462 //@@tmp = options->mg[ii] / world->param0[i];
1463 iitest += numpop;
1464 if (options->geo)
1465 {
1466 world->param0[iitest] = data->geo[iitest] * tmp;
1467 }
1468 else
1469 {
1470 world->param0[iitest] = tmp;
1471 }
1472 }
1473 }
1474 }
1475
1476 /// \brief sets the migration parameters to a start parameter using an FST value
1477 ///
1478 /// Sets the startparameter to a value from the FST calculation
set_mig_fststart(world_fmt * world,option_fmt * options,long locus)1479 void set_mig_fststart(world_fmt *world, option_fmt *options, long locus)
1480 {
1481 long i;
1482
1483 for (i = world->numpop; i < world->numpop2; i++)
1484 {
1485 if (world->fstparam[locus][i] > 0)
1486 {
1487 if (world->fstparam[locus][i] > 100)
1488 {
1489 world->param0[i] =
1490 1.0 / world->param0[(i - world->numpop) /
1491 (world->numpop)];
1492 if (world->param0[i] > 10000)
1493 world->param0[i] = 10000;
1494 }
1495 else
1496 world->param0[i] = world->fstparam[locus][i];
1497 }
1498 else
1499 {
1500 world->param0[i] =
1501 1.0 / world->param0[(i - world->numpop) / (world->numpop)];
1502 if (world->param0[i] > 10000)
1503 world->param0[i] = 10000;
1504 }
1505
1506 }
1507 }
1508
1509
1510 /// \brief set the starting parameters in main structure world from options
1511 ///
1512 /// Set the starting parameters in main structure world from options
1513 /// - Random starting parameters:\n
1514 /// will use an average and a standard deviation to generate a starting parameter
1515 /// - Own starting parameters
1516 /// - Starting parameters derived from FST
1517 /// - check that values from FST are not ridicoulously off and that Bayes start is in min/max bounds
1518 /// Once all parameters are set they are synchronized using the custom-migration matrix
1519 /// that will set some values to zero or to the same value. The synchronization will
1520 /// override user-error, but does provide little control about the outcome
1521 /// error control is only done later when user can check the values in the logfile
1522 void
set_param(world_fmt * world,data_fmt * data,option_fmt * options,long locus)1523 set_param (world_fmt * world, data_fmt * data, option_fmt * options,
1524 long locus)
1525 {
1526 // set THETA
1527 switch (options->thetaguess)
1528 {
1529 case NRANDOMESTIMATE:
1530 set_theta_nrandomstart(world, options);
1531 break;
1532
1533 case URANDOMESTIMATE:
1534 set_theta_urandomstart(world, options);
1535 break;
1536
1537 case OWN:
1538 set_theta_ownstart(world,options);
1539 break;
1540
1541 case FST:
1542 default:
1543 set_theta_fststart(world,options, locus);
1544 break;
1545 }
1546 // set MIGRATION
1547 switch (options->migrguess)
1548 {
1549 case NRANDOMESTIMATE:
1550 set_mig_nrandomstart(world,options);
1551 break;
1552
1553 case URANDOMESTIMATE:
1554 set_mig_urandomstart(world,options);
1555 break;
1556
1557 case OWN:
1558 set_mig_ownstart(world,options, data);
1559 break;
1560 case SLATKIN:
1561 case FST:
1562 default:
1563 set_mig_fststart(world,options, locus);
1564 break;
1565 }
1566 if(options->bayes_infer)
1567 {
1568 bayes_check_and_fix_param(world,options);
1569 }
1570 synchronize_param (world, options);
1571 if (options->gamma)
1572 world->options->alphavalue = options->alphavalue;
1573 //perhaps to come
1574 //if(options->thetaguess==PARAMGRID && options->migrguess==PARAMGRID)
1575 //set_grid_param(world,options->gridpoints);
1576 }
1577
1578 /// \todo not yet implemented, is this a good feature?
1579 void
set_grid_param(world_fmt * world,long gridpoints)1580 set_grid_param (world_fmt * world, long gridpoints)
1581 {
1582 static long which = 0;
1583 static long z = 0;
1584
1585 static MYREAL level = 0.1;
1586 static MYREAL bottom = -2.302585093; // => level = 0.1
1587 MYREAL top = 2.302585093; // => level = 10
1588 MYREAL len = top - bottom;
1589 MYREAL diff = len / (gridpoints - 1.);
1590
1591 if (z >= gridpoints)
1592 {
1593 z = 0;
1594 which++;
1595 }
1596 level = EXP (bottom + z * diff);
1597 world->param0[which] = world->param0[which] * level;
1598 z++;
1599
1600 }
1601
1602 /// \brief synchronize parameters
1603 void
synchronize_param(world_fmt * world,option_fmt * options)1604 synchronize_param (world_fmt * world, option_fmt * options)
1605 {
1606 char type;
1607 boolean found = FALSE;
1608 long i, z = 0, zz = 0, zzz = 0, len;
1609 long ii, jj;
1610 long ns = 0, ss = 0, ss2 = 0, xs = 0, migm = 0;
1611 boolean allsymmig = FALSE;
1612 boolean allsamesize = FALSE;
1613 boolean partmig = FALSE;
1614 MYREAL summ;
1615 MYREAL nsum = 0;
1616 len = (long) strlen (world->options->custm2);
1617 // not needed!
1618 // world->options->custm2 =
1619 // (char *) myrealloc (world->options->custm2,
1620 // sizeof (char) * (world->numpop2 + 2));
1621 if(world->options->bayes_infer)
1622 world->bayes->custm2 = world->options->custm2;
1623 if (len < world->numpop2)
1624 {
1625 fillup_custm (len, world, options);
1626 }
1627 migm = scan_connect (world->options->custm2, 0L, world->numpop, 'm');
1628 world->options->tmn = migm; // are there mean theta values?
1629 migm = scan_connect (world->options->custm2, world->numpop, world->numpop2, 'm');
1630 world->options->mmn = migm;
1631 for (i = 0; i < world->numpop2; i++)
1632 {
1633 if (!(allsymmig && allsamesize))
1634 {
1635 type = world->options->custm2[i];
1636 switch (type)
1637 {
1638 case '*':
1639 xs++;
1640 break;
1641 case 's': // M is symmetric
1642 if (i >= world->numpop)
1643 {
1644 z = i;
1645 m2mm (z, world->numpop, &ii, &jj);
1646 zz = mm2m (jj, ii, world->numpop);
1647 world->options->symparam = (twin_fmt *) myrealloc
1648 (world->options->symparam, sizeof (twin_fmt) * (ss + 2));
1649 world->options->symparam[ss][0] = zz;
1650 world->options->symparam[ss++][1] = z;
1651 world->options->symn = ss;
1652 summ = (world->param0[z] + world->param0[zz]) / 2.;
1653 world->param0[zz] = world->param0[z] = summ;
1654 }
1655 break;
1656 case 'S': // 4Nm is symmetric, not completely
1657 // implemented yet, -> derivatives.c
1658 if (i >= world->numpop)
1659 {
1660 z = i;
1661 m2mm (z, world->numpop, &ii, &jj);
1662 zz = mm2m (jj, ii, world->numpop);
1663 zzz = 0;
1664 found = FALSE;
1665 while (zzz < ss2)
1666 {
1667 if (world->options->sym2param[zzz][1] == zz)
1668 found = TRUE;
1669 zzz++;
1670 }
1671 if (found)
1672 break;
1673 world->options->sym2param = (quad_fmt *)
1674 myrealloc (world->options->sym2param,
1675 sizeof (quad_fmt) * (ss2 + 2));
1676 world->options->sym2param[ss2][0] = zz;
1677 world->options->sym2param[ss2][1] = z;
1678 world->options->sym2param[ss2][2] = ii;
1679 world->options->sym2param[ss2++][3] = jj;
1680 world->options->sym2n = ss2;
1681 summ = (world->param0[z] * world->param0[jj] +
1682 world->param0[zz] * world->param0[ii]) / 2.;
1683 world->param0[z] = summ / world->param0[jj];
1684 world->param0[zz] = summ / world->param0[ii];
1685 }
1686 break;
1687 case 'C':
1688 case 'c':
1689 world->options->constparam = (long *) myrealloc
1690 (world->options->constparam, sizeof (long) * (ns + 2));
1691 world->options->constparam[ns++] = i;
1692 world->options->constn = ns;
1693 break;
1694 case '0':
1695 z = i;
1696 world->param0[z] = 0;
1697 world->options->zeroparam = (long *) myrealloc
1698 (world->options->zeroparam, sizeof (long) * (ns + 2));
1699 world->options->zeroparam[ns++] = i;
1700 world->options->zeron = ns;
1701 break;
1702 case 'm':
1703 summ = 0;
1704 if (i < world->numpop) /*theta */
1705 {
1706 if (!allsamesize)
1707 {
1708 nsum = 0;
1709 allsamesize = TRUE;
1710 for (z = 0; z < world->numpop; z++)
1711 {
1712 if (world->options->custm2[z] == 'm')
1713 {
1714 nsum++;
1715 summ += world->param0[z];
1716 }
1717 }
1718 summ /= nsum;
1719 for (z = 0; z < world->numpop; z++)
1720 if (world->options->custm2[z] == 'm')
1721 world->param0[z] = summ;
1722 }
1723 }
1724 else /* migration */
1725 {
1726 summ = 0;
1727 if (!partmig)
1728 {
1729 nsum = 0;
1730 partmig = TRUE;
1731 for (z = world->numpop; z < world->numpop2; z++)
1732 {
1733 if (world->options->custm2[z] == 'm')
1734 {
1735 nsum++;
1736 summ += world->param0[z];
1737 }
1738 }
1739 summ /= nsum;
1740 for (z = world->numpop; z < world->numpop2; z++)
1741 if (world->options->custm2[z] == 'm')
1742 world->param0[z] = summ;
1743 }
1744 }
1745 break;
1746 case 'M':
1747 summ = 0;
1748 if (i < world->numpop) /*theta */
1749 {
1750 if (!allsamesize)
1751 {
1752 nsum = 0;
1753 allsamesize = TRUE;
1754 for (z = 0; z < world->numpop; z++)
1755 {
1756 if (world->options->custm2[z] == 'm')
1757 {
1758 nsum++;
1759 summ += world->param0[z];
1760 }
1761 }
1762 summ /= nsum;
1763 for (z = 0; z < world->numpop; z++)
1764 if (world->options->custm2[z] == 'm')
1765 world->param0[z] = summ;
1766 }
1767 }
1768 else /* migration */
1769 {
1770 summ = 0;
1771 if (!partmig)
1772 {
1773 nsum = 0;
1774 partmig = TRUE;
1775 for (z = world->numpop; z < world->numpop2; z++)
1776 {
1777 if (world->options->custm2[z] == 'M')
1778 {
1779 nsum++;
1780 m2mm (z, world->numpop, &ii, &jj);
1781 summ += world->param0[z] * world->param0[jj];
1782 }
1783 }
1784 summ /= nsum;
1785 for (z = world->numpop; z < world->numpop2; z++)
1786 {
1787 if (world->options->custm2[z] == 'M')
1788 {
1789 m2mm (z, world->numpop, &ii, &jj);
1790 world->param0[z] = summ/world->param0[jj];
1791 }
1792 }
1793 }
1794 }
1795 break;
1796 default:
1797 warning ("The migration connection matrix is misspecified\nOnly these are allowed\n* s S m M 0 (=zero)\nSupplied value was: %c\ncustom-migration=%s\n", type, world->options->custm);
1798 error("[Program exits]");
1799 }
1800 }
1801 }
1802 //--------gamma stuff
1803 if (world->options->gamma)
1804 {
1805 if (world->options->custm2[world->numpop2] == 'c')
1806 {
1807 world->options->constparam = (long *) myrealloc
1808 (world->options->constparam, sizeof (long) * (ns + 2));
1809 world->options->constparam[ns++] = i;
1810 world->options->constn = ns;
1811 }
1812 else
1813 {
1814 world->options->custm2[world->numpop2] = '*';
1815 world->options->custm2[world->numpop2 + 1] = '\0';
1816 }
1817 }
1818 }
1819
1820
1821 void
resynchronize_param(world_fmt * world)1822 resynchronize_param (world_fmt * world)
1823 {
1824 char type;
1825 long i, j = 0, z = 0, zz = 0;
1826 //long len;
1827 long ii, jj;
1828 long ns = 0, ss = 0, ss2 = 0, xs = 0, migm = 0;
1829 boolean allsymmig = FALSE;
1830 boolean allsamesize = FALSE;
1831 boolean partmig = FALSE;
1832 MYREAL summ;
1833 MYREAL nsum = 0;
1834 //xcode len = (long) strlen (world->options->custm2);
1835 world->options->symn = 0;
1836 world->options->sym2n = 0;
1837 world->options->zeron = 0;
1838 world->options->constn = 0;
1839 world->options->mmn = 0;
1840 migm = scan_connect (world->options->custm2, 0L, world->numpop, 'm');
1841 world->options->tmn = migm; // are there mean theta values?
1842 migm = scan_connect (world->options->custm2, world->numpop, world->numpop2, 'm');
1843 world->options->mmn = migm;
1844 for (i = 0; i < world->numpop2; i++)
1845 {
1846 if (!(allsymmig && allsamesize))
1847 {
1848 type = world->options->custm2[i];
1849 switch (type)
1850 {
1851 case '*':
1852 xs++;
1853 break;
1854 case 's': // M is symmetric
1855 if (i >= world->numpop)
1856 {
1857 z = i;
1858 m2mm (z, world->numpop, &ii, &jj);
1859 zz = mm2m (jj, ii, world->numpop);
1860 world->options->symparam = (twin_fmt *) myrealloc
1861 (world->options->symparam, sizeof (twin_fmt) * (ss + 2));
1862 world->options->symparam[ss][0] = zz;
1863 world->options->symparam[ss++][1] = z;
1864 world->options->symn = ss;
1865 summ = (world->param0[z] + world->param0[zz]) / 2.;
1866 world->param0[zz] = world->param0[z] = summ;
1867 }
1868 break;
1869 case 'S': // 4Nm is symmetric, not completely
1870 // implemented yet, -> derivatives.c
1871 if (i >= world->numpop)
1872 {
1873 z = i;
1874 m2mm (z, world->numpop, &ii, &jj);
1875 zz = mm2m (jj, ii, world->numpop);
1876 world->options->sym2param = (quad_fmt *)
1877 myrealloc (world->options->sym2param,
1878 sizeof (quad_fmt) * (ss2 + 2));
1879 world->options->sym2param[ss2][0] = zz;
1880 world->options->sym2param[ss2][1] = z;
1881 world->options->sym2param[ss2][2] = i;
1882 world->options->sym2param[ss2++][3] = j;
1883 world->options->sym2n = ss2;
1884 summ = (world->param0[z] * world->param0[i] +
1885 world->param0[zz] * world->param0[j]) / 2.;
1886 world->param0[z] = summ / world->param0[i];
1887 world->param0[zz] = summ / world->param0[j];
1888 }
1889 break;
1890 case 'C':
1891 case 'c':
1892 world->options->constparam = (long *) myrealloc
1893 (world->options->constparam, sizeof (long) * (ns + 2));
1894 world->options->constparam[ns++] = i;
1895 world->options->constn = ns;
1896 break;
1897 case '0':
1898 z = i;
1899 world->param0[z] = 0;
1900 world->options->zeroparam = (long *) myrealloc
1901 (world->options->zeroparam, sizeof (long) * (ns + 2));
1902 world->options->zeroparam[ns++] = i;
1903 world->options->zeron = ns;
1904 break;
1905 case 'm':
1906 summ = 0;
1907 if (i < world->numpop) /*theta */
1908 {
1909 if (!allsamesize)
1910 {
1911 allsamesize = TRUE;
1912 nsum = 0;
1913 for (z = 0; z < world->numpop; z++)
1914 {
1915 if (world->options->custm2[z] == 'm')
1916 {
1917 nsum++;
1918 summ += world->param0[z];
1919 }
1920 }
1921 summ /= nsum;
1922 for (z = 0; z < world->numpop; z++)
1923 if (world->options->custm2[z] == 'm')
1924 world->param0[z] = summ;
1925 }
1926 }
1927 else /* migration */
1928 {
1929 summ = 0;
1930 if (!partmig)
1931 {
1932 nsum = 0;
1933 partmig = TRUE;
1934 for (z = world->numpop; z < world->numpop2; z++)
1935 {
1936 if (world->options->custm2[z] == 'm')
1937 {
1938 nsum++;
1939 summ += world->param0[z];
1940 }
1941 }
1942 summ /= nsum;
1943 for (z = world->numpop; z < world->numpop2; z++)
1944 if (world->options->custm2[z] == 'm')
1945 world->param0[z] = summ;
1946 }
1947 }
1948 break;
1949 default:
1950 error ("no defaults allowed in resynchronize_param()\n");
1951 }
1952 }
1953 }
1954 //--------gamma stuff
1955 if (world->options->gamma)
1956 {
1957 if (world->options->custm2[world->numpop2] == 'c')
1958 {
1959 world->options->constparam = (long *) myrealloc
1960 (world->options->constparam, sizeof (long) * (ns + 2));
1961 world->options->constparam[ns++] = i;
1962 world->options->constn = ns;
1963 }
1964 else
1965 {
1966 world->options->custm2[world->numpop2] = '*';
1967 world->options->custm2[world->numpop2 + 1] = '\0';
1968 }
1969 }
1970 }
1971
1972
1973 long
scan_connect(char * custm2,long start,long stop,int check)1974 scan_connect (char *custm2, long start, long stop, int check)
1975 {
1976 long i, summ = 0;
1977 for (i = start; i < stop; i++)
1978 {
1979
1980 if (check == custm2[i])
1981 summ++;
1982 }
1983 return summ;
1984 }
1985
1986 void
set_partmean_mig(long ** mmparam,MYREAL * param,char * custm2,long migm,long numpop2)1987 set_partmean_mig (long **mmparam, MYREAL *param, char *custm2, long migm,
1988 long numpop2)
1989 {
1990 long i, z = 0;
1991 MYREAL summ = 0;
1992 long start = (long) sqrt ((MYREAL) numpop2);
1993 (*mmparam) = (long *) myrealloc ((*mmparam), sizeof (long) * (migm + 2));
1994
1995 for (i = start; i < numpop2; i++)
1996 {
1997 if (custm2[i] == 'm')
1998 {
1999 summ += param[i];
2000 (*mmparam)[z++] = i;
2001 }
2002 }
2003 summ /= migm;
2004 for (i = start; i < numpop2; i++)
2005 {
2006 if (custm2[i] == 'm')
2007 param[i] = summ;
2008 }
2009 }
2010
free_options_filenames(option_fmt * options)2011 void free_options_filenames(option_fmt * options)
2012 {
2013 myfree( options->parmfilename);
2014 myfree( options->infilename);
2015 myfree( options->outfilename);
2016 #ifdef PRETTY
2017 myfree( options->pdfoutfilename);
2018 #endif
2019 myfree( options->logfilename);
2020 myfree( options->mathfilename);
2021 myfree( options->sumfilename);
2022 myfree( options->treefilename);
2023 myfree( options->utreefilename);
2024 myfree( options->catfilename);
2025 myfree( options->weightfilename);
2026 myfree( options->mighistfilename);
2027 myfree( options->skylinefilename);
2028 myfree( options->distfilename);
2029 myfree( options->geofilename);
2030 myfree( options->divfilename);
2031 myfree( options->bootfilename);
2032 myfree( options->aicfilename);
2033 myfree( options->seedfilename);
2034
2035 #ifdef UEP
2036 myfree( options->uepfilename);
2037 #endif
2038 myfree( options->bayesfilename);
2039 myfree( options->bayesmdimfilename);
2040 myfree( options->datefilename);
2041 }
2042
2043 void
destroy_options(option_fmt * options)2044 destroy_options (option_fmt * options)
2045 {
2046 myfree(options->thetag);
2047 myfree(options->mg);
2048 myfree(options->mu_rates);
2049 myfree(options->inheritance_scalars);
2050 myfree(options->newpops);
2051 myfree(options->ttratio);
2052 myfree(options->rate);
2053 myfree(options->probcat);
2054 myfree(options->rrate);
2055
2056 while (--options->lratio->alloccounter >= 0)
2057 {
2058 myfree(options->lratio->data[options->lratio->alloccounter].value1);
2059 myfree(options->lratio->data[options->lratio->alloccounter].value2);
2060 myfree(options->lratio->data[options->lratio->alloccounter].connect);
2061 }
2062 myfree(options->lratio->data);
2063 myfree(options->lratio);
2064 myfree(options->custm);
2065 myfree(options->custm2);
2066 if (options->plot)
2067 {
2068 myfree(options->plotxvalues);
2069 myfree(options->plotyvalues);
2070 }
2071 myfree(options->bayespriortheta);
2072 myfree(options->mutationrate_year);
2073 free_options_filenames(options);
2074 myfree(options);
2075 }
2076
2077 void
decide_plot(worldoption_fmt * options,long chain,long chains,char type)2078 decide_plot (worldoption_fmt * options, long chain, long chains, char type)
2079 {
2080 if (options->plot && (chain >= chains - 1) && (type == 'l'))
2081 options->plotnow = TRUE;
2082 else
2083 options->plotnow = FALSE;
2084 }
2085
2086 void
set_plot(option_fmt * options)2087 set_plot (option_fmt * options)
2088 {
2089 long intervals = options->plotintervals;
2090 MYREAL prangex[2];
2091 MYREAL prangey[2];
2092 if (!options->plot)
2093 return;
2094 prangex[0] = options->plotrange[0];
2095 prangex[1] = options->plotrange[1];
2096 prangey[0] = options->plotrange[2];
2097 prangey[1] = options->plotrange[3];
2098
2099
2100 options->plotxvalues = (MYREAL *) mycalloc (1, sizeof (MYREAL) * intervals);
2101 options->plotyvalues = (MYREAL *) mycalloc (1, sizeof (MYREAL) * intervals);
2102 set_plot_values (&options->plotxvalues, prangex, options->plotintervals,
2103 options->plotscale);
2104 set_plot_values (&options->plotyvalues, prangey, options->plotintervals,
2105 options->plotscale);
2106 }
2107
2108 void
set_plot_values(MYREAL ** values,MYREAL plotrange[],long intervals,int type)2109 set_plot_values (MYREAL **values, MYREAL plotrange[], long intervals,
2110 int type)
2111 {
2112 long i;
2113 MYREAL diff = 0;
2114 MYREAL logstart = 0;
2115 (*values)[0] = plotrange[0];
2116 (*values)[intervals - 1] = plotrange[1];
2117 if (type == PLOTSCALELOG)
2118 {
2119 logstart = log10 ((*values)[0]);
2120 diff =
2121 (log10 ((*values)[intervals - 1]) - logstart) / (MYREAL) (intervals -
2122 1);
2123 for (i = 1; i < intervals - 1; i++)
2124 {
2125 (*values)[i] = pow (10., (logstart + i * diff));
2126 }
2127 }
2128 else
2129 {
2130 diff =
2131 ((*values)[intervals - 1] - (*values)[0]) / (MYREAL) (intervals - 1);
2132 for (i = 1; i < intervals - 1; i++)
2133 {
2134 (*values)[i] = (*values)[i - 1] + diff;
2135 }
2136 }
2137 }
2138
2139 ///
2140 /// save all the options into a file (default) called parmfile;
2141 /// uses save_options_buffer()
save_parmfile(option_fmt * options,world_fmt * world,data_fmt * data)2142 long save_parmfile (option_fmt * options, world_fmt * world, data_fmt *data)
2143 {
2144 FILE *fp;
2145 long bufsize;
2146 long allocbufsize = LONGLINESIZE;
2147 char *sbuffer;
2148 char *buffer;
2149 sbuffer = (char *) mycalloc (allocbufsize, sizeof (char));
2150 fp = options->parmfile;
2151 if (fp)
2152 {
2153 fclose (fp);
2154 openfile (&fp, options->parmfilename, "w", NULL);
2155 }
2156 else
2157 {
2158 openfile (&fp, options->parmfilename, "w", NULL);
2159 }
2160 options->parmfile = fp;
2161 bufsize = save_options_buffer (&sbuffer, &allocbufsize, options, data);
2162 buffer = sbuffer;
2163 fprintf (fp, "%s", buffer);
2164 fflush (fp);
2165 printf ("\n\n+++ Options were written to file %s in current directory +++\n\n", options->parmfilename);
2166 myfree(sbuffer);
2167 return bufsize;
2168 }
2169
2170
2171
2172 /// prints delimiters between sections in the parmfile
print_parm_delimiter(long * bufsize,char ** buffer,long * allocbufsize)2173 void print_parm_delimiter(long *bufsize, char **buffer, long *allocbufsize)
2174 {
2175 add_to_buffer("################################################################################\n",bufsize,buffer, allocbufsize);
2176 }
2177
2178 /// prints delimiters between sections in the parmfile
print_parm_smalldelimiter(long * bufsize,char ** buffer,long * allocbufsize)2179 void print_parm_smalldelimiter(long *bufsize, char **buffer, long *allocbufsize)
2180 {
2181 add_to_buffer("#-------------------------------------------------------------------------------\n", bufsize,buffer, allocbufsize);
2182 }
2183
2184 /// prints delimiters between sections in the parmfile
print_parm_br(long * bufsize,char ** buffer,long * allocbufsize)2185 void print_parm_br(long *bufsize, char **buffer, long *allocbufsize)
2186 {
2187 add_to_buffer("#\n", bufsize,buffer, allocbufsize);
2188 }
2189 /// prints parmfile single comment line
print_parm_comment(long * bufsize,char ** buffer,long * allocbufsize,char message[])2190 void print_parm_comment(long *bufsize, char **buffer, long *allocbufsize, char message[])
2191 {
2192 char fp[LINESIZE];
2193 sprintf(fp,"# %s\n",message);
2194 add_to_buffer(fp, bufsize,buffer, allocbufsize);
2195 }
2196
2197 /// prints parmfile mutable comment line
print_parm_mutable_comment(long * bufsize,char ** buffer,long * allocbufsize,char string[],...)2198 void print_parm_mutable_comment(long *bufsize, char **buffer, long *allocbufsize, char string[], ...)
2199 {
2200 char message[LINESIZE];
2201 char fp[LINESIZE];
2202 va_list args;
2203 va_start (args, string);
2204 vsprintf (message, string, args);
2205 va_end (args);
2206 sprintf(fp,"# %s\n",message);
2207 add_to_buffer(fp, bufsize,buffer, allocbufsize);
2208 }
2209
2210 /// prints parmfile mutable option line
print_parm_mutable(long * bufsize,char ** buffer,long * allocbufsize,char string[],...)2211 void print_parm_mutable(long *bufsize, char **buffer, long *allocbufsize, char string[], ...)
2212 {
2213 char message[LINESIZE];
2214 char fp[LINESIZE];
2215 va_list args;
2216 va_start (args, string);
2217 vsprintf (message, string, args);
2218 va_end (args);
2219 sprintf(fp,"%s\n",message);
2220 add_to_buffer(fp, bufsize,buffer, allocbufsize);
2221 }
2222
2223 /// prints parmfile fixed option line
print_parm(long * bufsize,char ** buffer,long * allocbufsize,char string[])2224 void print_parm(long *bufsize, char **buffer, long *allocbufsize, char string[])
2225 {
2226 char fp[LINESIZE];
2227 sprintf(fp,"%s\n",string);
2228 add_to_buffer(fp, bufsize,buffer, allocbufsize);
2229 }
2230 /// prints a title of section in the parmfile
print_parm_title(long * bufsize,char ** buffer,long * allocbufsize,char message[])2231 void print_parm_title(long *bufsize, char **buffer, long *allocbufsize, char message[])
2232 {
2233 print_parm_delimiter(bufsize, buffer, allocbufsize);
2234 print_parm_comment(bufsize, buffer, allocbufsize, message);
2235 print_parm_delimiter(bufsize, buffer, allocbufsize);
2236 }
2237
2238
2239
2240 /// print the ttratio into the parmfile buffer
print_parm_ttratio(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2241 void print_parm_ttratio(long *bufsize, char **buffer, long *allocbufsize, option_fmt *options)
2242 {
2243 long i=1;
2244 char fp[LINESIZE];
2245 sprintf(fp,"ttratio=%f ", options->ttratio[0]);
2246 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2247 while (options->ttratio[i] != 0.0)
2248 {
2249 sprintf (fp, "%f ", options->ttratio[i++]);
2250 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2251 }
2252 sprintf (fp, "\n");
2253 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2254 }
2255
2256 /// print base frequency into parmfile buffer
print_parm_freqfrom(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2257 void print_parm_freqfrom(long *bufsize, char **buffer, long *allocbufsize, option_fmt * options)
2258 {
2259 if (options->freqsfrom)
2260 print_parm(bufsize, buffer, allocbufsize, "freqs-from-data=YES");
2261 else
2262 print_parm_mutable(bufsize,buffer, allocbufsize, "freqs-from-data=NO:%f,%f, %f, %f\n", options->freqa,
2263 options->freqc, options->freqg, options->freqt);
2264 }
2265
2266 /// print whether one uses categories and in what file they are into parmfile buffer
print_parm_categs(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2267 void print_parm_categs(long *bufsize, char **buffer, long *allocbufsize, option_fmt * options)
2268 {
2269 if (options->categs>1)
2270 print_parm_mutable(bufsize, buffer, allocbufsize, "categories=%li:%s",options->categs, options->catfilename);
2271 else
2272 print_parm(bufsize, buffer, allocbufsize, "categories=1 #no categories file specified");
2273 }
2274
2275 /// print whether one uses weights and in what file they are into parmfile buffer
print_parm_weights(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2276 void print_parm_weights(long *bufsize, char **buffer, long *allocbufsize, option_fmt * options)
2277 {
2278 if (options->weights)
2279 print_parm_mutable(bufsize, buffer, allocbufsize, "weights=YES:%s", options->weightfilename);
2280 else
2281 print_parm(bufsize, buffer, allocbufsize, "weights=NO");
2282 }
2283
2284 /// print whether one uses rate categories, autocorrelation and in what file they are into parmfile buffer
print_parm_rates(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2285 void print_parm_rates(long *bufsize, char **buffer, long *allocbufsize, option_fmt * options)
2286 {
2287 long i;
2288 char fp[LINESIZE];
2289 long lbufsize = 0L;
2290 char *lbuffer;
2291 long alloclbufsize;
2292 // rates
2293 lbuffer = (char *) mycalloc(LINESIZE,sizeof(char));
2294 alloclbufsize = LINESIZE;
2295 for (i = 0; i < options->rcategs; i++)
2296 {
2297 sprintf (fp, "%f ", options->rrate[i]);
2298 add_to_buffer(fp,&lbufsize, &lbuffer, &alloclbufsize);
2299 }
2300 print_parm_mutable(bufsize, buffer, allocbufsize, "rates=%li: %s",options->rcategs, lbuffer);
2301 lbufsize = 0L;
2302 //reset the lbuffer
2303 lbuffer[0] = '\0';
2304 // probablities
2305 for (i = 0; i < options->rcategs; i++)
2306 {
2307 sprintf (fp, "%f ", options->probcat[i]);
2308 add_to_buffer(fp,&lbufsize,&lbuffer, &alloclbufsize);
2309 }
2310 print_parm_mutable(bufsize, buffer, allocbufsize, "prob-rates=%li: %s",options->rcategs, lbuffer);
2311 myfree(lbuffer);
2312 // autocorrelation
2313 if (!options->autocorr)
2314 print_parm(bufsize, buffer, allocbufsize, "autocorrelation=NO");
2315 else
2316 print_parm_mutable(bufsize, buffer, allocbufsize, "autocorrelation=YES:%f", 1. / options->lambda);
2317 }
2318
2319
2320 /// print data-option into parmfile buffer
print_parm_datatype(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2321 void print_parm_datatype(long *bufsize, char **buffer, long * allocbufsize, option_fmt *options)
2322 {
2323 switch (options->datatype)
2324 {
2325 case 'a':
2326 print_parm(bufsize, buffer, allocbufsize, "datatype=AllelicData");
2327 print_parm_mutable(bufsize, buffer, allocbufsize,"include-unknown=%s", options->include_unknown ? "YES" : "NO");
2328 break;
2329 case 'b':
2330 print_parm(bufsize, buffer, allocbufsize, "datatype=BrownianMicrosatelliteData");
2331 print_parm_mutable(bufsize, buffer, allocbufsize, "include-unknown=%s", options->include_unknown ? "YES" : "NO");
2332 break;
2333 case 'm':
2334 print_parm(bufsize, buffer, allocbufsize, "datatype=MicrosatelliteData");
2335 if(options->msat_option == SINGLESTEP)
2336 print_parm_mutable(bufsize, buffer, allocbufsize, "micro-submodel=%li", options->msat_option);
2337 else
2338 print_parm_mutable(bufsize, buffer, allocbufsize, "micro-submodel=%li:{%f, %f}", options->msat_option, options->msat_tuning[0], options->msat_tuning[1]);
2339 print_parm_mutable(bufsize, buffer, allocbufsize, "micro-threshold=%li", options->micro_threshold);
2340 print_parm_mutable(bufsize, buffer, allocbufsize, "include-unknown=%s", options->include_unknown ? "YES" : "NO");
2341 break;
2342 case 's':
2343 print_parm(bufsize, buffer, allocbufsize, "datatype=SequenceData");
2344 print_parm_ttratio(bufsize, buffer, allocbufsize, options);
2345 print_parm_freqfrom(bufsize, buffer, allocbufsize, options);
2346 print_parm_mutable(bufsize, buffer, allocbufsize, "seqerror-rate=%f", options->seqerror);
2347 print_parm_categs(bufsize, buffer, allocbufsize, options);
2348 print_parm_rates(bufsize, buffer, allocbufsize, options);
2349 print_parm_weights(bufsize, buffer, allocbufsize, options);
2350 print_parm_mutable(bufsize, buffer, allocbufsize, "interleaved=%s", options->interleaved ? "YES" : "NO");
2351 print_parm_mutable(bufsize, buffer, allocbufsize, "fast-likelihood=%s", options->fastlike ? "YES" : "NO");
2352 break;
2353 case 'h':
2354 print_parm(bufsize, buffer, allocbufsize, "datatype=HapmapSNPfrequencydata\n");
2355 print_parm_ttratio(bufsize, buffer, allocbufsize, options);
2356 print_parm_freqfrom(bufsize, buffer, allocbufsize, options);
2357 print_parm_mutable(bufsize, buffer, allocbufsize, "seqerror-rate=%f", options->seqerror);
2358 print_parm_categs(bufsize, buffer, allocbufsize, options);
2359 print_parm_rates(bufsize, buffer, allocbufsize, options);
2360 print_parm_weights(bufsize, buffer, allocbufsize, options);
2361 print_parm_mutable(bufsize, buffer, allocbufsize, "interleaved=%s", options->interleaved ? "YES" : "NO");
2362 print_parm_mutable(bufsize, buffer, allocbufsize, "fast-likelihood=%s", options->fastlike ? "YES" : "NO");
2363 break;
2364 case 'n':
2365 print_parm(bufsize, buffer, allocbufsize,"datatype=NucleotidePolymorphismData");
2366 print_parm_ttratio(bufsize, buffer, allocbufsize, options);
2367 print_parm_freqfrom(bufsize, buffer, allocbufsize, options);
2368 print_parm_mutable(bufsize, buffer, allocbufsize, "seqerror-rate=%f", options->seqerror);
2369 print_parm_categs(bufsize, buffer, allocbufsize, options);
2370 print_parm_rates(bufsize, buffer, allocbufsize, options);
2371 print_parm_weights(bufsize, buffer, allocbufsize, options);
2372 print_parm_mutable(bufsize, buffer, allocbufsize, "interleaved=%s", options->interleaved ? "YES" : "NO");
2373 print_parm_mutable(bufsize, buffer, allocbufsize, "fast-likelihood=%s", options->fastlike ? "YES" : "NO");
2374 break;
2375 case 'u':
2376 print_parm(bufsize, buffer, allocbufsize,"datatype=UnlinkedSNPData");
2377 print_parm_ttratio(bufsize, buffer, allocbufsize, options);
2378 print_parm_freqfrom(bufsize, buffer, allocbufsize, options);
2379 print_parm_mutable(bufsize, buffer, allocbufsize, "seqerror-rate=%f", options->seqerror);
2380 print_parm_categs(bufsize, buffer, allocbufsize, options);
2381 print_parm_rates(bufsize, buffer, allocbufsize, options);
2382 print_parm_weights(bufsize, buffer, allocbufsize, options);
2383 print_parm_mutable(bufsize, buffer, allocbufsize, "interleaved=%s", options->interleaved ? "YES" : "NO");
2384 print_parm_mutable(bufsize, buffer, allocbufsize, "fast-likelihood=%s", options->fastlike ? "YES" : "NO");
2385 break;
2386 case 'f':
2387 print_parm(bufsize, buffer, allocbufsize,"datatype=F-Ancestral state method");
2388 print_parm_ttratio(bufsize, buffer, allocbufsize, options);
2389 print_parm_freqfrom(bufsize, buffer, allocbufsize, options);
2390 print_parm_mutable(bufsize, buffer, allocbufsize, "seqerror-rate=%f", options->seqerror);
2391 print_parm_categs(bufsize, buffer, allocbufsize, options);
2392 print_parm_rates(bufsize, buffer, allocbufsize, options);
2393 print_parm_weights(bufsize, buffer, allocbufsize, options);
2394 print_parm_mutable(bufsize, buffer, allocbufsize, "interleaved=%s", options->interleaved ? "YES" : "NO");
2395 print_parm_mutable(bufsize, buffer, allocbufsize, "fast-likelihood=%s", options->fastlike ? "YES" : "NO");
2396 break;
2397 case 'g':
2398 print_parm(bufsize, buffer, allocbufsize,"datatype=GenealogySummaryOlderRun");
2399 break;
2400 default:
2401 error ("the parmfile-writer contains an error");
2402 }
2403 }
2404
print_parm_tipdate(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options,data_fmt * data)2405 void print_parm_tipdate(long *bufsize, char **buffer, long *allocbufsize, option_fmt *options, data_fmt *data)
2406 {
2407 long pos;
2408 long locus;
2409 char *input;
2410 if(options->has_datefile)
2411 {
2412 print_parm_mutable(bufsize, buffer, allocbufsize, "tipdate-file=YES:%s", options->datefilename);
2413 print_parm_mutable(bufsize, buffer, allocbufsize, "generation-per-year=%f", options->generation_year);
2414 input = (char *) mycalloc(options->mutationrate_year_numalloc * 50, sizeof(char));
2415 pos = sprintf(input, "{%20.20f", options->mutationrate_year[0]);
2416 for(locus=1; locus < options->mutationrate_year_numalloc; locus++)
2417 {
2418 pos += sprintf(input + pos,", %20.20f", options->mutationrate_year[locus]);
2419 }
2420 //xcode pos +=
2421 sprintf(input + pos,"}");
2422 print_parm_mutable(bufsize, buffer, allocbufsize, "mutationrate-per-year=%s", input);
2423 myfree(input);
2424 }
2425 }
2426
2427 ///
2428 /// print the parmfile entry for the inheritance scalar
print_parm_inheritence(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options,data_fmt * data)2429 void print_parm_inheritence(long *bufsize, char **buffer, long *allocbufsize, option_fmt *options, data_fmt *data)
2430 {
2431 long pos;
2432 long locus;
2433 char *input;
2434 input = (char *) mycalloc(options->inheritance_scalars_numalloc * 50, sizeof(char));
2435 pos = sprintf(input, "inheritance-scalars={%20.20f", options->inheritance_scalars[0]);
2436 for(locus=1; locus < options->inheritance_scalars_numalloc; locus++)
2437 {
2438 pos += sprintf(input + pos,", %20.20f", options->inheritance_scalars[locus]);
2439 }
2440 //xcode pos +=
2441 sprintf(input + pos,"}");
2442 print_parm_mutable(bufsize, buffer, allocbufsize, "%s", input);
2443 myfree(input);
2444 }
2445 ///
2446 /// print the parmfile entry for the population relabeling
print_parm_newpops(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options,data_fmt * data)2447 void print_parm_newpops(long *bufsize, char **buffer, long *allocbufsize, option_fmt *options, data_fmt *data)
2448 {
2449 long pos;
2450 long i;
2451 char *input;
2452 input = (char *) mycalloc(options->newpops_numalloc * 50, sizeof(char));
2453 pos = sprintf(input, "population-relabel={%li", options->newpops[0]);
2454 for(i=1; i < options->newpops_numalloc; i++)
2455 {
2456 pos += sprintf(input + pos,", %li", options->newpops[i]);
2457 }
2458 //xcode pos +=
2459 sprintf(input + pos,"}");
2460 print_parm_mutable(bufsize, buffer, allocbufsize, "%s", input);
2461 myfree(input);
2462 }
2463
print_parm_randomsubset(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2464 void print_parm_randomsubset(long * bufsize, char **buffer, long *allocbufsize, option_fmt *options)
2465 {
2466 if(options->randomsubset>0)
2467 {
2468 if(options->randomsubsetseed>0)
2469 print_parm_mutable(bufsize, buffer, allocbufsize, "random-subset=%li:%li", options->randomsubset,options->randomsubsetseed);
2470 else
2471 print_parm_mutable(bufsize, buffer, allocbufsize, "random-subset=%li", options->randomsubset);
2472 return;
2473 }
2474 }
2475
print_parm_usertree(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2476 void print_parm_usertree(long * bufsize, char **buffer, long *allocbufsize, option_fmt *options)
2477 {
2478 if(options->usertree)
2479 {
2480 print_parm_mutable(bufsize, buffer, allocbufsize, "usertree=TREE:%s", options->utreefilename);
2481 return;
2482 }
2483 if (options->randomtree)
2484 {
2485 print_parm(bufsize, buffer, allocbufsize, "usertree=RANDOMTREE");
2486 return;
2487 }
2488 if(options->dist)
2489 {
2490 print_parm_mutable(bufsize, buffer, allocbufsize, "usertree=DISTANCE:%s", options->distfilename);
2491 return;
2492 }
2493 print_parm(bufsize, buffer, allocbufsize, "usertree=AUTOMATIC");
2494 }
2495
2496 /// print the theta starting parameters
print_parm_theta(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2497 void print_parm_theta(long *bufsize, char ** buffer, long *allocbufsize, option_fmt * options)
2498 {
2499 long i;
2500 char fp[LINESIZE];
2501 char *lbuffer;
2502 long lbufsize = 1;
2503 long alloclbufsize = LINESIZE;
2504 lbuffer = (char *) mycalloc(alloclbufsize,sizeof(char));
2505
2506 switch (options->numthetag)
2507 {
2508 case 0:
2509 if (strchr ("snupf", options->datatype))
2510 {
2511 sprintf (fp, "theta=Own:0.01");
2512 add_to_buffer(fp, &lbufsize, &lbuffer, &alloclbufsize);
2513 }
2514 else
2515 {
2516 sprintf (fp, "theta=Own:1.0");
2517 add_to_buffer(fp, &lbufsize, &lbuffer, &alloclbufsize);
2518 }
2519 break;
2520 case 1:
2521 sprintf (fp, "theta=Own:%f",options->thetag[0]);
2522 add_to_buffer(fp, &lbufsize, &lbuffer, &alloclbufsize);
2523 break;
2524 default:
2525 sprintf (fp, "theta=Own:{");
2526 add_to_buffer(fp,&lbufsize, &lbuffer, &alloclbufsize);
2527 for (i = 0; i < options->numthetag - 1; i++)
2528 {
2529 sprintf (fp, "%f ", options->thetag[i]);
2530 add_to_buffer(fp,&lbufsize,&lbuffer, &alloclbufsize);
2531 }
2532 sprintf (fp, "%f}", options->thetag[i]);
2533 add_to_buffer(fp,&lbufsize, &lbuffer, &alloclbufsize);
2534 break;
2535 }
2536 print_parm_mutable(bufsize, buffer, allocbufsize, "%s", lbuffer);
2537 }
2538
2539 /// print M starting parameters to the parmfile buffer
print_parm_m(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2540 void print_parm_m(long *bufsize, char **buffer, long *allocbufsize, option_fmt * options)
2541 {
2542 char fp[LINESIZE];
2543 long i, j, z, num;
2544 switch (options->nummg)
2545 {
2546 case 0:
2547 sprintf (fp, "migration=Own:1\n");
2548 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2549 break;
2550 case 1:
2551 sprintf (fp, "migration=Own:%f\n", options->mg[0]);
2552 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2553 break;
2554 default:
2555 sprintf (fp, "migration=Own:{ ");
2556 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2557 z = 0;
2558 num = (long) (1. + sqrt (4. * (MYREAL) options->nummg + 1.) / 2.);
2559 for (i = 0; i < num; i++)
2560 {
2561 for (j = 0; j < num; j++)
2562 {
2563 if (i == j)
2564 {
2565 sprintf (fp, "- ");
2566 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2567 }
2568 else
2569 {
2570 sprintf (fp, "%f ", options->mg[z++]);
2571 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2572 }
2573 }
2574 }
2575 sprintf (fp, "}\n");
2576 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2577 }
2578 }
2579
print_parm_heating(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2580 void print_parm_heating(long *bufsize, char **buffer, long *allocbufsize, option_fmt *options)
2581 {
2582 long i;
2583 char fp[LINESIZE];
2584 sprintf (fp, "heating=%s", options->heating ? (options->adaptiveheat!=NOTADAPTIVE ? (options->adaptiveheat==STANDARD ? "ADAPTIVE_standard" : "Bounded_adaptive") : "YES") : "NO\n");
2585 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2586
2587 if (options->heating)
2588 {
2589 sprintf (fp, ":%li:{%f,", options->heating_interval, options->heat[0]);
2590 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2591
2592 for (i = 1; i < options->heated_chains - 1; i++)
2593 {
2594 sprintf (fp, "%f,", options->heat[i]);
2595 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2596 }
2597 sprintf (fp, "%f}\nheated-swap=%s\n", options->heat[i], (options->heatedswap_off ? "NO" : "YES") );
2598 add_to_buffer(fp,bufsize,buffer, allocbufsize);
2599 }
2600 }
2601
2602 ///
2603 /// returns a TRUE when the option is set and sets a filename, returns FALSE when the
2604 /// option is not used and then of course does not set the filename
set_filename(char * value,char comparison[],char ** filename)2605 boolean set_filename(char *value, char comparison[], char ** filename)
2606 {
2607 unsigned long len=1;
2608 char *newvalue;
2609 char *temp;
2610 temp = (char *) mycalloc(LINESIZE,sizeof(char));
2611 upper(value, &temp);
2612 len = strlen(comparison);
2613 if (strncmp (temp, comparison, len)==0)
2614 {
2615 newvalue = strchr(value,':');
2616 if(newvalue!=NULL)
2617 get_filename (filename, newvalue+1);
2618 myfree(temp);
2619 return TRUE;
2620 }
2621 myfree(temp);
2622 return FALSE;
2623 }
2624
2625 /// assumes that an options will be used and takes the filename after the ":"
set_filename_only(boolean check,char * value,char ** filename)2626 void set_filename_only(boolean check, char *value, char ** filename)
2627 {
2628 char *newvalue;
2629 if(check)
2630 {
2631 newvalue = strchr(value,':');
2632 if(newvalue != NULL)
2633 get_filename(filename, newvalue + 1);
2634 }
2635 }
2636
2637
set_parm_prior_values(int priortype,prior_fmt * prior,char * mytext)2638 void set_parm_prior_values(int priortype, prior_fmt * prior, char * mytext)
2639 {
2640 char tmp1[LINESIZE];
2641 char tmp2[LINESIZE];
2642 char tmp3[LINESIZE];
2643 char tmp4[LINESIZE];
2644 mytext[0]='\0';
2645 show_priormin(tmp1, prior, priortype);
2646 strcat(mytext,tmp1);
2647 switch(priortype)
2648 {
2649 // case SLICE:
2650 //show_priormax(tmp3, prior, priortype);
2651 //strcat(mytext,tmp3);
2652 //break;
2653 // case MULTPRIOR:
2654 // show_priormax(tmp3, prior, priortype);
2655 // show_priordelta(tmp4, prior,priortype);
2656 // strcat(mytext,tmp3);
2657 // strcat(mytext,tmp4);
2658 // break;
2659 case EXPPRIOR:
2660 show_priormean(tmp2,prior, priortype);
2661 show_priormax(tmp3, prior,priortype);
2662 strcat(mytext,tmp2);
2663 strcat(mytext,tmp3);
2664 break;
2665 case WEXPPRIOR:
2666 show_priormean(tmp2, prior, priortype);
2667 show_priormax(tmp3, prior, priortype);
2668 show_priordelta(tmp4, prior, priortype);
2669 strcat(mytext,tmp2);
2670 strcat(mytext,tmp3);
2671 strcat(mytext,tmp4);
2672 break;
2673 case GAMMAPRIOR:
2674 show_priormean(tmp2, prior, priortype);
2675 show_priormax(tmp3, prior, priortype);
2676 show_prioralpha(tmp4, prior, priortype);
2677 strcat(mytext,tmp2);
2678 strcat(mytext,tmp3);
2679 strcat(mytext,tmp4);
2680 break;
2681 case UNIFORMPRIOR:
2682 default:
2683 show_priormax(tmp3, prior, priortype);
2684 show_priordelta(tmp4, prior, priortype);
2685 strcat(mytext,tmp3);
2686 strcat(mytext,tmp4);
2687 break;
2688 }
2689 }
2690
2691
2692 /// \brief returns proposal type string
2693 /// returns a string that shows what proposal type is set
show_proposaltype(boolean priorset)2694 char * show_proposaltype(boolean priorset)
2695 {
2696 if(priorset)
2697 {
2698 return "SLICE Sampler" ;
2699 }
2700 else
2701 {
2702 return "METROPOLIS-HASTINGS Sampler";
2703 }
2704 }
2705
2706 /// \brief returns priortype sting
2707 /// returns a string that shows what prior distribution is set
show_parmpriortype(int priorset)2708 char * show_parmpriortype(int priorset)
2709 {
2710 switch(priorset)
2711 {
2712 // case MULTPRIOR: return "MULTPRIOR" ;
2713 case EXPPRIOR: return "EXPPRIOR" ;
2714 case WEXPPRIOR: return "WEXPPRIOR";
2715 case GAMMAPRIOR: return "GAMMAPRIOR";
2716 case UNIFORMPRIOR: return "UNIFORMPRIOR";
2717 default: return "UNIFORMPRIOR";
2718 }
2719 }
2720
2721 ///
2722 /// print prior values to buffer for parmfile
print_parm_proposal(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2723 void print_parm_proposal(long *bufsize, char **buffer, long *allocbufsize, option_fmt *options)
2724 {
2725 print_parm_mutable(bufsize, buffer, allocbufsize, "bayes-proposals= THETA %s",
2726 show_proposaltype(options->slice_sampling[THETAPRIOR]));
2727
2728 print_parm_mutable(bufsize, buffer, allocbufsize, "bayes-proposals= MIG %s",
2729 show_proposaltype(options->slice_sampling[MIGPRIOR]));
2730 if(options->bayesmurates)
2731 {
2732 print_parm_mutable(bufsize, buffer, allocbufsize, "bayes-proposals= RATE %s",
2733 show_proposaltype(options->slice_sampling[RATEPRIOR]));
2734 }
2735 }
2736
2737 ///
2738 /// print prior values to buffer for parmfile
print_parm_prior(long * bufsize,char ** buffer,long * allocbufsize,option_fmt * options)2739 void print_parm_prior(long *bufsize, char **buffer, long *allocbufsize, option_fmt *options)
2740 {
2741 char mytext[LINESIZE]="";
2742
2743 set_parm_prior_values(options->bayesprior[THETAPRIOR], options->bayespriortheta, mytext);
2744 print_parm_mutable(bufsize, buffer, allocbufsize, "bayes-priors= THETA %s: %s",
2745 show_parmpriortype(options->bayesprior[THETAPRIOR]), mytext);
2746
2747 set_parm_prior_values(options->bayesprior[MIGPRIOR], options->bayespriorm, mytext);
2748 print_parm_mutable(bufsize, buffer, allocbufsize, "bayes-priors= MIG %s: %s",
2749 show_parmpriortype(options->bayesprior[MIGPRIOR]), mytext);
2750 if(options->bayesmurates)
2751 {
2752 set_parm_prior_values(options->bayesprior[RATEPRIOR], options->bayespriorrate, mytext);
2753 print_parm_mutable(bufsize, buffer, allocbufsize, "bayes-priors= RATE %s: %s",
2754 show_parmpriortype(options->bayesprior[RATEPRIOR]), mytext);
2755 }
2756 }
2757
2758 ///copy mu_rates into world->options
set_meanmu(worldoption_fmt * wopt,option_fmt * options,long loci)2759 void set_meanmu(worldoption_fmt * wopt, option_fmt * options, long loci)
2760 {
2761 long i;
2762 long n = 0;
2763 MYREAL last;
2764
2765 if(options->mutationrate_year_numalloc == 0)
2766 {
2767 last = 1.0;
2768 }
2769 else
2770 {
2771 last = 0.0;
2772 }
2773
2774 wopt->meanmu = (MYREAL *) mycalloc(loci, sizeof(MYREAL));
2775
2776 for(i=0;i< options->mutationrate_year_numalloc;i++)
2777 {
2778 n++;
2779 wopt->meanmu[i] = options->mutationrate_year[i];
2780 last += (wopt->meanmu[i] - last ) / n;
2781 }
2782
2783 for(i = options->mutationrate_year_numalloc; i < loci;i++)
2784 {
2785 wopt->meanmu[i] = last;
2786 }
2787 }
2788
2789 ///save option->mu_rates into a buffer
save_mu_rates_buffer(char ** buffer,long * allocbufsize,option_fmt * options)2790 long save_mu_rates_buffer (char **buffer, long *allocbufsize, option_fmt * options)
2791 {
2792 long i;
2793 long bufsize = 0;
2794 print_parm_mutable(&bufsize, buffer, allocbufsize, "%li ",options->muloci);
2795 for(i=0; i < options->muloci; i++)
2796 {
2797 print_parm_mutable(&bufsize, buffer, allocbufsize, "%f ",options->mu_rates[i]);
2798 }
2799 // fprintf(stdout,"muloci=%li, muratebuffer:>%s<\n",options->muloci, *buffer, allocbufsize);
2800 return bufsize;
2801 }
2802
2803 ///save option->inheritance_scalars into a buffer
save_inheritance_scalars_buffer(char ** buffer,long * allocbufsize,option_fmt * options)2804 long save_inheritance_scalars_buffer (char **buffer, long *allocbufsize, option_fmt * options)
2805 {
2806 long i;
2807 long bufsize = 0;
2808 print_parm_mutable(&bufsize, buffer, allocbufsize, "%li ",options->inheritance_scalars_numalloc);
2809 for(i=0; i < options->inheritance_scalars_numalloc; i++)
2810 {
2811 print_parm_mutable(&bufsize, buffer, allocbufsize, "%f ",options->inheritance_scalars[i]);
2812 }
2813 return bufsize;
2814 }
2815 ///save option->newpops into a buffer
save_newpops_buffer(char ** buffer,long * allocbufsize,option_fmt * options)2816 long save_newpops_buffer (char **buffer, long *allocbufsize, option_fmt * options)
2817 {
2818 long i;
2819 long bufsize = 0;
2820 print_parm_mutable(&bufsize, buffer, allocbufsize, "%li ",options->newpops_numalloc);
2821 for(i=0; i < options->newpops_numalloc; i++)
2822 {
2823 print_parm_mutable(&bufsize, buffer, allocbufsize, "%f ",options->newpops[i]);
2824 }
2825 return bufsize;
2826 }
2827
2828 ///save options into a buffer
save_options_buffer(char ** buffer,long * allocbufsize,option_fmt * options,data_fmt * data)2829 long save_options_buffer (char **buffer, long *allocbufsize, option_fmt * options, data_fmt *data)
2830 {
2831 long i;
2832 char nowstr[LINESIZE] = "----";
2833 char fp[LINESIZE];
2834 long bufsize = 0;
2835 get_time (nowstr, "%c");
2836 if(options->bayes_infer)
2837 options->schains=0;
2838 // header for parmfile
2839 print_parm_delimiter(&bufsize, buffer, allocbufsize);
2840 print_parm_mutable_comment(&bufsize, buffer, allocbufsize, "Parmfile for Migrate %s-%s [do not remove these first TWO lines]", MIGRATEVERSION,MIGRATESUBVERSION);
2841 print_parm_comment(&bufsize, buffer, allocbufsize, "generated automatically on");
2842 print_parm_mutable_comment(&bufsize, buffer, allocbufsize, "%s", nowstr);
2843 print_parm_br(&bufsize, buffer, allocbufsize);
2844 print_parm_comment(&bufsize, buffer, allocbufsize, "please report problems to Peter Beerli");
2845 print_parm_comment(&bufsize, buffer, allocbufsize, " email: beerli@fsu.edu");
2846 print_parm_comment(&bufsize, buffer, allocbufsize," http://popgen.sc.fsu.edu/migrate.html");
2847 print_parm_delimiter(&bufsize, buffer, allocbufsize);
2848 // general options
2849 print_parm_br(&bufsize, buffer, allocbufsize);
2850 print_parm_title(&bufsize, buffer, allocbufsize, "General options");
2851 print_parm_br(&bufsize, buffer, allocbufsize);
2852 print_parm_comment(&bufsize, buffer, allocbufsize,"Interactive or batch job usage");
2853 print_parm_comment(&bufsize, buffer, allocbufsize," Syntax: menu= < YES | NO > ");
2854 print_parm_comment(&bufsize, buffer, allocbufsize,"For batch runs it needs to be set to NO");
2855 print_parm_mutable(&bufsize, buffer, allocbufsize,"menu=%s", options->menu ? "YES" : "NO ");
2856 print_parm_br(&bufsize, buffer, allocbufsize);
2857 print_parm_comment(&bufsize, buffer, allocbufsize,"Specification of length of names of indiviudals");
2858 print_parm_comment(&bufsize, buffer, allocbufsize," Syntax: nmlength=<INTEGER between 0 .. 30>");
2859 print_parm_mutable(&bufsize, buffer, allocbufsize,"nmlength=%li", options->nmlength);
2860 print_parm_br(&bufsize, buffer, allocbufsize);
2861
2862 print_parm_br(&bufsize, buffer, allocbufsize);
2863 print_parm_title(&bufsize, buffer, allocbufsize, "Data options");
2864 print_parm_br(&bufsize, buffer, allocbufsize);
2865 print_parm_comment(&bufsize, buffer, allocbufsize,"Several different main datatypes are possible:");
2866 print_parm_comment(&bufsize, buffer, allocbufsize,"INFINITE ALLELE: usable for electrophoretic markers,");
2867 print_parm_comment(&bufsize, buffer, allocbufsize," other markers with unknown mutation model");
2868 print_parm_comment(&bufsize, buffer, allocbufsize,"STEPWISE MUTATION: usable for microsatellite data or");
2869 print_parm_comment(&bufsize, buffer, allocbufsize," other markers with stepwise change");
2870 print_parm_comment(&bufsize, buffer, allocbufsize," from one allele to another");
2871 print_parm_comment(&bufsize, buffer, allocbufsize," [singlestep versus multistep model, see micro-submodel option]");
2872 print_parm_comment(&bufsize, buffer, allocbufsize,"FINITE SITES MUTATION: standard DNA/RNA sequence mutation");
2873 print_parm_comment(&bufsize, buffer, allocbufsize," model, usable for DNA or RNA contiguous");
2874 print_parm_comment(&bufsize, buffer, allocbufsize," sequences or varialbe sites only (SNP)");
2875 print_parm_comment(&bufsize, buffer, allocbufsize,"GENEALOGY SUMMARY: reanalyzing an old migrate run");
2876 print_parm_br(&bufsize, buffer, allocbufsize);
2877 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
2878 print_parm_comment(&bufsize, buffer, allocbufsize,"INFINITE ALLELE");
2879 print_parm_comment(&bufsize, buffer, allocbufsize," Syntax: datatype=ALLELICDATA ");
2880 print_parm_comment(&bufsize, buffer, allocbufsize," include-unknown=<YES | NO> with YES unknown alleles");
2881 print_parm_comment(&bufsize, buffer, allocbufsize," are included into analysis, NO is the default");
2882 print_parm_br(&bufsize, buffer, allocbufsize);
2883 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
2884 print_parm_br(&bufsize, buffer, allocbufsize);
2885 print_parm_comment(&bufsize, buffer, allocbufsize,"STEPWISE MUTATION");
2886 print_parm_comment(&bufsize, buffer, allocbufsize," Syntax: datatype=<MICROSATELLITEDATA | BROWNIANDATA");
2887 print_parm_comment(&bufsize, buffer, allocbufsize," MICRO specifies the standard stepwise mutation");
2888 print_parm_comment(&bufsize, buffer, allocbufsize," model, the BROWNIAN is an approximation to this");
2889 print_parm_comment(&bufsize, buffer, allocbufsize," micro-submodel=<1|2:{tune,pinc}>");
2890 print_parm_comment(&bufsize, buffer, allocbufsize," 1 means singlestep mutation model (this is the default and the standard");
2891 print_parm_comment(&bufsize, buffer, allocbufsize," 2 is the Multistep model (see Watkins 2007 TPB, section 4.2) it needs");
2892 print_parm_comment(&bufsize, buffer, allocbufsize," two parameters: tune specifies how close the model is to a singlestep model");
2893 print_parm_comment(&bufsize, buffer, allocbufsize," so tune=0 --> singlestep, tune=1 --> infinite allele model;");
2894 print_parm_comment(&bufsize, buffer, allocbufsize," the second parameter defines the probability that the repeat number");
2895 print_parm_comment(&bufsize, buffer, allocbufsize," is increasing, this value cannot be larger than 0.666, I suggest 0.5.");
2896 print_parm_comment(&bufsize, buffer, allocbufsize," Example: micro-submodel=2:{0.5,0.5}");
2897 print_parm_comment(&bufsize, buffer, allocbufsize," micro-threshold=<INTEGER> Default is 10 [MICRO only, NEEDS TO BE EVEN!],");
2898 print_parm_comment(&bufsize, buffer, allocbufsize," smaller values speed up analysis, but might also");
2899 print_parm_comment(&bufsize, buffer, allocbufsize," crash, large values slow down analysis considerably.");
2900 print_parm_comment(&bufsize, buffer, allocbufsize," Change this value only when you suspect that your");
2901 print_parm_comment(&bufsize, buffer, allocbufsize," data has huge gaps in repeat length.");
2902 print_parm_comment(&bufsize, buffer, allocbufsize," include-unknown=<YES | NO> with YES unknown alleles");
2903 print_parm_comment(&bufsize, buffer, allocbufsize," are included into analysis, NO is the default");
2904 print_parm_br(&bufsize, buffer, allocbufsize);
2905 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
2906 print_parm_br(&bufsize, buffer, allocbufsize);
2907 print_parm_comment(&bufsize, buffer, allocbufsize,"FINITE SITES MUTATION");
2908 print_parm_comment(&bufsize, buffer, allocbufsize," Syntax: datatype=<SEQUENCEDATA | NUCLEOTIDE | UNLINKEDSNPS | ANCESTRAL");
2909 print_parm_comment(&bufsize, buffer, allocbufsize," SEQENCEDATA: typical linked stretches of DNA, for example mtDNA");
2910 print_parm_comment(&bufsize, buffer, allocbufsize," NUCLEOTIDE: linked DNA stretches, all invariable sites removed");
2911 print_parm_comment(&bufsize, buffer, allocbufsize," UNLINKEDSNPS: each variable site is a locus, DO NOT USE THIS YET");
2912 print_parm_comment(&bufsize, buffer, allocbufsize," ANCESTRAL: instead taking into account all posible states, use");
2913 print_parm_comment(&bufsize, buffer, allocbufsize," use only the most likely state probability, DON'T USE THIS YET");
2914 print_parm_br(&bufsize, buffer, allocbufsize);
2915 print_parm_comment(&bufsize, buffer, allocbufsize," freqs-from-data=<YES | NO: freq(A), freq(C), freq(G), freq(T)>");
2916 print_parm_comment(&bufsize, buffer, allocbufsize, " calculate the prior base frequencies from the data,");
2917 print_parm_comment(&bufsize, buffer, allocbufsize, " or specify the frequencies");
2918 print_parm_comment(&bufsize, buffer, allocbufsize, " ttratio=<RATIO1 RATIO2 ....> Default is 2.0,");
2919 print_parm_comment(&bufsize, buffer, allocbufsize, " ratio between transitions and transversions.");
2920 print_parm_comment(&bufsize, buffer, allocbufsize, " seq-error=<VALUE> Default is 0.0, typical values for ABI 3700 ");
2921 print_parm_comment(&bufsize, buffer, allocbufsize, " sequencers after base calling are around 0.001 (1/650)");
2922 print_parm_comment(&bufsize, buffer, allocbufsize, " categories=<VALUE:CATFILE> The categories are integers or letters");
2923 print_parm_comment(&bufsize, buffer, allocbufsize, " specified in file called CATFILE, this assumes that all");
2924 print_parm_comment(&bufsize, buffer, allocbufsize, " sites belong to known categories, this can be used to");
2925 print_parm_comment(&bufsize, buffer, allocbufsize, " weight third positions etc.");
2926 print_parm_comment(&bufsize, buffer, allocbufsize, " rates=<VALUE1 VALUE2 ...> the rates are specified arbitrarily or");
2927 print_parm_comment(&bufsize, buffer, allocbufsize, " then are from a Gamma distribution with alpha=x, currently");
2928 print_parm_comment(&bufsize, buffer, allocbufsize, " the alpha value gets lost and is not recorded in the parmfile");
2929 print_parm_comment(&bufsize, buffer, allocbufsize, " prob-rates=<RATE2 RATE1 ... > These rates can be arbitrary or ");
2930 print_parm_comment(&bufsize, buffer, allocbufsize, " generated with gamma-deviated rates and then are derived");
2931 print_parm_comment(&bufsize, buffer, allocbufsize, " using Laguerre's quadrature, this should get better");
2932 print_parm_comment(&bufsize, buffer, allocbufsize, " results than equal probability methods.");
2933 print_parm_comment(&bufsize, buffer, allocbufsize, " autocorrelation=<NO | YES:VALUE> Default is NO");
2934 print_parm_comment(&bufsize, buffer, allocbufsize, " autocorrelation makes only sense with rates,");
2935 print_parm_comment(&bufsize, buffer, allocbufsize, " VALUE should be >1.0");
2936 print_parm_comment(&bufsize, buffer, allocbufsize, " weights=<NO | YES:WEIGHTFILE> The weights are specified");
2937 print_parm_comment(&bufsize, buffer, allocbufsize, " in file called WEIGHTFILE, this assumes that all sites");
2938 print_parm_comment(&bufsize, buffer, allocbufsize, " belong to known weights, this can be used to weight");
2939 print_parm_comment(&bufsize, buffer, allocbufsize, " portions of the sequence etc.");
2940 print_parm_comment(&bufsize, buffer, allocbufsize, " interleaved=<YES | NO> Use either an interleaved or ");
2941 print_parm_comment(&bufsize, buffer, allocbufsize, " non-interleaved format. Default is NO,");
2942 print_parm_comment(&bufsize, buffer, allocbufsize, " interleaved=YES is discouraged");
2943 print_parm_comment(&bufsize, buffer, allocbufsize, " fast-likelihood=<YES | NO> Default is YES, use NO when you");
2944 print_parm_comment(&bufsize, buffer, allocbufsize, " have many hundred individuals and get strange errors");
2945 print_parm_comment(&bufsize, buffer, allocbufsize, " during a run, NO is scaling the conditional likelihood");
2946 print_parm_comment(&bufsize, buffer, allocbufsize, " so that very small values are >0.00000");
2947 print_parm_comment(&bufsize, buffer, allocbufsize, " inheritance-scalars={values for each locus}");
2948 print_parm_comment(&bufsize, buffer, allocbufsize, " these values are multiplied with Theta, for example having");
2949 print_parm_comment(&bufsize, buffer, allocbufsize, " two autosomal and a locus on X- and one on Y-chromosome we would give ");
2950 print_parm_comment(&bufsize, buffer, allocbufsize, " inheritance-scalars={1 1 0.75 0.25}");
2951 print_parm_comment(&bufsize, buffer, allocbufsize, " [if all loci have the same scalar, just use {1}, even for many loci]]");
2952 print_parm_comment(&bufsize, buffer, allocbufsize, " population-relabel={assignment for each location in the infile}");
2953 print_parm_comment(&bufsize, buffer, allocbufsize, " example is population-relabel={1 2 2}");
2954 print_parm_comment(&bufsize, buffer, allocbufsize, " random-subset=number<:seed>");
2955 print_parm_comment(&bufsize, buffer, allocbufsize, " allows to subset the dataset randomly, if number > sample in population");
2956 print_parm_comment(&bufsize, buffer, allocbufsize, " all samples are taken, if number is smaller then the pop sample is shuffled and");
2957 print_parm_comment(&bufsize, buffer, allocbufsize, " and the first number samples are taken.");
2958 print_parm_comment(&bufsize, buffer, allocbufsize, " the random number seed guarantees that the");
2959 print_parm_comment(&bufsize, buffer, allocbufsize, " same subset is chosen in different runs");
2960 print_parm_comment(&bufsize, buffer, allocbufsize, " usertree=<NO | UPGMA | AUTOMATIC | TREE:TREEFILE | DISTANCE:DISTFILE | RANDOM>");
2961 print_parm_comment(&bufsize, buffer, allocbufsize, " Default is RANDOM, NO delivers a UPGMA tree using the data");
2962 print_parm_comment(&bufsize, buffer, allocbufsize, " with TREE and DISTANCE the user needs to ");
2963 print_parm_comment(&bufsize, buffer, allocbufsize, " give a usertreefile or a pairwise distance file, with RANDOM");
2964 print_parm_comment(&bufsize, buffer, allocbufsize, " a random tree will be the starting tree");
2965 print_parm_br(&bufsize, buffer, allocbufsize);
2966 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
2967 print_parm_br(&bufsize, buffer, allocbufsize);
2968 print_parm_br(&bufsize, buffer, allocbufsize);
2969
2970
2971 print_parm_datatype(&bufsize, buffer, allocbufsize, options);
2972 print_parm_tipdate(&bufsize, buffer, allocbufsize, options, data);
2973 print_parm_inheritence(&bufsize, buffer, allocbufsize, options, data);
2974 print_parm_newpops(&bufsize, buffer, allocbufsize, options, data);
2975 print_parm_randomsubset(&bufsize, buffer, allocbufsize, options);
2976 print_parm_usertree(&bufsize, buffer, allocbufsize, options);
2977
2978 #ifdef UEP
2979 // unique event polymorphisms
2980 print_parm_br(&bufsize, buffer, allocbufsize);
2981 print_parm_title(&bufsize, buffer, allocbufsize, "Unique event polymorphism options");
2982 print_parm_br(&bufsize, buffer, allocbufsize);
2983 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax uep=<NO | YES:UEPFILE >");
2984 print_parm_comment(&bufsize, buffer, allocbufsize, " Default is NO, with YES the user needs to ");
2985 print_parm_comment(&bufsize, buffer, allocbufsize, " give a file for each individual (same order as in datafile");
2986 print_parm_comment(&bufsize, buffer, allocbufsize, " the value: indiviudal<10 characters> uep-state<0|1>");
2987 print_parm_comment(&bufsize, buffer, allocbufsize, " uep-rates= mu[0->1] : nu[0->1] mutation rate between the two UEP states");
2988 print_parm_comment(&bufsize, buffer, allocbufsize, " uep-bases= FREQ1 : FREQ2 prior frequency of the two alleles");
2989 print_parm_br(&bufsize, buffer, allocbufsize);
2990
2991 if (options->uep)
2992 {
2993 print_parm_mutable(&bufsize, buffer, allocbufsize, "uep=YES:%s", options->uepfilename);
2994 print_parm_mutable(&bufsize, buffer, allocbufsize, "uep-rates=%f:%f", options->uepmu, options->uepnu);
2995 print_parm_mutable(&bufsize, buffer, allocbufsize, "uep-bases=%f:%f", options->uepfreq0, options->uepfreq1);
2996 }
2997 else
2998 {
2999 print_parm(&bufsize, buffer, allocbufsize, "uep=NO");
3000 }
3001 #endif
3002 //input and output options
3003 print_parm_br(&bufsize, buffer, allocbufsize);
3004 print_parm_title(&bufsize, buffer, allocbufsize, "Input options");
3005 print_parm_br(&bufsize, buffer, allocbufsize);
3006
3007 print_parm_comment(&bufsize, buffer, allocbufsize, "input file location");
3008 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax infile=FILEPATH");
3009 print_parm_mutable(&bufsize, buffer, allocbufsize, "infile=%s", options->infilename);
3010 print_parm_br(&bufsize, buffer, allocbufsize);
3011 if(options->prioralone)
3012 {
3013 print_parm_mutable(&bufsize, buffer, allocbufsize, "NODATA=Yes");
3014 print_parm_br(&bufsize, buffer, allocbufsize);
3015 }
3016 print_parm_comment(&bufsize, buffer, allocbufsize, "Random number seed specification");
3017 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax random-seed=<AUTO | OWN:< seedfile | value >");
3018 print_parm_comment(&bufsize, buffer, allocbufsize, " AUTO uses computer system clock to generate seed");
3019 print_parm_comment(&bufsize, buffer, allocbufsize, " OWN:seedfile uses file seedfile with random number seed");
3020 print_parm_comment(&bufsize, buffer, allocbufsize, " OWN:value uses number value for seed");
3021 switch (options->autoseed)
3022 {
3023 case NOAUTO:
3024 print_parm_mutable(&bufsize, buffer, allocbufsize, "random-seed=OWN:%s",options->seedfilename);
3025 break;
3026 case AUTO:
3027 print_parm_mutable(&bufsize, buffer, allocbufsize, "random-seed=AUTO #OWN:%li", options->inseed);
3028 break;
3029 case NOAUTOSELF:
3030 print_parm_mutable(&bufsize, buffer, allocbufsize, "random-seed=OWN:%li", options->inseed);
3031 break;
3032 default:
3033 error ("error in wrinting parmfile start seed method unknown");
3034 }
3035 print_parm_br(&bufsize, buffer, allocbufsize);
3036
3037 print_parm_comment(&bufsize, buffer, allocbufsize, "Specify the title of the run, will be overridden by title in datafile");
3038 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: title=title text [up to 80 characters]");
3039 print_parm_mutable(&bufsize, buffer, allocbufsize, "title=%s", options->title);
3040 print_parm_br(&bufsize, buffer, allocbufsize);
3041 print_parm_br(&bufsize, buffer, allocbufsize);
3042 print_parm_title(&bufsize, buffer, allocbufsize, "Output options");
3043 print_parm_br(&bufsize, buffer, allocbufsize);
3044 print_parm_comment(&bufsize, buffer, allocbufsize, "Progress report to the window where the program was started");
3045 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: progress=<NO | YES | VERBOSE>");
3046 print_parm_comment(&bufsize, buffer, allocbufsize, " NO nothing is printed to the console");
3047 print_parm_comment(&bufsize, buffer, allocbufsize, " YES some messages about progress are reported [default]");
3048 print_parm_comment(&bufsize, buffer, allocbufsize, " VERBOSE more messages are reported to console");
3049 print_parm_mutable(&bufsize, buffer, allocbufsize, "progress=%s",
3050 options->progress ? (options->verbose ? "VERBOSE" : "YES") : "NO ");
3051 print_parm_br(&bufsize, buffer, allocbufsize);
3052 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3053 print_parm_br(&bufsize, buffer, allocbufsize);
3054
3055 print_parm_comment(&bufsize, buffer, allocbufsize, "Recording messages to screen into logfile");
3056 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax logfile=<NO | YES:logfilename>");
3057 print_parm_comment(&bufsize, buffer, allocbufsize, " NONE no recording of progress");
3058 print_parm_comment(&bufsize, buffer, allocbufsize, " logfilename path to logfile");
3059 if(options->writelog)
3060 print_parm_mutable(&bufsize, buffer, allocbufsize, "logfile=YES:%s", options->logfilename);
3061 else
3062 print_parm(&bufsize, buffer, allocbufsize, "logfile=NO");
3063 print_parm_br(&bufsize, buffer, allocbufsize);
3064 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3065 print_parm_br(&bufsize, buffer, allocbufsize);
3066
3067 print_parm_comment(&bufsize, buffer, allocbufsize, "Print the data as read into the program");
3068 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax print-data=<NO | YES>");
3069 print_parm_mutable(&bufsize, buffer, allocbufsize, "print-data=%s", options->printdata ? "YES" : "NO");
3070 print_parm_br(&bufsize, buffer, allocbufsize);
3071 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3072 print_parm_br(&bufsize, buffer, allocbufsize);
3073
3074 print_parm_comment(&bufsize, buffer, allocbufsize, "Print output to file [default is outfile]");
3075 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax outfile=outfilename");
3076 print_parm_mutable(&bufsize, buffer, allocbufsize, "outfile=%s", options->outfilename);
3077 print_parm_br(&bufsize, buffer, allocbufsize);
3078 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3079 print_parm_br(&bufsize, buffer, allocbufsize);
3080
3081 #ifdef PRETTY
3082 print_parm_comment(&bufsize, buffer, allocbufsize, "Print output to a PDF file [default is outfile.pdf]");
3083 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax pdf-outfile=outfilename.pdf");
3084 print_parm_mutable(&bufsize, buffer, allocbufsize, "pdf-outfile=%s", options->pdfoutfilename);
3085 print_parm_br(&bufsize, buffer, allocbufsize);
3086 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3087 print_parm_br(&bufsize, buffer, allocbufsize);
3088 #endif
3089
3090 print_parm_comment(&bufsize, buffer, allocbufsize, "Report M (=migration rate/mutation rate) instead of 4Nm or 2 Nm or Nm");
3091 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax use-M=<NO | YES> Default is YES, the name 4Nm is ambiguous");
3092 print_parm_comment(&bufsize, buffer, allocbufsize, " for non-diploid data");
3093 print_parm_mutable(&bufsize, buffer, allocbufsize, "use-M=%s", options->usem ? "YES" : "NO");
3094 print_parm_br(&bufsize, buffer, allocbufsize);
3095 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3096 print_parm_br(&bufsize, buffer, allocbufsize);
3097
3098 print_parm_comment(&bufsize, buffer, allocbufsize, "Plotting parameters: migration versus population size, such that Theta1 x immigration_.1");
3099 print_parm_comment(&bufsize, buffer, allocbufsize, "this shows the sum of all imigrations int a population");
3100 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax plot=<NO | YES:<BOTH | OUTFILE>:<LOG | STD>");
3101 print_parm_comment(&bufsize, buffer, allocbufsize, " {x-start, x-end, y-start, y-end}:<N | M>:interval>");
3102 print_parm_comment(&bufsize, buffer, allocbufsize, " NO do not show a plot");
3103 print_parm_comment(&bufsize, buffer, allocbufsize, " YES show plot with following specifications");
3104 print_parm_comment(&bufsize, buffer, allocbufsize, " BOTH print raw coordinates into MATHFILE and plot to OUTFILE");
3105 print_parm_comment(&bufsize, buffer, allocbufsize, " OUTFILE plot only to OUTFILE");
3106 print_parm_comment(&bufsize, buffer, allocbufsize, " LOG scaling of both axes");
3107 print_parm_comment(&bufsize, buffer, allocbufsize, " STD non-log scaling");
3108 print_parm_comment(&bufsize, buffer, allocbufsize, " {...} plot range of both parameters");
3109 print_parm_comment(&bufsize, buffer, allocbufsize, " N use xNm to plot immigration, x=<1,2,3,4>");
3110 print_parm_comment(&bufsize, buffer, allocbufsize, " depending on the inheritance characteristic of the data");
3111 print_parm_comment(&bufsize, buffer, allocbufsize, " M plot migration rate/mutation rate as immigration axis");
3112 print_parm_comment(&bufsize, buffer, allocbufsize, " interval the plot range is broken up into interval intervals");
3113 if (options->plot)
3114 {
3115 print_parm_mutable(&bufsize, buffer, allocbufsize, "plot=YES:%s:%s:{%f,%f,%f,%f}:%1.1s%li",
3116 options->plotmethod == PLOTALL ? "BOTH" : "OUTFILE",
3117 options->plotscale == PLOTSCALELOG ? "LOG" : "STD",
3118 options->plotrange[0],options->plotrange[1], options->plotrange[2], options->plotrange[3],
3119 options->plotvar == PLOT4NM ? "N" : "M", options->plotintervals);
3120 }
3121 else
3122 {
3123 print_parm(&bufsize, buffer, allocbufsize, "plot=NO");
3124 }
3125 print_parm_br(&bufsize, buffer, allocbufsize);
3126 print_parm_comment(&bufsize, buffer, allocbufsize, "Print plot data into a file ");
3127 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: mathtfile=mathfile the values are printed in a mathematica readable way");
3128 print_parm_mutable(&bufsize, buffer, allocbufsize, "mathfile=%s", options->mathfilename);
3129 print_parm_br(&bufsize, buffer, allocbufsize);
3130 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3131 print_parm_br(&bufsize, buffer, allocbufsize);
3132
3133 print_parm_comment(&bufsize, buffer, allocbufsize, "Profile likelihood for each estimated parameter");
3134 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax profile=<NONE | <ALL | TABLES | SUMMARY>:");
3135 print_parm_comment(&bufsize, buffer, allocbufsize, " <PRECISE | DISCRETE | QUICK | FAST> >");
3136 print_parm_comment(&bufsize, buffer, allocbufsize, " NONE do not calculate profile likelihoods");
3137 print_parm_comment(&bufsize, buffer, allocbufsize, " ALL print individual profile tables and summary [default]");
3138 print_parm_comment(&bufsize, buffer, allocbufsize, " TABLES show only tables and no summary");
3139 print_parm_comment(&bufsize, buffer, allocbufsize, " SUMMARY show only summary");
3140 print_parm_comment(&bufsize, buffer, allocbufsize, " PRECISE evaluate profile likelihood at percentiles [Default]");
3141 // print_parm_comment(&bufsize, buffer, allocbufsize, " SPLINE use splines to calculate precentile values [DOES NOT WORK!]");
3142 print_parm_comment(&bufsize, buffer, allocbufsize, " QUICK assumes that there is no interaction of parameters");
3143 print_parm_comment(&bufsize, buffer, allocbufsize, " FAST same as QUICK except in last calculation cycle assumes interaction");
3144 print_parm_comment(&bufsize, buffer, allocbufsize, " DISCRETE uses fixed mutipliers: 0.02,0.1,0.2,0.5,1,2,5,10,50");
3145 switch (options->profilemethod)
3146 {
3147 case 's':
3148 sprintf (fp, ":SPLINE");
3149 break;
3150 case 'd':
3151 sprintf (fp, ":DISCRETE");
3152 break;
3153 case 'q':
3154 sprintf (fp, ":QUICKANDDIRTY");
3155 break;
3156 case 'f':
3157 sprintf (fp, ":FAST");
3158 break;
3159 case 'p':
3160 default:
3161 sprintf (fp, ":PRECISE");
3162 break;
3163 }
3164
3165 switch (options->profile)
3166 {
3167 case ALL:
3168 print_parm_mutable(&bufsize, buffer, allocbufsize, "profile=ALL%s",fp);
3169 break;
3170 case TABLES:
3171 print_parm_mutable(&bufsize, buffer, allocbufsize, "profile=TABLES%s",fp);
3172 break;
3173 case SUMMARY:
3174 print_parm_mutable(&bufsize, buffer, allocbufsize, "profile=SUMMARY%s",fp);
3175 break;
3176 case _NONE:
3177 default:
3178 print_parm_mutable(&bufsize, buffer, allocbufsize, "profile=NONE");
3179 break;
3180 }
3181 print_parm_br(&bufsize, buffer, allocbufsize);
3182 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3183 print_parm_br(&bufsize, buffer, allocbufsize);
3184
3185
3186 print_parm_comment(&bufsize, buffer, allocbufsize, "Print tree into treefile");
3187 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax print-tree=< NONE | <ALL | BEST | LASTCHAIN:Increment>:treefile >");
3188 print_parm_comment(&bufsize, buffer, allocbufsize, " NONE no tree printed [Default, and only choice using parallel");
3189 print_parm_comment(&bufsize, buffer, allocbufsize, " ALL print all visited genealogies [careful this will be huge]");
3190 print_parm_comment(&bufsize, buffer, allocbufsize, " BEST print only the best tree visited");
3191 print_parm_comment(&bufsize, buffer, allocbufsize, " LASTCHAIN print all trees in last chain");
3192 print_parm_comment(&bufsize, buffer, allocbufsize, " with increment INCREMENT");
3193 switch (options->treeprint)
3194 {
3195 case _NONE:
3196 print_parm(&bufsize, buffer, allocbufsize, "print-tree=NONE");
3197 break;
3198 case ALL:
3199 print_parm_mutable(&bufsize, buffer, allocbufsize, "print-tree=ALL:%s",options->treefilename);
3200 break;
3201 case BEST:
3202 print_parm_mutable(&bufsize, buffer, allocbufsize, "print-tree=BEST:%s",options->treefilename);
3203 break;
3204 case LASTCHAIN:
3205 print_parm_mutable(&bufsize, buffer, allocbufsize, "print-tree=LASTCHAIN:%li:%s", options->treeinc, options->treefilename);
3206 break;
3207 default:
3208 print_parm(&bufsize, buffer, allocbufsize, "print-tree=NONE");
3209 break;
3210 }
3211 print_parm_br(&bufsize, buffer, allocbufsize);
3212 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3213 print_parm_br(&bufsize, buffer, allocbufsize);
3214
3215
3216 print_parm_comment(&bufsize, buffer, allocbufsize, "write intermediate minimal statistics into a file for later use");
3217 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax write-summary=<NO | YES:SUMFILE >");
3218 print_parm_comment(&bufsize, buffer, allocbufsize, " Default is NO, with YES the user needs to ");
3219 print_parm_comment(&bufsize, buffer, allocbufsize, " give a file to record the summary statistics");
3220 if(options->writesum)
3221 print_parm_mutable(&bufsize, buffer, allocbufsize, "write-summary=YES:%s",options->sumfilename);
3222 else
3223 print_parm(&bufsize, buffer, allocbufsize, "write-summary=NO");
3224
3225 print_parm_br(&bufsize, buffer, allocbufsize);
3226 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3227 print_parm_br(&bufsize, buffer, allocbufsize);
3228
3229 print_parm_comment(&bufsize, buffer, allocbufsize, "Likelihood ratio test");
3230 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax l-ratio=<NO | YES:values_to_test>");
3231 print_parm_comment(&bufsize, buffer, allocbufsize, " Values_to_test are compared to the values generated in the run");
3232 print_parm_comment(&bufsize, buffer, allocbufsize, " values_to_test={ab..bbab..ba ... a}");
3233 print_parm_comment(&bufsize, buffer, allocbufsize, " the {} is a square matrix with values for the population sizes");
3234 print_parm_comment(&bufsize, buffer, allocbufsize, " on the diagonal and migration rates off-diagonal");
3235 print_parm_comment(&bufsize, buffer, allocbufsize, " the values a for the diagonal can be any of these:");
3236 print_parm_comment(&bufsize, buffer, allocbufsize, " number constant, the value is for example 0.002");
3237 print_parm_comment(&bufsize, buffer, allocbufsize, " * free to vary, the default is * for every parameter");
3238 print_parm_comment(&bufsize, buffer, allocbufsize, " m mean of theta, this can be a subgroup of all thetas");
3239 print_parm_comment(&bufsize, buffer, allocbufsize, " for example the theta 1-3 are averaged and thetas 4,5 are estimated");
3240 print_parm_comment(&bufsize, buffer, allocbufsize, " the values b for the migration rates can be any of these:");
3241 print_parm_comment(&bufsize, buffer, allocbufsize, " number constant, the value is for example 45.0 or 0.0");
3242 print_parm_comment(&bufsize, buffer, allocbufsize, " * free to vary, the default is * for every parameter");
3243 print_parm_comment(&bufsize, buffer, allocbufsize, " m mean of M_ij, this can be a subgroup of migration rates");
3244 print_parm_comment(&bufsize, buffer, allocbufsize, " for example the M_1-3i are averaged and M_4,5i are estimated");
3245 print_parm_comment(&bufsize, buffer, allocbufsize, " M means of 4Nm (diploid), 2Nm (haploid), Nm (mtDNA, Y-chromosome)");
3246 print_parm_comment(&bufsize, buffer, allocbufsize, " s symmetric migration rates M");
3247 print_parm_comment(&bufsize, buffer, allocbufsize, " S symmetric migrants 4Nm");
3248 print_parm_comment(&bufsize, buffer, allocbufsize, " an example for 5 populations could look like this:");
3249 print_parm_comment(&bufsize, buffer, allocbufsize, " l-ratio=YES:{*s00s s*s00 0s*s0 00s*s s00s*");
3250 print_parm_comment(&bufsize, buffer, allocbufsize, " this describes a circular stepping stone model with 5 symmetric rates");
3251 print_parm_comment(&bufsize, buffer, allocbufsize, " and independent sizes, a very basic stepping stone with 2 parameters would");
3252 print_parm_comment(&bufsize, buffer, allocbufsize, " look like this l-ratio=YES:{mm00m mmm00 0mmm0 00mmm m00mm}");
3253 print_parm_comment(&bufsize, buffer, allocbufsize, " [The L-RATIO statement can be repeated]");
3254 print_parm_comment(&bufsize, buffer, allocbufsize, " Default: l-ratio=NO");
3255 for (i = 0; i < options->lratio->counter; i++)
3256 {
3257 if (options->lratio->data[i].type == MLE)
3258 print_parm_mutable(&bufsize, buffer, allocbufsize, "l-ratio=%s:%s", "YES",
3259 options->lratio->data[i].value1);
3260 else
3261 print_parm_mutable(&bufsize, buffer, allocbufsize, "l-ratio=%s","NO");
3262
3263 /* this code is unused, but this needs careful checking else where
3264 else
3265 print_parm_mutable(&bufsize, buffer, allocbufsize, "l-ratio=%s:%s:%s", "ARBITRARY",
3266 options->lratio->data[i].value1,
3267 options->lratio->data[i].value2);
3268 */
3269 }
3270 print_parm_br(&bufsize, buffer, allocbufsize);
3271 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3272 print_parm_br(&bufsize, buffer, allocbufsize);
3273
3274 print_parm_comment(&bufsize, buffer, allocbufsize, "AIC model selection [do not use yet, will come in Summer 2004]");
3275 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax aic-modeltest=<NO | YES:<FAST | EXHAUSTIVE>>");
3276 print_parm_comment(&bufsize, buffer, allocbufsize, " FAST [do not use yet]");
3277 print_parm_comment(&bufsize, buffer, allocbufsize, " EXHAUSTIVE [do not use yet]");
3278 if (options->aic)
3279 {
3280 if (options->fast_aic)
3281 print_parm(&bufsize, buffer, allocbufsize, "aic-modeltest=YES:FAST");
3282 else
3283 print_parm(&bufsize, buffer, allocbufsize, "aic-modeltest=YES");
3284 }
3285 else
3286 print_parm(&bufsize, buffer, allocbufsize, "aic-modeltest=NO");
3287 print_parm_br(&bufsize, buffer, allocbufsize);
3288
3289 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3290 print_parm_br(&bufsize, buffer, allocbufsize);
3291 print_parm_comment(&bufsize, buffer, allocbufsize, "Print a histogram of the time of migration events for each M(i->j)");
3292 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax mig-histogram=<NO | <ALL | MIGRATIONEVENTSONLY>:binsize:mighistfile >");
3293 print_parm_comment(&bufsize, buffer, allocbufsize, " NO do not record any events");
3294 print_parm_comment(&bufsize, buffer, allocbufsize, " ALL record migration and coalescence event");
3295 print_parm_comment(&bufsize, buffer, allocbufsize, " MIGRATIONEVENTSONLY record only migration events");
3296 print_parm_comment(&bufsize, buffer, allocbufsize, " binsize has to be in mutation units, with an average Theta=0.01 try 0.001");
3297 print_parm_comment(&bufsize, buffer, allocbufsize, "Print a histogram of the parameters through time (skyline plot)");
3298 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax skyline=<NO | YES>:binsize:skylinefile >");
3299 print_parm_comment(&bufsize, buffer, allocbufsize, " NO do not calculate parameter estimates through time");
3300 print_parm_comment(&bufsize, buffer, allocbufsize, " YES calculate parameters through time");
3301 print_parm_comment(&bufsize, buffer, allocbufsize, " binsize has to be in mutation units, with an average Theta=0.01 try 0.001");
3302 print_parm_comment(&bufsize, buffer, allocbufsize, " If the interval is too fine the output will be very noisy");
3303
3304 if(options->mighist)
3305 {
3306 if(options->mighist_all)
3307 print_parm_mutable(&bufsize, buffer, allocbufsize, "mig-histogram=ALL:%f:%s", options->eventbinsize,options->mighistfilename);
3308 else
3309 print_parm_mutable(&bufsize, buffer, allocbufsize, "mig-histogram=YES:%f:%s", options->eventbinsize, options->mighistfilename);
3310 if(options->skyline)
3311 print_parm_mutable(&bufsize, buffer, allocbufsize, "skyline=YES:%f:%s #needs mig-histogram", options->eventbinsize, options->skylinefilename);
3312 }
3313 else
3314 {
3315 print_parm(&bufsize, buffer, allocbufsize, "mig-histogram=NO");
3316 print_parm(&bufsize, buffer, allocbufsize, "skyline=NO #needs mig-histogram=ALL:...");
3317 }
3318 print_parm_br(&bufsize, buffer, allocbufsize);
3319
3320
3321 print_parm_br(&bufsize, buffer, allocbufsize);
3322 print_parm_title(&bufsize, buffer, allocbufsize, "Parameter start settings");
3323 print_parm_br(&bufsize, buffer, allocbufsize);
3324 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: theta=<FST | OWN:<{value} | {value1, value2, ...., valuen} | NRANDOM:{mean std} | URANDOM{min,max}>");
3325 print_parm_comment(&bufsize, buffer, allocbufsize, " migrationt=<FST | OWN:<{value} | {value1, value2, ...., valuen} | NRANDOM:{mean std} | URANDOM{min,max}>");
3326 print_parm_comment(&bufsize, buffer, allocbufsize, " FST starting parameter are derived from");
3327 print_parm_comment(&bufsize, buffer, allocbufsize, " an FST-like calculation (Beerli&Felsenstein 1999");
3328 print_parm_comment(&bufsize, buffer, allocbufsize, " OWN starting values are supplied by user");
3329 print_parm_comment(&bufsize, buffer, allocbufsize, " {value} if only one value is supplied then all population");
3330 print_parm_comment(&bufsize, buffer, allocbufsize, " have the same starting value");
3331 print_parm_comment(&bufsize, buffer, allocbufsize, " {value1, value2, ..., valuen} each population has its");
3332 print_parm_comment(&bufsize, buffer, allocbufsize, " own starting value, if the number of values is");
3333 print_parm_comment(&bufsize, buffer, allocbufsize, " insuffient, then the last value is the template");
3334 print_parm_comment(&bufsize, buffer, allocbufsize, " for the remaining populations");
3335 print_parm_comment(&bufsize, buffer, allocbufsize, " NRANDOM starting parameter is drawn randomely from a Normal distribution");
3336 print_parm_comment(&bufsize, buffer, allocbufsize, " {mean std} with mean and standard deviation");
3337 print_parm_comment(&bufsize, buffer, allocbufsize, " URANDOM starting parameter is drawn randomely from a Uniform distribution");
3338 print_parm_comment(&bufsize, buffer, allocbufsize, " {min max} with minimum and maximum values");
3339 switch(options->thetaguess)
3340 {
3341 case FST:
3342 print_parm(&bufsize, buffer, allocbufsize, "theta=FST");
3343 break;
3344 case NRANDOMESTIMATE:
3345 print_parm_mutable(&bufsize, buffer, allocbufsize, "theta=NRANDOM: {%f %f}",options->thetag[0],options->thetag[1]);
3346 break;
3347 case URANDOMESTIMATE:
3348 print_parm_mutable(&bufsize, buffer, allocbufsize, "theta=URANDOM: {%f %f}",options->thetag[0],options->thetag[1]);
3349 break;
3350
3351 default:
3352 print_parm_theta(&bufsize, buffer, allocbufsize, options);
3353 break;
3354 }
3355
3356 switch (options->migrguess)
3357 {
3358 case FST:
3359 print_parm(&bufsize, buffer, allocbufsize, "migration=FST");
3360 break;
3361 case NRANDOMESTIMATE:
3362 print_parm_mutable(&bufsize, buffer, allocbufsize, "migration=NRANDOM: {%f %f}",options->mg[0],options->mg[1]);
3363 break;
3364 case URANDOMESTIMATE:
3365 print_parm_mutable(&bufsize, buffer, allocbufsize, "migration=URANDOM: {%f %f}",options->mg[0],options->mg[1]);
3366 break;
3367 default:
3368 print_parm_m(&bufsize, buffer, allocbufsize, options);
3369 break;
3370 }
3371 print_parm_br(&bufsize, buffer, allocbufsize);
3372
3373 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3374 print_parm_comment(&bufsize, buffer, allocbufsize, "Mutation rate modifiers");
3375 print_parm_br(&bufsize, buffer, allocbufsize);
3376 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: mutation=<NOGAMMA | CONSTANT | ESTIMATE | GAMMA:alpha | OWN:loci: rate1 rate2 ... rate_loci>");
3377 print_parm_comment(&bufsize, buffer, allocbufsize, " NOGAMMA all loci have same mutation rate");
3378 print_parm_comment(&bufsize, buffer, allocbufsize, " CONSTANT all loci have same mutation rate");
3379 print_parm_comment(&bufsize, buffer, allocbufsize, " ESTIMATE BAYESIAN estimate: mutation rate is drawn from prior");
3380 print_parm_comment(&bufsize, buffer, allocbufsize, " GAMMA:alpha ML estimate: mutation rate has Gamma distribution with alpha");
3381 print_parm_comment(&bufsize, buffer, allocbufsize, " OWN mutation rate is different for every locus, but fixed");
3382 print_parm_comment(&bufsize, buffer, allocbufsize, " :loci: rate1, ... number of loci, rate of locus 1, locus 2 etc.");
3383 print_parm_comment(&bufsize, buffer, allocbufsize, " DATA mutation rate modifier is deducted from loci in the data");
3384 print_parm_comment(&bufsize, buffer, allocbufsize, " using Watterson's Theta and then scaling all rates Theta_locus/mean(Theta_loci");
3385
3386 if (options->gamma)
3387 {
3388 print_parm_mutable(&bufsize, buffer, allocbufsize, "mutation=%s:%f", "GAMMA", options->alphavalue);
3389 }
3390 else
3391 {
3392 if (options->murates)
3393 {
3394 if(options->murates_fromdata)
3395 {
3396 sprintf (fp, "mutation=DATA\n");
3397 add_to_buffer(fp,&bufsize,buffer, allocbufsize);
3398 }
3399 else
3400 {
3401 sprintf (fp, "mutation=OWN:%li: ", options->muloci);
3402 add_to_buffer(fp,&bufsize,buffer, allocbufsize);
3403 for (i = 0; i < options->muloci; i++)
3404 {
3405 sprintf (fp, "%f ", options->mu_rates[i]);
3406 add_to_buffer(fp,&bufsize,buffer, allocbufsize);
3407 }
3408 sprintf (fp, " \n");
3409 add_to_buffer(fp,&bufsize,buffer, allocbufsize);
3410 }
3411 }
3412 else
3413 {
3414 if(options->bayesmurates)
3415 {
3416 print_parm(&bufsize, buffer, allocbufsize, "mutation=ESTIMATE");
3417 }
3418 else
3419 {
3420 print_parm(&bufsize, buffer, allocbufsize, "mutation=CONSTANT");
3421 }
3422 }
3423 }
3424 print_parm_br(&bufsize, buffer, allocbufsize);
3425 print_parm_comment(&bufsize, buffer, allocbufsize, "Treatment of inviariant sequence loci");
3426 print_parm_comment(&bufsize, buffer, allocbufsize, "Syntax: analyze-loci=<A | F | V>");
3427 print_parm_comment(&bufsize, buffer, allocbufsize, " A = analyze all loci (Default!)");
3428 print_parm_comment(&bufsize, buffer, allocbufsize, " F = analyze all variable loci and ONE invariant and extrapolate");
3429 print_parm_comment(&bufsize, buffer, allocbufsize, " V = analyze only variable loci");
3430 print_parm_mutable(&bufsize, buffer, allocbufsize, "#analyze-loci=%c",( options->onlyvariable ? 'V' :(options->has_variableandone ? 'F' : 'A')));
3431 print_parm_br(&bufsize, buffer, allocbufsize);
3432 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3433 print_parm_comment(&bufsize, buffer, allocbufsize, "FST model");
3434 print_parm_br(&bufsize, buffer, allocbufsize);
3435 print_parm_mutable(&bufsize, buffer, allocbufsize, "fst-type=%s", options->fsttype ? "THETA" : "MIGRATION");
3436 print_parm_br(&bufsize, buffer, allocbufsize);
3437
3438 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3439 print_parm_comment(&bufsize, buffer, allocbufsize, "Custom migration model");
3440 print_parm_br(&bufsize, buffer, allocbufsize);
3441 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: custom-migration={ab..bbab..ba ... a}");
3442 print_parm_comment(&bufsize, buffer, allocbufsize, " the {} is a square matrix with values for the population sizes");
3443 print_parm_comment(&bufsize, buffer, allocbufsize, " on the diagonal and migration rates off-diagonal");
3444 print_parm_comment(&bufsize, buffer, allocbufsize, " the values a for the diagonal can be any of these:");
3445 print_parm_comment(&bufsize, buffer, allocbufsize, " c constant, the value needs to be defined in the theta option");
3446 print_parm_comment(&bufsize, buffer, allocbufsize, " * free to vary, the default is * for every parameter");
3447 print_parm_comment(&bufsize, buffer, allocbufsize, " m mean of theta, this can be a subgroup of all thetas");
3448 print_parm_comment(&bufsize, buffer, allocbufsize, " for example the theta 1-3 are averaged and thetas 4,5 are estimated");
3449 print_parm_comment(&bufsize, buffer, allocbufsize, " the values b for the migration rates can be any of these:");
3450 print_parm_comment(&bufsize, buffer, allocbufsize, " c constant, the value needs to be defined in the migration option");
3451 print_parm_comment(&bufsize, buffer, allocbufsize, " * free to vary, the default is * for every parameter");
3452 print_parm_comment(&bufsize, buffer, allocbufsize, " m mean of M_ij, this can be a subgroup of migration rates");
3453 print_parm_comment(&bufsize, buffer, allocbufsize, " for example the M_1-3i are averaged and M_4,5i are estimated");
3454 print_parm_comment(&bufsize, buffer, allocbufsize, " M means of 4Nm (diploid), 2Nm (haploid), Nm (mtDNA, Y-chromosome)");
3455 print_parm_comment(&bufsize, buffer, allocbufsize, " s symmetric migration rates M");
3456 print_parm_comment(&bufsize, buffer, allocbufsize, " S symmetric migrants 4Nm");
3457 print_parm_comment(&bufsize, buffer, allocbufsize, " an example for 5 populations could look like this:");
3458 print_parm_comment(&bufsize, buffer, allocbufsize, " custom-migration={*s00s s*s00 0s*s0 00s*s s00s*");
3459 print_parm_comment(&bufsize, buffer, allocbufsize, " this describes a circular stepping stone model with 5 symmetric rates");
3460 print_parm_comment(&bufsize, buffer, allocbufsize, " and independent sizes, a very basic stepping stone with 2 parameters would");
3461 print_parm_comment(&bufsize, buffer, allocbufsize, " look like this custom-migration={mm00m mmm00 0mmm0 00mmm m00mm}");
3462 print_parm_mutable(&bufsize, buffer, allocbufsize, "custom-migration={%s}", options->custm);
3463 print_parm_br(&bufsize, buffer, allocbufsize);
3464 print_parm_comment(&bufsize, buffer, allocbufsize, "Influence of geography on migration rate");
3465 print_parm_comment(&bufsize, buffer, allocbufsize, "a distance matrix between populations changes the migration rate matrix so that");
3466 print_parm_comment(&bufsize, buffer, allocbufsize, "(genetic?) migration rates = inferred migration rate / distance ~ a dispersion coefficient");
3467 print_parm_comment(&bufsize, buffer, allocbufsize, "the geofile contains a number of populations, names for populations (10 characters), they");
3468 print_parm_comment(&bufsize, buffer, allocbufsize, "need to be in order of the dataset. And the distances between the populations, they do not");
3469 print_parm_comment(&bufsize, buffer, allocbufsize, "need to be symmetric");
3470 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: geo:<NO | YES:filename>");
3471 print_parm_comment(&bufsize, buffer, allocbufsize, " NO distances among populations are considered to be 1 [all equal]");
3472 print_parm_comment(&bufsize, buffer, allocbufsize, " YES distances are read from a file");
3473
3474 if (options->geo)
3475 {
3476 print_parm_mutable(&bufsize, buffer, allocbufsize, "geo=YES:%s", options->geofilename);
3477 }
3478 else
3479 {
3480 print_parm(&bufsize, buffer, allocbufsize, "geo=NO");
3481 }
3482 print_parm_br(&bufsize, buffer, allocbufsize);
3483 print_parm_br(&bufsize, buffer, allocbufsize);
3484 // SEARCH STRATEGIES
3485 print_parm_title(&bufsize, buffer, allocbufsize, "Search strategies");
3486 print_parm_br(&bufsize, buffer, allocbufsize);
3487
3488 print_parm_comment(&bufsize, buffer, allocbufsize, "MCMC Strategy method");
3489 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: bayes-update=< NO | YES>");
3490 print_parm_comment(&bufsize, buffer, allocbufsize, " NO maximum likelihood method");
3491 print_parm_comment(&bufsize, buffer, allocbufsize, " YES Bayesian method");
3492 print_parm_comment(&bufsize, buffer, allocbufsize, "Some of the options are only available in one or other mode");
3493 print_parm_comment(&bufsize, buffer, allocbufsize, "BAYESIAN OPTIONS");
3494 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-updatefreq=VALUE ");
3495 print_parm_comment(&bufsize, buffer, allocbufsize, " VALUE is a ratio between 0 and 1");
3496 print_parm_comment(&bufsize, buffer, allocbufsize, " ratio of how many times the genealogy is updated compared to the parameters");
3497 print_parm_comment(&bufsize, buffer, allocbufsize, " If the value is 0.4 in a 2 population scenario and with 1000000 steps");
3498 print_parm_comment(&bufsize, buffer, allocbufsize, " The tree will be evaluated 400000 times, Theta_1, Theta_2, M_21, and M_12");
3499 print_parm_comment(&bufsize, buffer, allocbufsize, " will be each evaluated 125000 times.");
3500 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-posteriorbins=VALUE VALUE");
3501 print_parm_comment(&bufsize, buffer, allocbufsize, " VALUE is the number of bins in the psterior distribution histogram for Theta or M");
3502 #ifdef PRETTY
3503 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-posteriormaxtype=< ALL | P99 | MAXP99 | P100 >");
3504 print_parm_comment(&bufsize, buffer, allocbufsize, " ALL plots the WHOLE prior-parameter range\n");
3505 print_parm_comment(&bufsize, buffer, allocbufsize, " P99 plots from the minimum prior range value to\n");
3506 print_parm_comment(&bufsize, buffer, allocbufsize, " the 99% percentile value of EACH parameter\n");
3507 print_parm_comment(&bufsize, buffer, allocbufsize, " MAXP99 sets all axes from minimum to the maximal\n");
3508 print_parm_comment(&bufsize, buffer, allocbufsize, " 99% percentile value of ALL parameter\n");
3509 print_parm_comment(&bufsize, buffer, allocbufsize, " P100 plots from the minimum prior range value to\n");
3510 print_parm_comment(&bufsize, buffer, allocbufsize, " the 100% percentile value of EACH parameter\n");
3511 #endif
3512 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-file=<YES:FILENAME|NO>");
3513 print_parm_comment(&bufsize, buffer, allocbufsize, " FILENAME is the name of the file that will contain");
3514 print_parm_comment(&bufsize, buffer, allocbufsize, " the results for the posterior distribution");
3515 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-allfile=<<YES|TEMP>:INTERVAL:FILENAME|NO>");
3516 print_parm_comment(&bufsize, buffer, allocbufsize, " FILENAME is the name of the file that will contain");
3517 print_parm_comment(&bufsize, buffer, allocbufsize, " all parameters of the posterior distribution [HUGE]");
3518 //OLD print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-allfileinterval=INTERVAL");
3519 print_parm_comment(&bufsize, buffer, allocbufsize, " INTERVAL is the interval at which all parameters are written to file\n");
3520 print_parm_comment(&bufsize, buffer, allocbufsize, " ");
3521 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-proposals= THETA < SLICE | METROPOLIS >");
3522 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-proposals= MIG < SLICE | METROPOLIS >");
3523 // print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-proposal= RATE < SLICE | METROPOLIS >");
3524 print_parm_comment(&bufsize, buffer, allocbufsize, " SLICE uses the slice sampler to propose new parameter values");
3525 print_parm_comment(&bufsize, buffer, allocbufsize, " METROPOLIS uses the Metropolis-Hastings sampler");
3526 print_parm_comment(&bufsize, buffer, allocbufsize, " (this is done for each parameter group: THETA or MIGration)");
3527 print_parm_comment(&bufsize, buffer, allocbufsize, " ");
3528 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-priors= THETA <UNIFORM unipriorvalues | EXP exppriorvalues | WINDOWEXP wexppriorvalues ");//| MULT multpriorvalues | FIXED>");
3529 print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-priors= MIG <UNIFORM unipriorvalues | EXP exppriorvalues | WINDOWEXP wexppriorvalues ");//| MULT multpriorvalues | FIXED>");
3530 // print_parm_comment(&bufsize, buffer, allocbufsize, " bayes-priors= RATE <UNIFORM unipriorvalues | EXP exppriorvalues | WINDOWEXP wexppriorvalues | MULT multpriorvalues | FIXED>");
3531 print_parm_comment(&bufsize, buffer, allocbufsize, " unipriorvalues: min max delta");
3532 print_parm_comment(&bufsize, buffer, allocbufsize, " exppriorvalues: min mean max");
3533 print_parm_comment(&bufsize, buffer, allocbufsize, " wexppriorvalues: min mean max delta");
3534 // print_parm_comment(&bufsize, buffer, allocbufsize, " multpriorvalues: min max mult");
3535 //print_parm_comment(&bufsize, buffer, allocbufsize, " fixed = no value is changed during the run");
3536 print_parm_br(&bufsize, buffer, allocbufsize);
3537
3538 print_parm_comment(&bufsize, buffer, allocbufsize, "Maximum likelihood OPTIONS");
3539 print_parm_comment(&bufsize, buffer, allocbufsize, " short-chains=VALUE VALUE is 1..n [Default is 10]");
3540 print_parm_comment(&bufsize, buffer, allocbufsize, " short-inc=VALUE VALUE is the number of updates that are not recorded");
3541 print_parm_comment(&bufsize, buffer, allocbufsize, " short-sample=VALUE VALUE is the number of sampled updates");
3542 print_parm_br(&bufsize, buffer, allocbufsize);
3543
3544 print_parm_comment(&bufsize, buffer, allocbufsize, "Search OPTIONS for both strategies");
3545 print_parm_comment(&bufsize, buffer, allocbufsize, " long-chains=VALUE VALUE is 1..n [Default is 3]");
3546 print_parm_comment(&bufsize, buffer, allocbufsize, " long-inc=VALUE VALUE is the number of updates that are not recorded");
3547 print_parm_comment(&bufsize, buffer, allocbufsize, " long-sample=VALUE VALUE is the number of sampled updates");
3548 print_parm_comment(&bufsize, buffer, allocbufsize, " burn-in=VALUE VALUE is the number of updates to discard at the beginning");
3549 print_parm_comment(&bufsize, buffer, allocbufsize, " auto-tune=YES:VALUE VALUE the the target acceptance ratio");
3550 print_parm_comment(&bufsize, buffer, allocbufsize, " if value is missing, it is set to 0.44");
3551 print_parm_br(&bufsize, buffer, allocbufsize);
3552
3553 if(options->bayes_infer)
3554 {
3555 print_parm(&bufsize, buffer, allocbufsize, "bayes-update=YES");
3556 print_parm_mutable(&bufsize, buffer, allocbufsize, "bayes-updatefreq=%f", options->updateratio);
3557 if(options->bayesmurates)
3558 print_parm_mutable(&bufsize, buffer, allocbufsize, "bayes-posteriorbins=%li %li %li",
3559 options->bayespriortheta->bins,
3560 options->bayespriorm->bins,
3561 options->bayespriorrate->bins);
3562 else
3563 print_parm_mutable(&bufsize, buffer, allocbufsize, "bayes-posteriorbins=%li %li",
3564 options->bayespriortheta->bins,
3565 options->bayespriorm->bins);
3566 print_parm_mutable(&bufsize, buffer, allocbufsize, "bayes-posteriormaxtype=%s",
3567 (options->bayespretty == PRETTY_P99 ? "P99" :
3568 (options->bayespretty == PRETTY_MAX ? "ALL" :
3569 (options->bayespretty == PRETTY_P100 ? "TOTAL" : "TOTAL"))));
3570
3571 if(options->has_bayesfile)
3572 print_parm_mutable(&bufsize, buffer, allocbufsize, "bayes-file=YES:%s", options->bayesfilename);
3573 else
3574 print_parm_mutable(&bufsize, buffer, allocbufsize, "bayes-file=NO");
3575
3576 if(options->has_bayesmdimfile)
3577 {
3578 print_parm_mutable(&bufsize, buffer, allocbufsize, "bayes-allfile=%s:%li:%s",
3579 (options->mdimdelete ? "TEMP" : "YES"),
3580 options->bayesmdiminterval, options->bayesmdimfilename);
3581 }
3582 else
3583 print_parm_mutable(&bufsize, buffer, allocbufsize, "bayes-allfile=NO");
3584
3585 print_parm_proposal(&bufsize, buffer, allocbufsize, options);
3586 print_parm_prior(&bufsize, buffer, allocbufsize, options);
3587 }
3588 else
3589 {
3590 print_parm(&bufsize, buffer, allocbufsize, "bayes-update=NO");
3591 print_parm_mutable(&bufsize, buffer, allocbufsize, "short-chains=%li", options->schains);
3592 print_parm_mutable(&bufsize, buffer, allocbufsize, "short-inc=%li", options->sincrement);
3593 print_parm_mutable(&bufsize, buffer, allocbufsize, "short-sample=%li", options->ssteps);
3594 }
3595 print_parm_br(&bufsize, buffer, allocbufsize);
3596 print_parm_mutable(&bufsize, buffer, allocbufsize, "long-chains=%li", options->lchains);
3597 print_parm_mutable(&bufsize, buffer, allocbufsize, "long-inc=%li", options->lincrement);
3598 print_parm_mutable(&bufsize, buffer, allocbufsize, "long-sample=%li", options->lsteps);
3599 print_parm_mutable(&bufsize, buffer, allocbufsize, "burn-in=%li %c", options->burn_in, options->burnin_autostop);
3600 if (options->has_autotune)
3601 print_parm_mutable(&bufsize, buffer, allocbufsize, "auto-tune=YES:%f", options->autotune);
3602 else
3603 print_parm_mutable(&bufsize, buffer, allocbufsize, "auto-tune=NO");
3604
3605
3606 print_parm_br(&bufsize, buffer, allocbufsize);
3607
3608 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3609 print_parm_comment(&bufsize, buffer, allocbufsize, "Schemes to improve MCMC searching and/or thermodynamic integration");
3610 print_parm_br(&bufsize, buffer, allocbufsize);
3611 print_parm_comment(&bufsize, buffer, allocbufsize, "Heating schemes {MCMCMC = MC cubed}");
3612 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: heating=< NO | <YES | ADAPTIVE>:SKIP:TEMPERATURES");
3613 print_parm_comment(&bufsize, buffer, allocbufsize, " NO No heating");
3614 print_parm_comment(&bufsize, buffer, allocbufsize, " YES heating using TEMPERATURES");
3615 print_parm_comment(&bufsize, buffer, allocbufsize, " ADAPTIVE adaptive heating using start TEMPERATURES [fails sometimes!!]");
3616 print_parm_comment(&bufsize, buffer, allocbufsize, " SKIP skip that many comparisons, this lengthens the run by SKIP");
3617 print_parm_comment(&bufsize, buffer, allocbufsize, " TEMPERATURES { 1.0, temp1, temp2, temp3 .. tempn}");
3618 print_parm_comment(&bufsize, buffer, allocbufsize, " Example: heating=YES:1:{1.0, 1.2, 3.0,1000000.0}");
3619 print_parm_comment(&bufsize, buffer, allocbufsize, "Heating: swapping chains");
3620 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: heated-swap=< YES | NO >");
3621 print_parm_comment(&bufsize, buffer, allocbufsize, " YES swapping of chains enabled [DEFAULT]");
3622 print_parm_comment(&bufsize, buffer, allocbufsize, " NO swapping of chains disabled");
3623 print_parm_comment(&bufsize, buffer, allocbufsize, " Example: heated-swap=YES");
3624 print_parm_heating(&bufsize, buffer, allocbufsize, options);
3625
3626 print_parm_br(&bufsize, buffer, allocbufsize);
3627 print_parm_comment(&bufsize, buffer, allocbufsize, "Lengthening chain schemes");
3628 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: moving-steps=< NO | YES:VALUE>");
3629 print_parm_comment(&bufsize, buffer, allocbufsize, " VALUE frequency is between 0..1");
3630
3631 if (options->movingsteps)
3632 {
3633 print_parm_mutable(&bufsize, buffer, allocbufsize, "moving-steps=YES:%f", options->acceptfreq);
3634 }
3635 else
3636 {
3637 print_parm(&bufsize, buffer, allocbufsize, "moving-steps=NO");
3638 }
3639 print_parm_br(&bufsize, buffer, allocbufsize);
3640 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: long-chain-epsilon=VALUE");
3641 print_parm_comment(&bufsize, buffer, allocbufsize, " VALUE is between 0..INFINITY");
3642 print_parm_comment(&bufsize, buffer, allocbufsize, " the VALUE is the likelihood ratio between the old and thew chain");
3643 print_parm_comment(&bufsize, buffer, allocbufsize, " the VALUE depends on the number of parameters: with 1 values of 0.5 are great");
3644 print_parm_comment(&bufsize, buffer, allocbufsize, " but with many parameters values and bad data >20 is more reasonable");
3645 if (options->lcepsilon < LONGCHAINEPSILON)
3646 print_parm_mutable(&bufsize, buffer, allocbufsize, "long-chain-epsilon=%f", options->lcepsilon);
3647 else
3648 print_parm_mutable(&bufsize, buffer, allocbufsize, "long-chain-epsilon=INFINITY");
3649
3650
3651 print_parm_br(&bufsize, buffer, allocbufsize);
3652 print_parm_comment(&bufsize, buffer, allocbufsize, " Convergence statistic [Gelman and Rubin]");
3653 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: gelman-convergence=< YES:Pairs|Summary | NO >");
3654 print_parm_comment(&bufsize, buffer, allocbufsize, " NO do not use Gelman's convergence criterium");
3655 print_parm_comment(&bufsize, buffer, allocbufsize, " YES use Gelman's convergence criteria between chain i, and i-1");
3656 print_parm_comment(&bufsize, buffer, allocbufsize, " PAIRS reports all replicate pairs");
3657 print_parm_comment(&bufsize, buffer, allocbufsize, " SUM reports only mean and maxima");
3658 print_parm_mutable(&bufsize, buffer, allocbufsize, "gelman-convergence=%s", options->gelman ? (options->gelmanpairs ? "Yes:Pairs" : "Yes:Sum" ) : "No");
3659 print_parm_br(&bufsize, buffer, allocbufsize);
3660 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax: replicate=< NO | YES:<VALUE | LastChains> >");
3661 print_parm_comment(&bufsize, buffer, allocbufsize, " NO no replication of run");
3662 print_parm_comment(&bufsize, buffer, allocbufsize, " YES replicate run");
3663 print_parm_comment(&bufsize, buffer, allocbufsize, " VALUE number between 2 and many, complete replicates");
3664 print_parm_comment(&bufsize, buffer, allocbufsize, " LastChains replications over last chains");
3665 if(!options->replicate)
3666 {
3667 print_parm(&bufsize, buffer, allocbufsize, "replicate=NO");
3668 }
3669 else
3670 {
3671 if(options->replicatenum==0)
3672 print_parm(&bufsize, buffer, allocbufsize, "replicate=YES:LastChains");
3673 else
3674 print_parm_mutable(&bufsize, buffer, allocbufsize, "replicate=YES:%li", options->replicatenum);
3675 }
3676 print_parm_br(&bufsize, buffer, allocbufsize);
3677
3678 /// \todo remove this cpu section
3679 // if (options->cpu > 1)
3680 // {
3681 // sprintf (fp, "# number of cpus available [not really used right now]cpu=%i", options->cpu);
3682 // add_to_buffer(fp,&bufsize,buffer, allocbufsize);
3683 // }
3684
3685 print_parm_comment(&bufsize, buffer, allocbufsize, "Migration rates are attracted to zero (fatal attraction)");
3686 print_parm_comment(&bufsize, buffer, allocbufsize, "Resistance is the lowest migration value for all but the last chain");
3687 print_parm_comment(&bufsize, buffer, allocbufsize, " Syntax resistance=VALUE");
3688 print_parm_comment(&bufsize, buffer, allocbufsize, " VALUE is the lowest migration rate value allowed during all but the last chain");
3689 print_parm_comment(&bufsize, buffer, allocbufsize, " typical values are 0.01 or _lower_ for data with sequences and 0.0001 or _lower_ for other data");
3690 print_parm_mutable(&bufsize, buffer, allocbufsize, "resistance=%f", options->minmigsumstat);
3691 print_parm_br(&bufsize, buffer, allocbufsize);
3692
3693 print_parm_smalldelimiter(&bufsize, buffer, allocbufsize);
3694 print_parm_br(&bufsize, buffer, allocbufsize);
3695 print_parm(&bufsize, buffer, allocbufsize, "end");
3696 return bufsize;
3697 }
3698
3699 /*private functions============================================= */
3700 ///
3701 /// fill the prior information of bayes theta
set_bayes_options(char * value,option_fmt * options)3702 void set_bayes_options(char *value, option_fmt *options)
3703 {
3704 int ptype;
3705 char paramtype[LINESIZE], priortype[LINESIZE];
3706 prior_fmt *prior=NULL;
3707 sscanf(value,"%s%s", paramtype, priortype);
3708 switch(uppercase(paramtype[0]))
3709 {
3710 case 'T'/* THETA*/: prior = options->bayespriortheta; ptype = THETAPRIOR; break;
3711 case 'M'/* MIG */: prior = options->bayespriorm; ptype = MIGPRIOR; break;
3712 case 'R'/* RATE */: prior = options->bayespriorrate; ptype = RATEPRIOR; break;
3713 default: return;
3714 }
3715 # ifdef USE_MYREAL_FLOAT
3716 switch(uppercase(priortype[0]))
3717 {
3718 case 'S'/*slice sampler with uniform prior*/:
3719 keys = sscanf(value,"%s%s%f%f", paramtype, priortype, &prior->min, &prior->max);
3720 options->bayesprior[ptype] = SLICE;
3721 options->slice_sampling[ptype] = TRUE;
3722 prior_consistency(options, ptype, uppercase(priortype[0]));
3723 break;
3724 case 'M'/*multprior */:
3725 keys = sscanf(value,"%s%s%f%f%f", paramtype, priortype, &prior->min, &prior->max, &prior->delta);
3726 options->bayesprior[ptype] = MULTPRIOR;
3727 prior_consistency(options, ptype, uppercase(priortype[0]));
3728 break;
3729 case 'E'/*expprior */:
3730 keys = sscanf(value,"%s%s%f%f%f", paramtype, priortype, &prior->min, &prior->mean, &prior->max);
3731 options->bayesprior[ptype] = EXPPRIOR;
3732 prior_consistency(options, ptype, uppercase(priortype[0]));
3733 break;
3734 case 'W'/*wexpprior */:
3735 keys = sscanf(value,"%s%s%f%f%f%f", paramtype, priortype, &prior->min, &prior->mean, &prior->max, &prior->delta);
3736 options->bayesprior[ptype] = WEXPPRIOR;
3737 prior_consistency(options, ptype, uppercase(priortype[0]));
3738 break;
3739 case 'G'/*gammaprior */:
3740 keys = sscanf(value,"%s%s%f%f%f%f", paramtype, priortype, &prior->min, &prior->mean, &prior->max, &prior->alpha);
3741 options->bayesprior[ptype] = GAMMAPRIOR;
3742 prior->delta = (prior->max + prior->min) / 10.;
3743 prior_consistency(options, ptype, uppercase(priortype[0]));
3744 break;
3745 case 'U'/*uniformprior*/:
3746 keys = sscanf(value,"%s%s%f%f%f", paramtype, priortype, &prior->min, &prior->max,&prior->delta);
3747 options->bayesprior[ptype] = UNIFORMPRIOR;
3748 prior_consistency(options, ptype, uppercase(priortype[0]));
3749 prior->mean = (prior->max + prior->min) / 2.;
3750 //prior->delta = (prior->max - prior->min) / 10.;
3751 break;
3752 }
3753 #else
3754 switch(uppercase(priortype[0]))
3755 {
3756 case 'S'/* slice sampler with uniform prior */:
3757 sscanf(value,"%s%s%lf%lf%lf", paramtype, priortype, &prior->min, &prior->max,&prior->delta);
3758 options->bayesprior[ptype] = SLICE;
3759 options->slice_sampling[ptype] = TRUE;
3760 prior->mean = (prior->max + prior->min) / 2.;
3761 prior_consistency(options, ptype, uppercase(priortype[0]));
3762 //prior->delta = (prior->max - prior->min) / 10.;
3763 break;
3764 case 'M'/*multprior */:
3765 sscanf(value,"%s%s%lf%lf%lf", paramtype, priortype, &prior->min, &prior->max, &prior->delta);
3766 options->bayesprior[ptype] = MULTPRIOR;
3767 prior_consistency(options, ptype, uppercase(priortype[0]));
3768 break;
3769 case 'E'/*expprior */: ;
3770 sscanf(value,"%s%s%lf%lf%lf", paramtype, priortype, &prior->min, &prior->mean, &prior->max);
3771 options->bayesprior[ptype] = EXPPRIOR;
3772 prior_consistency(options, ptype, uppercase(priortype[0]));
3773 break;
3774 case 'W'/*wexpprior */:
3775 sscanf(value,"%s%s%lf%lf%lf%lf", paramtype, priortype, &prior->min, &prior->mean, &prior->max, &prior->delta);
3776 options->bayesprior[ptype] = WEXPPRIOR;
3777 prior_consistency(options, ptype, uppercase(priortype[0]));
3778 break;
3779 case 'G'/*gammaprior */:
3780 sscanf(value,"%s%s%lf%lf%lf%lf", paramtype, priortype, &prior->min, &prior->mean, &prior->max, &prior->alpha);
3781 options->bayesprior[ptype] = GAMMAPRIOR;
3782 prior_consistency(options, ptype, uppercase(priortype[0]));
3783 break;
3784 case 'U'/*uniformprior*/:
3785 sscanf(value,"%s%s%lf%lf%lf", paramtype, priortype, &prior->min, &prior->max,&prior->delta);
3786 options->bayesprior[ptype] = UNIFORMPRIOR;
3787 prior_consistency(options, ptype, uppercase(priortype[0]));
3788 prior->mean = (prior->max + prior->min) / 2.;
3789 // prior->delta = (prior->max - prior->min) / 10.;
3790 break;
3791 }
3792 #endif
3793 }
3794
3795
3796 ///
3797 /// checks prior consistency
prior_consistency(option_fmt * options,int paramtype,char priortype)3798 void prior_consistency(option_fmt *options, int paramtype, char priortype)
3799 {
3800 MYREAL tmp=0.0;
3801 MYREAL minim = 0.;
3802 MYREAL maxim = HUGE;
3803 prior_fmt *prior=NULL;
3804 switch(paramtype)
3805 {
3806 case THETAPRIOR:
3807 minim = SMALLEST_THETA;
3808 maxim = BIGGEST_THETA;
3809 prior = options->bayespriortheta;
3810 break;
3811 case MIGPRIOR:
3812 minim = SMALLEST_MIGRATION;
3813 maxim = BIGGEST_MIGRATION;
3814 prior = options->bayespriorm;
3815 break;
3816 case RATEPRIOR:
3817 minim = SMALLEST_RATE;
3818 maxim = BIGGEST_RATE;
3819 prior = options->bayespriorrate;
3820 break;
3821 }
3822 if(prior->min < minim)
3823 prior->min = minim;
3824 if(prior->max > maxim)
3825 prior->max = maxim;
3826 if(prior->min > prior->max)
3827 {
3828 warning("maximum of prior is smaller (%f) than minimum (%f)! correct this problem!",prior->max, prior->min);
3829 tmp = prior->max;
3830 prior->max = prior->min;
3831 prior->min = tmp;
3832 }
3833 if( 2.0 * prior->delta > prior->max-prior->min)
3834 {
3835 warning("Your proposal window (delta for the prior) covers more than 50% of the prior range, reduced to 25%\n");
3836 prior->delta = 0.25 * (prior->max-prior->min);
3837 }
3838 switch(priortype)
3839 {
3840 case 'S': break;
3841 case 'M': break;
3842 case 'E': break;
3843 case 'W': break;
3844 case 'G': break;
3845 if(prior->alpha==0.0)
3846 {
3847 prior->alpha=1.0;
3848 warning("Gamma Prior: shape parameter alpha was reset from zero to 1.0! Fix your parmfile or menu_input\n");
3849 }
3850 case 'U': break;
3851 }
3852
3853 }
3854
3855 ///
3856 /// fills the prior information into the option structure old version for compatibility
set_bayes_options_oldstyle(char * value,option_fmt * options,int priortype)3857 void set_bayes_options_oldstyle(char *value, option_fmt *options, int priortype)
3858 {
3859 long keys = 0;
3860 char *tmp;
3861 char *keeptmp;
3862 char priort[LINESIZE], part[LINESIZE], parm[LINESIZE];
3863
3864 error("please change the bayes prior description, there is a more modern way to do that!");
3865
3866 tmp = (char *) mycalloc (LINESIZE, sizeof (char));
3867 keeptmp = tmp;
3868
3869 options->bayesprior[THETAPRIOR] = priortype;
3870 options->bayesprior[MIGPRIOR] = priortype;
3871 options->bayesprior[RATEPRIOR] = priortype;
3872 switch(priortype)
3873 {
3874 case SLICE:
3875 # ifdef USE_MYREAL_FLOAT
3876 keys = sscanf(value,"%s%s%f%f%s%f%f", priort, part, &options->bayespriortheta->min, &options->bayespriortheta->max,
3877 parm, &options->bayespriorm->min, &options->bayespriorm->max);
3878 #else
3879 keys = sscanf(value,"%s%s%lf%lf%s%lf%lf", priort, part, &options->bayespriortheta->min, &options->bayespriortheta->max,
3880 parm, &options->bayespriorm->min, &options->bayespriorm->max);
3881 #endif
3882 if(keys!=7)
3883 error("abort due to problem with bayes settings in parmfile: error in SLICE sampler specification\n");
3884 prior_consistency(options, THETAPRIOR,'S');
3885 prior_consistency(options, MIGPRIOR,'S');
3886 options->bayespriortheta->mean = (options->bayespriortheta->max + options->bayespriortheta->min)/2.;
3887 options->bayespriortheta->delta = (options->bayespriortheta->max - options->bayespriortheta->min)/10.;//about a tenth of the possible range
3888 options->bayespriorm->mean = (options->bayespriorm->max + options->bayespriorm->min)/2.;
3889 options->bayespriorm->delta = (options->bayespriorm->max - options->bayespriorm->min)/10.;
3890 break;
3891 case UNIFORMPRIOR:
3892 # ifdef USE_MYREAL_FLOAT
3893 keys = sscanf(value,"%s%s%f%f%s%f%f", priort, part, &options->bayespriortheta->min, &options->bayespriortheta->max,
3894 parm, &options->bayespriorm->min, &options->bayespriorm->max);
3895 #else
3896 keys = sscanf(value,"%s%s%lf%lf%s%lf%lf", priort, part, &options->bayespriortheta->min, &options->bayespriortheta->max,
3897 parm, &options->bayespriorm->min, &options->bayespriorm->max);
3898 #endif
3899 if(keys!=7)
3900 error("abort due to problem with bayes settings in parmfile: error in UNIFORM prior specification\n");
3901 prior_consistency(options, THETAPRIOR, 'U');
3902 prior_consistency(options, MIGPRIOR, 'U');
3903
3904 options->bayespriortheta->mean = (options->bayespriortheta->max + options->bayespriortheta->min)/2.;
3905 options->bayespriortheta->delta = (options->bayespriortheta->max - options->bayespriortheta->min)/10.;//about a tenth of the possible range
3906 options->bayespriorm->mean = (options->bayespriorm->max + options->bayespriorm->min)/2.;
3907 options->bayespriorm->delta = (options->bayespriorm->max - options->bayespriorm->min)/10.;
3908 break;
3909 case EXPPRIOR:
3910 # ifdef USE_MYREAL_FLOAT
3911 keys = sscanf(value,"%s%s%f%f%f%s%f%f%f", priort, part,
3912 &options->bayespriortheta->min,
3913 &options->bayespriortheta->mean,
3914 &options->bayespriortheta->max,
3915 parm, &options->bayespriorm->min,
3916 &options->bayespriorm->mean,
3917 &options->bayespriorm->max);
3918 #else
3919 keys = sscanf(value,"%s%s%lf%lf%lf%s%lf%lf%lf", priort, part,
3920 &options->bayespriortheta->min,
3921 &options->bayespriortheta->mean,
3922 &options->bayespriortheta->max,
3923 parm, &options->bayespriorm->min,
3924 &options->bayespriorm->mean,
3925 &options->bayespriorm->max);
3926 #endif
3927 if(keys!=9)
3928 error("abort due to problem with bayes settings in parmfile: error in EXP prior specification\n");
3929 prior_consistency(options, THETAPRIOR, 'E');
3930 prior_consistency(options, MIGPRIOR, 'E');
3931
3932 options->bayespriortheta->delta = (options->bayespriortheta->max - options->bayespriortheta->min)/10.;//about a tenth of the possible range
3933 options->bayespriorm->delta = (options->bayespriorm->max - options->bayespriorm->min)/10.;
3934 break;
3935 case WEXPPRIOR:
3936 # ifdef USE_MYREAL_FLOAT
3937 keys = sscanf(value,"%s%s%f%f%f%f%s%f%f%f%f", priort, part,
3938 &options->bayespriortheta->min,
3939 &options->bayespriortheta->mean,
3940 &options->bayespriortheta->max,
3941 &options->bayespriortheta->delta,
3942 parm, &options->bayespriorm->min,
3943 &options->bayespriorm->mean,
3944 &options->bayespriorm->max,
3945 &options->bayespriorm->delta);
3946 #else
3947 keys = sscanf(value,"%s%s%lf%lf%lf%lf%s%lf%lf%lf%lf", priort, part,
3948 &options->bayespriortheta->min,
3949 &options->bayespriortheta->mean,
3950 &options->bayespriortheta->max,
3951 &options->bayespriortheta->delta,
3952 parm, &options->bayespriorm->min,
3953 &options->bayespriorm->mean,
3954 &options->bayespriorm->max,
3955 &options->bayespriorm->delta);
3956 #endif
3957 if(keys!=11)
3958 error("abort due to problem with bayes settings in parmfile: error in WEXP prior specification\n");
3959 prior_consistency(options, THETAPRIOR, 'W');
3960 prior_consistency(options, MIGPRIOR, 'W');
3961
3962 break;
3963 case MULTPRIOR:
3964 # ifdef USE_MYREAL_FLOAT
3965 keys = sscanf(value,"%s%s%f%f%f%f%s%f%f%f%f", priort, part,
3966 &options->bayespriortheta->min,
3967 &options->bayespriortheta->mean,
3968 &options->bayespriortheta->max,
3969 &options->bayespriortheta->delta,
3970 parm, &options->bayespriorm->min,
3971 &options->bayespriorm->mean,
3972 &options->bayespriorm->max,
3973 &options->bayespriorm->delta);
3974 #else
3975 keys = sscanf(value,"%s%s%lf%lf%lf%s%lf%lf%lf", priort, part,
3976 &options->bayespriortheta->min,
3977 &options->bayespriortheta->max,
3978 &options->bayespriortheta->delta,
3979 parm, &options->bayespriorm->min,
3980 &options->bayespriorm->max,
3981 &options->bayespriorm->delta);
3982 #endif
3983 if(keys!=9)
3984 error("abort due to problem with bayes settings in parmfile: error in MULT prior specification\n");
3985 prior_consistency(options, THETAPRIOR, 'M');
3986 prior_consistency(options, MIGPRIOR, 'M');
3987
3988 break;
3989 }
3990 myfree(keeptmp);
3991 }
3992 ///
3993 /// checks all elements in the options file (parmfile) that involve boolean yes no type options
3994 /// typical format is option=<NO | YES<: filename or parameters etc>
3995 long
boolcheck(char ch)3996 boolcheck (char ch)
3997 {
3998 char c = uppercase (ch);
3999 if ((c == 'F') || (c == 'N'))
4000 return 0;
4001 else if ((c == 'T') || (c == 'Y'))
4002 return 1;
4003 else
4004 return -1;
4005 } /* boolcheck */
4006
4007 boolean
booleancheck(option_fmt * options,char * var,char * value)4008 booleancheck (option_fmt * options, char *var, char *value)
4009 {
4010 long i, check;
4011 char *booltokens[NUMBOOL] = BOOLTOKENS;
4012 char *tmp;
4013 long ltemp;
4014 char *extension;
4015
4016 check = boolcheck (value[0]);
4017 if (check == -1)
4018 {
4019 return FALSE;
4020 }
4021 i = 0;
4022 while (i < NUMBOOL && strcmp (var, booltokens[i]))
4023 i++;
4024 switch ((short) i)
4025 {
4026 case 0: /*menu = <yes | no> */
4027 options->menu = (boolean) (check);
4028 break;
4029 case 1: /*interleaved =<yes | no> */
4030 options->interleaved = (boolean) (check);
4031 break;
4032 case 2: /*print-data = <yes | no> */
4033 options->printdata = (boolean) (check);
4034 break;
4035 case 3: /* mixplot=<yes:mixfilename | no> */
4036 options->mixplot = set_filename(value, "YES", &options->mixfilename);
4037 break;
4038 case 4: /* moving-steps = <yes | no> */
4039 options->movingsteps = (boolean) (check);
4040 if (options->movingsteps)
4041 {
4042 strtok (value, ":");
4043 tmp = strtok (NULL, " ,\n");
4044 if (tmp != NULL)
4045 options->acceptfreq = atof ((char *) tmp);
4046 else
4047 options->acceptfreq = 0.1;
4048 }
4049 break;
4050 case 5: /* freqs-from-data = <yes | no> */
4051 options->freqsfrom = (boolean) (check);
4052 if (!options->freqsfrom)
4053 {
4054 strtok (value, ":");
4055 tmp = strtok (NULL, " ,");
4056 if (tmp != NULL)
4057 options->freqa = atof ((char *) tmp);
4058 tmp = strtok (NULL, " ,");
4059 if (tmp != NULL)
4060 options->freqc = atof ((char *) tmp);
4061 tmp = strtok (NULL, " ,");
4062 if (tmp != NULL)
4063 options->freqg = atof ((char *) tmp);
4064 tmp = strtok (NULL, "[,\n");
4065 if (tmp != NULL)
4066 {
4067 options->freqt = atof ((char *) tmp);
4068 // if we encounter a bracket here then we asssume the user
4069 // gave us a number of total sites examined S, this will then
4070 // be used to calculate weight for the invariant sites,
4071 // freq * (S-N)
4072 tmp = strtok (NULL, "\n");
4073 if (tmp != NULL)
4074 {
4075 options->totalsites = atol ((char *) tmp);
4076 }
4077 }
4078 }
4079 break;
4080 case 6: /* useroldtree = < NO | YES OBSOLETE use usetree*/
4081 options->usertree = set_filename(value, "YES", &options->utreefilename); //USERTREE
4082 break;
4083 case 7:
4084 { /* autocorrelation=<YES:value | NO> */
4085 options->autocorr = (boolean) (check);
4086 if (options->autocorr)
4087 {
4088 strtok (value, ":");
4089 tmp = strtok (NULL, " ;\n");
4090 if (tmp != NULL)
4091 options->lambda = 1.0 / atof ((char *) tmp);
4092 }
4093 break;
4094 }
4095 case 8: /* simulation = <yes | no> */
4096 options->simulation = (boolean) check;
4097 break;
4098 case 9: /* plot = <not | yes:<outfile | both><:std|:log><:{xs,xe,ys,ye}<:N|:M><#of_intervals>>> */
4099 options->plot = (boolean) check;
4100 if (options->plot)
4101 {
4102 strtok (value, ":");
4103 if (uppercase (value[0]) == 'Y' || uppercase (value[0]) == 'Y')
4104 {
4105 tmp = strtok (NULL, ":;\n");
4106 if (tmp == NULL)
4107 options->plotmethod = PLOTALL;
4108 else
4109 {
4110 switch (lowercase (tmp[0]))
4111 {
4112 case 'o':
4113 options->plotmethod = PLOTOUTFILE;
4114 break;
4115 case 'b':
4116 options->plotmethod = PLOTALL;
4117 break;
4118 default:
4119 options->plotmethod = PLOTALL;
4120 break;
4121 }
4122 tmp = strtok (NULL, ":;\n");
4123 if (tmp != NULL)
4124 {
4125 switch (lowercase (tmp[0]))
4126 {
4127 case 'l':
4128 options->plotscale = PLOTSCALELOG;
4129 break;
4130 case 's':
4131 options->plotscale = PLOTSCALESTD;
4132 break;
4133 default:
4134 options->plotscale = PLOTSCALELOG;
4135 }
4136 tmp = strtok (NULL, ":;\n");
4137 #ifdef USE_MYREAL_FLOAT
4138 if (4 !=
4139 sscanf (tmp, "{%f,%f,%f,%f",
4140 &options->plotrange[0],
4141 &options->plotrange[1],
4142 &options->plotrange[2],
4143 &options->plotrange[3]))
4144 sscanf (tmp, "{%f%f%f%f", &options->plotrange[0],
4145 &options->plotrange[1],
4146 &options->plotrange[2],
4147 &options->plotrange[3]);
4148 #else
4149 if (4 !=
4150 sscanf (tmp, "{%lf,%lf,%lf,%lf",
4151 &options->plotrange[0],
4152 &options->plotrange[1],
4153 &options->plotrange[2],
4154 &options->plotrange[3]))
4155 sscanf (tmp, "{%lf%lf%lf%lf", &options->plotrange[0],
4156 &options->plotrange[1],
4157 &options->plotrange[2],
4158 &options->plotrange[3]);
4159 #endif
4160 tmp = strtok (NULL, ":;\n");
4161 if (tmp != NULL)
4162 {
4163 switch (lowercase (tmp[0]))
4164 {
4165 case 'm':
4166 options->plotvar = 1;
4167 while (!isdigit (*tmp) && *tmp != '\0')
4168 tmp++;
4169 if ((ltemp =
4170 strtol (tmp, (char **) NULL, 10)) > 0)
4171 options->plotintervals = ltemp;
4172 break;
4173 case 'n':
4174 default:
4175 options->plotvar = PLOT4NM;
4176 while (!isdigit (*tmp) && *tmp != '\0')
4177 tmp++;
4178 if ((ltemp =
4179 strtol (tmp, (char **) NULL, 10)) > 0)
4180 options->plotintervals = ltemp;
4181
4182 break;
4183 }
4184 }
4185 }
4186 }
4187 }
4188 }
4189 break;
4190 case 10: /* weights = <yes | no> */
4191 options->weights = (boolean) check;
4192 set_filename(value, "YES", &options->weightfilename);
4193 break;
4194 case 11: /* read-summary <yes | no> */
4195 options->readsum = (boolean) check;
4196 options->datatype = 'g';
4197 break;
4198 case 12: /* write-summary = <yes | no> */
4199 options->writesum = (boolean) check;
4200 set_filename(value, "YES", &options->sumfilename);
4201 break;
4202 case 13: /* NODATA = <yes | no > */
4203 options->prioralone = (boolean) check;
4204 break;
4205 case 14: /* include-unknown=<yes | no> */
4206 options->include_unknown = (boolean) check;
4207 break;
4208 // old case 14(heating) moved to numbercheck
4209 case 15: /* print-fst = <yes | no> */
4210 options->printfst = (boolean) check;
4211 break;
4212 case 16: /* distfile = <yes | no> */
4213 warning("OBSOLETE option distance=YES|NO ===> use usertree=DISTANCE:distfile or usertree=NO\n");
4214 options->dist = (boolean) check;
4215 break;
4216 case 17: /* geofile = <yes:geofile | no> */
4217 options->geo = (boolean) check;
4218 set_filename(value, "YES", &options->geofilename);
4219 break;
4220 case 18: /* gelman-convergence = <yes:pairs|sum | no> */
4221 options->gelman = (boolean) check;
4222 options->gelmanpairs = FALSE;
4223 if (options->gelman)
4224 {
4225 get_next_word(&value, ":\n", &tmp);
4226 if (value != NULL)
4227 {
4228 if (uppercase (value[0]) == 'P')
4229 {
4230 options->gelmanpairs = TRUE;
4231 }
4232 }
4233 }
4234 break;
4235 case 19: /* randomtree start = <yes | no> */
4236 warning("OBSOLETE option randomtree=YES|NO ===> use usertree=RANDOM or usertree=NO\n");
4237 options->randomtree = (boolean) check;
4238 break;
4239 case 20: /* fast-likelihood calculator = <yes | no> */
4240 options->fastlike = (boolean) check;
4241 break;
4242 case 21: /*aic-modeltest=yes/no */
4243 options->aic = (boolean) check;
4244 if (options->aic)
4245 {
4246 strtok (value, ":\n");
4247 tmp = strtok (NULL, ":");
4248 if (tmp != NULL)
4249 {
4250 if (uppercase (tmp[0]) == 'F')
4251 {
4252 options->fast_aic = TRUE;
4253 tmp = strtok (NULL, ":");
4254 if (tmp != NULL)
4255 {
4256 options->aicmod = atof (tmp);
4257 }
4258 else
4259 options->aicmod = 2.;
4260 }
4261 else
4262 options->aicmod = atof (tmp);
4263 }
4264 else
4265 options->aicmod = 2.;
4266 }
4267 else
4268 options->aicmod = 2.;
4269 break;
4270 case 22: //use-M=true/false
4271 options->usem = (boolean) check;
4272 set_usem_related (options);
4273 break;
4274 case 23: //alternative to above is use-4Nm=false/true
4275 options->usem = !(boolean) check;
4276 set_usem_related (options);
4277 break;
4278 case 24: //bayes-update
4279 options->bayes_infer = (boolean) check;
4280 if(options->bayes_infer)
4281 options->lchains = 1;
4282 break;
4283 case 25: // bayes-allfile: write all bayes parameter estimates to a file
4284 options->has_bayesmdimfile = (boolean) check;
4285 if(options->has_bayesmdimfile)
4286 {
4287 get_next_word(&value,":",&tmp);// extract the word YES, value is now interval:filename
4288 options->mdimdelete=FALSE;
4289 if (!strcmp(tmp,"TEMP"))
4290 {
4291 options->mdimdelete=TRUE;
4292 }
4293 get_next_word(&value,":",&tmp);// extract the interval, value is now filename
4294 if(value==NULL)
4295 {
4296 strncpy (options->bayesmdimfilename, tmp, 255);
4297 }
4298 else
4299 {
4300 options->bayesmdiminterval = atol(tmp);
4301 strncpy (options->bayesmdimfilename, value, 255);
4302 }
4303 unpad(options->bayesmdimfilename," ");
4304 extension = strrchr(options->bayesmdimfilename,'.');
4305 #ifdef ZNZ
4306 if(extension!=NULL && !strncmp(extension,".gz",3))
4307 {
4308 options->use_compressed = 1;
4309 }
4310 else
4311 {
4312 options->use_compressed = 0;
4313 }
4314 #else
4315 options->use_compressed = 0;
4316 #endif
4317 }
4318 break;
4319 case 26: // bayes-file: write marginal parameter posterior distribution
4320 options->has_bayesfile = (boolean) check;
4321 set_filename(value, "YES", &options->bayesfilename);
4322 break;
4323 #ifdef NEWVERSION /* divfile datefile */
4324 case 27: /* divfile = <yes:divfile | no> */
4325 options->div = (boolean) check;
4326 set_filename(value, "YES", &options->divfilename);
4327 break;
4328 #endif
4329 case 28: /* tipdate-file: <yes:datefile | no > */
4330 options->has_datefile = (boolean) check;
4331 set_filename(value, "YES", &options->datefilename);
4332 break;
4333 case 29: /* heated-swap: <yes | no > */
4334 options->heatedswap_off = !((boolean) check);
4335 break;
4336 case 30: // auto-tune=<YES:acceptance-ratio | NO >
4337 options->has_autotune = (boolean) check;
4338 if(options->has_autotune)
4339 {
4340 get_next_word(&value,":",&tmp);// extract the word YES, value is now interval:filename
4341 get_next_word(&value,":",&tmp);// extract the interval, value is now filename
4342 if(value==NULL)
4343 {
4344 options->autotune = 0.44;// using a value suggested by Roberts and Rosenthal (2009)
4345 }
4346 else
4347 {
4348 options->autotune = atof(value);
4349 }
4350 }
4351 break;
4352 default:
4353 return FALSE;
4354 }
4355 return TRUE;
4356 } /* booleancheck */
4357
4358 void
set_usem_related(option_fmt * options)4359 set_usem_related (option_fmt * options)
4360 {
4361 if (options->usem)
4362 {
4363 options->plotvar = PLOTM;
4364 options->migvar = PLOTM;
4365 options->profileparamtype = PLOTM;
4366
4367 }
4368 else
4369 {
4370 options->plotvar = PLOT4NM;
4371 options->migvar = PLOT4NM;
4372 options->profileparamtype = PLOT4NM;
4373 }
4374
4375 }
4376
4377 boolean
numbercheck(option_fmt * options,char * var,char * value)4378 numbercheck (option_fmt * options, char *var, char *value)
4379 {
4380 int retval;
4381 MYREAL musum = 0., lastrate = 1.;
4382 long i = 0, z, cc = 0;
4383 char *tmp, *temp, *temp2, *keeptmp;
4384 char *numbertokens[NUMNUMBER] = NUMBERTOKENS ;
4385 tmp = (char *) mycalloc (LINESIZE, sizeof (char));
4386 keeptmp = tmp;
4387 while (i < NUMNUMBER && strcmp (var, numbertokens[i]))
4388 i++;
4389 switch ((short) i)
4390 {
4391 case 0: /*ttratio = value */
4392 z = 0;
4393 temp = strtok (value, " ,;\n\0");
4394 while (temp != NULL)
4395 {
4396 options->ttratio[z++] = atof (temp);
4397 options->ttratio =
4398 (MYREAL *) myrealloc (options->ttratio, sizeof (MYREAL) * (z + 1));
4399 options->ttratio[z] = 0.0;
4400
4401 temp = strtok (NULL, " ,;\n\0");
4402 }
4403 break;
4404 case 1: /*short-chains = value */
4405 options->schains = atol (value);
4406 break;
4407 case 2: /*short-steps = value */
4408 case 32: /* short-sample = value */
4409 options->ssteps = atol (value);
4410 break;
4411 case 3: /*short-increment = value */
4412 options->sincrement = atol (value);
4413 break;
4414 case 4: /*long-chains = value */
4415 options->lchains = atol (value);
4416 break;
4417 case 5: /*long-steps = value */
4418 case 33: /*long-sample = value */
4419 options->lsteps = atol (value);
4420 break;
4421 case 6: /*long-increment = value */
4422 options->lincrement = atol (value);
4423 break;
4424 case 7:
4425 break; /* theta: already handled in read_theta() */
4426 case 8: /*nmlength = value */
4427 options->nmlength = atoi (value); //strtol (value, (char **) NULL, 10);
4428 break;
4429 case 9: /* seed = <Auto | seedfile | Own:value> */
4430 switch (value[0])
4431 {
4432 case 'A':
4433 case 'a':
4434 case '0':
4435 options->autoseed = AUTO;
4436 options->inseed = (long) time (0) / 4 + 1;
4437 break;
4438 case 'S':
4439 case 's':
4440 case '1':
4441 options->autoseed = NOAUTO;
4442 openfile(&options->seedfile, options->seedfilename, "r", NULL);
4443 if (options->seedfile == NULL)
4444 {
4445 usererror ("cannot find seedfile\n");
4446 }
4447 retval = fscanf (options->seedfile, "%ld%*[^\n]", &options->inseed);
4448 fclose (options->seedfile);
4449 break;
4450 case 'O':
4451 case 'o':
4452 case '2':
4453 options->autoseed = NOAUTOSELF;
4454 strtok (value, ":");
4455 tmp = strtok (NULL, " ;\n");
4456 if (tmp != NULL)
4457 options->inseed = atol ((char *) tmp);
4458 if (options->inseed > 0)
4459 break;
4460 default:
4461 options->autoseed = AUTO;
4462 options->inseed = (long) time (0) / 4 + 1;
4463 usererror ("Failure to read seed method, should be\n \
4464 random-seed=auto or random-seed=seedfile or random-seed=own:value\nwhere value is a positive integer\nUsing AUTOMATIC seed=%li\n", options->inseed);
4465 break;
4466 }
4467 break;
4468 case 10:
4469 break; /*"migration" fake: this is already handled in read_migrate */
4470 case 11: /*mutation= <auto=gamma | nogamma | constant | estimate | own | data> */
4471 switch (value[0])
4472 {
4473 case 'A': /*automatic */
4474 case 'a':
4475 case 'E': /*automatic */
4476 case 'e':
4477 options->bayesmurates=TRUE;
4478 options->murates = FALSE;
4479 options->murates_fromdata = FALSE;
4480 options->gamma = FALSE;
4481 break;
4482 case 'G':
4483 case 'g':
4484 options->gamma = TRUE;
4485 options->bayesmurates=TRUE;
4486 options->murates = FALSE;
4487 options->murates_fromdata = FALSE;
4488 (void) strtok (value, " :");
4489 temp = strtok (NULL, " ,;\n");
4490 if (temp != NULL)
4491 options->alphavalue = atof (temp);
4492 else
4493 options->alphavalue = START_ALPHA;
4494 break;
4495 case 'O':
4496 case 'o':
4497 options->murates = TRUE;
4498 options->gamma = FALSE;
4499 options->murates_fromdata = FALSE;
4500 options->bayesmurates = FALSE;
4501 (void) strtok (value, " :");
4502 temp = strtok (NULL, " ,;\n");
4503 if (temp != NULL)
4504 {
4505 options->muloci = atol (temp);
4506 options->mu_rates =
4507 (MYREAL *) mycalloc (options->muloci, sizeof (MYREAL));
4508 musum = 0.;
4509 for (i = 0; i < options->muloci; i++)
4510 {
4511 temp = strtok (NULL, " ,;\n");
4512 if (temp == NULL)
4513 {
4514 while (i < options->muloci)
4515 {
4516 options->mu_rates[i] = lastrate;
4517 musum += options->mu_rates[i];
4518 i++;
4519 }
4520 }
4521 lastrate = options->mu_rates[i] = atof (temp);
4522 musum += options->mu_rates[i];
4523 }
4524 // mean must be 1.
4525 musum /= options->muloci;
4526 for (i = 0; i < options->muloci; i++)
4527 {
4528 options->mu_rates[i] /= musum;
4529 }
4530 }
4531 break;
4532 case 'D': /* mutaton is varying and calculated from the data*/
4533 case 'd':
4534 options->gamma = FALSE;
4535 options->murates = TRUE;
4536 options->murates_fromdata = TRUE;
4537 options->bayesmurates = FALSE;
4538 break;
4539 case 'N': /*nogamma, none, all loci have same mu */
4540 case 'n':
4541 case 'C':
4542 case 'c':
4543 default:
4544 options->murates = FALSE;
4545 options->gamma = FALSE;
4546 options->murates_fromdata = FALSE;
4547 options->bayesmurates = FALSE;
4548 break;
4549 }
4550 break;
4551 case 12: /*datatype=<allele|microsatellite|brownian|sequence|f-ancestral states|genealogies> */
4552 switch (value[0])
4553 {
4554 case 'a':
4555 case 'A':
4556 options->datatype = 'a';
4557 break;
4558 case 'm':
4559 case 'M':
4560 options->datatype = 'm';
4561 break;
4562 case 'b':
4563 case 'B':
4564 options->datatype = 'b';
4565 break;
4566 case 's':
4567 case 'S':
4568 options->datatype = 's';
4569 break;
4570 case 'n':
4571 case 'N':
4572 options->datatype = 'n';
4573 break;
4574 case 'h':
4575 case 'H':
4576 options->datatype = 'h';
4577 options->fastlike=FALSE;
4578 break;
4579 case 'u':
4580 case 'U':
4581 options->datatype = 'u';
4582 break;
4583 case 'f':
4584 case 'F':
4585 options->datatype = 'f';
4586 break;
4587 case 'g':
4588 case 'G':
4589 options->datatype = 'g';
4590 options->readsum = TRUE;
4591 break;
4592 default:
4593 options->datatype = 's';
4594 break;
4595 }
4596 break;
4597 case 13: /* categories=<None | value> */
4598 if (uppercase (value[0] == 'N'))
4599 {
4600 options->categs = ONECATEG;
4601 break;
4602 }
4603 else
4604 {
4605 options->categs = strtol (value, (char **) NULL, 10);
4606 /* needs to read auxilliary file catfile */
4607 sprintf(tmp,"%li",options->categs);
4608 set_filename(value, tmp, &options->catfilename);
4609 }
4610 break;
4611 case 14: /*create rates=value:list of rates */
4612 strncpy (tmp, value, strcspn (value, ":"));
4613 if (strtol (tmp, (char **) NULL, 10) /*;atoi (tmp) */ > 1)
4614 { /* rate categories */
4615 options->rcategs = strtol (tmp, (char **) NULL, 10);
4616 options->rrate =
4617 (MYREAL *) myrealloc (options->rrate,
4618 sizeof (MYREAL) * (options->rcategs + 1));
4619 (void) strtok (value, " :");
4620 temp = strtok (NULL, " ,;\n");
4621 z = 0;
4622 while (temp != NULL)
4623 {
4624 if (z > options->rcategs)
4625 usererror ("check parmfile-option rates, missing rate\n");
4626 options->rrate[z++] = atof (temp);
4627 temp = strtok (NULL, " ,;\n");
4628 }
4629 }
4630 break;
4631 case 15: /* probabilities for each rate category */
4632 strncpy (tmp, value, strcspn (value, ":"));
4633 if (strtol (tmp, (char **) NULL, 10) > 1)
4634 { /* probabilities for each rate category */
4635 options->rcategs = strtol (tmp, (char **) NULL, 10);
4636 options->probcat =
4637 (MYREAL *) myrealloc (options->probcat,
4638 sizeof (MYREAL) * (options->rcategs + 1));
4639 (void) strtok (value, " :");
4640 temp = strtok (NULL, " ,;\n");
4641 z = 0;
4642 while (temp != NULL)
4643 {
4644 if (z > options->rcategs)
4645 usererror
4646 ("check parmfile prob-rates, missing rate probability\n");
4647 options->probcat[z++] = atof (temp);
4648 temp = strtok (NULL, " ,;\n");
4649 }
4650 }
4651 break;
4652 case 16: /*micro-stepmax */
4653 // options->micro_stepnum = strtol (value, (char **) NULL, 10);
4654
4655 break;
4656 case 17: /*micro-threshold */
4657 options->micro_threshold = atol(value);
4658 if(options->micro_threshold % 2 != 0)
4659 options->micro_threshold += 1;
4660 break;
4661 case 18: /*delimiter */
4662 options->dlm = value[0];
4663 break;
4664 case 19: /*burn-in */
4665 options->burn_in = atol (value);
4666 if(strchr(value,'a')!=NULL)
4667 options->burnin_autostop = 'a';
4668 else
4669 {
4670
4671 if(strchr(value,'t')!=NULL)
4672 {
4673 options->burnin_autostop = 't';
4674 }
4675 else
4676
4677 if(strchr(value,'e')!=NULL)
4678 {
4679 options->burnin_autostop = 'e';
4680 }
4681 else
4682 options->burnin_autostop = ' ';
4683 }
4684 break;
4685 case 20: /*infilename */
4686 get_filename (&options->infilename, value);
4687 break;
4688 case 21: /*outfilename */
4689 get_filename(&options->outfilename,value);
4690 break;
4691 case 22: /*mathfilename */
4692 get_filename (&options->mathfilename, value);
4693 break;
4694 case 23: /*title */
4695 strncpy (options->title, value, 80);
4696 break;
4697 case 24: /*long-chain-epsilon */
4698 options->lcepsilon = atof (value);
4699 if (options->lcepsilon <= 0)
4700 options->lcepsilon = LONGCHAINEPSILON;
4701 break;
4702 case 25: /* print tree options */
4703 switch (uppercase (value[0]))
4704 {
4705 case 'N':
4706 options->treeprint = _NONE;
4707 break;
4708 case 'A':
4709 options->treeprint = ALL;
4710 options->treeinc = 1;
4711 set_filename(value, "ALL", &options->treefilename);
4712 break;
4713 case 'B':
4714 options->treeprint = BEST;
4715 set_filename(value, "BEST", &options->treefilename);
4716 break;
4717 case 'L':
4718 options->treeprint = LASTCHAIN;
4719 get_next_word(&value,":",&tmp);// the word LASTCHAIN
4720 get_next_word(&value,":",&tmp);
4721 options->treeinc = atol(tmp);
4722 set_filename(value, "LASTCHAIN", &options->treefilename);
4723 break;
4724 default:
4725 options->treeprint = _NONE;
4726 break;
4727 }
4728 break;
4729 case 26: /* progress: No, Yes, Verbose */
4730 switch (uppercase (value[0]))
4731 {
4732 case 'F':
4733 case 'N':
4734 options->progress = FALSE;
4735 options->verbose = FALSE;
4736 break;
4737 case 'T':
4738 case 'Y':
4739 options->progress = TRUE;
4740 options->verbose = FALSE;
4741 break;
4742 case 'V':
4743 options->progress = TRUE;
4744 options->verbose = TRUE;
4745 break;
4746 }
4747 break;
4748 case 27: /* l-ratio: <NO | YES>:val1,val2,val3,val4,val5 */
4749 cc = options->lratio->counter;
4750 switch (uppercase (value[0]))
4751 {
4752 case 'Y':
4753 options->lratio->data[cc].type = MLE;
4754 break;
4755 case 'N':
4756 default:
4757 myfree(keeptmp);
4758 return FALSE;
4759 }
4760 (void) strtok (value, ":");
4761 temp = strtok (NULL, "\n");
4762 if (temp != NULL)
4763 {
4764 temp2 = strchr (temp, ':');
4765 if (temp2 != NULL)
4766 {
4767 strcpy (options->lratio->data[cc].value2, temp2);
4768 *temp2 = '\0';
4769 strcpy (options->lratio->data[cc].value1, temp);
4770 }
4771 else
4772 strcpy (options->lratio->data[cc].value1, temp);
4773 }
4774 if (cc + 1 == options->lratio->alloccounter)
4775 {
4776 options->lratio->alloccounter += 2;
4777 options->lratio->data =
4778 (lr_data_fmt *) myrealloc (options->lratio->data,
4779 sizeof (lr_data_fmt) *
4780 (options->lratio->alloccounter+1));
4781 for (i = cc + 1; i < options->lratio->alloccounter; i++)
4782 {
4783 options->lratio->data[i].elem = 0;
4784 options->lratio->data[i].value1 =
4785 (char *) mycalloc (1, sizeof (char) * LONGLINESIZE);
4786 options->lratio->data[i].elem = 0;
4787 options->lratio->data[i].value2 =
4788 (char *) mycalloc (1, sizeof (char) * LONGLINESIZE);
4789 options->lratio->data[i].connect =
4790 (char *) mycalloc (1, sizeof (char) * LONGLINESIZE);
4791
4792 }
4793 }
4794 options->lratio->counter += 1;
4795 break;
4796 case 28: /* fst-type: <Theta | Migration> */
4797 switch (uppercase (value[0]))
4798 {
4799 case 'T':
4800 options->fsttype = 'T';
4801 break;
4802 case 'M':
4803 default:
4804 options->fsttype = 'M';
4805 break;
4806 }
4807 fst_type (options->fsttype);
4808 break;
4809 case 29: /*profile=<NO| NONE | YES | ALL | TABLES | SUMMARY>><: <FAST | */
4810 switch (uppercase (value[0]))
4811 {
4812 case 'S':
4813 options->profile = SUMMARY;
4814 break;
4815 case 'Y':
4816 case 'A':
4817 options->profile = ALL;
4818 break;
4819 case 'N':
4820 options->profile = _NONE;
4821 break;
4822 case 'T':
4823 options->profile = TABLES;
4824 break;
4825 default: /*A */
4826 options->profile = ALL;
4827 break;
4828 }
4829 (void) strtok (value, ":;\n");
4830 temp = strtok (NULL, ":;\n");
4831 if (temp != NULL)
4832 {
4833 switch (lowercase (temp[0]))
4834 {
4835 case 'p': /*precise percentiles */
4836 options->profilemethod = 'p';
4837 break;
4838 case 'd': /*discrete steps see at start of file */
4839 options->profilemethod = 'd';
4840 break;
4841 //case 's':
4842 //options->profilemethod = 's';
4843 //break;
4844 case 'x': /* x-rated */
4845 case 'u': /* uncorrelated */
4846 case 'q': /* quick and dirty */
4847 options->profilemethod = 'q';
4848 break;
4849 case 'f': /* quick and exact mixture */
4850 options->profilemethod = 'f';
4851 break;
4852 default:
4853 options->profilemethod = 'f';
4854 options->printprofsummary = TRUE;
4855 break;
4856 }
4857 temp = strtok (NULL, ":;\n");
4858 if (temp != NULL)
4859 {
4860 switch (lowercase (temp[0]))
4861 {
4862 case 'm':
4863 options->profileparamtype = 1;
4864 break;
4865 default:
4866 options->profileparamtype = PLOT4NM;
4867 }
4868 }
4869 }
4870 set_profile_options (options);
4871 break;
4872 case 30: /* custom-migration:<{> migration matrix and theta on
4873 diagonal:
4874 0 means not estimated,
4875 x means estimated, s means symmetrically
4876 estimated, m means all are the same <}> */
4877 if (myID == MASTER)
4878 read_custom_migration (options->parmfile, options, value,
4879 options->numpop);
4880 #ifdef MPI
4881
4882 else
4883 read_custom_migration_worker (options->buffer, options, value,
4884 options->numpop);
4885 #endif
4886
4887 break;
4888 case 31: /*sumfilename */
4889 strcpy (options->sumfilename, value);
4890 break;
4891 /*case 32 and case 33 are fallthroughs to 2 and 3 */
4892 case 34: /*replicate */
4893 switch (uppercase (value[0]))
4894 {
4895 case 'T':
4896 case 'Y':
4897 options->replicate = TRUE;
4898 temp = strtok (value, ":;\n");
4899 if (temp != NULL)
4900 {
4901 temp = strtok (NULL, ":;\n");
4902 if (uppercase (temp[0]) == 'L')
4903 options->replicatenum = 0;
4904 else
4905 options->replicatenum = strtol (temp, (char **) NULL, 10);
4906 }
4907 else
4908 options->replicatenum = 0;
4909 break;
4910 default:
4911 options->replicate = FALSE;
4912 options->replicatenum = 0;
4913 }
4914 break;
4915 case 35: /* cpu number */
4916 options->cpu = (short) ATOI (value);
4917 break;
4918 case 36: /* logfile=<YES:logfile | NO> do we write a logfile or not */
4919 options->writelog = set_filename(value,"YES",&options->logfilename);
4920 break;
4921 case 37: /* sequencing error */
4922 options->seqerror = atof (value);
4923 if (options->seqerror < 0.0)
4924 {
4925 warning
4926 ("Sequencing error was misspecified in parmfile, reset to 0.0\n");
4927 options->seqerror = 0.0;
4928 }
4929 break;
4930 #ifdef UEP
4931
4932 case 38: /* do we have a uep file or not, function returns TRUE if uep=YES:filename*/
4933 options->uep = set_filename(value, "YES", &options->uepfilename);
4934 break;
4935 case 39: /* do we have uep-rates */
4936 temp = strtok (value, ":;\n");
4937 if (temp != NULL)
4938 {
4939 options->uepmu = atof (temp);
4940 temp = strtok (NULL, ":;\n");
4941 if (temp != NULL)
4942 {
4943 options->uepnu = atof (temp);
4944 }
4945 options->ueprate = options->uepmu;
4946 }
4947 break;
4948 case 40: /* do we have uep-bases */
4949 temp = strtok (value, ":; \n");
4950 if (temp != NULL)
4951 {
4952 options->uepfreq1 = atof (temp);
4953 temp = strtok (NULL, " :;\n");
4954 if (temp != NULL)
4955 {
4956 options->uepfreq0 = atof (temp);
4957 }
4958 else
4959 options->uepfreq0=1. - options->uepfreq1;
4960 }
4961 break;
4962 #endif
4963
4964 case 41: // mu-rates??????
4965 break;
4966 case 42: /* heating=<no | <yes | adaptive | bounded>:numintervals:{temperatures}> */
4967 switch (uppercase (value[0]))
4968 {
4969 case 'A': //adaptive heating on
4970 options->heating = 1;
4971 options->adaptiveheat = STANDARD;
4972 break;
4973 case 'B': //adaptive heating on
4974 options->heating = 1;
4975 options->adaptiveheat = BOUNDED;
4976 break;
4977 case 'Y':
4978 case 'P':
4979 options->heating = 1;
4980 options->adaptiveheat = NOTADAPTIVE;
4981 break;
4982 case 'N':
4983 default:
4984 options->heating=0;
4985 options->adaptiveheat=NOTADAPTIVE;
4986 break;
4987 }
4988 if (options->heating == 1)
4989 {
4990 strtok (value, ":\n");
4991 tmp = strtok (NULL, ": ");
4992 if (tmp != NULL)
4993 {
4994 options->heating_interval = atol (tmp);
4995 tmp = strtok (NULL, "{, ");
4996 if (tmp != NULL)
4997 {
4998 z = 0;
4999 while (1)
5000 {
5001 options->heat[z++] = atof (tmp);
5002 tmp = strtok (NULL, ", :\n");
5003 if (tmp == NULL || z >= 1000)
5004 {
5005 if(z>=1000)
5006 warning("Ignored the all temperature above 1000 heated chains\n");
5007 break;
5008 }
5009 }
5010 options->heated_chains = z;
5011 }
5012 }
5013 }
5014 break;
5015 #ifdef LONGSUM
5016
5017 case 43: /* fluctuate=<no | <yes>:{rate_1today, time1today,rate_1middle,time1middle, rate_1past, time1past,rate2_today,time2_today,...}> */
5018 switch (uppercase (value[0]))
5019 {
5020 case 'Y':
5021 case 'P':
5022 options->fluctuate = TRUE;
5023 break;
5024 case 'N':
5025 default:
5026 options->fluctuate = FALSE;
5027 break;
5028 }
5029 if (options->fluctuate)
5030 {
5031 strtok (value, ":\n");
5032 tmp = strtok (NULL, ":{,");
5033 if (tmp != NULL)
5034 {
5035 z = 0;
5036 while (1)
5037 {
5038 options->flucrates = myrealloc(options->flucrates, sizeof(MYREAL) * (z+1));
5039 options->flucrates[z++] = atof (tmp);
5040 tmp = strtok (NULL, ", {}:\n");
5041 if (tmp == NULL)
5042 break;
5043 }
5044 }
5045 }
5046 if(z%3!=0)
5047 error("Fluctuating rates and times need to be 3 rates and 3 times per population");
5048 break;
5049 #endif /*LONGSUM*/
5050 case 44: /*resistance = value [to fatal attraction to zero*/
5051 //xcode z = 0;
5052 temp = strtok (value, " ,;\n\0");
5053 if (temp != NULL)
5054 options->minmigsumstat = atof (temp);
5055 else
5056 options->minmigsumstat = MINMIGSUMSTAT;
5057 break;
5058 case 45: /* bayes-updatefreq=[0..1] */
5059 temp = strtok (value, " ,;\n\0");
5060 if (temp != NULL)
5061 options->updateratio = atof (temp);
5062 else
5063 options->updateratio= HALF;
5064 break;
5065 case 49: /*bayes-posteriorbins*/
5066 temp = strtok (value, " ,;\n\0");
5067 if (temp != NULL)
5068 {
5069 options->bayespriortheta->bins = atol (temp);
5070 temp = strtok (NULL, " ,;\n\0");
5071 if (temp != NULL)
5072 options->bayespriorm->bins = atol (temp);
5073 else
5074 options->bayespriorm->bins = BAYESNUMBIN;
5075 }
5076 else
5077 {
5078 options->bayespriortheta->bins = BAYESNUMBIN;
5079 options->bayespriorm->bins = BAYESNUMBIN;
5080 }
5081 break;
5082
5083 case 46: /*bayesfile=FILENAME set bayesfile and bayes analysis? */
5084 strncpy (options->bayesfilename, value, 255); //this is legacy code and got replaced
5085 // by bayes-file=<YES | NO>:FILENAME
5086 options->has_bayesfile = TRUE;
5087 break;
5088
5089 case 47: /*set bayes-prior*/
5090 switch(toupper(value[0]))
5091 {
5092 //case 'S':
5093 //set_bayes_options_oldstyle(value, options,SLICE);
5094 //break;
5095 case 'M':
5096 set_bayes_options_oldstyle(value, options,MULTPRIOR);
5097 break;
5098 case 'E':
5099 set_bayes_options_oldstyle(value, options,EXPPRIOR);
5100 break;
5101 case 'W':
5102 set_bayes_options_oldstyle(value, options,WEXPPRIOR);
5103 break;
5104 case 'U':
5105 default:
5106 set_bayes_options_oldstyle(value, options,UNIFORMPRIOR);
5107 break;
5108 }
5109 break;
5110 case 48: /* usertree = < NO | UPGMA | AUTOMATIC | TREE:intreefilename | RANDOM | DISTANCE:distfilename */
5111 options->usertree = set_filename(value, "T", &options->utreefilename); //USERTREE
5112 options->randomtree = (value[0] =='R') ? TRUE : FALSE ; //RANDOMTREE
5113 options->dist = set_filename(value, "D", &options->distfilename); //DISTANCE
5114 break;
5115 /* case 49 is further back*/
5116 case 50: /* mig-histogram=<<yes|all>:histogram-binsize:filename | no> */
5117 get_next_word(&value,":",&tmp);
5118 options->mighist_all = ((uppercase(tmp[0]) == 'A') ? TRUE : FALSE);
5119 options->mighist = ((uppercase(tmp[0]) == 'Y') ? TRUE : FALSE);
5120 if(options->mighist_all || options->mighist)
5121 {
5122 options->mighist = TRUE;
5123 get_next_word(&value,":",&tmp);
5124 if(value!=NULL)
5125 {
5126 options->eventbinsize = (float) atof(tmp);
5127 strncpy (options->mighistfilename, value, 255);
5128 }
5129 else
5130 {
5131 strncpy (options->mighistfilename, tmp, 255);
5132 }
5133 }
5134 break;
5135 case 51:
5136 switch(toupper(value[0]))
5137 {
5138 case 'A':
5139 options->bayespretty = PRETTY_MAX;
5140 break;
5141 case 'P':
5142 options->bayespretty = PRETTY_P99;
5143 break;
5144 case 'M':
5145 options->bayespretty = PRETTY_P99MAX;
5146 break;
5147 case 'T':
5148 default:
5149 options->bayespretty = PRETTY_P100;
5150 break;
5151 }
5152 break;
5153 #ifdef PRETTY
5154 case 52: /*PDF-outfile name */
5155 strcpy (options->pdfoutfilename, value);
5156 break;
5157 #endif
5158 case 53: /*bayes interval for writing all parameters to file */
5159 options->bayesmdiminterval = atol(value);
5160 break;
5161 case 54: /*set bayes-priorS*/
5162 // new parmfile variable so that the old stype setting still work
5163 set_bayes_options(value, options);
5164 break;
5165 case 55: /*set skyline-histogram*/
5166 get_next_word(&value,":",&tmp);
5167 options->skyline = tmp[0] == 'Y' ? TRUE : FALSE;
5168 if(options->skyline)
5169 {
5170 get_next_word(&value,":",&tmp);
5171 options->eventbinsize = (float) atof(tmp);
5172 strncpy (options->skylinefilename, value, 255);
5173 }
5174 break;
5175 case 56: /* rates-gamma */
5176 get_next_word(&value,":,; ",&tmp);
5177 options->seqrate_gamma_num = 0;
5178 while(tmp!=NULL)
5179 {
5180 options->seqrate_gamma[options->seqrate_gamma_num++] = atof(tmp);
5181 }
5182 break;
5183 case 57: /* Bayes-proposals*/
5184 get_next_word(&value,":,; ",&tmp);
5185 if(tmp != NULL)
5186 {
5187 switch(uppercase(tmp[0]))
5188 {
5189 case 'T': get_next_word(&value,":,; ",&tmp);
5190 if(tmp != NULL)
5191 {
5192 if(uppercase(tmp[0])=='S')
5193 options->slice_sampling[THETAPRIOR] = TRUE;
5194 else
5195 options->slice_sampling[THETAPRIOR] = FALSE;
5196 }
5197 break;
5198 case 'M': get_next_word(&value,":,; ",&tmp);
5199 if(tmp != NULL)
5200 {
5201 if(uppercase(tmp[0])=='S')
5202 options->slice_sampling[MIGPRIOR] = TRUE;
5203 else
5204 options->slice_sampling[MIGPRIOR] = FALSE;
5205 }
5206 break;
5207 case 'R': get_next_word(&value,":,; ",&tmp);
5208 if(tmp != NULL)
5209 {
5210 if(uppercase(tmp[0])=='S')
5211 options->slice_sampling[RATEPRIOR] = TRUE;
5212 else
5213 options->slice_sampling[RATEPRIOR] = FALSE;
5214 }
5215 break;
5216 }
5217 }
5218 break;
5219 case 58: /*generation-per-year*/
5220 get_next_word(&value,":,; ",&tmp);
5221 options->generation_year = atof(tmp);
5222 break;
5223 case 59: /*mutationrate-per-year absolute mutation rate per year for each locus*/
5224 options->mutationrate_year[0] = 1.0;
5225 get_next_word(&value,"{}:,; ",&tmp);
5226 z=0;
5227 while(tmp != NULL)
5228 {
5229 if(z >= options->mutationrate_year_numalloc)
5230 {
5231 options->mutationrate_year_numalloc = z+1;
5232 options->mutationrate_year = (MYREAL*) myrealloc(options->mutationrate_year, options->mutationrate_year_numalloc * sizeof(MYREAL));
5233 }
5234 options->mutationrate_year[z] = atof(tmp);
5235 printf("%i> mutationrate_year[%li]=%g %s\n",myID,z, options->mutationrate_year[z],tmp);
5236 z++;
5237 get_next_word(&value,"{}:,; ",&tmp);
5238 }
5239 break;
5240 case 60: /*inheritance-scalars defines the scalars so that all reference*/
5241 /* is made to 4 Ne mu when the scalar is used, otherwise the scalar is*/
5242 /* assume to be 1*/
5243 options->inheritance_scalars[0] = 1.0;
5244 get_next_word(&value,"{}:,; ",&tmp);
5245 z=0;
5246 while(tmp != NULL)
5247 {
5248 if(z >= options->inheritance_scalars_numalloc)
5249 {
5250 options->inheritance_scalars_numalloc = z+1;
5251 options->inheritance_scalars = (MYREAL*) myrealloc(options->inheritance_scalars, options->inheritance_scalars_numalloc * sizeof(MYREAL));
5252 }
5253 options->inheritance_scalars[z] = atof(tmp);
5254 //printf("%i> inheritance scalar[%li]=%g %s\n",myID,z, options->inheritance_scalars[z],tmp);
5255 z++;
5256 get_next_word(&value,"{}:,; ",&tmp);
5257 }
5258 break;
5259 case 61: /* micro-submodel */
5260 get_next_word(&value,":",&tmp);
5261 if(tmp==NULL)
5262 options->msat_option = atol(value);
5263 else
5264 {
5265 options->msat_option = atol(tmp);
5266 if(value==NULL)
5267 break;
5268 get_next_word(&value,", ",&tmp);
5269 if(tmp[0]=='{')
5270 options->msat_tuning[0] = atof(tmp+1);
5271 else
5272 options->msat_tuning[0] = atof(tmp);
5273 get_next_word(&value," }",&tmp);
5274 options->msat_tuning[1] = atof(tmp);
5275 }
5276 break;
5277
5278 case 62: /*random-subset*/
5279 get_next_word(&value,":,; ",&tmp);
5280 options->randomsubset = atol(tmp);
5281 get_next_word(&value,":,; ",&tmp);
5282 if (tmp != NULL)
5283 options->randomsubsetseed = atol(tmp);
5284 break;
5285
5286 case 63: /* pooling of populations an array of numbers indicates the pooling*/
5287 set_localities(&value, &tmp, options);
5288 break;
5289 case 64:/* analyze-loci=<All | Variable | Fast > */
5290 if (value!=NULL)
5291 {
5292 switch(value[0])
5293 {
5294 case 'V': // analyze only variable loci
5295 options->onlyvariable = TRUE;
5296 break;
5297 case 'F': // analyze all variable and one invariant weighted to speed up run
5298 options->onlyvariable = FALSE;
5299 options->has_variableandone = TRUE;
5300 break;
5301 case 'A': //analyze all, default
5302 default:
5303 options->onlyvariable = FALSE;
5304 break;
5305 }
5306 }
5307 break;
5308 default:
5309 myfree(keeptmp);
5310 return FALSE;
5311 }
5312 myfree(keeptmp);
5313
5314 return TRUE;
5315 } /* numbercheck */
5316
5317 /*void
5318 reset_oneline (option_fmt * options, long position)
5319 {
5320 fseek (options->parmfile, position, SEEK_SET);
5321 }*/
5322
5323 void
reset_oneline(option_fmt * options,fpos_t * position)5324 reset_oneline (option_fmt * options, fpos_t * position)
5325 {
5326 fsetpos (options->parmfile, position);
5327 }
5328
read_random_theta(option_fmt * options,char ** buffer)5329 void read_random_theta(option_fmt *options, char ** buffer)
5330 {
5331 char *tempstr;
5332 char *tmp;
5333 char *otempstr;
5334 char *otmp;
5335 char *retval;
5336
5337 tempstr = (char *) mycalloc(LINESIZE,sizeof(char));
5338 tmp = (char *) mycalloc(LINESIZE,sizeof(char));
5339 otempstr = tempstr;
5340 otmp = tmp;
5341
5342 options->numthetag = 2;
5343 options->thetag =
5344 (MYREAL *) myrealloc (options->thetag, sizeof (MYREAL) * 3);
5345 if( myID == MASTER)
5346 {
5347 retval = fgets (tmp, LINESIZE, options->parmfile);
5348 }
5349 else
5350 {
5351 sgets (tmp, LINESIZE, buffer);
5352 }
5353 get_next_word(&tmp,",{} \n",&tempstr);
5354 while(strchr("{ ", (tempstr)[0]))
5355 {
5356 get_next_word(&tmp,",{} \n",&tempstr);
5357 }
5358 options->thetag[0] = atof(tempstr);
5359 get_next_word(&tmp,",{} \n",&tempstr);
5360 options->thetag[1] = atof(tempstr);
5361 myfree(otmp);
5362 myfree(otempstr);
5363 }
5364
5365
5366
read_random_mig(option_fmt * options,char ** buffer)5367 void read_random_mig(option_fmt *options, char ** buffer)
5368 {
5369 char *tempstr;
5370 char *tmp;
5371 char *otempstr;
5372 char *otmp;
5373 char *retval;
5374 tempstr = (char *) mycalloc(LINESIZE,sizeof(char));
5375 tmp = (char *) mycalloc(LINESIZE,sizeof(char));
5376 otempstr = tempstr;
5377 otmp = tmp;
5378 options->nummg = 2;
5379 options->mg = (MYREAL *) myrealloc (options->mg, sizeof (MYREAL) * 3);
5380 if(myID == MASTER)
5381 {
5382 retval = fgets (tmp, LINESIZE, options->parmfile);
5383 }
5384 else
5385 {
5386 sgets (tmp, LINESIZE, buffer);
5387 }
5388 get_next_word(&tmp,",{} \n",&tempstr);
5389 // while(strchr("{ ", tempstr[0]))
5390 // tempstr = strsep(&tmp,",{} \n");
5391 options->mg[0] = atof(tempstr);
5392 get_next_word(&tmp,",{} \n",&tempstr);
5393 options->mg[1] = atof(tempstr);
5394 myfree(otmp);
5395 myfree(otempstr);
5396 }
5397
5398 ///
5399 /// read the values for start theta from the list in the parmfile
5400 /// called through read_options_master()
read_theta(option_fmt * options)5401 void read_theta (option_fmt * options)
5402 {
5403 long i = 0;
5404
5405 char varvalue[LINESIZE];
5406 char ch;
5407
5408 char *tempstr;
5409 char *tmp;
5410 char *otempstr;
5411 char *otmp;
5412
5413 tempstr = (char *) mycalloc(LINESIZE,sizeof(char));
5414 tmp = (char *) mycalloc(LINESIZE,sizeof(char));
5415 otempstr = tempstr;
5416 otmp = tmp;
5417 //xcode removed ch out of while
5418 while ((getc (options->parmfile)) != '=')
5419 ;
5420 ch = getc (options->parmfile);
5421 while (!isspace ((int) ch) && ch != ':' && ch != '{')
5422 {
5423 varvalue[i++] = ch;
5424 ch = getc (options->parmfile);
5425 }
5426 varvalue[i]='\0';
5427 // printf("%s, %c\n",varvalue,ch);
5428 switch (uppercase (varvalue[0]))
5429 {
5430 case 'N':
5431 options->thetaguess = NRANDOMESTIMATE;
5432 read_random_theta(options, NULL);
5433 case 'U':
5434 options->thetaguess = URANDOMESTIMATE;
5435 read_random_theta(options, NULL);
5436 break;
5437 case 'F':
5438 case '_':
5439 options->thetaguess = FST;
5440 break;
5441 //case 'G':
5442 case 'O':
5443 case '0':
5444 options->thetaguess = OWN;
5445 ch = skip_space (options);
5446 if (ch == '\0')
5447 return;
5448 if (ch == '{')
5449 {
5450 ch = skip_space (options);
5451 while (ch != '}')
5452 {
5453 i = 0;
5454 if (ch == '\0')
5455 return;
5456 while (!strchr(" ,\t}", ch))
5457 {
5458 tmp[i++] = ch;
5459 ch = getc (options->parmfile);
5460 }
5461 tmp[i] = '\0';
5462 options->thetag[options->numthetag] = atof (tmp);
5463 options->numthetag += 1;
5464 options->thetag = (MYREAL *) myrealloc (options->thetag, sizeof (MYREAL) * (1 + options->numthetag));
5465 ch = skip_space (options);
5466 }
5467 }
5468 else
5469 {
5470 i = 0;
5471 while (!isspace (ch))
5472 {
5473 tmp[i++] = ch;
5474 ch = getc (options->parmfile);
5475 }
5476 tmp[i] = '\0';
5477 options->thetag[options->numthetag] = atof (tmp);
5478 options->numthetag += 1;
5479 options->thetag =
5480 (MYREAL *) myrealloc (options->thetag,
5481 sizeof (MYREAL) * (1 + options->numthetag));
5482 }
5483 options->numpop = options->numthetag;
5484 if (options->numthetag == 0)
5485 {
5486 warning ("You forgot to add your guess value:\n");
5487 warning ("Theta=Own:{pop1,pop2, ...}\n");
5488 warning ("or Theta=Own:guess_pop (same value for all)\n");
5489 }
5490 break;
5491 default:
5492 #ifdef MPI
5493 fprintf(stderr,"%i> Problem with read_theta()\n",myID);
5494 #endif
5495 warning
5496 ("Failure to read start theta method, should be\ntheta=FST or theta=Own:x.x\n or theta=Own:{x.x, x.x , x.x, .....} or theta=NRANDOM:{mean,std} or URANDOM:{min,max}; evasive action: use FST start method instead");
5497 options->thetaguess = FST;
5498 break;
5499 }
5500 myfree(otmp);
5501 myfree(otempstr);
5502 }
5503
5504
5505 void
read_mig(option_fmt * options)5506 read_mig (option_fmt * options)
5507 {
5508 // char parmvar[LINESIZE];
5509 long test = 0;
5510 long i = 0;
5511
5512 char varvalue[LINESIZE];
5513 char ch;
5514
5515 char *tempstr;
5516 char *tmp;
5517 char *otempstr;
5518 char *otmp;
5519 int choint;
5520 char choice;
5521 varvalue[0]='#';// to keep the static analyzer happy
5522
5523 tempstr = (char *) mycalloc(LINESIZE,sizeof(char));
5524 tmp = (char *) mycalloc(LINESIZE,sizeof(char));
5525 otempstr = tempstr;
5526 otmp = tmp;
5527
5528 /* 1st example: 1.0 (n-island model)
5529 2nd example: {1.0} (migration matrix model, all the same start values
5530 3rd example: the dashes on the diagonal are NECESSARY, {} are facultativ
5531 - 1.0 0.1
5532 1.0 - 2.0
5533 0.9 1.2 -
5534 to specify real 0.0 you need to use the custom-migration settings.
5535 0.0 in the table will be change to SMALLES_MIGRATION
5536 */
5537 //xcode removed ch out of while
5538 while ((getc (options->parmfile)) != '=')
5539 ;
5540 // {
5541 // parmvar[i++] = ch;
5542 // }
5543 i = 0;
5544 ch = getc (options->parmfile);
5545 while (!isspace ((int) ch) && ch != ':' && ch != '{')
5546 {
5547 varvalue[i++] = ch;
5548 ch = getc (options->parmfile);
5549 }
5550 choint = (int) varvalue[0];
5551 choice = uppercase (choint);
5552 switch (choice)
5553 {
5554 case 'N':
5555 options->migrguess = NRANDOMESTIMATE;
5556 read_random_mig(options,NULL);
5557 break;
5558 case 'U':
5559 options->migrguess = URANDOMESTIMATE;
5560 read_random_mig(options,NULL);
5561 break;
5562 case 'F':
5563 case '_':
5564 options->migrguess = FST;
5565 break;
5566 // case 'G':
5567 case 'O':
5568 case '0':
5569 // if(varvalue[0]=='G')
5570 // options->migrguess = PARAMGRID;
5571 // else
5572 options->migrguess = OWN;
5573 ch = skip_space (options);
5574 if (ch == '\0')
5575 return;
5576 if (ch == '{')
5577 {
5578 options->migration_model = MATRIX;
5579 ch = skip_space (options);
5580 while (ch != '}')
5581 {
5582 if ((ch == '\0') || (ch == '#'))
5583 return;
5584 i = 0;
5585 while (!strchr(" ,\t}",ch))
5586 {
5587 tmp[i++] = ch;
5588 ch = getc (options->parmfile);
5589 if(ch == '#')
5590 {
5591 while(ch!='\n' || ch != '\0')
5592 ch = getc(options->parmfile);
5593 return;
5594 }
5595 }
5596 tmp[i] = '\0';
5597 if (strcmp (tmp, "-"))
5598 {
5599 options->mg[options->nummg] = atof (tmp);
5600 options->nummg += 1;
5601 options->mg =
5602 (MYREAL *) myrealloc (options->mg,
5603 sizeof (MYREAL) * (1 +
5604 options->nummg));
5605 }
5606 else
5607 {
5608 test++;
5609 }
5610 ch = skip_space (options);
5611 }
5612 options->numpop = test;
5613 }
5614 else
5615 {
5616 options->migration_model = ISLAND;
5617 i = 0;
5618 options->numpop = 1;
5619 while (!isspace (ch))
5620 {
5621 tmp[i++] = ch;
5622 ch = getc (options->parmfile);
5623 }
5624 options->mg[options->nummg] = atof (tmp);
5625 options->nummg += 1;
5626 }
5627 if (options->nummg == 0)
5628 {
5629 warning ("You forgot to add your guess value, use either:\n");
5630 warning ("migration=FST\n");
5631 warning ("migration=Own:{migration matrix, diagonal is -}\n");
5632 usererror
5633 ("migration=Own:migration_value, all matrix elements have the same value\n");
5634 }
5635 break;
5636 default:
5637 warning ("Failure to read start migration method; evasive action: use FST instead\n");
5638 options->migrguess = FST;
5639 break;
5640 }
5641 myfree(otmp);
5642 myfree(otempstr);
5643 }
5644
5645 #ifdef MPI
5646 void
read_theta_worker(char ** buffer,option_fmt * options)5647 read_theta_worker (char **buffer, option_fmt * options)
5648 {
5649 long i = 0;
5650
5651 char varvalue[LINESIZE];
5652 char ch;
5653
5654 //char *tempstr;
5655 char *tmp;
5656 //char *otempstr;
5657 char *otmp;
5658
5659 //tempstr = (char *) mycalloc(LINESIZE,sizeof(char));
5660 tmp = (char *) mycalloc(LINESIZE,sizeof(char));
5661 //otempstr = tempstr;
5662 otmp = tmp;
5663
5664 while ((ch = sgetc (buffer)) != '=')
5665 ;
5666 ch = sgetc (buffer);
5667 while (!isspace ((int) ch) && ch != ':' && ch != '{')
5668 {
5669 varvalue[i++] = ch;
5670 ch = sgetc (buffer);
5671 }
5672 switch (uppercase (varvalue[0]))
5673 {
5674 case 'N':
5675 options->thetaguess = NRANDOMESTIMATE;
5676 read_random_theta(options, buffer);
5677 break;
5678 case 'U':
5679 options->thetaguess = URANDOMESTIMATE;
5680 read_random_theta(options, buffer);
5681 break;
5682 break;
5683 case 'F':
5684 case '_':
5685 options->thetaguess = FST;
5686 break;
5687 //case 'G':
5688 case 'O':
5689 case '0':
5690 //perhaps to come
5691 // if(varvalue[0]=='G')
5692 //options->thetaguess = PARAMGRID;
5693 //else
5694 options->thetaguess = OWN;
5695 ch = skip_sspace (buffer);
5696 if (ch == '\0')
5697 return;
5698 if (ch == '{')
5699 {
5700
5701 while (ch != '}')
5702 {
5703 i = 0;
5704 ch = skip_sspace (buffer);
5705 if (ch == '\0')
5706 return;
5707 while (ch != ' ' && ch != ',' && ch != '}')
5708 {
5709 tmp[i++] = ch;
5710 ch = sgetc (buffer);
5711 }
5712 tmp[i] = '\0';
5713 options->thetag[options->numthetag] = atof (tmp);
5714 options->numthetag += 1;
5715 options->thetag =
5716 (MYREAL *) myrealloc (options->thetag,
5717 sizeof (MYREAL) * (1 +
5718 options->numthetag));
5719
5720 }
5721 }
5722 else
5723 {
5724 i = 0;
5725 tmp[i++] = ch;
5726 while (!isspace (ch))
5727 {
5728 tmp[i++] = ch;
5729 ch = sgetc (buffer);
5730 }
5731 tmp[i] = '\0';
5732 options->thetag[options->numthetag] = atof (tmp);
5733 options->numthetag += 1;
5734 options->thetag =
5735 (MYREAL *) myrealloc (options->thetag,
5736 sizeof (MYREAL) * (1 + options->numthetag));
5737 }
5738 options->numpop = options->numthetag;
5739 if (options->numthetag == 0)
5740 {
5741 warning ("You forgot to add your guess value:\n");
5742 warning ("Theta=Own:{pop1,pop2, ...}\n");
5743 warning ("or Theta=Own:guess_pop (same value for all)\n");
5744 }
5745 break;
5746 default:
5747 usererror
5748 ("Failure to read start theta method, should be\ntheta=FST or theta=Own:x.x\n or theta=Own:{x.x, x.x , x.x, .....}");
5749 }
5750 myfree(otmp);
5751 //myfree(otempstr);
5752 }
5753
5754
5755 void
read_mig_worker(char ** buffer,option_fmt * options)5756 read_mig_worker (char **buffer, option_fmt * options)
5757 {
5758 char ch;
5759 long test = 0;
5760 long i = 0;
5761
5762 // char input[LINESIZE];
5763 char varvalue[LINESIZE];
5764
5765 //char *tempstr;
5766 char *tmp;
5767 //char *otempstr;
5768 char *otmp;
5769
5770 //tempstr = (char *) mycalloc(LINESIZE,sizeof(char));
5771 tmp = (char *) mycalloc(LINESIZE,sizeof(char));
5772 //otempstr = tempstr;
5773 otmp = tmp;
5774
5775 /* 1st example: 1.0 (n-island model)
5776 2nd example: {1.0} (migration matrix model, all the same start values
5777 3rd example: the dashes on the diagonal are NECESSARY, {} are facultativ
5778 - 1.0 0.1
5779 1.0 - 2.0
5780 0.9 1.2 -
5781 to specify real 0.0 you need to use the custom-migration settings.
5782 0.0 in the table will be change to SMALLES_MIGRATION
5783 */
5784
5785 while ((ch = sgetc (buffer)) != '=')
5786 {
5787 // parmvar[i++] = ch;
5788 }
5789 i = 0;
5790 ch = sgetc (buffer);
5791 while (!isspace ((int) ch) && ch != ':' && ch != '{')
5792 {
5793 varvalue[i++] = ch;
5794 ch = sgetc (buffer);
5795 }
5796 switch (uppercase (varvalue[0]))
5797 {
5798 case 'N':
5799 options->migrguess = NRANDOMESTIMATE;
5800 read_random_mig(options, buffer);
5801 case 'U':
5802 options->migrguess = URANDOMESTIMATE;
5803 read_random_mig(options, buffer);
5804 break;
5805 case 'F':
5806 case '_':
5807 options->migrguess = FST;
5808 break;
5809 // case 'G':
5810 case 'O':
5811 case '0':
5812 // if(varvalue[0]=='G')
5813 // options->migrguess = PARAMGRID;
5814 // else
5815 options->migrguess = OWN;
5816 ch = skip_sspace (buffer);
5817 if (ch == '\0')
5818 return;
5819 if (ch == '{')
5820 {
5821 options->migration_model = MATRIX;
5822 while (ch != '}')
5823 {
5824 ch = skip_sspace (buffer);
5825 if ((ch == '\0') || (ch == '}'))
5826 return;
5827 i = 0;
5828 while (!isspace (ch) && ch != ',' && ch != '}')
5829 {
5830 tmp[i++] = ch;
5831 ch = sgetc (buffer);
5832 }
5833 tmp[i] = '\0';
5834 if (strcmp (tmp, "-"))
5835 {
5836 options->mg[options->nummg] = atof (tmp);
5837 options->nummg += 1;
5838 options->mg =
5839 (MYREAL *) myrealloc (options->mg,
5840 sizeof (MYREAL) * (1 +
5841 options->nummg));
5842 }
5843 else
5844 {
5845 test++;
5846 }
5847 }
5848 options->numpop = test;
5849 }
5850 else
5851 {
5852 options->migration_model = ISLAND;
5853 i = 0;
5854 options->numpop = 1;
5855 tmp[i++] = ch;
5856 while (!isspace (ch))
5857 {
5858 tmp[i++] = ch;
5859 ch = sgetc (buffer);
5860 }
5861 options->mg[options->nummg] = atof (tmp);
5862 options->nummg += 1;
5863 }
5864 if (options->nummg == 0)
5865 {
5866 warning ("You forgot to add your guess value, use either:\n");
5867 warning ("migration=FST\n");
5868 warning ("migration=Own:{migration matrix, diagonal is -}\n");
5869 usererror
5870 ("migration=Own:migration_value, all matrix elements have the same value\n");
5871 }
5872 break;
5873 default:
5874 usererror ("Failure to read start migration method\n");
5875 }
5876 myfree(otmp);
5877 //myfree(otempstr);
5878 }
5879 #endif
5880
5881 char
skip_space(option_fmt * options)5882 skip_space (option_fmt * options)
5883 {
5884 char ch = getc (options->parmfile);
5885 while (isspace ((int) ch) || ch == ',')
5886 {
5887 ch = getc (options->parmfile);
5888 }
5889 if (isalpha (ch))
5890 {
5891 ungetc (ch, options->parmfile);
5892 ch = '\0';
5893 }
5894 return ch;
5895 }
5896
5897 #ifdef MPI
5898 char
skip_sspace(char ** buffer)5899 skip_sspace (char **buffer)
5900 {
5901 char ch = sgetc (buffer);
5902 while (isspace ((int) ch) || ch == ',')
5903 {
5904 ch = sgetc (buffer);
5905 if (isalpha (ch))
5906 return ch;
5907 }
5908 return ch;
5909 }
5910 #endif
5911
5912 void
set_profile_options(option_fmt * options)5913 set_profile_options (option_fmt * options)
5914 {
5915 switch (options->profile)
5916 {
5917 case _NONE:
5918 options->printprofsummary = options->printprofile = FALSE;
5919 break;
5920 case ALL:
5921 options->printprofsummary = options->printprofile = TRUE;
5922 break;
5923 case TABLES:
5924 options->printprofsummary = FALSE;
5925 options->printprofile = TRUE;
5926 break;
5927 case SUMMARY:
5928 options->printprofsummary = TRUE;
5929 options->printprofile = FALSE;
5930 break;
5931 }
5932 if (options->profilemethod == 'd')
5933 options->printprofsummary = FALSE;
5934 }
5935
5936
5937
5938 /* custom-migration:<{> migration matrix and theta on
5939 diagonal:
5940 0 means not estimated,
5941 x means estimated, s means symmetrically
5942 estimated, m means all are the same,
5943 c means remains constant at start value <}>
5944 example:
5945 {* * s
5946 * c *
5947 s 0 *}
5948 */
5949 void
read_custom_migration(FILE * file,option_fmt * options,char * value,long customnumpop)5950 read_custom_migration (FILE * file, option_fmt * options, char *value,
5951 long customnumpop)
5952 {
5953 long cc = customnumpop;
5954 long zz = 0, z = 0;
5955 char ch = '\0';
5956 long lc, numpop, i, j, ii;
5957 // long position = 0;
5958 // fpos_t thePos;
5959 if (cc == 0)
5960 cc = 1000000;
5961 else
5962 cc *= cc;
5963
5964 z = 0;
5965 zz = 0;
5966 while (ch != '}' && zz < cc + options->gamma)
5967 {
5968 ch = value[z];
5969 switch (ch)
5970 {
5971 case '}':
5972 case '{':
5973 case ' ':
5974 case '\t':
5975 z++;
5976 break;
5977 case '\0':
5978 case '\n':
5979 z = 0;
5980 if (file == stdin)
5981 printf ("Enter the next value or list of values\n");
5982 FGETS (value, LINESIZE, file);
5983 break;
5984 default:
5985 options->custm =
5986 (char *) myrealloc (options->custm, sizeof (char) * (zz + 2));
5987 options->custm2 =
5988 (char *) myrealloc (options->custm2, sizeof (char) * (zz + 2));
5989 switch (ch)
5990 {
5991 case 'S':
5992 options->custm[zz++] = ch;
5993 break;
5994 case 'M': //do we have code for this?
5995 options->custm[zz++] = ch;
5996 break;
5997 case 'x':
5998 case 'X':
5999 options->custm[zz++] = '*';
6000 break;
6001 case 'c':
6002 case 'C':
6003 options->custm[zz++] = 'c';
6004 break;
6005 default:
6006 options->custm[zz++] = tolower (ch);
6007 break;
6008 }
6009 z++;
6010 }
6011 }
6012 options->custm[zz] = '\0';
6013 lc = (long) strlen (options->custm);
6014 numpop = (long) sqrt ((MYREAL) lc);
6015 z = numpop;
6016 for (i = 0; i < numpop; i++)
6017 {
6018 for (j = 0; j < numpop; j++)
6019 {
6020 ii = i * numpop + j;
6021 if (i == j)
6022 options->custm2[i] = options->custm[ii];
6023 else
6024 options->custm2[z++] = options->custm[ii];
6025 }
6026 }
6027 options->custm2[z] = '\0';
6028 specify_migration_type (options);
6029 if (file != stdin)
6030 {
6031 //position = ftell (options->parmfile);
6032 fgetpos(options->parmfile, &thePos);
6033 }
6034 while (file != stdin && !(strstr (value, "end") || strchr (value, '=')))
6035 {
6036 // position = ftell (options->parmfile);
6037 fgetpos(options->parmfile, &thePos);
6038 FGETS (value, LINESIZE, file);
6039 }
6040 if (file != stdin)
6041 {
6042 //reset_oneline (options, position);
6043 reset_oneline (options, &thePos);
6044 }
6045 //printf("%li %li %li : %s\n", zz, z, lc, options->custm);
6046 }
6047
6048 #ifdef MPI
6049 void
read_custom_migration_worker(char ** buffer,option_fmt * options,char * value,long customnumpop)6050 read_custom_migration_worker (char **buffer, option_fmt * options,
6051 char *value, long customnumpop)
6052 {
6053 long cc = customnumpop;
6054 long zz = 0, z = 0;
6055 char ch = '\0';
6056 long lc, numpop, i, j, ii;
6057 char *position;
6058
6059 if (cc == 0)
6060 cc = 1000000;
6061 else
6062 cc *= cc;
6063
6064 z = 0;
6065 zz = 0;
6066 while (ch != '}' && zz < cc)
6067 {
6068 ch = value[z];
6069 switch (ch)
6070 {
6071 case '}':
6072 case '{':
6073 case ' ':
6074 case '\t':
6075 z++;
6076 break;
6077 case '\0':
6078 case '\n':
6079 z = 0;
6080 sgets (value, LINESIZE, buffer);
6081 break;
6082 default:
6083 options->custm =
6084 (char *) myrealloc (options->custm, sizeof (char) * (zz + 2));
6085 options->custm2 =
6086 (char *) myrealloc (options->custm2, sizeof (char) * (zz + 2));
6087 switch (ch)
6088 {
6089 case 'S':
6090 options->custm[zz++] = ch;
6091 break;
6092 case 'M': //do we have code for this?
6093 options->custm[zz++] = ch;
6094 break;
6095 case 'x':
6096 case 'X':
6097 options->custm[zz++] = '*';
6098 break;
6099 default:
6100 options->custm[zz++] = tolower (ch);
6101 break;
6102 }
6103 z++;
6104 }
6105 }
6106 options->custm[zz] = '\0';
6107 lc = (long) strlen (options->custm);
6108 numpop = (long) sqrt ((MYREAL) lc);
6109 z = numpop;
6110 for (i = 0; i < numpop; i++)
6111 {
6112 for (j = 0; j < numpop; j++)
6113 {
6114 ii = i * numpop + j;
6115 if (i == j)
6116 options->custm2[i] = options->custm[ii];
6117 else
6118 options->custm2[z++] = options->custm[ii];
6119 }
6120 }
6121 options->custm2[z] = '\0';
6122 specify_migration_type (options);
6123 position = *buffer;
6124 while (!(strstr (value, "end") || strchr (value, '=')))
6125 {
6126 position = *buffer;
6127 sgets (value, LINESIZE, buffer);
6128 }
6129 *buffer = position;
6130 ;
6131 }
6132 #endif /*MPI*/
6133 void
specify_migration_type(option_fmt * options)6134 specify_migration_type (option_fmt * options)
6135 {
6136 long len = (long) strlen (options->custm);
6137 long ms = 0, xs = 0, ns = 0, ss = 0, len2, i;
6138 char *p;
6139 p = options->custm;
6140 while (*p != '\0')
6141 {
6142 switch (*p)
6143 {
6144 case 'm':
6145 ms++;
6146 break;
6147 case 'x':
6148 case '*':
6149 xs++;
6150 break;
6151 case '0':
6152 ns++;
6153 break;
6154 case 'S':
6155 case 's':
6156 ss++;
6157 break;
6158 case 'c':
6159 break;
6160 }
6161 p++;
6162 }
6163 if (ms >= len)
6164 {
6165 options->migration_model = ISLAND;
6166 return;
6167 }
6168 if (xs >= len)
6169 {
6170 options->migration_model = MATRIX;
6171 return;
6172 }
6173 if (ns >= len)
6174 {
6175 usererror ("Custom migration matrix was completely set to zero?!\n");
6176 return;
6177 }
6178 len2 = (long) sqrt ((MYREAL) len);
6179 if (ms == len2 && xs == len - len2)
6180 {
6181 for (i = 0; i < len2; i++)
6182 {
6183 if (options->custm[i * len2 + i] != 'm')
6184 {
6185 options->migration_model = MATRIX;
6186 return;
6187 }
6188 }
6189 options->migration_model = MATRIX_SAMETHETA;
6190 return;
6191 }
6192 if (xs == len2 && ms == len - len2)
6193 {
6194 for (i = 0; i < len2; i++)
6195 {
6196 if (options->custm[i * len2 + i] != '*')
6197 {
6198 options->migration_model = MATRIX;
6199 return;
6200 }
6201 }
6202 options->migration_model = ISLAND_VARTHETA;
6203 return;
6204 }
6205 options->migration_model = MATRIX_ARBITRARY;
6206 }
6207
6208
6209
6210 void
fillup_custm(long len,world_fmt * world,option_fmt * options)6211 fillup_custm (long len, world_fmt * world, option_fmt * options)
6212 {
6213 long i, j, ii, z;
6214 char *tmp;
6215 len = (long) strlen (options->custm);
6216 tmp = (char *) mycalloc (1, sizeof (char) * (world->numpop2 + 2));
6217 options->custm =
6218 (char *) myrealloc (options->custm, sizeof (char) * (world->numpop2 + 2));
6219
6220 options->custm2 =
6221 (char *) myrealloc (options->custm2, sizeof (char) * (world->numpop2 + 2));
6222 strncpy (tmp, options->custm, world->numpop2 + options->gamma);
6223 z = world->numpop;
6224 for (i = 0; i < world->numpop; i++)
6225 {
6226 for (j = 0; j < world->numpop; j++)
6227 {
6228 ii = i * world->numpop + j;
6229 if (ii < len)
6230 options->custm[ii] = tmp[ii];
6231 else
6232 options->custm[ii] = '*';
6233 if (i == j)
6234 options->custm2[i] = options->custm[ii];
6235 else
6236 options->custm2[z++] = options->custm[ii];
6237 }
6238 }
6239 if (options->gamma)
6240 {
6241 if (len <= world->numpop2)
6242 {
6243 options->custm[world->numpop2] = '*';
6244 options->custm2[world->numpop2] = '*';
6245 }
6246 else
6247 {
6248 options->custm[world->numpop2] = tmp[world->numpop2];
6249 options->custm2[world->numpop2] = tmp[world->numpop2];
6250 }
6251 options->custm[world->numpop2 + 1] = '\0';
6252 options->custm2[world->numpop2 + 1] = '\0';
6253 }
6254 else
6255 {
6256 options->custm[world->numpop2] = '\0';
6257 options->custm2[world->numpop2] = '\0';
6258 }
6259 strcpy (world->options->custm, options->custm);
6260 strcpy (world->options->custm2, options->custm2);
6261 myfree(tmp);
6262 }
6263
6264 void
print_arbitrary_migration_table(FILE * file,world_fmt * world,option_fmt * options,data_fmt * data)6265 print_arbitrary_migration_table (FILE * file, world_fmt * world, option_fmt *options,
6266 data_fmt * data)
6267 {
6268 long i,j,z;
6269 char mytext[LONGLINESIZE];
6270 switch (world->options->migration_model)
6271 {
6272 case ISLAND:
6273 strcpy (mytext, "N-Island migration model");
6274 fprintf (file, "Migration model:\n %-44.44s\n", mytext);
6275 break;
6276 case ISLAND_VARTHETA:
6277 strcpy (mytext, "N-Island migration model with variable Theta");
6278 fprintf (file, "Migration model:\n %-44.44s\n", mytext);
6279 break;
6280 case MATRIX:
6281 strcpy (mytext, "Migration matrix model with variable Theta ");
6282 fprintf (file, "Migration model:\n %-44.44s\n", mytext);
6283 break;
6284 case MATRIX_SAMETHETA:
6285 strcpy (mytext, "Migration matrix model with same Theta");
6286 fprintf (file, "Migration model:\n %-44.44s\n", mytext);
6287 break;
6288 case MATRIX_ARBITRARY:
6289 default:
6290 strcpy (mytext, "Arbitrary migration matrix model");
6291 fprintf (file, "Migration model: %-44.44s\n", mytext);
6292 fprintf (file,
6293 "[Legend: m = average (average over a group of Thetas or M]\n");
6294 fprintf (file,
6295 "[s = symmetric M, S = symmetric 4Nm,\n 0 = zero, and not estimated, ]\n");
6296 fprintf (file, "[* = free to vary, Thetas are on diagonal]\n");
6297 if (world->options->migration_model == MATRIX_ARBITRARY)
6298 {
6299 for (i = 0; i < data->numpop; i++)
6300 {
6301 fprintf (file, "%10.10s ", data->popnames[i]);
6302 for (j = 0; j < data->numpop; j++)
6303 {
6304 z = world->numpop * (options->newpops[i]-1) + (options->newpops[j]-1);
6305 fprintf (file, "%c ", world->options->custm[z]);
6306 }
6307 fprintf (file, "\n");
6308 }
6309 fprintf (file, "\n\n");
6310 }
6311 break;
6312 }
6313 }
6314
6315 void
print_distance_table(FILE * file,world_fmt * world,option_fmt * options,data_fmt * data)6316 print_distance_table (FILE * file, world_fmt * world, option_fmt * options,
6317 data_fmt * data)
6318 {
6319 long i, j;
6320
6321 if (options->geo)
6322 {
6323 fprintf (file,
6324 " Geographic distance matrix between locations\n ");
6325 for (i = 0; i < world->numpop; i++)
6326 {
6327 fprintf (file, "%-10.10s ", data->popnames[i]);
6328 for (j = 0; j < world->numpop; j++)
6329 {
6330 if (i == j)
6331 fprintf (file, " - ");
6332 else
6333 fprintf (file, "%6.3f ", data->ogeo[i][j]);
6334 }
6335 fprintf (file, "\n ");
6336 }
6337 fprintf (file, "\n");
6338 }
6339 }
6340
reorder_populations(world_fmt * world,option_fmt * options,data_fmt * data)6341 void reorder_populations(world_fmt *world, option_fmt *options, data_fmt *data)
6342 {
6343 // long newnumpop=0;
6344 long pop;
6345 long lim = data->numpop;
6346 long i;
6347 for(pop=0;pop<lim;pop++)
6348 {
6349 for(i=0; i<world->sumtips;i++)
6350 {
6351 if(world->nodep[i]->actualpop == pop)
6352 {
6353 world->nodep[i]->actualpop = world->nodep[i]->pop = options->newpops[pop]-1;
6354 }
6355 }
6356 }
6357 }
6358
set_localities(char ** value,char ** tmp,option_fmt * options)6359 void set_localities(char **value, char **tmp, option_fmt *options)
6360 {
6361 long z=0;
6362 options->newpops[0] = 1;
6363 get_next_word(value,"{}:,; ",tmp);
6364 while(*tmp != NULL)
6365 {
6366 if(z >= options->newpops_numalloc)
6367 {
6368 options->newpops_numalloc = z+1;
6369 options->newpops = (long*) myrealloc(options->newpops, options->newpops_numalloc * sizeof(long));
6370 }
6371 options->newpops[z] = atol(*tmp);
6372 //printf("%i> population relabel [%li]=%li (input: |%s|)\n",myID,z, options->newpops[z],*tmp);
6373 z++;
6374 get_next_word(value,"{}:,; ",tmp);
6375 }
6376 }
6377