1 /* ! \file menu.c */
2 /*------------------------------------------------------
3 Maximum likelihood estimation
4 of migration rate and effectice population size
5 using a Metropolis-Hastings Monte Carlo algorithm
6 -------------------------------------------------------
7 M E N U R O U T I N E S
8
9 presents the menu and its submenus.
10 Peter Beerli 1996, Seattle
11 beerli@fsu.edu
12
13 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
14 Copyright 2003-2007 Peter Beerli, Tallahassee FL
15
16 Permission is hereby granted, free of charge, to any person obtaining a copy
17 of this software and associated documentation files (the "Software"), to deal
18 in the Software without restriction, including without limitation the rights
19 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
20 of the Software, and to permit persons to whom the Software is furnished to do
21 so, subject to the following conditions:
22
23 The above copyright notice and this permission notice shall be included in all copies
24 or substantial portions of the Software.
25
26 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
27 INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
28 PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
29 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
30 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
31 OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
32
33
34 $Id: menu.c 2172 2013-10-30 17:50:08Z beerli $
35
36 -------------------------------------------------------*/
37 #include "migration.h"
38 #include "sighandler.h"
39 #include "tools.h"
40 #include "sequence.h"
41 #include "fst.h"
42 #include "options.h"
43 #include "migrate_mpi.h"
44 #include "menu.h"
45 #include "pretty.h"
46
47 #ifdef DMALLOC_FUNC_CHECK
48 #include <dmalloc.h>
49 #endif
50 /* prototypes ------------------------------------------- */
51 //void print_menu_title(FILE * file, option_fmt * options);
52 //void print_menu_accratio(long a, long b, world_fmt * world);
53 //long print_title(world_fmt * world, option_fmt * options);
54 //void get_menu(option_fmt * options);
55 /* private functions */
56 void setup_datatype(char *datatype, option_fmt * options);
57 void menuData(option_fmt * options, char datatype[]);
58 void menuInput(option_fmt * options);
59 void menuParameters(option_fmt * options);
60 void menuStrategy(option_fmt * options);
61 void menuHeat(option_fmt * options, char *input);
62 void display_ml_mcmc(option_fmt * options);
63 void display_bayes_mcmc(option_fmt * options);
64 //void menuSequences(option_fmt * options);
65 void read_custom_menu_migration(option_fmt * options);
66 void read_custom_menu_lratio(option_fmt * options);
67 char *custom_migration_type(long type);
68 void read_heatvalues(option_fmt * options);
69 void menuRandom(MYREAL * param, char type);
70 void start_tree_method(option_fmt * options);
71 void start_data_method(option_fmt * options);
72 void how_many_pop(long *numpop);
73 void set_menu_localities(option_fmt *options);
74 void set_localities_string(char *loc, option_fmt *options);
75 char dialog_lrt(long numpop, lratio_fmt * lratio);
76 char menuread_lrt_paramvalue(char *value, long *counter, long numpop2);
77 void get_plotmenu(option_fmt * options);
78 boolean menuStrategy_ml(option_fmt * options);
79 boolean menuStrategy_bayes(option_fmt * options);
80 long get_prior(char *input);
81 void set_prior(char *output, int * prior, boolean without_rate);
82 void set_proposal(char *output, boolean *proposal, boolean without_rate);
83 char * submodeltype(int type);
84
85 #ifdef POPMODEL
86 #include "popmodel.h"
87 #endif
88 #define MI_INFILE 1
89 #define MI_RAND 2
90 #define MI_TITLE 3
91 #define MI_SUMREAD 4
92 #define MI_PROGRESS 5
93 #define MI_PRINTDATA 6
94 #define MI_OUTFILE 7
95 #define MI_PLOT 8
96 #define MI_PROFILE 9
97 #define MI_LRT 10
98 #define MI_AIC 11
99 #define MI_TREES 12
100 #define MI_PLOTCOORD 13
101 #define MI_SUMFILE 14
102 #define MI_LOGFILE 15
103 #define MI_UEPFILE 16
104 #define MI_UEPRATE 17
105 #define MI_UEPFREQ 18
106 #define MI_MIGHISTOGRAM 19
107 #define MI_SKYLINE 20
108
109 ///
110 ///Enumeration of possible population models
111 enum popmodel_menu
112 {
113 WRIGHTFISHER, CANNING, MORAN
114 };
115
116 ///
117 ///Enumeration of possible choices for ML menu
118 enum ml_menu
119 {
120 MLSTRATEGY, MLSHORTCHAINS, MLSHORTSKIP, MLSHORTSAMPLES, MLLONGCHAINS, MLLONGSKIP, MLLONGSAMPLES,
121 MLBURNIN, MLREPLICATE, MLHEAT, MLMOVINGSTEPS, MLEPSILON, MLGELMAN
122 };
123
124 ///
125 ///Enumeration of possible choices for Datatype menu
126 enum dtseq_menu
127 {
128 DTSEQZERO, DTSEQTYPE, DTSEQTRATIO, DTSEQFREQ, DTSEQSITECATEGS, DTSEQRATES, DTSEQCORR, DTSEQWEIGHT, DTSEQINTERLEAVED, DTSEQERROR, DTSEQFAST, DTSEQTREE, DTSEQINHERITANCE, DTSEQRANDOMSUBSET, DTSEQTIPDATE, DTSEQMUTRATE, DTSEQGENERATION
129 };
130
131 enum dtmsat_menu {
132 DTMSATZERO, DTMSATTYPE, DTMSATMISS, DTMSATTRESH, DTMSATTREE, DTMSATINHERITANCE, DTMSATRANDOMSUBSET, DTMSATTIPDATE, DTMSATMUTRATE, DTMSATGENERATION
133 };
134
135 enum dtbrown_menu
136 {
137 DTBROWNZERO, DTBROWNTYPE, DTBROWNMISS, DTBROWNTREE, DTBROWNINHERITANCE, DTBROWNRANDOMSUBSET, DTBROWNTIPDATE, DTBROWNMUTRATE, DTBROWNGENERATION
138 };
139 enum dtep_menu {
140 DTEPZERO, DTEPTYPE, DTEPMISS, DTEPTREE, DTEPINHERITANCE, DTEPRANDOMSUBSET, DTEPTIPDATE, DTEPMUTRATE, DTEPGENERATION
141 };
142
143 ///
144 ///Enumeration of possible choices for Bayes menu
145 enum bayesian_menu {
146 BAYESSTRATEGY, BAYESOUT, BAYESMDIMOUT, BAYESBINNING, BAYESPRETTY, BAYESFREQ, BAYESPROPOSAL, BAYESPRIOR, BAYESLCHAINS, BAYESSKIP, BAYESSAMPLES, BAYESBURNIN,
147 BAYESREPLICATE, BAYESHEAT, BAYESMOVINGSTEPS, BAYESGELMAN, BAYESPRIORALONE
148 };
149
150 ///
151 ///get user supplied filenames or fill the filenames with default values
152 void
menu_get_filename(char message[],char thedefault[],char * filename)153 menu_get_filename(char message[], char thedefault[], char *filename)
154 {
155 char input[LINESIZE];
156 printf("%s\n[Default: %s]\n===> ", message, thedefault);
157 fflush(stdout); FGETS(input, LINESIZE, stdin);
158 if (input[0] == '\0')
159 strcpy(filename, thedefault);
160 else
161 strcpy(filename, input);
162 }
163
164 ///
165 ///sets the plotoptions default is NO PLOT
166 void
get_plotoptions(option_fmt * options)167 get_plotoptions(option_fmt * options)
168 {
169 char input[LINESIZE];
170 do {
171 printf(" Plot Likelihood surface:\n");
172 printf
173 (" (B)oth to outfile and mathfile, (O)utfile only, (N)o plot\n===> ");
174 fflush(stdout); FGETS(input, LINESIZE, stdin);
175 }
176 while (strchr("BON", (int) uppercase(input[0])) == NULL);
177 switch (uppercase(input[0])) {
178 case 'B':
179 options->plotmethod = PLOTALL;
180 break;
181 case 'O':
182 options->plotmethod = PLOTOUTFILE;
183 break;
184 case 'N':
185 default:
186 options->plot = FALSE;
187 return;
188 }
189 get_plotmenu(options);
190 }
191
192 ///
193 ///allows the setting of the tickmarks and ranges for the plots
194 ///[remember these produce only ASCII plots
get_plotmenu(option_fmt * options)195 void get_plotmenu(option_fmt * options) {
196 char input[LINESIZE];
197 char answer = 'N';
198 int count;
199 do {
200 printf(" Current plot settings are:\n");
201 printf
202 (" Parameter: %s, Scale: %s, Intervals: %li\n",
203 options->plotvar == PLOT4NM ? "{Theta, xNm}" : "{Theta, M}",
204 options->plotscale == PLOTSCALELOG ? "Log10" : "Standard",
205 options->plotintervals);
206 printf(" Ranges: X-%5.5s: %f - %f\n",
207 options->plotvar == PLOT4NM ? "xNm" : "M", options->plotrange[0], options->plotrange[1]);
208 printf(" Ranges: Y-%5.5s: %f - %f\n",
209 "Theta", options->plotrange[2],
210 options->plotrange[3]);
211 printf(" Is this OK? [YES or No]\n===> ");
212 fflush(stdout); FGETS(input, LINESIZE, stdin);
213 answer = uppercase(input[0]);
214 if (answer != 'Y') {
215 printf("1 Use as X-axis: xNm Y-axis: Theta\n");
216 printf("2 Use as X-axis: M Y-axis: Theta\n\n===> ");
217 fflush(stdout); FGETS(input, LINESIZE, stdin);
218 if (input[0] == '2')
219 options->plotvar = PLOTM;
220 else
221 options->plotvar = PLOT4NM;
222 printf("1 Scale is standard \n");
223 printf("2 Scale is in Log10\n\n===> ");
224 fflush(stdout); FGETS(input, LINESIZE, stdin);
225 if (input[0] == '2')
226 options->plotscale = PLOTSCALELOG;
227 else
228 options->plotscale = PLOTSCALESTD;
229 do {
230 printf("How many plot positions?\n");
231 printf
232 ("[Default is 36, this changes only in the mathfile]\n\n===> ");
233 fflush(stdout); FGETS(input, LINESIZE, stdin);
234 if (strlen(input) > 1)
235 options->plotintervals = ATOL(input);
236 }
237 while (options->plotintervals < 2);
238 printf("Change the ranges of the axes\n");
239 printf(" Ranges: X-%5.5s: %f - %f\n",
240 options->plotvar == PLOT4NM ? "xNm" : "M",
241 options->plotrange[0],
242 options->plotrange[1]);
243 printf(" Ranges: Y-%5.5s: %f - %f\n",
244 "Theta", options->plotrange[2],
245 options->plotrange[3]);
246 printf
247 ("Give lowest and highest value for X axis:\n===> ");
248 fflush(stdout); FGETS(input, LINESIZE, stdin);
249 if (strlen(input) > 1) {
250 #ifdef USE_MYREAL_FLOAT
251 count = sscanf(input, "%f%f", &options->plotrange[0],
252 &options->plotrange[1]);
253 #else
254 count = sscanf(input, "%lf%lf", &options->plotrange[0],
255 &options->plotrange[1]);
256 #endif
257 }
258 printf
259 ("Give lowest and highest value for Y axis:\n===> ");
260 fflush(stdout); FGETS(input, LINESIZE, stdin);
261 if (strlen(input) > 1) {
262 #ifdef USE_MYREAL_FLOAT
263 count = sscanf(input, "%f%f", &options->plotrange[2],
264 &options->plotrange[3]);
265 #else
266 count = sscanf(input, "%lf%lf", &options->plotrange[2],
267 &options->plotrange[3]);
268 #endif
269 }
270 }
271 }
272 while (answer != 'Y');
273 }
274
275
276 void
print_menu_title(FILE * file,option_fmt * options)277 print_menu_title(FILE * file, option_fmt * options)
278 {
279 char nowstr[LINESIZE];
280 if (options->menu || options->progress) {
281 FPRINTF(file, " =============================================\n");
282 FPRINTF(file, " MIGRATION RATE AND POPULATION SIZE ESTIMATION\n");
283 FPRINTF(file, " using Markov Chain Monte Carlo simulation\n");
284 FPRINTF(file, " =============================================\n");
285 #ifdef MPI
286 FPRINTF(file, " Compiled for a PARALLEL COMPUTER ARCHITECTURE\n");
287 FPRINTF(file, " One master and %i compute nodes are available.\n", numcpu-1);
288 #endif
289 #ifdef THREAD
290 FPRINTF(file, " Compiled for a SYMMETRIC MULTIPROCESSORS\n");
291 #endif
292 #ifdef ALTIVEC
293 FPRINTF(file, " ALTIVEC enabled\n");
294 #endif
295 #ifdef FAST_EXP
296 FPRINTF(file, " Fast approximation to Exp() and Log() used\n");
297 #endif
298 #ifdef PRETTY
299 FPRINTF(file, " PDF output enabled [%s]\n",
300 #ifdef A4PAPER
301 "A4-size"
302 #else
303 "Letter-size"
304 #endif
305 );
306 #endif
307
308 FPRINTF(file, " Version %s ", MIGRATEVERSION);
309 if(strlen(MIGRATESUBVERSION)>0)
310 FPRINTF(file, " [%s]\n", MIGRATESUBVERSION);
311 else
312 FPRINTF(file, "\n");
313 get_time(nowstr, " %c");
314 if (nowstr[0] != '\0')
315 FPRINTF(file, " Program started at %s\n", nowstr);
316 FPRINTF(file, "\n\n");
317 }
318 }
319
320 void
print_menu_accratio(long a,long b,world_fmt * world)321 print_menu_accratio(long a, long b, world_fmt * world)
322 {
323 char buffer[STRSIZE];
324 boolean writelog = world->options->writelog;
325 boolean progress = world->options->progress;
326
327 if (writelog || progress)
328 {
329 sprintf(buffer, " Acceptance-ratio = %li/%li (%f)\n\n", a, b,
330 (MYREAL) a / (MYREAL) b);
331 if (progress)
332 FPRINTF(stdout, "%s", buffer);
333 if (writelog)
334 FPRINTF(world->options->logfile, "%s", buffer);
335 }
336 }
337
338 long
print_title(world_fmt * world,option_fmt * options)339 print_title(world_fmt * world, option_fmt * options)
340 {
341 char nowstr[LINESIZE];
342 long len = 45, i, filepos = -1;
343 if (!world->options->simulation)
344 {
345 FPRINTF(world->outfile, " ");
346 if (options->title[0] != '\0')
347 {
348 len = (long) MAX(strlen(options->title), 45);
349 for (i = 0; i < len; i++)
350 FPRINTF(world->outfile, "=");
351 FPRINTF(world->outfile, "\n %-*s\n ", (int) len,options->title);
352 }
353 for (i = 0; i < len; i++)
354 FPRINTF(world->outfile, "=");
355
356 FPRINTF(world->outfile,
357 "\n MIGRATION RATE AND POPULATION SIZE ESTIMATION\n");
358 FPRINTF(world->outfile,
359 " using Markov Chain Monte Carlo simulation\n ");
360
361 for (i = 0; i < len; i++)
362 FPRINTF(world->outfile, "=");
363 FPRINTF(world->outfile, "\n Version %s\n", MIGRATEVERSION);
364 get_time(nowstr, "%c");
365 if (nowstr[0] != '\0') {
366 FPRINTF(world->outfile, "\n Program started at %s\n", nowstr);
367 filepos = ftell(world->outfile);
368 /*
369 * DO NOT REMOVE the blanks in the following print
370 * statement, they are overwritten at the program end
371 * with the timestamp
372 */
373 FPRINTF(world->outfile,
374 " \n");
375 }
376 FPRINTF(world->outfile, "\n\n");
377 }
378 return filepos;
379 }
380
381 ///
382 /// main menu display, if counter is bigger than 100 the program will exit, to
383 /// make sure that batch jobs are not filling up the available harddisk space
384 void
get_menu(option_fmt * options,world_fmt * world,data_fmt * data)385 get_menu(option_fmt * options, world_fmt *world, data_fmt *data)
386 {
387 long counter=0L;
388 char input[LINESIZE];
389 char *datatype;
390 datatype = (char *) mycalloc(LINESIZE, sizeof(char));
391 if (options->menu) {
392 setup_datatype(datatype, options);
393 do {
394 printf(" Settings for this run:\n");
395 printf(" D Data type currently set to: %-30.30s\n", datatype);
396 printf(" I Input/Output formats\n");
397 printf(" P Parameters [start, migration model]\n");
398 printf(" S Search strategy\n");
399 printf(" W Write a parmfile\n");
400 printf(" Q Quit the program\n");
401 printf("\n\n");
402 printf(" To change the settings type the letter for the menu to change\n");
403 printf(" Start the program with typing Yes or Y\n===> ");fflush(stdout);
404 fflush(stdout); FGETS(input, LINESIZE, stdin);
405 printf("\n");
406 switch (uppercase(input[0])) {
407 case 'D':
408 menuData(options, datatype);
409 break;
410 case 'I':
411 menuInput(options);
412 break;
413 case 'P':
414 menuParameters(options);
415 break;
416 case 'S':
417 menuStrategy(options);
418 break;
419 case 'W':
420 save_parmfile(options, world, data);
421 break;
422 case 'Q':
423 #ifdef MPI
424 MPI_Finalize ();
425 #endif /*MPI*/
426 exit(0);
427 default:
428 break;
429 }
430 if(counter++ > 100)
431 {
432 warning("Batch job uses menu ON option, and gets non-sensical input, set menu=NO in parmfile\nProgram continues but check!!!!!");
433 input[0]='Y';
434 break;
435 }
436 }
437 while (uppercase(input[0]) != 'Y');
438 if (options->usertree && !strchr(SEQUENCETYPES, options->datatype))
439 options->usertree = FALSE;
440 if (options->dist && !strchr(SEQUENCETYPES, options->datatype))
441 options->dist = FALSE;
442 //prevents incompatability
443 // of tree and distance designation, only for
444 //sequence data the datalines and tree - tips can match
445 // msats or EP data is ambiguos concerning the tree.
446 }
447 if (options->datatype == 'g')
448 //prevents overwriting sumfile
449 options->writesum = FALSE;
450
451 #ifdef UEP
452
453 if (options->uep)
454 options->randomtree = TRUE;
455 #endif
456 myfree(datatype);
457 }
458 /*
459 * private
460 * functions---------------------------------------------------
461 */
462
setup_categs(option_fmt * options)463 boolean setup_categs(option_fmt * options) {
464 boolean didchangecat = FALSE;
465 if (options->categs > 1)
466 didchangecat = TRUE;
467 else
468 didchangecat = FALSE;
469 return didchangecat;
470 }
471
setup_rcategs(option_fmt * options)472 boolean setup_rcategs(option_fmt * options) {
473 boolean didchangercat = FALSE;
474 if (options->rcategs > 1)
475 didchangercat = TRUE;
476 else
477 didchangercat = FALSE;
478
479 return didchangercat;
480 }
481
setup_datatype(char * datatype,option_fmt * options)482 void setup_datatype(char *datatype, option_fmt * options) {
483 switch (options->datatype) {
484 case 'a':
485 strcpy(datatype, "'infinite' allele model");
486 break;
487 case 'm':
488 if(options->msat_option == SINGLESTEP)
489 strcpy(datatype, "Singlestep mutation model");
490 else
491 strcpy(datatype, "Multistep mutation model");
492 break;
493 case 'b':
494 strcpy(datatype, "Brownian motion model");
495 break;
496 case 's':
497 strcpy(datatype, "DNA sequence model");
498 break;
499 case 'n':
500 strcpy(datatype, "SNP model");
501 break;
502 case 'h':
503 strcpy(datatype, "SNP model (Hapmap data)");
504 break;
505 //case 'u':
506 //strcpy(datatype, "Panel-SNP (panel is population 1)");
507 //break;
508 //case 'f':
509 //strcpy(datatype, "ancestral state reconstruction method");
510 //break;
511 case 'g':
512 strcpy(datatype, "Genealogy summary");
513 options->readsum = TRUE;
514 break;
515 default:
516 if (options->readsum) {
517 options->datatype = 'g';
518 strcpy(datatype, "Genealogy summary");
519 } else {
520 options->datatype = 's';
521 strcpy(datatype, "DNA sequence model");
522 }
523 break;
524 }
525 }
526
527 ///
528 /// set up start tree
setup_starttree(char * starttree,option_fmt * options)529 void setup_starttree(char *starttree, option_fmt * options) {
530 if (options->usertree && strchr(SEQUENCETYPES, options->datatype))
531 sprintf(starttree, "is supplied in %s", options->utreefilename);
532 else {
533 if (options->dist && strchr(SEQUENCETYPES, options->datatype))
534 sprintf(starttree, "generates using %s", options->distfilename);
535 else {
536 if (options->randomtree)
537 strcpy(starttree, "random start genealogy");
538 else
539 strcpy(starttree, "UPGMA based start genealogy");
540 }
541 }
542 }
543
544
545 ///
546 /// checks whether the data type can use a start tree or not
547 /// and also allows to change the datatype
check_type_tree(char * input,option_fmt * options)548 void check_type_tree(char *input, option_fmt * options) {
549 switch (options->datatype)
550 {
551 case 'a':
552 switch (atol(input)) {
553 case DTEPTREE:
554 start_tree_method(options);
555 break;
556 case DTEPTYPE:
557 start_data_method(options);
558 break;
559 default:
560 break;
561 }
562 break;
563 case 'm':
564 switch (atol(input))
565 {
566 case DTMSATTREE:
567 start_tree_method(options);
568 break;
569 case DTMSATTYPE:
570 start_data_method(options);
571 break;
572 default:
573 break;
574 }
575 break;
576 case 'b':
577 switch (atol(input))
578 {
579 case DTBROWNTREE:
580 start_tree_method(options);
581 break;
582 case DTBROWNTYPE:
583 start_data_method(options);
584 break;
585 default:
586 break;
587 }
588 break;
589 case 's':
590 case 'n':
591 case 'h':
592 switch (atol(input))
593 {
594 case DTSEQTREE:
595 start_tree_method(options);
596 break;
597 case DTSEQTYPE:
598 start_data_method(options);
599 break;
600 default:
601 break;
602 }
603 break;
604 //case 'u':
605 //case 'f':
606 case 'g':
607 switch (atol(input)) {
608 case 1:
609 start_data_method(options);
610 break;
611 default:
612 break;
613 }
614 break;
615 default:
616 break;
617 }
618 }
619
620
621 ///
622 /// changing tipdate treatment and filename
change_tipdate(char * input,option_fmt * options)623 void change_tipdate(char *input, option_fmt * options)
624 {
625 printf(" Samples (Tips) where sampled at different dates? [No]\n===> ");
626 fflush(stdout); FGETS(input, LINESIZE, stdin);
627 if (uppercase(input[0]) == 'Y')
628 {
629 options->has_datefile = TRUE;
630 menu_get_filename(" What is the filename that contains the sample dates?", TIPDATEFILE, options->datefilename);
631 }
632 else
633 {
634 options->has_datefile = FALSE;
635 }
636 input[0] = 'X';
637 }
638
639 ///
640 /// changing mutation rate associated with tipdate
change_mutationrate(char * input,option_fmt * options)641 void change_mutationrate(char *input, option_fmt * options)
642 {
643 char *in;
644 char *tmp;
645 long count=0;
646 printf(" Give the number of loci or\n when all loci have the same mutation rate, the rate\nor a C to cancel\n===> ");
647 fflush(stdout); FGETS(input, LINESIZE, stdin);
648 if (uppercase(input[0]) == 'C')
649 return;
650 if (input[0] == '0')
651 {
652 options->mutationrate_year_numalloc = 1;
653 options->mutationrate_year[0] = atof(input);
654 }
655 else
656 {
657 options->mutationrate_year_numalloc = atol(input);
658 options->mutationrate_year = (MYREAL *) myrealloc(options->mutationrate_year,
659 sizeof(MYREAL)*options->mutationrate_year_numalloc);
660 while(count < options->mutationrate_year_numalloc)
661 {
662 printf(" Give the mutation rate per year for each locus\n[either in groups or singly use a space to separate]\n===> ");
663 fflush(stdout); FGETS(input, LINESIZE, stdin);
664 in = input;
665 tmp = strsep(&in,";, ");
666 while(tmp != NULL)
667 {
668 options->mutationrate_year[count++] = atof(tmp);
669 if(in==NULL)
670 break;
671 tmp = strsep(&in,";, ");
672 }
673 }
674 }
675 input[0] = 'X';
676 }
677
678 ///
679 /// changing generation time associated with tipdates
change_generationtime(char * input,option_fmt * options)680 void change_generationtime(char *input, option_fmt * options)
681 {
682 printf(" Give the number of generation per year\nor a C to cancel\n===> ");
683 fflush(stdout); FGETS(input, LINESIZE, stdin);
684 if (uppercase(input[0]) == 'C')
685 return;
686 options->generation_year = atof(input);
687 input[0] = 'X';
688 }
689
690 ///
691 /// changing hwo many individuals are used for the run, subsetting the data randomly
change_randomsubset(char * input,option_fmt * options)692 void change_randomsubset(char *input, option_fmt * options)
693 {
694 printf(" How may individuals should be used\nGive a number or A for all or a C to cancel\n===> ");
695 fflush(stdout); FGETS(input, LINESIZE, stdin);
696 if (uppercase(input[0]) == 'C')
697 return;
698 if (uppercase(input[0]) == 'A')
699 {
700 options->randomsubset = 0;
701 options->randomsubsetseed = 0;
702 }
703 else
704 {
705 options->randomsubset = (float) atof(input);
706 printf(" If you want to generate a subset that be regenerated\nin another run of migrate you need to give\na random number seed (a number between 1 and 2147483648),\nif you give the same numbers in other runs migrate will pick the same set\n===> ");
707 fflush(stdout); FGETS(input, LINESIZE, stdin);
708 options->randomsubsetseed=atol(input);
709 }
710 input[0] = 'X';
711 }
712
713 ///
714 /// changing inheritance scalar for each locus
change_inheritance(char * input,option_fmt * options)715 void change_inheritance(char *input, option_fmt * options)
716 {
717 char *in;
718 char *tmp;
719 long count=0;
720 printf(" Inheritance scalars: Give the number of loci or\nor a C to cancel\n===> ");
721 fflush(stdout); FGETS(input, LINESIZE, stdin);
722 if (uppercase(input[0]) == 'C')
723 return;
724
725 options->inheritance_scalars_numalloc = atol(input);
726 options->inheritance_scalars = (MYREAL *) myrealloc(options->inheritance_scalars,
727 sizeof(MYREAL)*options->inheritance_scalars_numalloc);
728 while(count < options->inheritance_scalars_numalloc)
729 {
730 printf(" Give the inheritance scaler for each locus\n[either in groups or singly use a space to separate]\n===> ");
731 fflush(stdout); FGETS(input, LINESIZE, stdin);
732 in = input;
733 tmp = strsep(&in,";, ");
734 while(tmp != NULL)
735 {
736 options->inheritance_scalars[count++] = atof(tmp);
737 if(in==NULL)
738 break;
739 tmp = strsep(&in,";, ");
740 }
741 }
742 input[0] = 'X';
743 }
744
745 ///
746 ///let the user set all data related options
747 void
menuData(option_fmt * options,char datatype[])748 menuData(option_fmt * options, char datatype[]) {
749 static boolean didchangecat, didchangercat;
750
751 long i;
752 long z = 0;
753 long numchar = 0;
754
755 char input[LINESIZE];
756 char starttree[LINESIZE];
757
758 char *ostr = (char *) mycalloc(LINESIZE, sizeof(char));
759
760 didchangecat = setup_categs(options);
761 didchangercat = setup_rcategs(options);
762
763 do {
764 setup_datatype(datatype, options);
765 setup_starttree(starttree, options);
766 printf(" DATATYPE AND DATA SPECIFIC OPTIONS\n\n");
767 if(options->datatype=='m')
768 printf(" %2i change Datatype, currently set to: Stepwise mutation model (%10.10s)\n", DTSEQTYPE,submodeltype(options->msat_option));
769 else
770 printf(" %2i change Datatype, currently set to:%29.29s\n", DTSEQTYPE, datatype);
771 switch (options->datatype) {
772 // allelic data --------------------------------------------------------------------
773 case 'a':
774 printf(" %2i Discard missing data %37.37s\n", DTEPMISS, options->include_unknown ? " NO" : "YES");
775 printf(" %2i Start genealogy %37.37s\n", DTEPTREE, starttree);
776
777 printf(" %2i Inheritance scalar set %3.3s\n", DTEPINHERITANCE,
778 options->inheritance_scalars_numalloc>1 ? "YES" : " NO");
779 if(options->randomsubset > 0)
780 {
781 if(options->randomsubsetseed>0)
782 printf(" %2i Pick random subset of indiv. per pop. with random seed %4li:%8li\n", DTEPRANDOMSUBSET, options->randomsubset,options->randomsubsetseed);
783 else
784 printf(" %2i Pick random subset of individuals per population %4li\n", DTEPRANDOMSUBSET, options->randomsubset);
785 }
786 else
787 printf(" %2i Pick random subset per population of individuals NO\n", DTEPRANDOMSUBSET);
788
789 if(options->has_datefile)
790 {
791 printf(" %2i Tip date file %49.49s\n", DTEPTIPDATE, options->datefilename);
792 if(options->mutationrate_year_numalloc > 1)
793 sprintf(ostr,"multiple rates");
794 else
795 sprintf(ostr,"%.12f",options->mutationrate_year[0]);
796 printf(" %2i Mutation rate per locus and year %30.30s\n", DTEPMUTRATE, ostr);
797 sprintf(ostr,"%10.4f",options->generation_year);
798 printf(" %2i How many generations per year %30.30s\n", DTEPGENERATION, ostr);
799 }
800 else
801 {
802 printf(" %2i Tip date file None, all tips a contemporary\n", DTEPTIPDATE);
803 }
804 break;
805 // brownian motion data --------------------------------------------------------------------
806 case 'b':
807 printf(" %2i Discard missing data: %37.37s\n", DTBROWNMISS, options->include_unknown ? " NO" : "YES");
808 printf(" %2i Start genealogy %37.37s\n", DTBROWNTREE, starttree);
809
810 printf(" %2i Inheritance scalar set %3.3s\n", DTBROWNINHERITANCE,
811 options->inheritance_scalars_numalloc>1 ? "YES" : " NO");
812 if(options->randomsubset > 0)
813 {
814 if(options->randomsubsetseed>0)
815 printf(" %2i Pick random subset of indiv. per pop. with random seed %4li:%8li\n", DTBROWNRANDOMSUBSET, options->randomsubset,options->randomsubsetseed);
816 else
817 printf(" %2i Pick random subset of individuals per population %4li\n", DTBROWNRANDOMSUBSET, options->randomsubset);
818 }
819 else
820 printf(" %2i Pick random subset per population of individuals NO\n", DTBROWNRANDOMSUBSET);
821 if(options->has_datefile)
822 {
823 printf(" %2i Tip date file %37.37s\n", DTBROWNTIPDATE, options->datefilename);
824 if(options->mutationrate_year_numalloc > 1)
825 sprintf(ostr,"multiple rates");
826 else
827 sprintf(ostr,"%.12f",options->mutationrate_year[0]);
828 printf(" %2i Mutation rate per locus and year %30.30s\n", DTBROWNMUTRATE, ostr);
829 sprintf(ostr,"%10.4f",options->generation_year);
830 printf(" %2i Number of generations per year %30.30s\n",DTBROWNGENERATION, ostr);
831 }
832 else
833 {
834 printf(" %2i Tip date file None, all tips a contemporary\n", DTBROWNTIPDATE);
835 }
836 break;
837 case 'm':
838 printf(" %2i Discard missing data: %37.37s\n", DTMSATMISS, options->include_unknown ? " NO" : "YES");
839 printf(" %2i Threshold value: %10li\n",
840 DTMSATTRESH, options->micro_threshold);
841 printf(" %2i Start genealogy %37.37s\n", DTMSATTREE, starttree);
842
843 printf(" %2i Inheritance scalar set %3.3s\n", DTMSATINHERITANCE,
844 options->inheritance_scalars_numalloc>1 ? "YES" : " NO");
845 if(options->randomsubset > 0)
846 {
847 if(options->randomsubsetseed>0)
848 printf(" %2i Pick random subset of indiv. per pop. with random seed %4li:%8li\n", DTMSATRANDOMSUBSET, options->randomsubset,options->randomsubsetseed);
849 else
850 printf(" %2i Pick random subset of individuals per population %4li\n", DTMSATRANDOMSUBSET, options->randomsubset);
851 }
852 else
853 printf(" %2i Pick random subset per population of individuals NO\n", DTMSATRANDOMSUBSET);
854
855
856 if(options->has_datefile)
857 {
858 printf(" %2i Tip date file %37.37s\n", DTMSATTIPDATE, options->datefilename);
859 if(options->mutationrate_year_numalloc > 1)
860 sprintf(ostr,"multiple rates");
861 else
862 sprintf(ostr,"%.12f",options->mutationrate_year[0]);
863 printf(" %2i Mutation rate per locus and year %30.30s\n", DTMSATMUTRATE, ostr);
864 sprintf(ostr,"%10.4f",options->generation_year);
865 printf(" %2i Number of generations per year %30.30s\n", DTMSATGENERATION, ostr);
866 }
867 else
868 {
869 printf(" %2i Tip date file None, all tips a contemorary\n", DTMSATTIPDATE);
870 }
871
872 break;
873 case 'u':
874 case 'n':
875 case 'h':
876 case 's':
877 case 'f':
878 z = 0;
879 numchar = 0;
880 while (options->ttratio[z] > 0.0)
881 numchar += sprintf(ostr + numchar, "%8.4f ", options->ttratio[z++]);
882 printf(" %2i Transition/transversion ratio: %31s\n", DTSEQTRATIO, ostr);
883 printf(" %2i Use empirical base frequencies? %30s\n", DTSEQFREQ, (options->freqsfrom ? "YES" : "NO"));
884 printf(" %2i Fixed categories for each site? %30.30s\n", DTSEQSITECATEGS, options->categs == ONECATEG ?
885 "One category" :
886 "Categories supplied in file");
887 printf(" %2i Site rate variation?", DTSEQRATES);
888 if (options->rcategs == 1)
889 printf(" YES\n");
890 else {
891 printf(" %4li categories of regions\n", options->rcategs);
892 printf(" %2i Rates at adjacent sites correlated?", DTSEQCORR);
893 if (!options->autocorr)
894 printf(" NO, they are independent\n");
895 else
896 printf("YES, mean block length =%6.1f\n",
897 1.0 / options->lambda);
898 }
899 printf(" %2i Sites weighted? %36.36s\n", DTSEQWEIGHT,
900 (options->weights ? "YES" : "NO"));
901 printf(" %2i Input sequences interleaved? %34.34s\n", DTSEQINTERLEAVED,
902 (options->interleaved ? "YES" : "NO, sequential"));
903 printf(" %2i Sequencing error rate? [0.0 = no error] %4.3f\n", DTSEQERROR,
904 options->seqerror);
905 printf(" %2i Slow but safer Data likelihood calculation %20.20s\n", DTSEQFAST,
906 options->fastlike ? "NO" : "YES");
907 printf(" %2i Start genealogy %37.37s\n", DTSEQTREE, starttree);
908
909 printf(" %2i Inheritance scalar set %-3.3s\n", DTSEQINHERITANCE,
910 options->inheritance_scalars_numalloc>1 ? "YES" : " NO");
911 if(options->randomsubset > 0)
912 {
913 if(options->randomsubsetseed>0)
914 printf(" %2i Pick random subset of indiv. per pop. with random seed %4li:%8li\n", DTSEQRANDOMSUBSET, options->randomsubset,options->randomsubsetseed);
915 else
916 printf(" %2i Pick random subset of individuals per population %4li\n", DTSEQRANDOMSUBSET, options->randomsubset);
917 }
918 else
919 printf(" %2i Pick random subset per population of individuals NO\n", DTSEQRANDOMSUBSET);
920
921
922 if(options->has_datefile)
923 {
924 printf(" %2i Tip date file %37.37s\n", DTSEQTIPDATE, options->datefilename);
925 if(options->mutationrate_year_numalloc > 1)
926 sprintf(ostr,"multiple rates");
927 else
928 sprintf(ostr,"%.12f",options->mutationrate_year[0]);
929 printf(" %2i Mutation rate per locus and year %30.30s\n", DTSEQMUTRATE, ostr);
930 sprintf(ostr,"%10.4f",options->generation_year);
931 printf(" %2i Number of generations per year %30.30s\n", DTSEQGENERATION, ostr);
932 }
933 else
934 {
935 printf(" %2i Tip date file None, all tips a contemorary\n", DTSEQTIPDATE);
936 }
937 break;
938
939 case 'g':
940 printf(" [Reanalyze an old run]\n");
941 break;
942 }
943 printf("\n\n");
944 printf(" Are the settings correct?\n");
945 printf(" (Type Y or the number of the entry to change)\n===> ");
946 input[0] = '\0';
947 fflush(stdout); FGETS(input, LINESIZE, stdin);
948 if (strlen(input) == 0)
949 continue;
950 check_type_tree(input, options);
951 //start tree option for all and switch datatype
952 if (strchr("ab", options->datatype))
953 // this only works if brownian and allele datatype have the same number of options
954 {
955 switch (atol(input))
956 {
957 case DTEPMISS: /* discard unknowns */
958 printf(" Discard missing Alleles (YES|NO)?\n");
959 printf(" [Default is YES, this is the best choice for most situations]\n===> ");
960 fflush(stdout); FGETS(input, LINESIZE, stdin);
961 if (strchr("Yy", input[0]))
962 options->include_unknown = FALSE;
963 else
964 options->include_unknown = TRUE;
965 input[0] = '\0';
966 break;
967 case DTEPINHERITANCE:
968 change_inheritance(input,options);
969 break;
970 case DTEPRANDOMSUBSET:
971 change_randomsubset(input,options);
972 break;
973 case DTEPTIPDATE:
974 change_tipdate(input,options);
975 break;
976 case DTEPMUTRATE:
977 change_mutationrate(input,options);
978 break;
979 case DTEPGENERATION:
980 change_generationtime(input,options);
981 break;
982 default:
983 break;
984 }
985 }
986 if (options->datatype == 'm')
987 {
988 switch (atol(input)) {
989 case DTMSATMISS: /* discard unknowns */
990 printf(" Discard missing Alleles?\n");
991 printf(" [Default is YES, this is the best choice for most situations]\n===> ");
992 fflush(stdout); FGETS(input, LINESIZE, stdin);
993 if (strchr("Yy", input[0]))
994 options->include_unknown = TRUE;
995 else
996 options->include_unknown = FALSE;
997 break;
998 case DTMSATTRESH: /* micro-threshold */
999 printf(" What is the threshold value? [needs to be even!]\n");
1000 printf
1001 (" E.g. if your allele is 24 and the threshold is 10\n");
1002 printf(" there is some probability that the allele 24 can\n");
1003 printf
1004 (" change to allele 14 (or 38), but there is a probability\n");
1005 printf(" of 0.0 (ZERO) to go to 13 (39),\n");
1006 printf
1007 (" if you choose this too small, than the program will fail\n");
1008 printf
1009 (" The default is set 100, if the biggest difference in the data is smaller the value\n");
1010 printf(" will be adjust to that maximal difference\n");
1011 printf
1012 (" [the bigger the longer the run and the more\n accurate the estimate]\n===> ");
1013 fflush(stdout); FGETS(input, LINESIZE, stdin);
1014 options->micro_threshold = atol(input);
1015 if(options->micro_threshold % 2 != 0)
1016 options->micro_threshold += 1;
1017 input[0] = '\0';
1018 break;
1019 case DTMSATTIPDATE:
1020 change_tipdate(input,options);
1021 break;
1022 case DTMSATMUTRATE:
1023 change_mutationrate(input,options);
1024 break;
1025 case DTMSATGENERATION:
1026 change_generationtime(input,options);
1027 break;
1028 case DTMSATINHERITANCE:
1029 change_inheritance(input,options);
1030 break;
1031 case DTMSATRANDOMSUBSET:
1032 change_randomsubset(input,options);
1033 break;
1034 default:
1035 break;
1036 }
1037 }
1038 if (strchr(SEQUENCETYPES, options->datatype)) {
1039 switch (atol(input)) {
1040 case DTSEQTRATIO:
1041 initratio(options);
1042 break;
1043 case DTSEQFREQ:
1044 options->freqsfrom = !options->freqsfrom;
1045 if (!options->freqsfrom) {
1046 initfreqs(&options->freqa, &options->freqc,
1047 &options->freqg, &options->freqt);
1048 }
1049 break;
1050 case DTSEQSITECATEGS:
1051 if (options->categs == ONECATEG) {
1052 options->categs = MANYCATEGS;
1053 didchangecat = TRUE;
1054 } else {
1055 options->categs = ONECATEG;
1056 didchangecat = FALSE;
1057 }
1058 break;
1059 case DTSEQRATES:
1060 if (didchangercat) {
1061 options->autocorr = FALSE;
1062 didchangercat = FALSE;
1063 options->rcategs = ONECATEG;
1064 } else {
1065 printf("\n Regional rates:\n");
1066 initcatn(&options->rcategs);
1067 options->probcat =
1068 (MYREAL *) myrealloc(options->probcat,
1069 options->rcategs * sizeof(MYREAL));
1070 options->rrate =
1071 (MYREAL *) myrealloc(options->rrate,
1072 options->rcategs * sizeof(MYREAL));
1073 didchangercat = TRUE;
1074 options->gammarates =
1075 initcategs(options->rcategs, options->rrate,
1076 options->probcat);
1077 if (!options->gammarates)
1078 {
1079 initprobcat(options->rcategs, &options->probsum,
1080 options->probcat);
1081 constrain_rates(options->rcategs, options->rrate,
1082 options->probcat);
1083 }
1084 FPRINTF(stdout,
1085 "\n\nRegion type Rate of change Probability\n");
1086 FPRINTF(stdout,
1087 "---------------------------------------------\n");
1088 for (i = 0; i < options->rcategs; i++)
1089 FPRINTF(stdout, "%9ld%16.3f%17.3f\n",
1090 i + 1, options->rrate[i], options->probcat[i]);
1091 FPRINTF(stdout, "\n");
1092 }
1093 break;
1094 case DTSEQCORR:
1095 options->autocorr = !options->autocorr;
1096 if (options->autocorr)
1097 initlambda(options);
1098 break;
1099 case DTSEQWEIGHT:
1100 options->weights = !options->weights;
1101 break;
1102 case DTSEQINTERLEAVED:
1103 options->interleaved = !options->interleaved;
1104 break;
1105 case DTSEQERROR:
1106 FPRINTF(stdout,
1107 "Enter the sequencing error rate per site\n[Good values are 0.0 (=no error), or 0.01 (1%% error):\n===> ");
1108 fflush(stdout); FGETS(input, LINESIZE, stdin);
1109 options->seqerror = atof(input);
1110 input[0] = '\0';
1111 break;
1112 case DTSEQFAST:
1113 options->fastlike = !options->fastlike;
1114 break;
1115 case DTSEQTIPDATE:
1116 change_tipdate(input,options);
1117 break;
1118 case DTSEQMUTRATE:
1119 change_mutationrate(input,options);
1120 break;
1121 case DTSEQGENERATION:
1122 change_generationtime(input,options);
1123 break;
1124 case DTSEQINHERITANCE:
1125 change_inheritance(input,options);
1126 break;
1127 case DTSEQRANDOMSUBSET:
1128 change_randomsubset(input,options);
1129 break;
1130 default:
1131 break;
1132 }
1133 if (!didchangercat) {
1134 options->probcat =
1135 (MYREAL *) myrealloc(options->probcat, sizeof(MYREAL) * 2);
1136 options->rrate =
1137 (MYREAL *) myrealloc(options->rrate, sizeof(MYREAL) * 2);
1138 options->rrate[0] = 1.0;
1139 options->probcat[0] = 1.0;
1140 }
1141 if (!didchangecat) {
1142 options->rate =
1143 (MYREAL *) myrealloc(options->rate, sizeof(MYREAL) * 2);
1144 options->rate[0] = 1.0;
1145 }
1146 }
1147 if (options->datatype != 'g')
1148 options->readsum = FALSE;
1149 }
1150 while (uppercase(input[0]) != 'Y');
1151
1152 myfree(ostr);
1153
1154 }
1155
1156 void
menuInput(option_fmt * options)1157 menuInput(option_fmt * options) {
1158 long test = 0;
1159 long len = (long) strlen(options->title);
1160 unsigned long timeseed;
1161 char input[LINESIZE];
1162 char treeinc[LINESIZE];
1163 char outputstring[LINESIZE];
1164 char *stringstep, *string2, *string3, *string4;
1165 char treestr[4][20] = {"None", "All", "Best", "Last chain"};
1166 char progress[3][8] = {"Verbose", "YES", "NO"};
1167 char *extension;
1168 int retval;
1169 stringstep = (char *) mymalloc(sizeof(char) * 128);
1170 string2 = (char *) mymalloc(sizeof(char) * 128);
1171 string3 = (char *) mymalloc(sizeof(char) * 128);
1172 string4 = (char *) mymalloc(sizeof(char) * 128);
1173
1174 do {
1175 printf(" INPUT/OUTPUT FORMATS (for %s)\n -------------------- [approach can be changed in SEARCH Strategy]\n\n INPUT:\n",options->bayes_infer ? "Bayesian approach" : "Maximum likelihood approach");
1176 printf(" %2i Datafile name is %46.46s\n", MI_INFILE,
1177 options->infilename);
1178
1179 sprintf(treeinc,"[every %li]",options->treeinc);
1180
1181 switch (options->autoseed) {
1182 case AUTO:
1183 sprintf(stringstep, "YES");
1184 break;
1185 case NOAUTO:
1186 sprintf(stringstep, "NO, use seedfile");
1187 break;
1188 case NOAUTOSELF:
1189 sprintf(stringstep, "NO, seed=%li ", options->inseed);
1190 break;
1191 default:
1192 options->autoseed = AUTO;
1193 sprintf(stringstep, "YES");
1194 break;
1195 }
1196
1197 printf(" %2i Use automatic seed for randomisation? %24.24s\n",
1198 MI_RAND, stringstep);
1199 printf(" %2i Title of the analysis is", MI_TITLE);
1200
1201 if (len == 0 || (len == 5 && strstr(options->title, "Migration analysis")))
1202 printf("%39.39s\n", "<no title given>");
1203 else
1204 printf("%39.39s\n", options->title);
1205
1206 if (options->readsum && !options->bayes_infer) {
1207 printf(" %2i Summary of genealogies are read from %s\n",
1208 MI_SUMREAD, options->sumfilename);
1209 }
1210 else
1211 {
1212 if (options->readsum && options->bayes_infer) {
1213 printf(" %2i Bayesian output data are read from %s\n",
1214 MI_SUMREAD, options->bayesmdimfilename);
1215 }
1216 }
1217 printf("\n OUTPUT:\n");
1218
1219 printf(" %2i Print indications of progress of run? %25.25s\n",
1220 MI_PROGRESS, options->progress ? (options->
1221 verbose ? progress[0] :
1222 progress[1]) : progress[2]);
1223 printf(" %2i Print the data?%48.48s\n",
1224 MI_PRINTDATA, options->printdata ? "YES" : "NO");
1225
1226 printf(" %2i Outputfile name is %43.43s\n", MI_OUTFILE,
1227 options->outfilename);
1228 #ifdef PRETTY
1229 printf(" %43.43s\n", options->pdfoutfilename);
1230 #endif
1231 if(!options->bayes_infer)
1232 {
1233 if (options->plot)
1234 {
1235 switch (options->plotmethod)
1236 {
1237 case PLOTALL:
1238 strcpy(string2, "YES, to outfile and mathfile");
1239 break;
1240 default:
1241 strcpy(string2, "YES, to outfile");
1242 break;
1243 }
1244 }
1245 else
1246 {
1247 strcpy(string2, "NO");
1248 }
1249 printf(" %2i Plot likelihood surface? %38.38s\n",
1250 MI_PLOT, string2);
1251 switch (options->profile)
1252 {
1253 case _NONE:
1254 strcpy(string3, "NO");
1255 break;
1256 case ALL:
1257 strcpy(string3, "YES, tables and summary");
1258 break;
1259 case TABLES:
1260 strcpy(string3, "YES, tables");
1261 break;
1262 case SUMMARY:
1263 strcpy(string3, "YES, summary");
1264 break;
1265 }
1266 printf(" %2i Profile-likelihood?%44.44s\n", MI_PROFILE, string3);
1267 if (options->profile != _NONE)
1268 {
1269 switch (options->profilemethod)
1270 {
1271 case 'p':
1272 strcpy(string4, "[Percentiles using exact Bisection method]");
1273 break;
1274 case 's':
1275 strcpy(string4,
1276 "[Percentiles using not so exact Spline method]");
1277 break;
1278 case 'd':
1279 strcpy(string4,
1280 "[Evaluation at preset fractions of ML estimate]");
1281 break;
1282 case 'q':
1283 strcpy(string4, "[Assuming that parameter are independent]");
1284 break;
1285 case 'f':
1286 strcpy(string4, "[Mixture of 'p' and 'q']");
1287 break;
1288 }
1289 printf("%70.70s\n", string4);
1290 }
1291 printf(" %2i Likelihood-Ratio tests?%40.40s\n",
1292 MI_LRT, options->lratio->counter > 0 ? " YES" : "NO");
1293
1294 printf(" %2i AIC model selection? %40.40s\n",
1295 MI_AIC, options->aic ? (options->fast_aic ? "YES:Fast" : "YES") : "NO");
1296 } /*if !bayes_infer*/
1297
1298 switch (options->treeprint)
1299 {
1300 case ALL:
1301 printf(" %2i Print genealogies? %40.40s\n",
1302 MI_TREES, treestr[1]);
1303 break;
1304 case BEST:
1305 printf(" %2i Print genealogies? %40.40s\n",
1306 MI_TREES, treestr[2]);
1307 break;
1308 case LASTCHAIN:
1309 printf(" %2i Print genealogies? %19.19s %20.20s\n",
1310 MI_TREES, treeinc ,treestr[3]);
1311 break;
1312 case _NONE:
1313 default:
1314 printf(" %2i Print genealogies? %40.40s\n",
1315 MI_TREES, treestr[0]);
1316 break;
1317 }
1318 if(!options->bayes_infer)
1319 {
1320 printf(" %2i Plot coordinates are saved in %33.33s\n",
1321 MI_PLOTCOORD, options->mathfilename);
1322 if (options->writesum || options->datatype == 'g')
1323 {
1324 if (!options->readsum && options->writesum)
1325 {
1326 printf(" %2i Summary of genealogies saved in %31.31s\n",
1327 MI_SUMFILE, options->sumfilename);
1328 }
1329 }
1330 else
1331 {
1332 printf(" %2i Summary of genealogies %39.39s\n",
1333 MI_SUMFILE, "will not be saved");
1334 }
1335 }
1336 if (options->writelog)
1337 printf(" %2i Save logging information into %33.33s\n", MI_LOGFILE,
1338 options->logfilename);
1339 else
1340 printf(" %2i Save logging information? %36.36s\n", MI_LOGFILE, "NO");
1341 #ifdef UEP
1342
1343 if (options->uep) {
1344 printf(" %2i Read unique polymorphism? From file %30.30s\n",
1345 MI_UEPFILE, options->uepfilename);
1346 printf(" %2i Mutation rate for UEP is %10.5f %s\n",
1347 MI_UEPRATE, options->ueprate, " x mutation rate");
1348 printf(" %2i Base frequency for UEP is \"0\"=%f, \"1\"=%f\n",
1349 MI_UEPFREQ, options->uepfreq0, options->uepfreq1);
1350 } else
1351 printf(" %2i Read unique polymorphism? %20.20s\n", MI_UEPFILE, "NO");
1352 #endif
1353 if (options->mighist)
1354 {
1355 sprintf(outputstring,"%s %s", options->mighistfilename,
1356 options->mighist_all ? "(all events)" : "(migration events)");
1357 printf(" %2i Show event statistics%42.42s\n", MI_MIGHISTOGRAM, outputstring);
1358 if(options->mighist_increment > 1)
1359 sprintf(outputstring,"every %li %s", options->mighist_increment,"sample steps");
1360 else
1361 sprintf(outputstring,"every %s", "sample step");
1362 printf(" Events are recorded every %32.32s\n",outputstring);
1363 sprintf(outputstring,"%f",options->eventbinsize);
1364 printf(" Histogram bin width %32.32s\n",outputstring);
1365
1366 }
1367 else
1368 {
1369 sprintf(outputstring,"NO");
1370 printf(" %2i Show event statistics %32.32s\n", MI_MIGHISTOGRAM, outputstring);
1371
1372 }
1373 if (options->skyline)
1374 {
1375 remove_trailing_blanks(&options->skylinefilename);
1376 sprintf(outputstring,"%s", options->skylinefilename);
1377 printf(" %2i Record parameter change through time?%26.26s\n", MI_SKYLINE, outputstring);
1378 sprintf(outputstring,"%f",options->eventbinsize);
1379 printf(" Histogram bin width %32.32s\n",outputstring);
1380 }
1381 else
1382 {
1383 sprintf(outputstring,"NO");
1384 printf(" %2i Record parameter change through time?%26.26s\n", MI_SKYLINE, outputstring);
1385 }
1386
1387 printf("\n\n Are the settings correct?\n");
1388 printf
1389 (" (type Y to go back to the main menu or the letter for the entry to change)\n===> ");
1390 fflush(stdout); FGETS(input, LINESIZE, stdin);
1391 test = atoi(input);
1392 switch (test)
1393 {
1394 case MI_INFILE:
1395 printf(" What is the datafile name?\n[Default: %s]\n===> ", options->infilename);
1396 fflush(stdout); FGETS(input, LINESIZE, stdin);
1397 if (input[0] != '\0')
1398 {
1399 strcpy(options->infilename, input);
1400 }
1401 break;
1402 case MI_RAND:
1403 do
1404 {
1405 printf(" (A)utomatic or (S)eedfile or (O)wn\n");
1406 printf(" Start value for Random-generator seed\n===> ");
1407 fflush(stdout); FGETS(input, LINESIZE, stdin);
1408 switch (uppercase(input[0]))
1409 {
1410 case 'A':
1411 options->autoseed = AUTO;
1412 #ifndef MAC
1413 timeseed = (unsigned long) time(NULL) / 4;
1414 #else
1415 timeseed = (unsigned long) clock() / 4;
1416 #endif
1417 options->inseed = (long) timeseed + 1;
1418 break;
1419 case 'S':
1420 menu_get_filename(" What is the filename that contains the random number seed?", SEEDFILE, options->seedfilename);
1421 openfile(&options->seedfile, options->seedfilename, "r", NULL);
1422 if (options->seedfile) {
1423 options->autoseed = NOAUTO;
1424 retval = fscanf(options->seedfile, "%ld%*[^\n]",
1425 &options->inseed);
1426 fclose(options->seedfile);
1427 } else
1428 printf("\n\n There is no seedfile present\n");
1429 break;
1430 case 'O':
1431 options->autoseed = NOAUTOSELF;
1432 printf(" Random number seed (best values are x/4 +1)?\n");
1433 retval = scanf("%ld%*[^\n]", &options->inseed);
1434 break;
1435 }
1436 }
1437 while (options->autoseed < AUTO || options->autoseed > NOAUTOSELF);
1438 break;
1439 case MI_TITLE:
1440 printf(" Enter a title? [max 80 Characters are used]\n===> ");
1441 fflush(stdout); FGETS(input, LINESIZE, stdin);
1442 if (input[0] == '\0')
1443 options->title[0] = '\0';
1444 else
1445 sprintf(options->title,"%80.80s", input);
1446 break;
1447 case MI_SUMREAD:
1448 if(!options->bayes_infer)
1449 {
1450 printf
1451 (" What is the filename for the summary of genealogies\n[Default: %s]\n===> ",
1452 SUMFILE);
1453 fflush(stdout); FGETS(input, LINESIZE, stdin);
1454 if (input[0] == '\0')
1455 strcpy(options->sumfilename, SUMFILE);
1456 else {
1457 strcpy(options->sumfilename, input);
1458 }
1459 }
1460 else
1461 {
1462 printf
1463 (" What is the filename of the recorded Bayes data\n[Default: %s]\n===> ",
1464 BAYESMDIMFILE);
1465 fflush(stdout); FGETS(input, LINESIZE, stdin);
1466 if (input[0] == '\0')
1467 strcpy(options->bayesmdimfilename, BAYESMDIMFILE);
1468 else {
1469 strcpy(options->bayesmdimfilename, input);
1470 }
1471 unpad(options->bayesmdimfilename, " ");
1472 extension = strrchr(options->bayesmdimfilename,'.');
1473 #ifdef ZNZ
1474 if(extension!=NULL && !strncmp(extension,".gz",3))
1475 {
1476 options->use_compressed = 1;
1477 }
1478 else
1479 {
1480 options->use_compressed = 0;
1481 }
1482 #else
1483 options->use_compressed = 0;
1484 #endif
1485 }
1486 break;
1487 case MI_PROGRESS:
1488 printf(" Progress report during the run? <YES | Verbose | NO>\n===> ");
1489 fflush(stdout); FGETS(input, LINESIZE, stdin);
1490 switch (tolower(input[0])) {
1491 case 'n':
1492 options->progress = FALSE;
1493 options->verbose = FALSE;
1494 break;
1495 case 'v':
1496 options->progress = TRUE;
1497 options->verbose = TRUE;
1498 break;
1499 case '\0':
1500 case 'y':
1501 default:
1502 options->progress = TRUE;
1503 options->verbose = FALSE;
1504 break;
1505 }
1506 input[0] = 'X';
1507 break;
1508 case MI_PRINTDATA:
1509 options->printdata = !options->printdata;
1510 break;
1511 case MI_OUTFILE:
1512 menu_get_filename(" What is the output filename?", OUTFILE, options->outfilename);
1513 #ifdef PRETTY
1514 strcpy(options->pdfoutfilename,options->outfilename);
1515 strcat(options->pdfoutfilename,".pdf");
1516 #endif
1517 break;
1518 case MI_PLOT:
1519 options->plot = !options->plot;
1520 if (options->plot)
1521 {
1522 get_plotoptions(options);
1523 input[0] = '\0';
1524 }
1525 break;
1526 case MI_PROFILE:
1527 do
1528 {
1529 printf(" Evaluate profile likelihoods:\n");
1530 printf
1531 (" (N)o, (A)all [Tables and Summary],\n (T)ables, (S)ummary\n===> ");
1532 fflush(stdout); FGETS(input, LINESIZE, stdin);
1533 }
1534 while (strchr("NATS", (int) uppercase(input[0])) == NULL);
1535 switch (uppercase(input[0]))
1536 {
1537 case 'N':
1538 options->profile = _NONE;
1539 break;
1540 case 'A':
1541 options->profile = ALL;
1542 break;
1543 case 'T':
1544 options->profile = TABLES;
1545 break;
1546 case 'S':
1547 options->profile = SUMMARY;
1548 break;
1549 }
1550 if (options->profile != _NONE)
1551 {
1552 do
1553 {
1554 printf(" Method to evaluate the profiles?\n");
1555 printf(" (P)ercentiles: exact evaluation (slow)\n");
1556 printf
1557 (" (F)ast : Assuming parameters are indepenent\n + one complete profile-maximization.\n");
1558 //printf
1559 // (" (S)plines : percentiles using splines\n (fails sometimes)\n");
1560 printf
1561 (" (D)iscrete : evaluation at (0.02, 0.05, 0.1, .. , 50)*ML estimate\n");
1562 printf
1563 (" (Q)uick : Assuming parameters are independent\n===> ");
1564 fflush(stdout); FGETS(input, LINESIZE, stdin);
1565 }
1566 while (strchr("pdqf", (int) lowercase(input[0])) == NULL);
1567 switch (lowercase(input[0]))
1568 //in strchr removed s
1569 {
1570 case 'p':
1571 options->profilemethod = 'p';
1572 break;
1573 case 's':
1574 options->profilemethod = 's';
1575 break;
1576 case 'd':
1577 options->profilemethod = 'd';
1578 break;
1579 case 'q':
1580 options->profilemethod = 'q';
1581 break;
1582 case 'f':
1583 options->profilemethod = 'f';
1584 break;
1585 default:
1586 options->profilemethod = 's';
1587 break;
1588 }
1589 }
1590 set_profile_options(options);
1591 break;
1592 case MI_LRT:
1593 if (options->lratio->counter > 0)
1594 {
1595 options->lratio->counter = 0;
1596 }
1597 else
1598 {
1599 FPRINTF(stdout, " Likelihood-Ratio tests:\n");
1600 FPRINTF(stdout, " -----------------------\n");
1601 read_custom_menu_lratio(options);
1602 }
1603 break;
1604 case MI_AIC:
1605 printf
1606 (" Use AIC to select minimal migration model?\n[Default: NO]\n===> ");
1607 fflush(stdout); FGETS(input, LINESIZE, stdin);
1608 switch (uppercase(input[0]))
1609 {
1610 case 'F':
1611 options->fast_aic = TRUE;
1612 options->aic = TRUE;
1613 break;
1614 case 'Y':
1615 options->fast_aic = FALSE;
1616 options->aic = TRUE;
1617 break;
1618 case 'N':
1619 default:
1620 options->aic = FALSE;
1621 options->fast_aic = FALSE;
1622 break;
1623 }
1624 input[0] = 'x';
1625 break;
1626 case MI_TREES:
1627 do
1628 {
1629 printf(" Print genealogies:\n");
1630 printf(" (N)one, (A)all [!], (B)est, (L)ast chain\n===> ");
1631 fflush(stdout); FGETS(input, LINESIZE, stdin);
1632 }
1633 while (strchr("NABL", (int) uppercase(input[0])) == NULL);
1634 switch (uppercase(input[0]))
1635 {
1636 case 'N':
1637 options->treeprint = _NONE;
1638 break;
1639 case 'A':
1640 options->treeprint = ALL;
1641 menu_get_filename(" What is the filename for storing recorded genealogies?", TREEFILE, options->treefilename);
1642 break;
1643 case 'B':
1644 options->treeprint = BEST;
1645 menu_get_filename(" What is the filename for storing recorded genealogies?", TREEFILE, options->treefilename);
1646 break;
1647 case 'L':
1648 options->treeprint = LASTCHAIN;
1649 menu_get_filename(" What is the filename for storing recorded genealogies?", TREEFILE, options->treefilename);
1650 printf("Give the increment to record genealogies\n===>");
1651 fflush(stdout); FGETS(input, LINESIZE, stdin);
1652 retval = sscanf(input,"%li",&options->treeinc);
1653 break;
1654 default:
1655 options->treeprint = _NONE;
1656 break;
1657 }
1658 break;
1659 case MI_SUMFILE:
1660 printf(" Save genealogy summaries? [NO]\n===> ");
1661 fflush(stdout); FGETS(input, LINESIZE, stdin);
1662 if (uppercase(input[0]) == 'Y')
1663 {
1664 options->writesum = TRUE;
1665 menu_get_filename(" What is the filename for the summary of genealogies?", SUMFILE, options->sumfilename);
1666 }
1667 else
1668 {
1669 options->writesum = FALSE;
1670 }
1671 input[0] = 'X';
1672 break;
1673 case MI_PLOTCOORD:
1674 menu_get_filename(" What is the plot coordinate filename?", MATHFILE, options->mathfilename);
1675 break;
1676 case MI_LOGFILE:
1677 printf(" Save logging information? [NO]\n===> ");
1678 fflush(stdout); FGETS(input, LINESIZE, stdin);
1679 if (uppercase(input[0]) == 'Y')
1680 {
1681 options->writelog = TRUE;
1682 menu_get_filename(" What is the filename for logging?", LOGFILE, options->logfilename);
1683 }
1684 else
1685 {
1686 options->writelog = FALSE;
1687 }
1688 input[0] = 'X';
1689 break;
1690 #ifdef UEP
1691 case MI_UEPFILE:
1692 printf(" UEP estimation? [NO]\n===> ");
1693 fflush(stdout); FGETS(input, LINESIZE, stdin);
1694 if (uppercase(input[0]) == 'Y')
1695 {
1696 options->uep = TRUE;
1697 menu_get_filename(" What is the filename of the UEP data?", UEPFILE, options->uepfilename);
1698 }
1699 else
1700 {
1701 options->uep = FALSE;
1702 }
1703 input[0] = 'X';
1704 break;
1705 case MI_UEPRATE:
1706 do
1707 {
1708 FPRINTF(stdout, "What is mutation rate for the UEP locus\n");
1709 FPRINTF(stdout,
1710 "[expressed as a ratio of the point mutation rate\n");
1711 FPRINTF(stdout, "Reasonable values are 0.0 =< x <=1.0]\n===> ");
1712 fflush(stdout); FGETS(input, LINESIZE, stdin);
1713 options->ueprate = atof(input);
1714 }
1715 while (options->ueprate < 0.0);
1716 break;
1717 case MI_UEPFREQ:
1718 do
1719 {
1720 FPRINTF(stdout, "What is the base frequency for the \"1\" allele\n");
1721 FPRINTF(stdout, "Reasonable values are 0.0 < x <1.0]\n===> ");
1722 fflush(stdout); FGETS(input, LINESIZE, stdin);
1723 options->uepfreq1 = atof(input);
1724 options->uepfreq0 = 1. - options->uepfreq1;
1725 }
1726 while (options->uepfreq1 < 0.0);
1727 break;
1728 #endif
1729 case MI_MIGHISTOGRAM:
1730 printf(" Display event statistic and save events into file?\n [options are: NO, ALL events, MIGration events]\n===> ");
1731 fflush(stdout); FGETS(input, LINESIZE, stdin);
1732 if (uppercase(input[0]) == 'N')
1733 {
1734 options->mighist = FALSE;
1735 }
1736 else
1737 {
1738 options->mighist = TRUE;
1739 if (uppercase(input[0]) == 'A')
1740 options->mighist_all = TRUE;
1741 else
1742 options->mighist_all = FALSE;
1743
1744 /* printf(" Thin the raw data and sample only every x sample steps,\n Enter increment [default is 1]:\n===> ");
1745 fflush(stdout); FGETS(input, LINESIZE, stdin);
1746 options->mighist_increment = atol(input);
1747 if(options->mighist_increment < 1)
1748 options->mighist_increment = 1;
1749 */
1750 menu_get_filename(" Specify a filename for the raw event data!", MIGHISTFILE, options->mighistfilename);
1751
1752 do
1753 {
1754 printf("Events are recorded per time intervals.\nWhat is the time interval size?\n");
1755 printf("[Current value: %f, try values that are about 10x smaller than the largest Theta]\n===> ",
1756 options->eventbinsize);
1757 fflush(stdout); FGETS(input, LINESIZE, stdin);
1758 if(input[0]=='\0' && options->eventbinsize > 0.0)
1759 break;
1760 options->eventbinsize = (float) atof(input);
1761 }
1762 while (options->eventbinsize <= 0.0);
1763
1764 }
1765 input[0] = 'X';
1766 break;
1767 case MI_SKYLINE:
1768 printf(" Display time dependent parameters and save them into file? [Choices are NO or YES]\n===> ");
1769 fflush(stdout); FGETS(input, LINESIZE, stdin);
1770 if (uppercase(input[0]) == 'N')
1771 {
1772 options->skyline = FALSE;
1773 }
1774 else
1775 {
1776 options->skyline = TRUE;
1777 do
1778 {
1779 printf("What is the time interval size?\n");
1780 printf("[Current value: %f, try values that are about 10x smaller than the largest Theta]\n===> ",
1781 options->eventbinsize);
1782 fflush(stdout); FGETS(input, LINESIZE, stdin);
1783 if(input[0]=='\0' && options->eventbinsize > 0.0)
1784 break;
1785 options->eventbinsize = (float) atof(input);
1786 }
1787 while (options->eventbinsize <= 0.0);
1788 menu_get_filename(" Specify a filename for the raw skyline output!", SKYLINEFILE, options->skylinefilename);
1789 if(!options->mighist)
1790 {
1791 printf("The skyline option needs the coalescence and migration event recording\n");
1792 menu_get_filename(" Specify a filename for the raw event data!", MIGHISTFILE, options->mighistfilename);
1793 options->mighist_all = TRUE;
1794 options->mighist_increment = 0;
1795 }
1796
1797 }
1798 input[0] = 'X';
1799 break;
1800 default:
1801 break;
1802 }
1803 }
1804 while (uppercase(input[0]) != 'Y');
1805 myfree(stringstep);
1806 myfree(string2);
1807 myfree(string3);
1808 myfree(string4);
1809 }
1810
1811 ///sets the start parameter options for theta
switch_theta(char * ostr,option_fmt * options)1812 void switch_theta(char *ostr, option_fmt * options) {
1813 long i;
1814 long numchar = 0;
1815 switch (options->thetaguess) {
1816 case FST:
1817 sprintf(ostr, "Estimate with FST (Fw/Fb) measure");
1818 break;
1819 case URANDOMESTIMATE:
1820 sprintf(ostr, "Random Theta from Uniform distribution (min=%f, max=%f)", options->thetag[0], options->thetag[1]);
1821 break;
1822 case NRANDOMESTIMATE:
1823 sprintf(ostr, "Random Theta from Normal distribution (mean=%f, std=%f)", options->thetag[0], options->thetag[1]);
1824 break;
1825 default:
1826 numchar += sprintf(ostr, "NO, initial Theta = {");
1827 for (i = 0; i < options->numthetag - 1; i++) {
1828 numchar += sprintf(ostr + numchar, "%.5f,", options->thetag[i]);
1829 }
1830 sprintf(ostr + numchar, "%.5f}", options->thetag[i]);
1831 }
1832 }
1833
switch_mig(char * ostr,option_fmt * options)1834 void switch_mig(char *ostr, option_fmt * options) {
1835 switch (options->migrguess) {
1836 case FST:
1837 sprintf(ostr, "Estimate with FST (Fw/Fb) measure");
1838 break;
1839 case URANDOMESTIMATE:
1840 sprintf(ostr, "Random migration rates from Uniform distribution (min=%f, max=%f)", options->mg[0], options->mg[1]);
1841 break;
1842 case NRANDOMESTIMATE:
1843 sprintf(ostr, "Random migration rates from Normal distribution (mean=%f, std=%f)", options->mg[0], options->mg[1]);
1844 break;
1845 default:
1846 sprintf(ostr, "NO, initial migration rates are given");
1847 }
1848 }
1849
1850
1851 ///
1852 ///allows to set all parameter related issues:start parameters, mutation type, migration model
1853 void
menuParameters(option_fmt * options)1854 menuParameters(option_fmt * options)
1855 {
1856 #ifdef POPMODEL
1857 char tmp[LINESIZE];
1858 #endif
1859 char input[LINESIZE];
1860 char custmexplain[LINESIZE];
1861 char localities[LINESIZE];
1862 char *temp;
1863 char *outputstring;
1864 long count = 0;
1865 long i, j, tt, check = 0, numpop = 0;
1866 MYREAL sum = 0.;
1867 double tmpdouble;
1868 int retval;
1869 outputstring = (char *) mycalloc(LINESIZE, sizeof(char));
1870
1871
1872 do {
1873 printf
1874 (" PARAMETERS\n ---------------------------\n Start parameters:\n");
1875 switch_theta(outputstring, options);
1876 printf(" 1 Use a simple estimate of theta as start?\n %64s\n", outputstring);
1877 switch_mig(outputstring, options);
1878 printf(" 2 Use a simple estimate of migration rate as start?\n %64s\n", outputstring);
1879 printf
1880 ("\n Gene flow parameter and Mutation rate variation among loci:\n");
1881 printf(" 3 Use M for the gene flow parameter %28.28s\n",
1882 options->usem ? "YES [M=m/mu]" : "NO [Theta * M]");
1883 if(options->gamma)
1884 {
1885 printf(" 4 Mutation rate is %46.46s%f\n","Estimated: Gamma distribution with alpha=",(float) options->alphavalue);
1886 }
1887 else
1888 {
1889 if(options->murates)
1890 {
1891 if(options->murates_fromdata)
1892 {
1893 printf(" 4 Mutation rate is %46.46s\n","Varying (from data)");
1894 }
1895 else
1896 {
1897 printf(" 4 Mutation rate is %46.46s\n","Varying (from user)");
1898 }
1899 }
1900 else
1901 {
1902 if(options->bayesmurates)
1903 {
1904 printf(" 4 Mutation rate is %46.46s\n","Estimated (see prior menu: Rate)");
1905 }
1906 else
1907 {
1908 printf(" 4 Mutation rate is %46.46s\n","Constant");
1909 }
1910 }
1911 }
1912 if (options->migrguess == FST || options->thetaguess == FST) {
1913 printf("\n FST-Calculation (for START value):\n");
1914 printf(" 5 Method: %56.56s\n",
1915 options->fsttype ==
1916 'T' ? "Variable Theta, M symmetric" :
1917 "Variable M, Theta is the same for all populations");
1918 printf(" 6 Print FST table:%48.48s\n", options->printfst ? "YES" : "NO");
1919 }
1920 strcpy(custmexplain, custom_migration_type(options->migration_model));
1921 set_localities_string(&localities[0],options);
1922 printf("\n Migration model and combination of localities:\n");
1923 printf(" 7 Sampling localities %43.43s\n", localities);
1924 printf(" 8 Model is set to %48.48s\n", custmexplain);
1925 if(options->geo)
1926 printf(" 9 Geographic distance matrix: YES:%34.34s\n", options->geofilename);
1927 else
1928 printf(" 9 Geographic distance matrix: %36.36s\n", "NO");
1929 #ifdef POPMODEL
1930 printf(" 10 Population model: %46.46s\n",
1931 which_popmodel(options,tmp));
1932 #endif /*POPMODEL*/
1933 printf("\n\n Are the settings correct?\n");
1934 printf
1935 (" (Type Y to go back to the main menu or the letter for an entry to change)\n===> ");
1936 fflush(stdout); FGETS(input, LINESIZE, stdin);
1937 switch (input[0]) {
1938 case '1':
1939 printf(" Which method? (F)st or (O)wn or (N)ormally or ((U)niform distributed random value\n===> ");
1940 fflush(stdout); FGETS(input, LINESIZE, stdin);
1941 switch (uppercase(input[0])) {
1942 case 'N':
1943 options->thetaguess = NRANDOMESTIMATE;
1944 options->thetag = (MYREAL *) myrealloc(options->thetag,
1945 sizeof(MYREAL) * 2);
1946 menuRandom(options->thetag,'N');
1947 options->numthetag = 2;
1948 break;
1949 case 'U':
1950 options->thetaguess = URANDOMESTIMATE;
1951 options->thetag = (MYREAL *) myrealloc(options->thetag,
1952 sizeof(MYREAL) * 2);
1953 menuRandom(options->thetag,'U');
1954 options->numthetag = 2;
1955 break;
1956 case 'F':
1957 options->thetaguess = FST;
1958 break;
1959 case 'O':
1960 options->thetaguess = OWN;
1961 if (numpop == 0) {
1962 printf(" I do not know yet how many populations the data set contains? Specify the number of populations > ");fflush(stdout);
1963 retval = scanf("%ld", &numpop);
1964 retval = scanf("%*[^\n]");
1965 }
1966 options->thetag =
1967 (MYREAL *) myrealloc(options->thetag,
1968 sizeof(MYREAL) * (numpop + 1));
1969 printf(" Give the initial Theta estimates:\n");
1970 for (i = 0; i < numpop; i++) {
1971 printf("Population %3li> ", i);fflush(stdout);
1972 #ifdef USE_MYREAL_FLOAT
1973 retval = scanf("%f", &options->thetag[i]);
1974 #else
1975 retval = scanf("%lf", &options->thetag[i]);
1976 #endif
1977 }
1978 options->numthetag = i;
1979 break;
1980 default:
1981 options->thetaguess = FST;
1982 break;
1983 }
1984
1985 break;
1986 case '2':
1987 printf("Which method? (F)st or (O)wn or (N)ormally or (U)niformly distributed random value\n===> ");
1988 fflush(stdout); FGETS(input, LINESIZE, stdin);
1989 switch (uppercase(input[0])) {
1990 case 'U':
1991 options->migrguess = URANDOMESTIMATE;
1992 options->mg = (MYREAL *) myrealloc(options->mg,
1993 sizeof(MYREAL) * 2);
1994 menuRandom(options->mg,'U');
1995 options->nummg = 2;
1996 break;
1997 case 'N':
1998 options->migrguess = NRANDOMESTIMATE;
1999 options->mg = (MYREAL *) myrealloc(options->mg,
2000 sizeof(MYREAL) * 2);
2001 menuRandom(options->mg,'N');
2002 options->nummg = 2;
2003 break;
2004 case 'F':
2005 options->migrguess = FST;
2006 break;
2007 case 'O':
2008 options->migrguess = OWN;
2009 if (numpop == 0) {
2010 printf(" I do not know yet how many populations the data set contains? Specify the number of populations > ");
2011 fflush(stdout);
2012 retval = scanf("%ld", &numpop);
2013 retval = scanf("%*[^\n]");
2014 }
2015 printf
2016 (" Initial migration rate estimate?\n[give 4Nm for diploid data\n 2Nm for haploid data\n Nm for mtDNA data]\n[If you use M (m/mu) instead of Nm then you need to specify M here !]");
2017 tt = 0;
2018 options->mg =
2019 (MYREAL *) myrealloc(options->mg,
2020 sizeof(MYREAL) * (numpop * numpop + 1));
2021 for (i = 0; i < numpop; i++) {
2022 for (j = 0; j < numpop; j++) {
2023 if (j == i) {
2024 printf("Population %3li > --\n",
2025 i + 1);
2026 continue;
2027 }
2028 printf("From population %-3li to %-3li> ", j + 1,
2029 i + 1);fflush(stdout);
2030 #ifdef USE_MYREAL_FLOAT
2031 check = scanf("%f", &options->mg[tt++]);
2032 #else
2033 check = scanf("%lf", &options->mg[tt++]);
2034 #endif
2035 retval = scanf("%*[^\n]");
2036 if (check == 0)
2037 break;
2038 }
2039 }
2040 fflush(stdout); FGETS(input, LINESIZE, stdin);
2041 options->nummg = tt;
2042 tt = 0;
2043 printf
2044 ("\n You typed in the following migration matrix\n\n ");
2045 for (i = 0; i < numpop; i++) {
2046 for (j = 0; j < numpop; j++) {
2047 if (i != j)
2048 printf("%5.2f ", options->mg[tt++]);
2049 else {
2050 printf("----- ");
2051 }
2052 }
2053 printf("\n ");
2054 }
2055 printf("\n [Press <Return> to continue]\n");
2056 printf("\n\n");
2057 getchar();
2058 break;
2059 default:
2060 options->migrguess = FST;
2061 break;
2062 }
2063 break;
2064 case '3':
2065 options->usem = !options->usem;
2066 if (options->usem) {
2067 options->plotvar = PLOTM;
2068 options->migvar = PLOTM;
2069 options->profileparamtype = PLOTM;
2070 } else {
2071 options->plotvar = PLOT4NM;
2072 options->migvar = PLOT4NM;
2073 options->profileparamtype = PLOT4NM;
2074 }
2075
2076 break;
2077 case '4':
2078 printf("Mutation rate among loci\n");
2079 printf
2080 ("(C)onstant All loci have the same mutation rate [default]\n");
2081 printf
2082 ("(E)stimate Mutation rate \n");
2083 if (options->gamma)
2084 {
2085 printf(" GAMMA with Shape parameter alpha=%f [only ML]\n",
2086 options->alphavalue);
2087 }
2088 printf("(V)arying Mutation rates are different among loci [user input]\n");
2089 printf("(R)elative Mutation rates estimated from data\n===> ");
2090 fflush(stdout); FGETS(input, LINESIZE, stdin);
2091 switch (uppercase(input[0])) {
2092 case 'E':
2093 if(!options->bayes_infer)
2094 {
2095 options->gamma = TRUE;
2096 printf("Enter the value of the shape parameter alpha to use\n===> ");
2097 fflush(stdout); FGETS(input, LINESIZE, stdin);
2098 count = sscanf(input,"%lf",&tmpdouble);
2099 options->alphavalue = (MYREAL) tmpdouble;
2100 options->bayesmurates = FALSE;
2101 }
2102 else
2103 {
2104 options->bayesmurates = TRUE;
2105 }
2106 options->murates = FALSE;
2107 options->murates_fromdata=FALSE;
2108
2109 break;
2110 case 'V':
2111 printf
2112 ("Enter the number of loci and a rate for each of them\n[the rates are normalized to average to 1.0]\n===> ");
2113 fflush(stdout); FGETS(input, LINESIZE, stdin);
2114 temp = strtok(input, " ;,;\n");
2115 if (temp == NULL) {
2116 options->murates = FALSE;
2117 options->gamma = FALSE;
2118 options->murates_fromdata=FALSE;
2119 options->bayesmurates = FALSE;
2120 } else {
2121 options->gamma = FALSE;
2122 options->murates = TRUE;
2123 options->murates_fromdata=FALSE;
2124 options->muloci = atol(temp);
2125 options->bayesmurates = FALSE;
2126 options->mu_rates =
2127 (MYREAL *) myrealloc(options->mu_rates,
2128 sizeof(MYREAL) * options->muloci);
2129 // printf("%i> menu.c: opt realloc murate size %li\n",myID, options->muloci * sizeof (MYREAL));
2130 sum = 0.;
2131 count = 0;
2132 while (temp != NULL) {
2133 temp = strtok(NULL, " ;,;\n");
2134 if (temp == NULL)
2135 break;
2136 options->mu_rates[count] = atof(temp);
2137 sum += options->mu_rates[count];
2138 count++;
2139 }
2140 for (i = count; i < options->muloci; i++) {
2141 options->mu_rates[i] = options->mu_rates[count - 1];
2142 sum += options->mu_rates[count];
2143 }
2144 sum /= options->muloci;
2145 for (i = 0; i < options->muloci; i++) {
2146 options->mu_rates[i] /= sum;
2147 }
2148 }
2149 break;
2150 case 'R':
2151 options->gamma = FALSE;
2152 options->murates = TRUE;
2153 options->murates_fromdata=TRUE;
2154 options->bayesmurates = FALSE;
2155 break;
2156 default:
2157 options->gamma = FALSE;
2158 options->murates = FALSE;
2159 options->murates_fromdata=FALSE;
2160 options->bayesmurates = FALSE;
2161 break;
2162 }
2163 break;
2164 case '5':
2165 printf("Which FST calculation method?\n");
2166 printf("(T)heta can be different for each population\n");
2167 printf(" and migration rates are symmetric.\n");
2168 printf(" (Number of populations >= 2)\n");
2169 printf("(M)igration rate can be asymmetric\n");
2170 printf(" and Theta is the same for both populations\n");
2171 printf(" (Number of populations = 2)\n===> ");
2172 fflush(stdout); FGETS(input, LINESIZE, stdin);
2173 switch (uppercase(input[0])) {
2174 case 'M':
2175 options->fsttype = 'M';
2176 fst_type('M');
2177 break;
2178 case 'T':
2179 options->fsttype = 'T';
2180 fst_type('T');
2181 break;
2182 default:
2183 options->fsttype = 'T';
2184 fst_type('T');
2185 break;
2186 }
2187 break;
2188 case '6':
2189 options->printfst = !options->printfst;
2190 break;
2191 case '7':
2192 set_menu_localities(options);
2193 break;
2194 case '8': /* fill in the custom migration
2195 * matrix */
2196 read_custom_menu_migration(options);
2197 break;
2198 case '9':
2199 options->geo = !options->geo;
2200 if(options->geo)
2201 {
2202 menu_get_filename(" What is the filename for the distance matrix file?",
2203 GEOFILE, options->geofilename);
2204 }
2205 break;
2206 #ifdef POPMODEL
2207 case '10':
2208 set_popmodel(options);
2209 break;
2210 #endif
2211 }
2212 }
2213 while (uppercase(input[0]) != 'Y');
2214 myfree(outputstring);
2215 }
2216
set_menu_localities(option_fmt * options)2217 void set_menu_localities(option_fmt *options)
2218 {
2219 char *input;
2220 char *tmp;
2221 int i;
2222 input = (char *) mycalloc(LINESIZE,sizeof(char));
2223 tmp = (char *) mycalloc(LINESIZE,sizeof(char));
2224 printf("Associate sampling locations with populations\n");
2225 printf("---------------------------------------------\n");
2226 printf("This menu allows to combine sample locations into populations\n");
2227 printf("For example there are 4 locations: 1, 2, 3, 4\n");
2228 printf("They can be combined into 2 populations\n");
2229 printf("by mapping the 4 positions 1, 2, 3, 4 to 1, 1, 1, 2\n");
2230 printf("Migrate will now combine the first three locations\n\n");
2231 printf("Give (1) the number of populations <return> then (2) the mapping\n\n");
2232 printf("How many localities are in the data set?\n");
2233 printf("[Default: every sampling location is a population]\n");
2234 printf("> ");fflush(stdout);
2235 FGETS(input,LINESIZE,stdin);
2236 if(input[0]!='\0')
2237 {
2238 options->newpops_numalloc = atoi(input);
2239 options->newpops = (long*) myrealloc(options->newpops, options->newpops_numalloc * sizeof(long));
2240 printf("Enter now the remappings (little checking is done with this, enter exactly the number of locations)\n");
2241 for(i=1;i<=options->newpops_numalloc;i++)
2242 {
2243 switch((int) log10((double) options->newpops_numalloc))
2244 {
2245 case 0:
2246 case 1:
2247 printf("%2i",i); break;
2248 case 2:
2249 printf("%3i",i); break;
2250 default:
2251 printf("%4i",i); break;
2252 }
2253 }
2254 printf("\n>");fflush(stdout);
2255 FGETS(input,LINESIZE,stdin);
2256 set_localities(&input,&tmp,options);
2257 }
2258 else
2259 {
2260 printf("All locations are treated as populations.");
2261 options->newpops = (long*) myrealloc(options->newpops, sizeof(long));
2262 options->newpops[0] = 1;
2263 options->newpops_numalloc = 1;
2264 }
2265 myfree(input);
2266 myfree(tmp);
2267 }
2268
2269 ///print the standard question "are the settings correct? ....." at the end of the menu
print_bottom_menu_part()2270 void print_bottom_menu_part()
2271 {
2272 printf("\n\n Are the settings correct?\n");
2273 printf(" (Type Y to go back to the main menu or the number for a menu to change)\n===>");
2274 }
2275
2276 ///
2277 ///displays and makes changes to the options that manage the analysis strategy and run condition
2278 ///
menuStrategy(option_fmt * options)2279 void menuStrategy(option_fmt * options) {
2280 boolean done = FALSE;
2281 do {
2282 printf(" SEARCH STRATEGY\n\n");
2283 if (!options->bayes_infer) {
2284 display_ml_mcmc(options);
2285 print_bottom_menu_part();
2286 done = menuStrategy_ml(options);
2287 } else {
2288 display_bayes_mcmc(options);
2289 print_bottom_menu_part();
2290 done = menuStrategy_bayes(options);
2291 }
2292 } while (done == FALSE);
2293 }
2294
2295 //-------------------------------------------------------------------------------------
2296 // prior menu functions
2297
2298 /// \brief Menu for multiplier proposal
2299 /// Menu for multiplier proposal
set_mult_prior(int paramgroupm,prior_fmt * prior)2300 void set_mult_prior(int paramgroupm, prior_fmt *prior)
2301 {
2302 char input[LINESIZE];
2303 long count = 0;
2304 do {
2305 printf("Specify the MINIMUM, MAXIMUM and maximal MULTIPLIER\n[Current setting: {%f %f %f}\n===> ",
2306 prior->min, prior->max, prior->delta);
2307 fflush(stdout); FGETS(input, LINESIZE, stdin);
2308 if (input[0] == '\0')
2309 break;
2310 #ifdef USE_MYREAL_FLOAT
2311 count = sscanf(input, "%f%f%f", &prior->min, &prior->max, &prior->delta);
2312 #else
2313 count = sscanf(input, "%lf%lf%lf", &prior->min, &prior->max, &prior->delta);
2314 #endif
2315 } while (count != 3 && count > 0);
2316 }
2317
2318 /// \brief Menu for exponential proposal
2319 /// Menu for exponential proposal
set_exp_prior(int paramgroupm,prior_fmt * prior)2320 void set_exp_prior(int paramgroupm, prior_fmt *prior)
2321 {
2322 char input[LINESIZE];
2323 long count = 0;
2324 do {
2325 printf("Specify the MINIMUM, MEAN, MAXIMUM\n[Current setting: {%f %f %f}\n===> ",
2326 prior->min, prior->mean, prior->max);
2327 fflush(stdout); FGETS(input, LINESIZE, stdin);
2328 if (input[0] == '\0')
2329 break;
2330 #ifdef USE_MYREAL_FLOAT
2331 count = sscanf(input, "%f%f%f", &prior->min, &prior->mean, &prior->max);
2332 #else
2333 count = sscanf(input, "%lf%lf%lf", &prior->min, &prior->mean, &prior->max);
2334 #endif
2335 } while (count != 3 && count > 0);
2336 }
2337
2338 /// \brief Menu for gamma proposal
2339 /// Menu for gamma proposal, mean = alpha * beta = 1 ->
set_gamma_prior(int paramgroupm,prior_fmt * prior)2340 void set_gamma_prior(int paramgroupm, prior_fmt *prior)
2341 {
2342 char input[LINESIZE];
2343 long count = 0;
2344 do {
2345 printf("Specify the MINIMUM, MEAN, MAXIMUM, ALPHA\n[Current setting: {%f %f %f %f}\n===> ",
2346 prior->min, prior->mean, prior->max, prior->alpha);
2347 fflush(stdout); FGETS(input, LINESIZE, stdin);
2348 if(input[0]=='\0')
2349 return;
2350 #ifdef USE_MYREAL_FLOAT
2351 count = sscanf(input, "%f%f%f%f", &prior->min, &prior->mean, &prior->max, &prior->alpha);
2352 #else
2353 count = sscanf(input, "%lf%lf%lf%lf", &prior->min, &prior->mean, &prior->max, &prior->alpha);
2354 #endif
2355 } while (count != 4 && count > 0);
2356 }
2357
2358 /// \brief Menu for exponential with window proposal
2359 /// Menu for exponential with window proposal
set_wexp_prior(int paramgroupm,prior_fmt * prior)2360 void set_wexp_prior(int paramgroupm, prior_fmt *prior)
2361 {
2362 char input[LINESIZE];
2363 long count = 0;
2364 do {
2365 printf("Specify the MINIMUM, MEAN, MAXIMUM and DELTA\nUse for DELTA about 1/10 of the MIN-MAX range\n[Current setting: {%f %f %f %f}\n===> ",
2366 prior->min, prior->mean, prior->max, prior->delta);
2367 fflush(stdout); FGETS(input, LINESIZE, stdin);
2368 if(input[0]=='\0')
2369 return;
2370 #ifdef USE_MYREAL_FLOAT
2371 count = sscanf(input, "%f%f%f%f", &prior->min, &prior->mean, &prior->max, &prior->delta);
2372 #else
2373 count = sscanf(input, "%lf%lf%lf%lf", &prior->min, &prior->mean, &prior->max, &prior->delta);
2374 #endif
2375 } while (count != 4 && count > 0);
2376 }
2377
2378 /// \brief Menu for multiplier proposal
2379 /// Menu for multiplier proposal
set_uni_prior(int paramgroupm,prior_fmt * prior)2380 void set_uni_prior(int paramgroupm, prior_fmt *prior)
2381 {
2382 char input[LINESIZE];
2383 long count = 0;
2384 do {
2385 printf("Specify the MINIMUM, MAXIMUM, DELTA\nUse for DELTA about 1/10 of the MIN-MAX range\n[Current setting: {%f %f %f}\n===> ",
2386 prior->min, prior->max, prior->delta);
2387 fflush(stdout); FGETS(input, LINESIZE, stdin);
2388 if(input[0]=='\0')
2389 return;
2390 #ifdef USE_MYREAL_FLOAT
2391 count = sscanf(input, "%f%f%f", &prior->min, &prior->max, &prior->delta);
2392 #else
2393 count = sscanf(input, "%lf%lf%lf", &prior->min, &prior->max, &prior->delta);
2394 #endif
2395 } while (count != 3 && count > 0);
2396 prior->mean = (prior->max + prior->min) / 2.;
2397 //prior->delta = (prior->max - prior->min) / 10.;
2398 }
2399
2400 /// \brief return "set" "-" depending on prior set
2401 /// returns a short string that shows whether this prior distribution is set
is_prior(int priortype,int priorset)2402 char * is_prior(int priortype, int priorset)
2403 {
2404 return priortype == priorset ? "Set" : "-";
2405 }
2406
2407 /// \brief returns priortype sting
2408 /// returns a string that shows what prior distribution is set
is_proposaltype(boolean proposalset)2409 char * is_proposaltype(boolean proposalset)
2410 {
2411 switch(proposalset)
2412 {
2413 case TRUE: return "Slice sampling" ;
2414 default: return "Metropolis sampling";
2415 }
2416 }
2417
2418
2419 /// \brief Menu to set all prior distributions
2420 /// menu to set proposal distribution, returns TRUE when standard course of action
2421 /// FALSE when a problem occured.
set_proposal_menu(int paramgroup,option_fmt * options)2422 boolean set_proposal_menu(int paramgroup, option_fmt *options)
2423 {
2424 char input[LINESIZE];
2425 boolean proposal;
2426
2427 switch(paramgroup)
2428 {
2429 case THETAPRIOR:
2430 printf("Proposal setting for all population sizes [Theta]:\n");
2431 proposal = options->slice_sampling[THETAPRIOR];
2432 break;
2433 case MIGPRIOR:
2434 printf("Proposal setting for all migration rates [%s]:\n",options->usem ? "M" : "Theta*M");
2435 proposal = options->slice_sampling[MIGPRIOR];
2436 break;
2437 case RATEPRIOR:
2438 printf("Proposal setting for mutation rate modifier [r]:\n");
2439 proposal = options->slice_sampling[RATEPRIOR];
2440 break;
2441 default:
2442 return FALSE;
2443 }
2444 input[0] = '\0';
2445 printf("Which proposal distribution? [current: %11s] (choices: Slice or Metropolis)\n===> ",
2446 is_proposaltype(proposal));
2447 fflush(stdout); FGETS(input, LINESIZE, stdin);
2448 switch (uppercase(input[0]))
2449 {
2450 case 'M':
2451 options->slice_sampling[paramgroup] = FALSE;
2452 break;
2453 case 'S':
2454 default:
2455 options->slice_sampling[paramgroup] = TRUE;
2456 break;
2457 }
2458 input[0]='\0';
2459 return TRUE;
2460 }
2461
2462
2463 /// \brief returns priortype sting
2464 /// returns a string that shows what prior distribution is set
is_priortype(int priorset)2465 char * is_priortype(int priorset)
2466 {
2467 switch(priorset)
2468 {
2469 // case SLICE: return "Uniform/Slice" ;
2470 case MULTPRIOR: return "Multiplier" ;
2471 case EXPPRIOR: return "Exponential" ;
2472 case WEXPPRIOR: return "Exp window";
2473 case GAMMAPRIOR: return "Gamma";
2474 case UNIFORMPRIOR: return "Uniform";
2475 default: return "-";
2476 }
2477 }
2478
2479
2480 /// \brief Menu to set all prior distributions
2481 /// menu to set all prior distributions, returns TRUE when standard course of action
2482 /// FALSE when a problem occured.
set_prior_menu(int paramgroup,option_fmt * options)2483 boolean set_prior_menu(int paramgroup, option_fmt *options)
2484 {
2485 int tmp;
2486 char input[LINESIZE];
2487 prior_fmt *prior;
2488
2489 switch(paramgroup)
2490 {
2491 case THETAPRIOR:
2492 printf("Prior setting for all population sizes [Theta]:\n");
2493 prior = options->bayespriortheta;
2494 break;
2495 case MIGPRIOR:
2496 printf("Prior setting for all migration rates [%s]:\n",options->usem ? "M" : "Theta*M");
2497 prior = options->bayespriorm;
2498 break;
2499 case RATEPRIOR:
2500 printf("Prior setting for mutation rate modifier [r]:\n");
2501 prior = options->bayespriorrate;
2502 break;
2503 default:
2504 return FALSE;
2505 }
2506 do{
2507 input[0] = '\0';
2508 // printf("\nSet the prior distribution and its parameters\n\n");
2509 //printf(" %i Set Slice sampler? %11s\n", SLICE,
2510 // is_prior(SLICE,options->bayesprior[paramgroup]));
2511 printf(" %i Set Multiplication prior distribution? %11s\n", MULTPRIOR,
2512 is_prior(MULTPRIOR,options->bayesprior[paramgroup]));
2513 printf(" %i Set Exponential prior distribution? %11s\n", EXPPRIOR,
2514 is_prior(EXPPRIOR,options->bayesprior[paramgroup]));
2515 printf(" %i Set Exponential prior with window distribution?%11s\n", WEXPPRIOR,
2516 is_prior(WEXPPRIOR,options->bayesprior[paramgroup]));
2517 printf(" %i Set Gamma prior distribution? %11s\n", GAMMAPRIOR,
2518 is_prior(GAMMAPRIOR,options->bayesprior[paramgroup]));
2519 printf(" %i Set Uniform prior distribution? %11s\n\n", UNIFORMPRIOR,
2520 is_prior(UNIFORMPRIOR,options->bayesprior[paramgroup]));
2521 printf(" (Type Y to go back to the main menu or the letter for an entry to change)\n===> ");
2522 fflush(stdout); FGETS(input, LINESIZE, stdin);
2523 if(uppercase(input[0]) == 'Y')
2524 return TRUE;
2525 input[1] = '\0';
2526 tmp = atoi(input);
2527 switch (tmp)
2528 {
2529 case UNIFORMPRIOR:
2530 set_uni_prior(paramgroup, prior);
2531 options->bayesprior[paramgroup] = UNIFORMPRIOR;
2532 break;
2533 case EXPPRIOR:
2534 set_exp_prior(paramgroup, prior);
2535 options->bayesprior[paramgroup] = EXPPRIOR;
2536 break;
2537 case WEXPPRIOR:
2538 set_wexp_prior(paramgroup, prior);
2539 options->bayesprior[paramgroup] = WEXPPRIOR;
2540 break;
2541 case MULTPRIOR:
2542 set_mult_prior(paramgroup, prior);
2543 options->bayesprior[paramgroup] = MULTPRIOR;
2544 break;
2545 case GAMMAPRIOR:
2546 set_gamma_prior(paramgroup, prior);
2547 options->bayesprior[paramgroup] = GAMMAPRIOR;
2548 break;
2549 default:
2550 set_uni_prior(paramgroup, prior);
2551 options->bayesprior[paramgroup] = UNIFORMPRIOR;
2552 break;
2553 }
2554 } while (uppercase(input[0]) != 'Y');
2555 input[0]='\0';
2556 return TRUE;
2557 }
2558
set_autotuning(option_fmt * options)2559 void set_autotuning(option_fmt *options)
2560 {
2561
2562 char input[LINESIZE];
2563 printf("Autotuning for Metropolis-Hastings sampling (YES or NO)\n");
2564 printf("===> ");
2565 fflush(stdout);
2566 FGETS(input, LINESIZE, stdin);
2567 if (input[0] != '\0')
2568 {
2569 if(input[0]=='Y' || input[0]=='y')
2570 options->has_autotune = TRUE;
2571 else
2572 options->has_autotune=FALSE;
2573 }
2574 if (options->has_autotune)
2575 {
2576 printf("Give a target acceptance ratio [suggested 0.44]\n");
2577 printf("===> ");
2578 fflush(stdout);
2579 FGETS(input, LINESIZE, stdin);
2580 if (input[0] != '\0')
2581 options->autotune = atof(input);
2582 else
2583 {
2584 if (options->autotune>1.0 || options->autotune<0.01)
2585 options->autotune=AUTOTUNEDEFAULT;
2586 }
2587 }
2588 }
2589
2590 ///
2591 /// sets proposal distribution
2592 void
menuProposal(option_fmt * options)2593 menuProposal(option_fmt * options) {
2594 char input[LINESIZE];
2595 do
2596 {
2597 input[0] = '\0';
2598 printf (" Proposal distribution setting [You still need to set the PRIOR distribution!]:\n");
2599 printf(" 1 Set proposal distribution for Theta? %11s\n",
2600 is_proposaltype(options->slice_sampling[THETAPRIOR]));
2601 printf(" 2 Set proposal distribution for migration? %11s\n",
2602 is_proposaltype(options->slice_sampling[MIGPRIOR]));
2603 if(options->bayesmurates)
2604 {
2605 printf(" 3 Set mutation rate modifier prior distribution?%11s\n",
2606 is_proposaltype(options->slice_sampling[RATEPRIOR]));
2607 if(options->has_autotune)
2608 printf(" 4 Autotuning of the acceptance ratio %7s %.2f\n","YES",options->autotune);
2609 else
2610 printf(" 4 Autotuning of the acceptance ratio %11s\n","NO");
2611 }
2612 else
2613 {
2614 if(options->has_autotune)
2615 printf(" 3 Autotuning of the acceptance ratio %7s %.2f\n","YES",options->autotune);
2616 else
2617 printf(" 3 Autotuning of the acceptance ratio %11s\n","NO");
2618 }
2619 printf("\n (Type Y to go back to the main menu or the letter for a menu to change)\n===> ");
2620 fflush(stdout); FGETS(input, LINESIZE, stdin);
2621 switch (input[0])
2622 {
2623 case '1':
2624 set_proposal_menu(THETAPRIOR, options);
2625 break;
2626 case '2':
2627 set_proposal_menu(MIGPRIOR, options);
2628 break;
2629 case '3':
2630 if(options->bayesmurates)
2631 {
2632 set_proposal_menu(RATEPRIOR, options);
2633 }
2634 else
2635 {
2636 set_autotuning(options);
2637 }
2638 break;
2639 case '4':
2640 set_autotuning(options);
2641 break;
2642 default:
2643 break;
2644 }
2645 } while (uppercase(input[0]) != 'Y');;
2646 }
2647 ///
2648 /// sets prior distribution parameters
2649 void
menuPrior(option_fmt * options)2650 menuPrior(option_fmt * options) {
2651 char input[LINESIZE];
2652 do
2653 {
2654 input[0] = '\0';
2655 printf (" Prior distribution setting:\n");
2656 printf(" 1 Set Theta prior distribution? %11s\n",
2657 is_priortype(options->bayesprior[THETAPRIOR]));
2658 printf(" 2 Set Migration prior distribution? %11s\n",
2659 is_priortype(options->bayesprior[MIGPRIOR]));
2660 if(options->bayesmurates)
2661 {
2662 printf(" 3 Set mutation rate modifier prior distribution?%11s\n",
2663 is_priortype(options->bayesprior[RATEPRIOR]));
2664 }
2665 printf("\n (Type Y to go back to the main menu or the letter for a menu to change)\n===> ");
2666 fflush(stdout); FGETS(input, LINESIZE, stdin);
2667 switch (input[0])
2668 {
2669 case '1':
2670 set_prior_menu(THETAPRIOR, options);
2671 break;
2672 case '2':
2673 set_prior_menu(MIGPRIOR, options);
2674 break;
2675 case '3':
2676 if(options->bayesmurates)
2677 {
2678 set_prior_menu(RATEPRIOR, options);
2679 }
2680 break;
2681 default:
2682 break;
2683 }
2684 } while (uppercase(input[0]) != 'Y');;
2685 }
2686
2687 ///
2688 ///reads bayes menu choices and sets options accordingly
menuStrategy_bayes(option_fmt * options)2689 boolean menuStrategy_bayes(option_fmt * options) {
2690 char input[LINESIZE];
2691 char *priorstring;
2692 long count = 0;
2693 #ifdef ZNZ
2694 char *extension=NULL;
2695 #endif
2696 //read user input
2697 fflush(stdout); FGETS(input, LINESIZE, stdin);
2698 //if user types yes or YES or YES exit and tell menuStrategy that we are done
2699 if (strchr("Yy", input[0]))
2700 return TRUE;
2701
2702 priorstring = (char *) mycalloc(LINESIZE, sizeof(char));
2703 //make changes to options
2704 switch (atoi(input)) {
2705 case BAYESSTRATEGY:
2706 //strategy change
2707
2708 if (strchr(input,'0'))
2709 {
2710 options->bayes_infer = !options->bayes_infer;
2711 if(options->bayes_infer)
2712 options->lchains=1;
2713 }
2714 input[0] = '\0';
2715 break;
2716 case BAYESOUT:
2717 printf(" Save posterior distribution (frequency histogram values)? [YES]\n===> ");
2718 fflush(stdout); FGETS(input, LINESIZE, stdin);
2719 if (uppercase(input[0]) != 'N')
2720 {
2721 options->has_bayesfile = TRUE;
2722 menu_get_filename(" What is the filename for the posterior distribution (frequency histogram)?",
2723 BAYESFILE, options->bayesfilename);
2724 }
2725 else
2726 {
2727 options->has_bayesfile = FALSE;
2728 }
2729 break;
2730 case BAYESMDIMOUT:
2731 printf(" Save posterior distribution (all raw parameter values)? [NO]\n===> ");
2732 fflush(stdout); FGETS(input, LINESIZE, stdin);
2733 if (uppercase(input[0]) == 'Y')
2734 {
2735 options->has_bayesmdimfile = TRUE;
2736 #ifdef ZNZ
2737 menu_get_filename("What is the filename for the complete posterior distribution (raw parameter values)?\n[if the extension of the file is .gz then the will be compressed]",
2738 BAYESMDIMFILE, options->bayesmdimfilename);
2739 unpad(options->bayesmdimfilename, " ");
2740 extension = strrchr(options->bayesmdimfilename,'.');
2741 if(extension!=NULL && !strncmp(extension,".gz",3))
2742 {
2743 options->use_compressed = 1;
2744 }
2745 else
2746 {
2747 options->use_compressed = 0;
2748 }
2749 #else
2750 menu_get_filename("What is the filename for the complete posterior distribution (raw parameter values)?",
2751 BAYESMDIMFILE, options->bayesmdimfilename);
2752 options->use_compressed = 0;
2753 #endif
2754 do
2755 {
2756 printf("Sampling interval for the raw parameter values\n");
2757 printf("[Small values can result in a HUGE file!!]\n");
2758 printf("[Default is the same as the sampling increment]\n");
2759 printf("[Examples: 1 --> saving all parameters every sampling increment]\n");
2760 printf("[ 2 --> saving every second sampling increment]\n");
2761 printf("[ 100 --> saving only very hundredth sample]\n");
2762 printf("[ if long-inc=50 and here we have here 100 then only]\n");
2763 printf("[ every 5,000th step is saved to file]\n");
2764 printf("[ with many loci and large run this is still a large file]\n");
2765 printf("Current setting: %li\n===> ",options->bayesmdiminterval);
2766 fflush(stdout); FGETS(input, LINESIZE, stdin);
2767 if(input[0] == '\0' && options->bayesmdiminterval > 0)
2768 break;
2769 else
2770 options->bayesmdiminterval = atol(input);
2771 }
2772 while(options->bayesmdiminterval<1);
2773 }
2774 else
2775 {
2776 options->has_bayesmdimfile = FALSE;
2777 }
2778 break;
2779 case BAYESBINNING:
2780 do {
2781 printf("Specify the number of bins for THETA and %s!\n", options->usem ? "M" : "Theta*M");
2782 printf("Current setting: Bins Interval width\n");
2783 printf(" Theta %li %f\n", options->bayespriortheta->bins,
2784 (options->bayespriortheta->max - options->bayespriortheta->min) / (float) options->bayespriortheta->bins);
2785 printf("%7.7s %li %f\n", options->usem ? "M" : "Theta*M", options->bayespriorm->bins,
2786 (options->bayespriorm->max - options->bayespriorm->min) / (float) options->bayespriorm->bins);
2787 if(options->bayesmurates)
2788 printf(" Rate %li %f\n", options->bayespriorrate->bins,
2789 (options->bayespriorrate->max - options->bayespriorrate->min) / (float) options->bayespriorrate->bins);
2790 printf("\n==> ");
2791 fflush(stdout); FGETS(input, LINESIZE, stdin);
2792 if (input[0] == '\0')
2793 break;
2794 if(options->bayesmurates)
2795 count = sscanf(input, "%li%li%li", &options->bayespriortheta->bins, &options->bayespriorm->bins, &options->bayespriorrate->bins);
2796 else
2797 count = sscanf(input, "%li%li", &options->bayespriortheta->bins, &options->bayespriorm->bins);
2798 } while (count != (options->bayesmurates ? 3 : 2) && count > 0);
2799 break;
2800 case BAYESPRETTY:
2801 printf("Specify the method for plotting the posterior distribution!\n");
2802 printf("Valid options are:\n");
2803 printf("ALL use the range of the prior distribution\n");
2804 printf("P99 do not plot the values higher than the 99%% percentile for each parameter\n");
2805 printf("MAXP99 plot up to the maximum 99%% percentile of all parameters.\n");
2806 printf("TOTAL plot up to the 100%% percentile of all parameters.\n");
2807 printf("Current setting: %s\n",
2808 (options->bayespretty == PRETTY_P99 ? "P99" :
2809 (options->bayespretty == PRETTY_MAX ? "ALL" :
2810 (options->bayespretty == PRETTY_P100 ? "TOTAL" : "MAXP99"))));
2811 printf("\n\n==> ");
2812 fflush(stdout); FGETS(input, LINESIZE, stdin);
2813 if (input[0] == '\0')
2814 break;
2815 switch(uppercase(input[0]))
2816 {
2817 case 'A':
2818 options->bayespretty = PRETTY_MAX;
2819 break;
2820 case 'P':
2821 options->bayespretty = PRETTY_P99;
2822 break;
2823 case 'T':
2824 options->bayespretty = PRETTY_P100;
2825 break;
2826 case 'M':
2827 default:
2828 options->bayespretty = PRETTY_P99MAX;
2829 break;
2830 }
2831 break;
2832 case BAYESFREQ:
2833 do {
2834 printf("What is the frequency of the tree updates vs the parameter updates?\n");
2835 printf("[A frequncy of 1.0 updates only trees, frequency of 0.0 updates only parameters]\n");
2836 printf("[Suggestion for small problems: 0.5, actual value: %f]\n===> ", options->updateratio);
2837 fflush(stdout); FGETS(input, LINESIZE, stdin);
2838 if (input[0] != '\0')
2839 options->updateratio = atof(input);
2840 } while (options->updateratio > 1);
2841 break;
2842 case BAYESPROPOSAL:
2843 menuProposal(options);
2844 break;
2845 case BAYESPRIOR:
2846 menuPrior(options);
2847 break;
2848 case BAYESLCHAINS:
2849 do {
2850 printf(" How many Long Chains?\n===> ");
2851 fflush(stdout); FGETS(input, LINESIZE, stdin);
2852 options->lchains = atoi(input);
2853 if (options->lchains < 0)
2854 printf(" Must be non-negative\n");
2855 }
2856 while (options->lchains < 0);
2857 break;
2858 case BAYESSKIP:
2859 do {
2860 printf(" How many steps (tree changes, parameter changes) to skip?\n===> ");
2861 fflush(stdout); FGETS(input, LINESIZE, stdin);
2862 options->lincrement = atoi(input);
2863 if (options->lincrement <= 0)
2864 printf(" Must be positive\n");
2865 }
2866 while (options->lincrement <= 0);
2867 break;
2868 case BAYESSAMPLES:
2869 do {
2870 printf(" How many trees and parameter steps to record?\n===> ");
2871 fflush(stdout); FGETS(input, LINESIZE, stdin);
2872 options->lsteps = atoi(input);
2873 if (options->lsteps <= 0)
2874 printf(" Must be a positive integer\n");
2875 }
2876 while (options->lsteps <= 0);
2877 break;
2878 case BAYESBURNIN:
2879 do {
2880 printf(" How many steps to discard? (burn-in)\n)");
2881 //printf(" For standard fixed burn-in simply give a number, for example 50000\n");
2882 //printf(" With Metropolis-Hastings sampling (ignored with Slice-sampling)\n");
2883 //printf(" append a space and a character (a, e, or t) to invoke autostopping\n");
2884 //printf(" for example like this: 100000 a\n");
2885 //printf(" a using variance ratios every 1000 steps to stop: abs(1-var(n)/var(n-1))<0.01\n");
2886 //printf(" e using the effective samples size to stop: min(ess(param_i)>%li\n",(long) ESSMINIMUM);
2887 //printf(" t using the average acceptance ratio: %f<avg(acceptance(param_i)<%f\n",(MYREAL) options->autotune-0.05, (MYREAL) options->autotune+0.05);
2888 printf("===> ");
2889 fflush(stdout);
2890 FGETS(input, LINESIZE, stdin);
2891 // a autostop using variance ratios every 1000 steps
2892 // e autostop using effective sample size of run > ESSMINIMUM (of worst parameter)
2893 // t autostop using when target acceptance ratio is reached (of the average)
2894 if(strchr(input,'a'))
2895 {
2896 options->burnin_autostop = 'a';
2897 }
2898 else
2899 {
2900 if(strchr(input,'t'))
2901 {
2902 options->burnin_autostop = 't';
2903 if(!options->has_autotune)
2904 {
2905 options->has_autotune=TRUE;
2906 options->autotune=AUTOTUNEDEFAULT;
2907 }
2908 }
2909 else
2910 {
2911 if(strchr(input,'e'))
2912 {
2913 options->burnin_autostop = 'e';
2914 }
2915 else
2916 {
2917 options->burnin_autostop = ' ';
2918 }
2919 }
2920 }
2921 options->burn_in = atoi(input);
2922 if (options->burn_in <= 0)
2923 printf(" Must be a positive integer or zero (0)\n");
2924 }
2925 while (options->burn_in < 0);
2926 break;
2927 case BAYESREPLICATE:
2928 do {
2929 printf(" Summarize over multiple chains? [YES | NO] \n===> ");
2930 fflush(stdout); FGETS(input, LINESIZE, stdin);
2931 if (uppercase(input[0]) == 'Y')
2932 {
2933 options->replicate = TRUE;
2934 printf(" How many independent runs\n===> ");
2935 fflush(stdout); FGETS(input, LINESIZE, stdin);
2936 options->replicatenum = ATOL(input);
2937 if (options->replicatenum < 1)
2938 printf(" Enter a number >= 1\n");
2939 }
2940 else
2941 {
2942 options->replicate = FALSE;
2943 options->replicatenum = 0;
2944 }
2945 }
2946 while (options->replicatenum < 0);
2947 break;
2948 case BAYESHEAT:
2949 printf
2950 (" Heating scheme? < NO | YES | STATIC | BOUNDED_ADAPTIVE >\n===> ");
2951 fflush(stdout); FGETS(input, LINESIZE, stdin);
2952 switch (tolower(input[0])) {
2953 case 'b':
2954 options->heating = 1;
2955 options->adaptiveheat = BOUNDED;
2956 options->heating_interval = 1;
2957 menuHeat(options, input);
2958 break;
2959 case 'y':
2960 case 's':
2961 options->heating = 1;
2962 options->adaptiveheat = NOTADAPTIVE;
2963 options->heating_interval = 1;
2964 menuHeat(options, input);
2965 break;
2966 case '\0':
2967 case 'n':
2968 default:
2969 options->heating = 0;
2970 options->adaptiveheat = NOTADAPTIVE;
2971 options->heated_chains = 1;
2972 break;
2973 }
2974 break;
2975 case BAYESMOVINGSTEPS:
2976 options->movingsteps = !options->movingsteps;
2977 if (options->movingsteps) {
2978 do {
2979 printf
2980 (" How big should the fraction of new genealogies\n");
2981 printf
2982 (" of the originally proposed number of samples be?\n[Use this option only when acceptance ratio is small (<0.01)\nRuntime can increase tremendeously]\n===> ");
2983 fflush(stdout); FGETS(input, LINESIZE, stdin);
2984 options->acceptfreq = atof(input);
2985 if (options->acceptfreq < 0 || options->acceptfreq > 1)
2986 printf(" Range should be between 0 - 1, and not %f\n",
2987 options->acceptfreq);
2988 }
2989 while (options->acceptfreq < 0 || options->acceptfreq > 1);
2990 }
2991 break;
2992 case BAYESGELMAN:
2993 printf(" Convergence statistic options: Pairs , Sum, NO\n");
2994 printf(" [Currently set to %15s]\n===> ",
2995 options->gelman ? (options->gelmanpairs ? "YES:Pairs" : "YES:Summary" ) : " NO");
2996 fflush(stdout); FGETS(input, LINESIZE, stdin);
2997 if (input[0] == '\0')
2998 break;
2999 switch(uppercase(input[0]))
3000 {
3001 case 'N':
3002 options->gelman = FALSE;
3003 options->gelmanpairs = FALSE;
3004 break;
3005 case 'P':
3006 options->gelman = TRUE;
3007 options->gelmanpairs = TRUE;
3008 break;
3009 default:
3010 options->gelman = TRUE;
3011 options->gelmanpairs = FALSE;
3012 break;
3013 }
3014 break;
3015 case BAYESPRIORALONE:
3016 printf
3017 (" Show only the prior distributions? < NO | YES >\nwith YES no data is analyzed!!!\n===> ");
3018 fflush(stdout); FGETS(input, LINESIZE, stdin);
3019 switch (tolower(input[0])) {
3020 case 'y':
3021 options->prioralone = TRUE;
3022 break;
3023 case 'n':
3024 default:
3025 options->prioralone = FALSE;
3026 }
3027 break;
3028 default:
3029 break;
3030 }
3031 myfree(priorstring);
3032 return FALSE;
3033 }
3034
3035
menuStrategy_ml(option_fmt * options)3036 boolean menuStrategy_ml(option_fmt * options)
3037 {
3038 char input[LINESIZE];
3039 //read user input
3040 fflush(stdout); FGETS(input, LINESIZE, stdin);
3041 //if user types yes or Yes or YES exit and tell menuStrategy that we are done
3042 if (strchr("Yy", input[0]))
3043 return TRUE;
3044
3045 //make changes to options
3046 switch (atoi(input)) {
3047 case MLSTRATEGY://strategy change
3048 if (strchr(input,'0'))
3049 options->bayes_infer = !options->bayes_infer;
3050 input[0] = '\0';
3051 break;
3052 case MLSHORTCHAINS:
3053 do {
3054 printf(" How many Short Chains?\n===> ");
3055 fflush(stdout); FGETS(input, LINESIZE, stdin);
3056 options->schains = atoi(input);
3057 if (options->schains < 0)
3058 printf(" Must be non-negative\n");
3059 }
3060 while (options->schains < 0);
3061 break;
3062 case MLSHORTSKIP:
3063 do {
3064 printf(" How many trees to skip?\n===> ");
3065 fflush(stdout); FGETS(input, LINESIZE, stdin);
3066 options->sincrement = atoi(input);
3067 if (options->sincrement <= 0)
3068 printf(" Must be positive\n");
3069 }
3070 while (options->sincrement <= 0);
3071 break;
3072 case MLSHORTSAMPLES:
3073 do {
3074 printf(" How many trees to record?\n===> ");
3075 fflush(stdout); FGETS(input, LINESIZE, stdin);
3076 options->ssteps = atoi(input);
3077 if (options->ssteps <= 0)
3078 printf(" Must be a positive integer\n");
3079 }
3080 while (options->ssteps <= 0);
3081 break;
3082 case MLLONGCHAINS:
3083 do {
3084 printf(" How many Long Chains?\n===> ");
3085 fflush(stdout); FGETS(input, LINESIZE, stdin);
3086 options->lchains = atoi(input);
3087 if (options->lchains < 0)
3088 printf(" Must be non-negative\n");
3089 }
3090 while (options->lchains < 0);
3091 break;
3092 case MLLONGSKIP:
3093 do {
3094 printf(" How many trees to skip?\n===> ");
3095 fflush(stdout); FGETS(input, LINESIZE, stdin);
3096 options->lincrement = atoi(input);
3097 if (options->lincrement <= 0)
3098 printf(" Must be positive\n");
3099 }
3100 while (options->lincrement <= 0);
3101 break;
3102 case MLLONGSAMPLES:
3103 do {
3104 printf(" How many trees to record?\n===> ");
3105 fflush(stdout); FGETS(input, LINESIZE, stdin);
3106 options->lsteps = atoi(input);
3107 if (options->lsteps <= 0)
3108 printf(" Must be a positive integer\n");
3109 }
3110 while (options->lsteps <= 0);
3111 break;
3112 case MLBURNIN:
3113 do {
3114 printf(" How many genealogies to discard?\n===> ");
3115 fflush(stdout); FGETS(input, LINESIZE, stdin);
3116 options->burn_in = atoi(input);
3117 if (options->burn_in <= 0)
3118 printf(" Must be a positive integer or zero (0)\n");
3119 }
3120 while (options->burn_in < 0);
3121 break;
3122 case MLREPLICATE:
3123 options->replicate = !options->replicate;
3124 if (options->replicate) {
3125 do {
3126 printf(" Summarize over (L)ong chains?\n");
3127 printf(" or (M)ultiple runs ?\n");
3128 printf(" [Default is (L)]\n===> ");
3129 fflush(stdout); FGETS(input, LINESIZE, stdin);
3130 if (uppercase(input[0]) == 'M') {
3131 printf(" How many independent runs\n===> ");
3132 fflush(stdout); FGETS(input, LINESIZE, stdin);
3133 options->replicatenum = ATOL(input);
3134 if (options->lcepsilon < 0)
3135 printf(" Enter a number >= 1\n");
3136 } else
3137 options->replicatenum = 0;
3138 }
3139 while (options->replicatenum < 0);
3140 } else {
3141 options->replicatenum = 0;
3142 }
3143 break;
3144 case MLHEAT:
3145 printf
3146 (" Heating scheme? < NO | YES | ADAPTIVE | BOUNDED_ADAPTIVE>\n===> ");
3147 fflush(stdout); FGETS(input, LINESIZE, stdin);
3148 switch (tolower(input[0])) {
3149 case 'a':
3150 options->heating = 1;
3151 options->adaptiveheat = STANDARD;
3152 options->heating_interval = 1;
3153 menuHeat(options, input);
3154 break;
3155 case 'b':
3156 options->heating = 1;
3157 options->adaptiveheat = BOUNDED;
3158 options->heating_interval = 1;
3159 menuHeat(options, input);
3160 break;
3161 case 'y':
3162 options->heating = 1;
3163 options->adaptiveheat = NOTADAPTIVE;
3164 options->heating_interval = 1;
3165 menuHeat(options, input);
3166 break;
3167 case '\0':
3168 case 'n':
3169 default:
3170 options->heating = 0;
3171 options->adaptiveheat = NOTADAPTIVE;
3172 break;
3173 }
3174 break;
3175 case MLMOVINGSTEPS:
3176 options->movingsteps = !options->movingsteps;
3177 if (options->movingsteps) {
3178 do {
3179 printf
3180 (" How big should the fraction of new genealogies\n");
3181 printf
3182 (" of the originally proposed number of samples be?\n===> ");
3183 fflush(stdout); FGETS(input, LINESIZE, stdin);
3184 options->acceptfreq = atof(input);
3185 if (options->acceptfreq < 0 || options->acceptfreq > 1)
3186 printf(" Range should be between 0 - 1, and not %f\n",
3187 options->acceptfreq);
3188 }
3189 while (options->acceptfreq < 0 || options->acceptfreq > 1);
3190 }
3191 break;
3192 case MLEPSILON:
3193 do {
3194 printf(" Parameter likelihood epsilon?\n[INF is Default]\n===> ");
3195 fflush(stdout); FGETS(input, LINESIZE, stdin);
3196 if (uppercase(input[0]) == 'I')
3197 options->lcepsilon = LONGCHAINEPSILON;
3198 else
3199 options->lcepsilon = atof(input);
3200 if (options->lcepsilon <= 0)
3201 printf
3202 (" Must be a positive value, be warned: too small values will run the program forever\n");
3203 }
3204 while (options->lcepsilon <= 0);
3205 break;
3206 case MLGELMAN:
3207 printf(" Convergence statistic options: Pairs , Sum, NO\n");
3208 printf(" [Currently set to %10s]\n===> ",
3209 options->gelman ? (options->gelmanpairs ? "YES:Pairs" : "YES:Summary" ) : " NO");
3210 fflush(stdout); FGETS(input, LINESIZE, stdin);
3211 if (input[0] == '\0')
3212 break;
3213 switch(uppercase(input[0]))
3214 {
3215 case 'N':
3216 options->gelman = FALSE;
3217 options->gelmanpairs = FALSE;
3218 break;
3219 case 'P':
3220 options->gelman = TRUE;
3221 options->gelmanpairs = TRUE;
3222 break;
3223 default:
3224 options->gelman = TRUE;
3225 options->gelmanpairs = FALSE;
3226 break;
3227 }
3228 break;
3229 }
3230 return FALSE;
3231 }
3232
menuHeat(option_fmt * options,char * input)3233 void menuHeat(option_fmt * options, char *input)
3234 {
3235 if (options->heating > 0) {
3236 printf("Enter the number of different \"heated\" chains.\nMinimum is 4\n===> ");
3237 fflush(stdout); FGETS(input, LINESIZE, stdin);
3238 options->heated_chains = atol(input);
3239 if (options->heated_chains < 4)
3240 options->heated_chains = HEATED_CHAIN_NUM;
3241 printf("Enter the interval between swapping trees\nEnter 0 (zero) for NO swapping\n[Current interval is %li]\n===> ",
3242 options->heating_interval);
3243 fflush(stdout); FGETS(input, LINESIZE, stdin);
3244 // if (input[0] == '\0')
3245 // return;
3246 if (atol(input) > 0) {
3247 options->heating_interval = atol(input);
3248 options->heating = 1;
3249 } else {
3250 options->heating = 1;
3251 options->heating_interval = 0;
3252 options->heatedswap_off=TRUE;
3253 }
3254 read_heatvalues(options);
3255 }
3256 }
3257
3258 void
display_ml_mcmc(option_fmt * options)3259 display_ml_mcmc(option_fmt * options) {
3260 char temp[LINESIZE];
3261
3262 printf(" 0 Strategy: %42.42s\n", options->bayes_infer ?
3263 "Bayesian Inference" :
3264 "Maximum Likelihood");
3265 printf(" 1 Number of short chains to run? %6ld\n",
3266 options->schains);
3267 if (options->schains > 0) {
3268 printf
3269 (" 2 Short sampling increment? %6ld\n",
3270 options->sincrement);
3271 printf
3272 (" 3 Number of recorded genealogies in short chain?%6ld\n",
3273 options->ssteps);
3274 }
3275 printf(" 4 Number of long chains to run? %6ld\n",
3276 options->lchains);
3277 if (options->lchains > 0) {
3278 printf
3279 (" 5 Long sampling increment? %6ld\n",
3280 options->lincrement);
3281 printf
3282 (" 6 Number of recorded genealogies in long chain? %6ld\n",
3283 options->lsteps);
3284 }
3285 switch (options->burnin_autostop)
3286 {
3287 case 'a':
3288 sprintf(temp,"(stopping crit: variance) %10li", options->burn_in);
3289 break;
3290 case 'e':
3291 sprintf(temp,"(stopping crit: ESS) %10li", options->burn_in);
3292 break;
3293 case ' ':
3294 default:
3295 sprintf(temp," %10li", options->burn_in);
3296 }
3297 printf(" 7 Burn-in for each chain:%30.30s\n",temp);
3298 if (!options->replicate)
3299 printf
3300 (" 8 Combine chains or runs for estimates? NO\n");
3301 else {
3302 if (options->replicatenum == 0)
3303 printf
3304 (" 8 Combine chains for estimates? YES, long chains\n");
3305 else
3306 printf
3307 (" 8 Combine chains for estimates? YES, over %3li runs\n",
3308 options->replicatenum);
3309 }
3310 if (options->heating == 0)
3311 printf
3312 (" 9 Heating: NO\n");
3313 else {
3314 if (options->adaptiveheat!=NOTADAPTIVE)
3315 sprintf(temp, "Adaptive (%s %3li chains, swap interval is %3li)\n",
3316 options->adaptiveheat == STANDARD ? "Standard" : "Bounded",
3317 options->heated_chains, options->heating_interval);
3318 else
3319 sprintf(temp, " YES (%3li chains, swap interval is %3li)\n",
3320 options->heated_chains, options->heating_interval);
3321 printf(" 9 Heating: %s", temp);
3322 }
3323 printf
3324 ("\n -------------------------------------------------------------\n");
3325 printf(" Obscure options (consult the documentation on these)\n\n");
3326 if (options->movingsteps)
3327 printf
3328 (" 10 Sample a fraction of %2.2f new genealogies? YES\n",
3329 options->acceptfreq);
3330 else
3331 printf
3332 (" 10 Sample at least a fraction of new genealogies? NO\n");
3333 if (options->lcepsilon < LONGCHAINEPSILON)
3334 sprintf(temp, "%17.2f", options->lcepsilon);
3335 else
3336 sprintf(temp, "%17.17s", "infinity");
3337 printf(" 11 Epsilon of parameter likelihood %s\n", temp);
3338 printf(" 12 Use Gelman's convergence criterium? %13s\n",
3339 options->gelman ? (options->gelmanpairs ? "YES:Pairs" : "YES:Summary" ) : " NO");
3340 }
3341
3342 void
display_bayes_mcmc(option_fmt * options)3343 display_bayes_mcmc(option_fmt * options) {
3344 char *outputstring;
3345 outputstring = (char *) mycalloc(LINESIZE, sizeof(char));
3346
3347 printf(" %2i Strategy:%54.54s\n", BAYESSTRATEGY, options->bayes_infer ?
3348 "Bayesian Inference" :
3349 "Maximum Likelihood");
3350 if(options->has_bayesfile)
3351 sprintf(outputstring,"%s%s", "YES:", options->bayesfilename);
3352 else
3353 sprintf(outputstring,"NO");
3354 printf( " %2i File for recording posterior distribution?%21s\n", BAYESOUT, outputstring);
3355 if(options->has_bayesmdimfile)
3356 sprintf(outputstring,"%s%s", "YES:", options->bayesmdimfilename);
3357 else
3358 sprintf(outputstring,"NO");
3359 printf( " %2i File for recording all parameter values? %21s\n", BAYESMDIMOUT, outputstring);
3360 if(options->has_bayesmdimfile)
3361 {
3362 sprintf(outputstring,"[save every %li x sample-increment = %li]", options->bayesmdiminterval,
3363 options->bayesmdiminterval*options->lincrement);
3364 printf( " %58.58s\n", outputstring);
3365 }
3366 sprintf(outputstring, "%li, %li", options->bayespriortheta->bins, options->bayespriorm->bins);
3367 printf(" %2i Number of bins of posterior [Theta,M]?%25s\n", BAYESBINNING, outputstring);
3368 printf(" %2i Plotting type of posterior distribution?%23s\n", BAYESPRETTY,
3369 (options->bayespretty == PRETTY_P99 ? "up to 99% percentile" :
3370 (options->bayespretty == PRETTY_MAX ? "prior distr. range" :
3371 (options->bayespretty == PRETTY_P100 ? "up to ~100% percentile" : "up to maximal 99%"))));
3372 printf(" %2i Frequency of tree updates vs. parameter updates? %6.2f\n",
3373 BAYESFREQ, options->updateratio);
3374 set_proposal(outputstring, options->slice_sampling, !options->bayesmurates);
3375 printf(" %2i Proposal distribution?%41.41s\n",
3376 BAYESPROPOSAL, outputstring);
3377 set_prior(outputstring, options->bayesprior, !options->bayesmurates);
3378 printf(" %2i Prior distribution?%44.44s\n",
3379 BAYESPRIOR, outputstring);
3380 printf(" %2i Number of long chains to run? %6ld\n",
3381 BAYESLCHAINS, options->lchains);
3382 printf(" %2i Sampling increment? %6ld\n",
3383 BAYESSKIP, options->lincrement);
3384 printf(" %2i Number of recorded steps in chain %12ld\n",
3385 BAYESSAMPLES, options->lsteps);
3386 switch (options->burnin_autostop)
3387 {
3388 case 'a':
3389 sprintf(outputstring,"(stopping crit: variance) %10li", options->burn_in);
3390 break;
3391 case 'e':
3392 sprintf(outputstring,"(stopping crit: ESS) %10li", options->burn_in);
3393 break;
3394 case ' ':
3395 default:
3396 sprintf(outputstring," %10li", options->burn_in);
3397 }
3398 printf(" %2i Burn-in for each chain:%38.38s\n",
3399 BAYESBURNIN,outputstring);
3400
3401 if(!options->replicate)
3402 sprintf(outputstring,"NO");
3403 else
3404 sprintf(outputstring,"YES (%li independent chains)",options->replicatenum);
3405 printf(" %2i Running multiple replicates: %33.33s\n", BAYESREPLICATE, outputstring);
3406
3407 if(options->heating == 0)
3408 printf(" %2i Heating: %43.43s\n", BAYESHEAT, "NO");
3409 else
3410 printf(" %2i Heating: %10s (%3li parallel chains)\n",
3411 BAYESHEAT, (options->adaptiveheat==STANDARD ? "ADAPTIVE" : (options->adaptiveheat==BOUNDED ? "BOUNDED" : "STATIC")) , options->heated_chains);
3412 sprintf(outputstring,"%f",options->acceptfreq);
3413 printf(" %2i Sampling at least fraction of new genealogies: %10.10s\n", BAYESMOVINGSTEPS, outputstring);
3414 sprintf(outputstring,"%s",options->gelman ? (options->gelmanpairs ? "YES:Pairs" : "YES:Summary") : "NO");
3415 printf(" %2i Convergence diagnostic for replicates: %13.13s\n", BAYESGELMAN, outputstring);
3416
3417 printf(" %2i Run analysis without data: %13.13s\n", BAYESPRIORALONE, options->prioralone ? "YES" : "NO");
3418
3419 myfree(outputstring);
3420 }
3421
3422 void
display_lratio_menutext(void)3423 display_lratio_menutext(void) {
3424 printf(" H0: the ML-estimates and your values are the same\n");
3425 printf(" H1: the ML-estimates and your values are different\n");
3426 printf(" You need to specify values for ALL parameters\n");
3427 printf(" Beware! If you restrict the migration model then some test\n");
3428 printf(" will not work. This test assumes that your hypothesis\n");
3429 printf(" is nested in the hypothesis you test against.\n");
3430 printf(" Currently you cannot specify \"<\" or \">\" comparisons AND\n");
3431 printf(" averages of migration rates are calculated and not\n");
3432 printf(" averages for number of migrants\n\n");
3433 }
3434
3435 void
display_lrt_syntax(void)3436 display_lrt_syntax(void) {
3437 printf(" Specify values as an {n x n} matrix\n");
3438 printf(" Theta values are on the diagonal, migration rates are\n");
3439 printf(" off-diagonal, newlines are allowed, and\n");
3440 printf(" values must be separated by space or commas or tabs.\n");
3441 printf("\n Syntax:\n");
3442 printf(" * = sames as the MLE, and not estimated again\n");
3443 printf(" s = averages of MIGRATION RATES from i->j and j->i\n");
3444 printf(" m = averages of THETAS or MIGRATION RATES [!!!]\n");
3445 printf(" value = specify your own value\n\n");
3446 printf(" [If you want to leave this editor, type \"quit\" or \"end\" and <return>\n");
3447 //To review the Syntax type \ "syntax\"\n");
3448 }
3449
3450 char
ask_lratio_type(lr_data_fmt * data,long counter)3451 ask_lratio_type(lr_data_fmt * data, long counter) {
3452 char answer = '\0';
3453 char input[LINESIZE];
3454 printf(" Is the test against the Maximum likelihood estimates or\n");
3455 printf(" or against arbitrary values? [(M)LE, (A)rbitrary]\n===> ");
3456 fflush(stdout); FGETS(input, 1024, stdin);
3457 switch (uppercase(input[0])) {
3458 case 'M':
3459 data[counter].type = MLE;
3460 answer = 'M';
3461 break;
3462 case 'A':
3463 answer = 'A';
3464 data[counter].type = ARBITRARY;
3465 break;
3466 case 'Q':
3467 case 'E':
3468 answer = 'Q';
3469 break;
3470 case 'S':
3471 answer = 'S';
3472 break;
3473 default:
3474 data[counter].type = MLE;
3475 break;
3476 }
3477 return answer;
3478 }
3479
3480 void
read_custom_menu_lratio(option_fmt * options)3481 read_custom_menu_lratio(option_fmt * options) {
3482 //static long numpop = 0, numpop2 = 0;
3483 //char input[LINESIZE];
3484 //long i, z = 0, pointer;
3485 //char *temp;
3486 //char *valpointer;
3487 char answer = '\0';
3488 //boolean done = FALSE;
3489 lratio_fmt *lratio = options->lratio;
3490 display_lratio_menutext();
3491 display_lrt_syntax();
3492 while (answer != 'Q') {
3493 //answer = ask_lratio_type(lratio->data, lratio->counter);
3494 answer = 'M';
3495 lratio->data[lratio->counter].type = MLE;
3496 switch (answer) {
3497 case 'S':
3498 display_lrt_syntax();
3499 break;
3500 case 'E':
3501 case 'Q':
3502 break;
3503 case 'A':
3504 case 'M':
3505 how_many_pop(&options->numthetag);
3506 answer = dialog_lrt(options->numthetag, lratio);
3507 break;
3508 }
3509 }
3510 }
3511
3512 void
how_many_pop(long * numpop)3513 how_many_pop(long *numpop) {
3514 char input[LINESIZE];
3515 if (*numpop == 0) {
3516 do {
3517 printf(" How many populations?\n===> ");
3518 fflush(stdout); FGETS(input, 1024, stdin);
3519 *numpop = atoi(input);
3520 }
3521 while (*numpop <= 0 && *numpop < 100);
3522 } else
3523 printf(" The data set contains %li populations\n", *numpop);
3524 }
3525
3526 void
check_lrt_allocation(lratio_fmt * lratio)3527 check_lrt_allocation(lratio_fmt * lratio) {
3528 long i;
3529 if (lratio->counter + 1 == lratio->alloccounter) {
3530 lratio->alloccounter += 2;
3531 lratio->data =
3532 (lr_data_fmt *) myrealloc(lratio->data, sizeof(lr_data_fmt) *
3533 lratio->alloccounter);
3534 for (i = lratio->counter + 1; i < lratio->alloccounter; i++) {
3535 lratio->data[i].elem = 0;
3536 lratio->data[i].value1 =
3537 (char *) mycalloc(1, sizeof(char) * LONGLINESIZE);
3538 lratio->data[i].value2 =
3539 (char *) mycalloc(1, sizeof(char) * LONGLINESIZE);
3540 lratio->data[i].connect =
3541 (char *) mycalloc(1, sizeof(char) * LONGLINESIZE);
3542
3543 }
3544 }
3545 }
3546
3547 char
menuread_lrt_paramvalue(char * value,long * counter,long numpop2)3548 menuread_lrt_paramvalue(char *value, long *counter, long numpop2) {
3549 char *valptr = value;
3550 long pointer = 0;
3551 long z = 0;
3552 char *temp;
3553 char input[LONGLINESIZE];
3554 while (z < numpop2) {
3555 fflush(stdout); FGETS(input, LONGLINESIZE, stdin);
3556 temp = strtok(input, ", \n\t");
3557 if (temp != NULL) {
3558 if (strchr("EQ", uppercase(temp[0]))) {
3559 //(*counter)--;
3560 return 'Q';
3561 }
3562 }
3563 while (temp != NULL) {
3564 sprintf(valptr + pointer, " %s,", temp);
3565 z++;
3566 pointer += (long) strlen(temp) + 2;
3567 temp = strtok(NULL, ", \n\t");
3568 }
3569 }
3570 return ' ';
3571 }
3572
3573 char
dialog_lrt(long numpop,lratio_fmt * lratio)3574 dialog_lrt(long numpop, lratio_fmt * lratio) {
3575 //char input[LINESIZE];
3576 //long z = 0;
3577 char answer = ' ';
3578 //char *valptr1, *valptr2;
3579 long numpop2 = numpop * numpop;
3580 check_lrt_allocation(lratio);
3581 printf(" %li. Likelihood ratio test\n", lratio->counter);
3582 printf(" Enter now the %li values for the FIRST parameter set:\n",
3583 numpop2);
3584 lratio->data[lratio->counter].type = MLE;
3585 answer = menuread_lrt_paramvalue(lratio->data[lratio->counter].value1,
3586 &lratio->counter, numpop2);
3587 if (answer == 'Q')
3588 return answer;
3589 if (lratio->data[lratio->counter].type == MLE) {
3590 lratio->counter++;
3591 return answer;
3592 }
3593 printf(" Enter now the values for the SECOND parameter set:\n");
3594 answer = menuread_lrt_paramvalue(lratio->data[lratio->counter].value1,
3595 &lratio->counter, numpop2);
3596 if (answer == 'Q')
3597 return answer;
3598 lratio->counter++;
3599 return answer;
3600 }
3601
3602 void
read_custom_menu_migration(option_fmt * options)3603 read_custom_menu_migration(option_fmt * options) {
3604 char input[LINESIZE];
3605 long z = 0, numpop = 0, numpop2;
3606 printf(" Specify the migration model as an {n x n} matrix\n");
3607 printf(" Theta values are on the diagonal, migration rates are\n");
3608 printf(" off-diagonal, spaces (\" \"), \"{\", \"}\", or newlines\n");
3609 printf(" are allowed, but not necessary.\n");
3610 printf("\n Syntax:\n");
3611 printf(" * = independent parameter\n");
3612 printf(" 0 = (zero) not estimated]\n");
3613 printf(" c = (constant) not estimated, taken from start-parameter]\n");
3614 printf(" s = symmetric migration rates (M=m/mu)\n");
3615 printf(" S = symmetric migration rates (xNm) \n");
3616 printf(" m = average of each label group [not k, or s]\n");
3617 // printf(" a,b,c,... = average of each label group [not k, or s]\n");
3618 // printf(" 1,2,3,... = if both migration rates are labeled then combine populations\n");
3619 //printf(" M = average, either ALL thetas and/or ALL migration rates\n");
3620
3621 if (options->numthetag > 0)
3622 numpop = options->numthetag;
3623 else {
3624 if (options->nummg > 0)
3625 numpop =
3626 (long) ((1. + sqrt(1. + (MYREAL) options->nummg * 4.)) / 2.);
3627 else {
3628 how_many_pop(&options->numthetag);
3629 numpop = options->numthetag;
3630 }
3631 }
3632 numpop2 = numpop * numpop;
3633 printf("\n You must give %li values\n===> ", numpop2);
3634 while (z < numpop2) {
3635 fflush(stdout); FGETS(input, 1024L, stdin);
3636 read_custom_migration(stdin, options, input, numpop);
3637 z = (long) strlen(options->custm);
3638 }
3639 }
3640
3641 char *
custom_migration_type(long type)3642 custom_migration_type(long type) {
3643 switch (type) {
3644 case MATRIX:
3645 return "Full migration matrix model";
3646 break;
3647 case MATRIX_SYMMETRIC:
3648 return "Symmetric migration matrix model";
3649 break;
3650 case MATRIX_SAMETHETA:
3651 return "Full migration matrix model (same Theta)";
3652 break;
3653 case MATRIX_ARBITRARY:
3654 return "User specified migration matrix model";
3655 break;
3656 case ISLAND:
3657 return "N-population island model";
3658 break;
3659 case ISLAND_VARTHETA:
3660 return "N-population island model (variable Theta)";
3661 break;
3662 case STEPSTONE:
3663 return "Stepping stone model";
3664 break;
3665 case CONTINUUM:
3666 return "Continuum model";
3667 break;
3668 case NEIGHBOR:
3669 return "Isolation by distance model";
3670 break;
3671 default:
3672 return "Illegal migration model";
3673 break;
3674 }
3675 }
3676
3677
3678 void
read_heatvalues(option_fmt * options)3679 read_heatvalues(option_fmt * options) {
3680 long i, z;
3681 char *tmp;
3682 char input[LINESIZE];
3683 MYREAL diff = 0.;
3684 MYREAL x;
3685 FPRINTF(stdout, " ");
3686 printf("Enter %li \"temperatures\"\n", options->heated_chains);
3687 printf("[The coldest temperature, which is the first, has to be 1\n");
3688 printf(" For example: 1.0 1.5 3.0 1000000.0]\n");
3689
3690 FPRINTF(stdout,
3691 "OR give a range of values [linear increase: 1 - 10]\n");
3692 FPRINTF(stdout,
3693 " [exponential increase: 1 @ 10]\n");
3694 FPRINTF(stdout,
3695 "or, most lazily, let me suggest a range [simply type a #]\n");
3696 FPRINTF(stdout,"@@@@@ For model comparison, the range of temperatures @@@@@\n");
3697 FPRINTF(stdout,"@@@@@ MUST include a very hot chain (>100000.0) @@@@@\n");
3698 printf("===> ");
3699 fflush(stdout); FGETS(input, LINESIZE, stdin);
3700 if (strstr(input, "-"))
3701 {
3702 tmp = strstr(input, "-");
3703 options->heat[options->heated_chains - 1] = fabs(atof(tmp + 1));
3704 options->heat[0] = 1.0;
3705 diff = options->heat[options->heated_chains - 1] - options->heat[0];
3706 diff /= options->heated_chains - 1.;
3707 for (i = 1; i < options->heated_chains; i++)
3708 options->heat[i] = options->heat[i - 1] + diff;
3709 }
3710 else
3711 {
3712 if (strstr(input, "#"))
3713 {
3714 diff = 1./(options->heated_chains - 1);
3715 x = 1 - diff;
3716 for (i = 1; i < options->heated_chains-1; i++)
3717 {
3718 options->heat[i] = 1./ x;
3719 x -= diff;
3720 }
3721 options->heat[i] = 1000000.;
3722 }
3723 else
3724 {
3725 if (strstr(input, "@"))
3726 {
3727 tmp = strstr(input, "@");
3728 options->heat[options->heated_chains - 1] = fabs(atof(tmp + 1));
3729 options->heat[0] = 1.0;
3730 diff = 2. * (options->heat[options->heated_chains - 1] -
3731 options->heat[0]);
3732 diff /= pow(2.0, (MYREAL) options->heated_chains) - 2.;
3733 for (i = 1; i < options->heated_chains; i++)
3734 options->heat[i] =
3735 options->heat[i - 1] + pow(2.0, i - 1.) * diff;
3736 }
3737 else
3738 {
3739 z = 0;
3740 tmp = strtok(input, " ,");
3741 if (tmp != NULL)
3742 {
3743 options->heat[z++] = atof(tmp);
3744 }
3745 for (;;)
3746 {
3747 tmp = strtok(NULL, " ,");
3748 if (tmp != NULL)
3749 {
3750 options->heat[z++] = atof(tmp);
3751 }
3752 else
3753 {
3754 if (z < options->heated_chains)
3755 {
3756 printf("%li more values are needed\n===> ",
3757 options->heated_chains - z);
3758 fflush(stdout); FGETS(input, LINESIZE, stdin);
3759 tmp = strtok(input, " ,");
3760 if (tmp != NULL)
3761 {
3762 options->heat[z++] = atof(tmp);
3763 }
3764 }
3765 else
3766 break;
3767 }
3768 }
3769 }
3770 }
3771 }
3772 printf("Chain Temperature\n");
3773 printf("-------------------------\n");
3774 for (i = options->heated_chains - 1; i >= 0; --i)
3775 printf("%4li %20.5f\n", i+1, options->heat[i]);
3776 printf("\n\nIs this correct? [YES or NO]\n===> ");
3777 fflush(stdout); FGETS(input, LINESIZE, stdin);
3778 if ('y' != tolower(input[0]))
3779 read_heatvalues(options);
3780 }
3781
3782 void
menuRandom(MYREAL * param,char type)3783 menuRandom(MYREAL * param, char type) {
3784 //long check;
3785 int retval;
3786 if(type == 'N')
3787 {
3788 printf
3789 ("Specify a MEAN and a STANDARD DEVIATION\nThis will be used to generate a random value from a Normal distribution\n");
3790 }
3791 else
3792 {
3793 printf
3794 ("Specify a MINIMUM and a MAXIMUM\nThis will be used to generate a random value from a Uniform distribution\n");
3795 }
3796 #ifdef USE_MYREAL_FLOAT
3797 retval = scanf("%f%f", ¶m[0], ¶m[1]);
3798 #else
3799 retval = scanf("%lf%lf", ¶m[0], ¶m[1]);
3800 #endif
3801 }
3802
3803 void
start_tree_method(option_fmt * options)3804 start_tree_method(option_fmt * options) {
3805 char input[LINESIZE];
3806 char compstring[LINESIZE];
3807 printf(" Start genealogy is created with\n");
3808 printf(" (a)utomatic\n");
3809 if (strchr(SEQUENCETYPES, options->datatype)) {
3810 printf(" (u)sertree\n");
3811 printf(" (d)istancematrix\n");
3812 }
3813 printf(" (r)andom\n\n===> ");
3814 strcpy(compstring, (strchr(SEQUENCETYPES, options->datatype) ?
3815 "daur" : "dar"));
3816 do {
3817 fflush(stdout); FGETS(input, LINESIZE, stdin);
3818 }
3819 while (strchr(compstring, (int) (lowercase(input[0]))) == NULL);
3820 switch (lowercase(input[0])) {
3821 case 'u':
3822 options->usertree = TRUE;
3823 options->dist = FALSE;
3824 break;
3825 case 'd':
3826 options->usertree = FALSE;
3827 options->dist = TRUE;
3828 options->randomtree = FALSE;
3829 break;
3830 case 'a':
3831 options->usertree = FALSE;
3832 options->dist = FALSE;
3833 options->randomtree = FALSE;
3834 break;
3835 case 'r':
3836 options->usertree = FALSE;
3837 options->dist = FALSE;
3838 options->randomtree = TRUE;
3839 break;
3840 }
3841 }
3842
submodeltype(int type)3843 char * submodeltype(int type)
3844 {
3845 switch(type)
3846 {
3847 case MULTISTEP:
3848 return "Multistep method";
3849 case SINGLESTEP:
3850 return "Singlestep method";
3851 default:
3852 return "Singlestep method";
3853 }
3854 }
3855
3856 /// \brief Menu for exponential proposal
3857 /// Menu for exponential proposal
menu_submodel(option_fmt * options)3858 void menu_submodel(option_fmt *options)
3859 {
3860 char input[LINESIZE];
3861 char text[LINESIZE];
3862 char val;
3863 double tune = 0. ;
3864 double pinc = 0.5;
3865 long count = 0;
3866 do {
3867 if(options->msat_option == MULTISTEP)
3868 {
3869 sprintf(text,"MULTISTEP (tune=%f, p_increase=%f)",options->msat_tuning[0], options->msat_tuning[1]);
3870 }
3871 else
3872 {
3873 sprintf(text,"SINGLESTEP");
3874 }
3875 printf("\nSingle step or Multi step model? [ENTER either S or M]\n");
3876 printf("Single Step model is the standard stepwise mutation model\n");
3877 printf("The multistep model is between the infinite model and and the\n");
3878 printf("singlestep model using two additional parameters:\n");
3879 printf("- chance of increasing repeat length\n");
3880 printf("- tuning between single step (tune=0) and infinite allele (tune=1)\n");
3881 printf("[Current value: %s]\n===> ",text);
3882 fflush(stdout); FGETS(input, LINESIZE, stdin);
3883 if (input[0] == '\0')
3884 break;
3885 count = sscanf(input, "%c", &val);
3886 if(uppercase(val) == 'S')
3887 options->msat_option = SINGLESTEP;
3888 else
3889 {
3890 options->msat_option = MULTISTEP;
3891 do
3892 {
3893 printf("MULTISTEP model\nGive two numbers separated by a space:\nFor the tune (range 0 to 1) and\nfor the probability of repeat number increase (0 to 2/3)\n===> ");
3894 fflush(stdout); FGETS(input, LINESIZE, stdin);
3895 count = sscanf(input, "%lf%lf", &tune, &pinc);
3896 options->msat_tuning[0] = (MYREAL) tune;
3897 options->msat_tuning[1] = (MYREAL) pinc;
3898 } while (count != 2 && count > 0);
3899 count = 1;
3900 }
3901 } while (count != 1 && count > 0);
3902 }
3903
3904
3905 void
start_data_method(option_fmt * options)3906 start_data_method(option_fmt * options) {
3907 char input[LINESIZE];
3908 do {
3909 printf(" (a)llele model\n");
3910 printf(" (m)icrosatellite model [SLOW! Ladder model: %s]\n", submodeltype(options->msat_option));
3911 printf(" (b)rownian microsatellite model [Brownian motion model]\n");
3912 printf(" (s)equence model\n");
3913
3914 printf(" (n)ucleotide polymorphism (SNP)\n");
3915 printf(" (h)apmap data format of SNPs\n");
3916 printf(" (g)enealogy summaries\n\n===> ");
3917 fflush(stdout);
3918 FGETS(input, LINESIZE, stdin);
3919 }
3920 while (strchr("ambnsgufh", (int) (lowercase(input[0]))) == NULL);
3921 options->datatype = input[0];
3922
3923 if(options->datatype == 'm')
3924 {
3925 menu_submodel(options);
3926 }
3927 if (!strchr(SEQUENCETYPES, options->datatype))
3928 options->usertree = FALSE;
3929 }
3930
get_prior(char * input)3931 long get_prior(char *input) {
3932 switch (toupper(input[0])) {
3933 //case 'S':
3934 //return SLICE;
3935 //break;
3936 case 'M':
3937 return MULTPRIOR;
3938 break;
3939 case 'E':
3940 return EXPPRIOR;
3941 break;
3942 case 'W':
3943 return WEXPPRIOR;
3944 break;
3945 case 'U':
3946 return UNIFORMPRIOR;
3947 default:
3948 return UNIFORMPRIOR;
3949 break;
3950 }
3951 }
3952
3953 ///
3954 /// prints type of proposal in the menu
set_proposal(char * output,boolean * proposal,boolean without_rate)3955 void set_proposal(char *output, boolean *proposal, boolean without_rate)
3956 {
3957 int i;
3958 int l=0;
3959 const char types[][LINESIZE]={"Theta","Mig","Rate"};
3960
3961 for(i=THETAPRIOR; i <= RATEPRIOR; i++)
3962 {
3963 if(without_rate && i==RATEPRIOR)
3964 continue;
3965 switch (proposal[i]) {
3966 case TRUE:
3967 l+=sprintf(output+l, " %s:Slice", types[i]);
3968 break;
3969 case FALSE:
3970 default:
3971 l+=sprintf(output+l, " %s:MH", types[i]);
3972 break;
3973 // default:
3974 //l+=sprintf(output+l, " %s:?", types[i]);
3975 //break;
3976 }
3977 }
3978 }
3979
3980 ///
3981 /// prints type of prior in the menu
set_prior(char * output,int * prior,boolean without_rate)3982 void set_prior(char *output, int *prior, boolean without_rate)
3983 {
3984 int i;
3985 int l=0;
3986 const char types[][LINESIZE]={"Theta","Mig","Rate"};
3987 for(i=THETAPRIOR; i <= RATEPRIOR; i++)
3988 {
3989 if(without_rate && i==RATEPRIOR)
3990 continue;
3991
3992 switch (prior[i]) {
3993 //case SLICE:
3994 //l+=sprintf(output+l, " Slice");
3995 //break;
3996 case MULTPRIOR:
3997 l+=sprintf(output+l, " %s:Mult.", types[i]);
3998 break;
3999 case EXPPRIOR:
4000 l+=sprintf(output+l, " %s:Exp.", types[i]);
4001 break;
4002 case WEXPPRIOR:
4003 l+=sprintf(output+l, " %s:WExp.", types[i]);
4004 break;
4005 case GAMMAPRIOR:
4006 l+=sprintf(output+l, " %s:Gamma", types[i]);
4007 break;
4008 case UNIFORMPRIOR:
4009 default:
4010 l+=sprintf(output+l, " %s:Unif.", types[i]);
4011 break;
4012 }
4013 }
4014 }
4015
set_localities_string(char * loc,option_fmt * options)4016 void set_localities_string(char *loc, option_fmt *options)
4017 {
4018 long z;
4019 long i;
4020 if(options->newpops_numalloc>1)
4021 {
4022 z=0;
4023 for(i=0; i<options->newpops_numalloc; i++)
4024 {
4025 z += sprintf(loc+z,"%li,",options->newpops[i]);
4026 }
4027 }
4028 else
4029 {
4030 sprintf(loc,"default");
4031 }
4032 }
4033