1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate and effectice population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6 H E L P E R R O U T I N E S
7
8 some math stuff and
9 string and file manipulation routines
10
11
12 Peter Beerli started 1996, Seattle
13 beerli@fsu.edu
14
15 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
16 Copyright 2003-2006 Peter Beerli, Tallahassee FL
17
18 Permission is hereby granted, free of charge, to any person obtaining
19 a copy of this software and associated documentation files (the
20 "Software"), to deal in the Software without restriction, including
21 without limitation the rights to use, copy, modify, merge, publish,
22 distribute, sublicense, and/or sell copies of the Software, and to
23 permit persons to whom the Software is furnished to do so, subject
24 to the following conditions:
25
26 The above copyright notice and this permission notice shall be
27 included in all copies or substantial portions of the Software.
28
29 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
30 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
31 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
32 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
33 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
34 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
35 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
36
37
38 $Id: tools.c 2158 2013-04-29 01:56:20Z beerli $
39 -------------------------------------------------------*/
40 /* \file tools.c */
41
42 #include <sys/types.h>
43 #include <sys/stat.h>
44
45 #include "migration.h"
46 #include "sighandler.h"
47 #include "data.h"
48 #include "world.h"
49 #include "tree.h"
50 #include "random.h"
51 #include "broyden.h"
52 #include "migrate_mpi.h"
53 #include "options.h"
54 #include "laguerre.h"
55 #include "pretty.h"
56
57 #include <stdio.h>
58 #include <math.h>
59 #ifdef DMALLOC_FUNC_CHECK
60 #include <dmalloc.h>
61 #endif
62
63 /* external globals*/
64 extern int myID;
65 #ifdef CAUTIOUS
66 extern boolean cautious;
67 #endif
68 /* prototypes ------------------------------------------- */
69 MYREAL lengthof (node * p);
70 MYINLINE node *crawlback (const node * theNode);
71 /*node *crawl(node * theNode); */
72 MYINLINE node *showtop (node * theNode);
73 void adjust_time (node * theNode, MYREAL tyme);
74 void adjust_time_all (node * theNode, MYREAL tyme);
75 void insert_migr_node (world_fmt * world, node * up, node * down,
76 migr_table_fmt * migr_table, long *migr_table_counter);
77 void children (node * mother, node ** brother, node ** sister);
78 /* math tools */
79 MYREAL incompletegamma (MYREAL tx, MYREAL talpha);
80 MYREAL polygamma (long n, MYREAL z);
81 void invert_matrix (MYREAL **a, long nsize);
82 boolean nrcheck (MYREAL **m, MYREAL **tm, MYREAL *v, long nrows, MYREAL *r1,
83 MYREAL *r2, boolean do_newton);
84 MYREAL rannor (MYREAL mean, MYREAL sd);
85 char lowercase (int c);
86 char uppercase (int c);
87 MYREAL calc_sum (MYREAL *vector, long n);
88
89 void gamma_rates (MYREAL *rate, MYREAL *probcat, long categs, char *input);
90 void calc_gamma (MYREAL alpha, MYREAL *gama, long categs);
91
92 MYREAL mylgamma (MYREAL z);
93
94 MYREAL logfac (long n);
95
96 void onepass_mean_std_start(MYREAL *mean, MYREAL *std, long *n);
97 void onepass_mean_std_calc(MYREAL *mean, MYREAL *std, long *n, MYREAL x);
98 void onepass_mean_std_end(MYREAL *mean, MYREAL *std, long *n);
99
100 //MYINLINE double fast_exp(double y) ;
101 //MYINLINE float fast_log (float vval);
102
103 MYREAL bezier(MYREAL x0, MYREAL y0, MYREAL x1, MYREAL y1);
104
105 /* vector initialization */
106 void doublevec1d (MYREAL **v, long size);
107 void doublevec2d (MYREAL ***v, long size1, long size2);
108 void floatvec2d (float ***v, long size1, long size2);
109 void charvec1d (char **v, long size1);
110 void charvec2d (char ***v, long size1, long size2);
111 void intvec2d (long ***v, long size1, long size2);
112 void free_doublevec2d (MYREAL **v);
113 void free_floatvec2d (float **v);
114 void free_charvec2d (char **v);
115 void free_intvec2d (long **v);
116 void setdoublevec1d (MYREAL **v, MYREAL *w, long size);
117 void add_vector (MYREAL *result, MYREAL *v, long size);
118
119 /*filemanipulation */
120 void init_files (world_fmt * world, data_fmt * data, option_fmt * options);
121 void exit_files (world_fmt * world, data_fmt * data, option_fmt * options);
122 void openfile (FILE ** fp, char *filename, char *mode, char *perm);
123 #ifdef ZNZ
124 void znzopenfile (znzFile * fp, char *filename, char *mode, int use_compressed);
125 #endif
126 long read_savesum (world_fmt * world, option_fmt * options, data_fmt * data);
127 void write_savesum (world_fmt * world);
128
129 /* string manipulation */
130 void translate (char *text, char from, char to);
131 long count_words (char *text);
132 void fprintf2(FILE *file, long filesize, const char *fmt, ...);
133 void print_line (FILE * outfile, char c, long nn, long flag);
134 void sprint_line (char *buffer, char c, long nn, long flag);
135 void add_to_buffer(char *fp, long *bufsize, char **buffer, long *allocbufsize);
136
137 /* time reporting */
138 void get_time (char *nowstr, char ts[]);
139 /*printing aid */
140 void print_llike (MYREAL llike, char *strllike);
141
142 /* searching and finding*/
143 boolean find (long i, long *list, long listlen);
144
145 /* conversion between the parameter schemes*/
146
147 long mstart (long pop, long numpop);
148 long mend (long pop, long numpop);
149 long mmstart (long pop, long numpop);
150 long mmend (long pop, long numpop);
151 long mm2m (long frompop, long topop, long numpop);
152 void m2mm (long i, long numpop, long *frompop, long *topop);
153 long m2mml (long i, long numpop);
154 long m2mml2 (long i, long topop, long numpop);
155
156 void set_paramstr(char * xsparamstr, long j, long numpop2, boolean usem);
157
158 /* private functions */
159 MYREAL alnorm (MYREAL x, int up);
160 void lu_decomp (MYREAL **m, long *indeks, long nrows);
161 void lu_substitution (MYREAL **m, long *indeks, MYREAL *v, long nrows);
162 MYREAL d1mach (long i);
163 long i1mach (long i);
164 int dpsifn (MYREAL *x, long *n, long kode, long m, MYREAL *ans, long *nz,
165 long *ierr);
166 MYREAL find_chi (long df, MYREAL prob);
167 MYREAL probchi (long df, double chi);
168 MYREAL chisquare (long df, MYREAL alpha);
169
170 MYREAL chiboundary (long zeros, long nonzeros, MYREAL alpha);
171
172 extern long calculate_newpop_numpop(option_fmt *options, data_fmt *data);
173
174 #ifdef MESS
175 long *check_collection;
176 long check_collection_count;
177 extern long unique_id_global;
178 #endif
179
180 /* vector initialization */
181 void
doublevec1d(MYREAL ** v,long size)182 doublevec1d (MYREAL **v, long size)
183 {
184 *v = (MYREAL *) mycalloc (size, sizeof (MYREAL));
185 }
186
187 /// allocate a 2D array with MYREALS
188 void
doublevec2d(MYREAL *** v,long size1,long size2)189 doublevec2d (MYREAL ***v, long size1, long size2)
190 {
191 long i;
192 *v = (MYREAL **) mycalloc (size1, sizeof (MYREAL *));
193 (*v)[0] = (MYREAL *) mycalloc (size1 * size2, sizeof (MYREAL));
194 for (i = 1; i < size1; i++)
195 {
196 (*v)[i] = (*v)[0] + size2 * i;
197 }
198 }
199
200 /// allocate a 2D array with floats
201 void
floatvec2d(float *** v,long size1,long size2)202 floatvec2d (float ***v, long size1, long size2)
203 {
204 long i;
205 *v = (float **) mycalloc (size1, sizeof (float *));
206 (*v)[0] = (float *) mycalloc (size1 * size2, sizeof (float));
207 for (i = 1; i < size1; i++)
208 {
209 (*v)[i] = (*v)[0] + size2 * i;
210 }
211 }
212
213 ///
214 /// allocates columns and rows for a matrix of strings
215 /// use freecharvec2d() to free this.
216 void
charvec2d(char *** v,long size1,long size2)217 charvec2d (char ***v, long size1, long size2)
218 {
219 long i;
220 *v = (char **) mycalloc (size1, sizeof (char *));
221 (*v)[0] = (char *) mycalloc (size1 * size2, sizeof (char));
222 for (i = 1; i < size1; i++)
223 {
224 (*v)[i] = (*v)[0] + size2 * i;
225 }
226 }
227 ///
228 /// allocates columns and rows for a matrix of strings
229 /// use free_intvec2d() to free this.
230 void
intvec2d(long *** v,long size1,long size2)231 intvec2d (long ***v, long size1, long size2)
232 {
233 long i;
234 *v = (long **) mycalloc (size1, sizeof (long *));
235 (*v)[0] = (long *) mycalloc (size1 * size2, sizeof (long));
236 for (i = 1; i < size1; i++)
237 {
238 (*v)[i] = (*v)[0] + size2 * i;
239 }
240 }
241
242 /// free 2D vector of MYREAL
243 void
free_doublevec2d(MYREAL ** v)244 free_doublevec2d (MYREAL **v)
245 {
246 myfree(v[0]);
247 myfree(v);
248 }
249
250 /// free 2D vector of floats
251 void
free_floatvec2d(float ** v)252 free_floatvec2d (float **v)
253 {
254 myfree(v[0]);
255 myfree(v);
256 }
257
258 void
free_charvec2d(char ** v)259 free_charvec2d (char **v)
260 {
261 myfree(v[0]);
262 myfree(v);
263 }
264 void
free_intvec2d(long ** v)265 free_intvec2d (long **v)
266 {
267 myfree(v[0]);
268 myfree(v);
269 }
270
271 void
setdoublevec1d(MYREAL ** v,MYREAL * w,long size)272 setdoublevec1d (MYREAL **v, MYREAL *w, long size)
273 {
274 doublevec1d (v, size);
275 memcpy (*v, w, sizeof (MYREAL) * size);
276 }
277
278 void
add_vector(MYREAL * result,MYREAL * v,long size)279 add_vector (MYREAL *result, MYREAL *v, long size)
280 {
281 long i;
282 for (i = 0; i < size; i++)
283 result[i] += v[i];
284 }
285
286
287 MYREAL
logfac(long n)288 logfac (long n)
289 {
290 /* log(n!) values were calculated with Mathematica
291 with a precision of 30 digits */
292 switch (n)
293 {
294 case 0:
295 return 0.;
296 case 1:
297 return 0.;
298 case 2:
299 return 0.693147180559945309417232121458;
300 case 3:
301 return 1.791759469228055000812477358381;
302 case 4:
303 return 3.1780538303479456196469416013;
304 case 5:
305 return 4.78749174278204599424770093452;
306 case 6:
307 return 6.5792512120101009950601782929;
308 case 7:
309 return 8.52516136106541430016553103635;
310 case 8:
311 return 10.60460290274525022841722740072;
312 case 9:
313 return 12.80182748008146961120771787457;
314 case 10:
315 return 15.10441257307551529522570932925;
316 case 11:
317 return 17.50230784587388583928765290722;
318 case 12:
319 return 19.98721449566188614951736238706;
320 case 13:
321 return 22.5521638531234228855708498286;
322 case 14:
323 return 25.1912211827386815000934346935;
324 case 15:
325 return 27.8992713838408915660894392637;
326 case 16:
327 return 30.6718601060806728037583677495;
328 case 17:
329 return 33.5050734501368888840079023674;
330 case 18:
331 return 36.3954452080330535762156249627;
332 case 19:
333 return 39.3398841871994940362246523946;
334 case 20:
335 return 42.3356164607534850296598759707;
336 case 21:
337 return 45.3801388984769080261604739511;
338 case 22:
339 return 48.4711813518352238796396496505;
340 case 23:
341 return 51.6066755677643735704464024823;
342 case 24:
343 return 54.7847293981123191900933440836;
344 case 25:
345 return 58.0036052229805199392948627501;
346 case 26:
347 return 61.2617017610020019847655823131;
348 case 27:
349 return 64.5575386270063310589513180238;
350 case 28:
351 return 67.8897431371815349828911350102;
352 case 29:
353 return 71.2570389671680090100744070426;
354 case 30:
355 return 74.6582363488301643854876437342;
356 default:
357 return LGAMMA (n + 1.);
358 }
359 }
360
361 ///
362 /// initialize the mean and standard deviation calculations
363 /// for the one-pass calculations that allows to throw away
364 /// the intermediate results
onepass_mean_std_start(MYREAL * mean,MYREAL * std,long * n)365 void onepass_mean_std_start(MYREAL *mean, MYREAL *std, long *n)
366 {
367 (*mean) = 0.0;
368 (*std) = 0.0;
369 (*n) = 0;
370 }
371
372 ///
373 /// calculates the mean and standard deviation using a one-pass filter
374 /// that allows to throw away the intermediate results
onepass_mean_std_calc(MYREAL * mean,MYREAL * std,long * n,MYREAL x)375 void onepass_mean_std_calc(MYREAL *mean, MYREAL *std, long *n, MYREAL x)
376 {
377
378 long nn = (*n);
379 MYREAL mm = (*mean);
380 MYREAL delta = x - mm;
381 mm = mm + delta / nn;
382 (*std) = (*std) + delta * (x - mm);
383 (*mean) = mm;
384 (*n) = nn;
385 }
386
387 ///
388 /// finishes one-pass mean and std calculations
onepass_mean_std_end(MYREAL * mean,MYREAL * std,long * n)389 void onepass_mean_std_end(MYREAL *mean, MYREAL *std, long *n)
390 {
391 (*std) = sqrt((*std) / ((*n) - 1));
392 }
393 #ifndef PRIORTEST
394 /*garbage collection and recycle*/
395 #ifdef MESS
find_in_list(long id)396 boolean find_in_list(long id)
397 {
398 long i;
399 for(i=0; i < check_collection_count; i++)
400 {
401 if(id == check_collection[i])
402 return TRUE;
403 }
404 return FALSE;
405 }
406
407 /// tags the lineage that was added the last time the tree was changed.
408 ///
traverse_check_id(node * theNode,long id)409 boolean traverse_check_id(node *theNode, long id)
410 {
411 boolean found = FALSE;
412 if(theNode->type=='r')
413 {
414 if(theNode->id == id || theNode->next->id == id || theNode->next->next->id == id)
415 return TRUE;
416 else
417 {
418 found = traverse_check_id(theNode->next->back, id);
419 return found;
420 }
421 }
422 if(theNode->type == 'i')
423 {
424 if(theNode->id == id || theNode->next->id == id || theNode->next->next->id == id)
425 return TRUE;
426 else
427 {
428 found = traverse_check_id(theNode->next->back, id);
429 if(found == TRUE)
430 return found;
431 found = traverse_check_id(theNode->next->next->back, id);
432 return found;
433 }
434 }
435 if(theNode->type == 'm')
436 {
437 if(theNode->id == id || theNode->next->id == id)
438 return TRUE;
439 else
440 {
441 found = traverse_check_id(theNode->next->back, id);
442 return found;
443 }
444 }
445 if(theNode->type == 't')
446 {
447 if(theNode->id == id)
448 return TRUE;
449 else
450 return FALSE;
451 }
452 error("node with no type");
453 return FALSE;
454 }
455
traverse_checker(node * root)456 void traverse_checker(node *root)
457 {
458 long i;
459 long x;
460 for(i=0;i<check_collection_count; i++)
461 {
462 x = check_collection[i];
463 if(!traverse_check_id(root,x))
464 printf ("Node id %li not found\n",x);
465 }
466 }
467 #endif /*MESS*/
468
start_node_collection(world_fmt * world)469 void start_node_collection(world_fmt *world)
470 {
471 long i;
472 #ifdef MESS
473 unique_id_global=0;
474 #endif
475 world->node_collection_allocated = 2 * HUNDRED + 2 * world->numpop * MIGRATION_LIMIT;
476 world->node_collection = (node **) calloc(world->node_collection_allocated,sizeof(node*));
477 #ifdef MESS
478 check_collection = (long *) calloc(world->node_collection_allocated, sizeof(long));
479 #endif
480 for(i = 0; i < world->node_collection_allocated; i++)
481 {
482 world->node_collection[i] = (node *) calloc(1, sizeof(node));
483 #ifdef MESS
484 (*world->node_collection[i]).id = unique_id_global++;
485 #endif
486 }
487 world->node_collection_count = world->node_collection_allocated;
488 }
489
stop_node_collection(world_fmt * world)490 void stop_node_collection(world_fmt *world)
491 {
492 long i;
493 for(i=0;i<world->node_collection_count;i++)
494 {
495 if(world->node_collection[i] != NULL)
496 {
497 free(world->node_collection[i]);
498 }
499 }
500 myfree(world->node_collection);
501 }
502
collect_nodelet(world_fmt * world,node * p)503 void collect_nodelet(world_fmt *world, node *p)
504 {
505 #ifdef MESS
506 long i;
507 #endif
508 long oldalloc;
509
510 if(world->node_collection_count < world->node_collection_allocated)
511 {
512 world->node_collection[world->node_collection_count] = p;
513 #ifdef MESS
514 for(i=0;i<check_collection_count;i++)
515 {
516 if(p->id == check_collection[i])
517 break;
518 }
519 if(i==check_collection_count)
520 {
521 warning("found node not accounted for: %li",p->id);
522 error("");
523 }
524 else
525 {
526 check_collection_count--;
527 check_collection[i] = check_collection[check_collection_count];
528 }
529 #endif
530 world->node_collection_count++;
531 // printf("d%li %c\n",p->id, p->type);
532 }
533 else
534 {
535 #ifdef DEBUG
536 warning("in node collection phase we should not need to allocate new space\n");
537 #endif
538 oldalloc = world->node_collection_allocated;
539 world->node_collection_allocated += 300;
540 world->node_collection = (node **) realloc(world->node_collection, world->node_collection_allocated * sizeof(node*));
541 memset(world->node_collection+oldalloc,0, 300 * sizeof(node*));
542 world->node_collection[world->node_collection_count] = p;
543 #ifdef MESS
544 for(i=0;i<check_collection_count;i++)
545 {
546 if(p->id == check_collection[i])
547 break;
548 }
549 if(i==check_collection_count)
550 {
551 warning("found node not accounted for: %li",p->id);
552 error("");
553 }
554 else
555 {
556 check_collection_count--;
557 check_collection[i] = check_collection[check_collection_count];
558 }
559 #endif
560 world->node_collection_count++;
561 #ifdef DEBUG
562 FPRINTF(stdout, "%i> Heat=%f: Collect_nodelet: recycler increased size to: %li/%li\n", myID, 1./world->heat, world->node_collection_count, world->node_collection_allocated);
563 #endif
564 }
565 //printf("\n");
566 }
567
568 ///
569 /// dispense memory for nodes from a list of avialable nodes, allocates more memory
570 /// no nodes are available in the list if nodes
dispense_nodelet(world_fmt * world)571 node * dispense_nodelet(world_fmt *world)
572 {
573 long i;
574 node *p = NULL;
575 if(world->node_collection_count > 0)
576 {
577 world->node_collection_count--;
578 p = world->node_collection[world->node_collection_count];
579 #ifdef MESS
580 check_collection[check_collection_count++] = p->id;
581 #endif
582 world->node_collection[world->node_collection_count] = NULL;
583 }
584 else
585 {
586 //all node memory is checked out, we
587 //increase available nodes in chunks of 300
588 world->node_collection_allocated += 300;
589 world->node_collection = (node **) myrealloc(world->node_collection,sizeof(node *)* world->node_collection_allocated);
590 // set new slots to NULL
591 memset(world->node_collection + world->node_collection_allocated - 300, 0,sizeof(node *)*300);
592 // we allocate memory of 300 nodes and start at the 0 element because all
593 // elements allocated earlier are checked out and the node_collection array
594 // should contain only NULLs
595 for(i = 0; i < 300; i++)
596 {
597 world->node_collection[i] = (node *) mycalloc(1, sizeof(node));
598 #ifdef MESS
599 world->node_collection[i]->id = unique_id_global++;
600 #endif
601 }
602 #ifdef MESS
603 check_collection = (long *) realloc(check_collection, world->node_collection_allocated * sizeof(long));
604 #endif
605 world->node_collection_count += 300;
606 world->node_collection_count--;
607 p = world->node_collection[world->node_collection_count]; // dispense
608 world->node_collection[world->node_collection_count] = NULL; // remove from list
609 #ifdef MESS
610 check_collection[check_collection_count++] = p->id;
611 #endif
612 }
613 return p;
614 }
615
616 ///
617 /// swaps the nodelet recycler between different temperatured chains
618 ///
swap_node_collection(world_fmt * tthis,world_fmt * that)619 void swap_node_collection(world_fmt * tthis, world_fmt * that)
620 {
621 node **temp;
622 //long tempalloc;
623 long tempcount;
624 long tempallocated;
625 temp = tthis->node_collection;
626 tempcount = tthis->node_collection_count;
627 tempallocated = tthis->node_collection_allocated;
628
629 tthis->node_collection = that->node_collection;
630 tthis->node_collection_count = that->node_collection_count;
631 tthis->node_collection_allocated = that->node_collection_allocated;
632
633 that->node_collection = temp;
634 that->node_collection_count = tempcount;
635 that->node_collection_allocated = tempallocated;
636 }
637
638 /*FILEMANIPULATION======================================================= */
639 void
init_files(world_fmt * world,data_fmt * data,option_fmt * options)640 init_files (world_fmt * world, data_fmt * data, option_fmt * options)
641 {
642 long *tempdb;
643 long maxfilenum = MYMAXFILENUM; //see definitions.h all files have a number and a filename;
644 #ifdef MPI
645 long i,j;
646 #ifdef PARALIO
647 long l;
648 boolean my_file_open_error;
649 char error_string[LINESIZE];
650 int length_of_error_string, error_class;
651 if(options->has_bayesmdimfile)
652 {
653 l = strlen(options->bayesmdimfilename);
654 sprintf(options->bayesmdimfilename+l,".%d",myID);
655 }
656 #endif
657 #endif
658
659 tempdb = (long *) mycalloc(maxfilenum,sizeof(long));
660
661 if (myID == MASTER)
662 {
663 #ifdef MPI
664 setup_filehandle_db((void *) stdout, world, options,data);
665 #endif
666 if (!options->readsum)
667 {
668 openfile (&data->infile, options->infilename, "r+", NULL);
669 if (options->usertree)
670 openfile (&data->utreefile, options->utreefilename, "r+", NULL);
671 if (options->weights)
672 openfile (&data->weightfile, options->weightfilename, "r+", NULL);
673 if (options->categs > 1)
674 openfile (&data->catfile, options->catfilename, "r+", NULL);
675 if (options->dist)
676 openfile (&data->distfile, options->distfilename, "r+", NULL);
677 if (options->writesum)
678 openfile (&data->sumfile, options->sumfilename, "w+", NULL);
679 }
680 else
681 {
682 if(!options->bayes_infer)
683 openfile (&data->sumfile, options->sumfilename, "r+", NULL);
684 }
685
686 openfile (&world->outfile, options->outfilename, "w+", NULL);
687 #ifdef MPI
688 setup_filehandle_db((void *) world->outfile, world, options,data);
689 #endif
690
691 if (options->treeprint > 0 && (!options->readsum))
692 {
693 openfile (&world->treefile, options->treefilename, "w+", NULL);
694 #ifdef MPI
695 setup_filehandle_db((void *) world->treefile, world, options,data);
696 #endif
697 }
698
699 if (options->mighist)
700 {
701 if(options->datatype != 'g')
702 {
703 openfile (&world->mighistfile, options->mighistfilename, "w+", NULL);
704 }
705 else
706 {
707 openfile (&world->mighistfile, options->mighistfilename, "r+", NULL);
708 }
709 #ifdef MPI
710 setup_filehandle_db((void *) world->mighistfile, world, options, data);
711 #endif
712 if(options->skyline)
713 {
714 if(options->datatype != 'g')
715 {
716 openfile (&world->skylinefile, options->skylinefilename, "w+", NULL);
717 }
718 else
719 {
720 openfile (&world->skylinefile, options->skylinefilename, "r+", NULL);
721 }
722 }
723 #ifdef MPI
724 setup_filehandle_db((void *) world->skylinefile, world, options, data);
725 #endif
726 }
727
728 if (options->writelog)
729 {
730 openfile (&options->logfile, options->logfilename, "w+", NULL);
731 #ifdef MPI
732 setup_filehandle_db((void *) options->logfile, world, options, data);
733 #endif
734 }
735
736 if (options->mixplot)
737 {
738 openfile (&options->mixfile, options->mixfilename, "w+", NULL);
739 #ifdef MPI
740 setup_filehandle_db((void *) options->mixfile, world, options, data);
741 #endif
742 }
743 if (options->aic)
744 {
745 openfile (&options->aicfile, options->aicfilename, "w+", NULL);
746 #ifdef MPI
747 setup_filehandle_db((void *) options->aicfile, world, options, data);
748 #endif
749 }
750 if (options->plot)
751 {
752 switch (options->plotmethod)
753 {
754 case PLOTALL:
755 openfile (&world->mathfile, options->mathfilename, "w+",
756 NULL);
757 #ifdef MPI
758 setup_filehandle_db((void *) world->mathfile, world, options, data);
759 #endif
760 break;
761 default: /*e.g. 0 this create just the plots in outfile */
762 break;
763 }
764 }
765 if (options->geo)
766 {
767 openfile (&data->geofile, options->geofilename, "r+", NULL);
768 }
769 if (options->has_datefile)
770 {
771 openfile (&data->datefile, options->datefilename, "r+", NULL);
772 }
773 #ifdef UEP
774
775 if (options->uep)
776 openfile (&data->uepfile, options->uepfilename, "r+", NULL);
777 #endif
778 if (options->bayes_infer)
779 {
780 if(options->has_bayesfile)
781 openfile (&world->bayesfile, options->bayesfilename, "w+", NULL);
782
783 if(options->has_bayesmdimfile)
784 {
785 if(options->datatype != 'g')
786 {
787 #ifdef MPI
788 #ifndef PARALIO
789 #ifdef ZNZ
790 znzopenfile (&world->bayesmdimfile, options->bayesmdimfilename, "w", options->use_compressed);
791 #else
792 openfile (&world->bayesmdimfile, options->bayesmdimfilename, "w+", NULL);
793 #endif /*znz*/
794 #endif /*not paralio*/
795 #else
796 // not MPI
797 #ifdef ZNZ
798 znzopenfile (&world->bayesmdimfile, options->bayesmdimfilename, "w", options->use_compressed);
799 #else
800 openfile (&world->bayesmdimfile, options->bayesmdimfilename, "w+", NULL);
801 #endif
802 #endif /*mpi*/
803 }
804 else
805 {
806 #ifdef ZNZ
807 znzopenfile (&world->bayesmdimfile, options->bayesmdimfilename, "r", options->use_compressed);
808 #else
809 openfile(&world->bayesmdimfile, options->bayesmdimfilename,"r+",NULL);
810 #endif
811 }
812 }
813 #ifdef MPI
814 if(options->has_bayesfile)
815 setup_filehandle_db((void *) world->bayesfile, world, options,data);
816 if(options->has_bayesmdimfile)
817 {
818 #ifdef PARALIO
819 fprintf(stdout,"%i> before open %s\n",myID,options->bayesmdimfilename);
820 my_file_open_error = MPI_File_open(MPI_COMM_SELF, options->bayesmdimfilename,
821 MPI_MODE_CREATE | MPI_MODE_WRONLY ,
822 MPI_INFO_NULL, &world->mpi_bayesmdimfile);
823 if (my_file_open_error != MPI_SUCCESS)
824 {
825 MPI_Error_class(my_file_open_error, &error_class);
826 MPI_Error_string(error_class, error_string, &length_of_error_string);
827 printf("%i> %s %s\n", myID, error_string, options->bayesmdimfilename);
828 MPI_Error_string(my_file_open_error, error_string,
829 &length_of_error_string);
830 printf("%i> %s\n", myID, error_string);
831 my_file_open_error = TRUE;
832 }
833 MPI_File_set_size(world->mpi_bayesmdimfile, 0);
834 #else
835 setup_filehandle_db((void *) world->bayesmdimfile, world, options,data);
836 #endif
837 } //has bayesmdimfile
838 #endif
839 }
840 #ifdef MPI
841 tempdb[0]=filenum;
842 for(i=0, j=1;i<filenum;i++)
843 {
844 tempdb[j++] = (long) filedb[i].file;
845 tempdb[j++] = (long) filedb[i].handle;
846 }
847 MYMPIBCAST (tempdb, maxfilenum, MPI_LONG, MASTER, comm_world);
848 #endif
849 }
850 else
851 {
852 // all non myID==MASTER
853 #ifdef MPI
854 #ifdef PARALIO
855 if(options->has_bayesmdimfile)
856 {
857 MPI_Info info;
858 fprintf(stdout,"%i> before open %s\n",myID,options->bayesmdimfilename);
859 my_file_open_error = MPI_File_open(MPI_COMM_SELF, options->bayesmdimfilename,
860 MPI_MODE_CREATE | MPI_MODE_WRONLY ,
861 MPI_INFO_NULL, &world->mpi_bayesmdimfile);
862 if (my_file_open_error != MPI_SUCCESS)
863 {
864 MPI_Error_class(my_file_open_error, &error_class);
865 MPI_Error_string(error_class, error_string, &length_of_error_string);
866 printf("%i> @%s %s\n", myID, error_string, options->bayesmdimfilename);
867 MPI_Error_string(my_file_open_error, error_string,
868 &length_of_error_string);
869 printf("%i> @%s\n", myID, error_string);
870 my_file_open_error = TRUE;
871 }
872 MPI_File_set_size(world->mpi_bayesmdimfile, 0);
873 }
874 #endif
875 MYMPIBCAST (tempdb, maxfilenum, MPI_LONG, MASTER, comm_world);
876 filenum = tempdb[0];
877 for(i=0, j=1;i<filenum;i++)
878 {
879 filedb[i].file = (FILE *) tempdb[j++];
880 filedb[i].handle = tempdb[j++];
881 }
882 //assumess same order as in the master
883 j=1; //zero is stdout
884 world->outfile = filedb[j++].file;
885 if (options->treeprint > 0 && (!options->readsum))
886 world->treefile = filedb[j++].file;
887 if (options->mighist)
888 world->mighistfile = filedb[j++].file;
889 if (options->skyline)
890 world->skylinefile = filedb[j++].file;
891 if (options->writelog)
892 options->logfile = filedb[j++].file;
893 if (options->aic)
894 options->aicfile = filedb[j++].file;
895 if (options->plot && options->plotmethod == PLOTALL)
896 world->mathfile = filedb[j++].file;
897 if (options->bayes_infer)
898 {
899 if(options->has_bayesfile)
900 world->bayesfile = filedb[j++].file;
901 // TEST: experiment with parallele I/O
902 if(options->has_bayesmdimfile)
903 {
904 #ifdef PARALIO
905 //has to be dealt with before the broadcast MPI_File_open(comm_world,options->bayesmdimfilename,MPI_MODE_CREATE,MPI_INFO_NULL, &world->mpi_bayesmdimfile);
906 #else
907 #ifdef ZNZ
908 world->bayesmdimfile = (znzFile) filedb[j++].file;
909 #else
910 world->bayesmdimfile = filedb[j++].file;
911 #endif
912 #endif
913 }
914 }
915 #endif /*MPI*/
916 }
917 myfree(tempdb);
918 }
919
920 void
exit_files(world_fmt * world,data_fmt * data,option_fmt * options)921 exit_files (world_fmt * world, data_fmt * data, option_fmt * options)
922 {
923 if(myID==MASTER)
924 {
925 if (!options->readsum)
926 {
927 FClose (data->infile);
928 // if (options->treeprint > 0)
929 // FClose (world->treefile);
930 if(options->usertree)
931 FClose(data->utreefile);
932 if (options->weights)
933 FClose (data->weightfile);
934 if (options->categs > 1)
935 FClose (data->catfile);
936 if (options->dist)
937 FClose (data->distfile);
938 }
939
940 if (options->writesum || options->readsum)
941 FClose (data->sumfile);
942
943 FClose (world->outfile);
944
945 if (options->writelog)
946 FClose (options->logfile);
947 if (options->mighist)
948 FClose (world->mighistfile);
949 if (options->geo)
950 FClose (data->geofile);
951 if (options->has_datefile)
952 FClose (data->datefile);
953 if (options->aic)
954 FClose (options->aicfile);
955 if (options->plot && options->plotmethod == PLOTALL)
956 FClose (world->mathfile);
957 #ifdef UEP
958
959 if (options->uep)
960 FClose (data->uepfile);
961 #endif
962 if (options->bayes_infer)
963 {
964 if(options->has_bayesfile)
965 FClose (world->bayesfile);
966 if(options->has_bayesmdimfile)
967 #ifdef ZNZ
968 znzClose (world->bayesmdimfile);
969 #else
970 FClose (world->bayesmdimfile);
971 #endif
972 }
973 }
974 }
975
976 /* string manipulation ================================== */
977 /* Converts any character from to character to in string text */
978 void
translate(char * text,char from,char to)979 translate (char *text, char from, char to)
980 {
981 int i, j, gap = 0;
982 while (text[gap] == from)
983 gap++;
984 for (i = gap, j = 0; text[i] != '\0'; i++)
985 {
986 if (text[i] != from)
987 {
988 text[j++] = text[i];
989 }
990 else
991 {
992 if (text[i - 1] != from)
993 {
994 text[j++] = to;
995 }
996 }
997 }
998 text[j] = '\0';
999 }
1000 /// unpad from the end of the word until there is none of the
1001 /// specified character left
1002 void
unpad(char * text,char removechars[])1003 unpad (char *text, char removechars[])
1004 {
1005 int gap = strlen(text)-1;
1006 while (gap > 0 && strchr(removechars,text[gap]))
1007 gap--;
1008 text[gap+1] = '\0';
1009 }
1010
1011 ///
1012 /// get next word from list with list of delimiters
1013 /// uses strsep
1014 void
get_next_word(char ** instring,char * delimiters,char ** nextword)1015 get_next_word(char **instring, char *delimiters, char **nextword)
1016 {
1017 *nextword = strsep(instring, delimiters);
1018 while(*nextword != NULL && strchr(delimiters,(*nextword)[0]))
1019 {
1020 // if((*instring)==NULL)
1021 // {
1022 // (*nextword) = NULL;
1023 // return;
1024 // }
1025 *nextword = strsep(instring,delimiters);
1026 }
1027 }
1028 /*===============================================
1029 count words in a string delimited by delimiter
1030 */
1031 long
count_words(char * text)1032 count_words (char *text)
1033 {
1034 long counts = 0;
1035 char *pt = text;
1036 while (isspace (*pt) && *pt != '\0')
1037 pt++;
1038 while (*pt != '\0')
1039 {
1040 while (!isspace (*pt) && *pt != '\0')
1041 pt++;
1042 while (isspace (*pt) && *pt != '\0')
1043 pt++;
1044 counts++;
1045 }
1046 return counts;
1047 }
1048
1049
1050 /*===============================================
1051 timer utility
1052
1053 ts = "%c" -> time + full date (see man strftime)
1054 = "%H:%M:%S" -> time hours:minutes:seconds */
1055
1056 void
get_time(char * nowstr,char ts[])1057 get_time (char *nowstr, char ts[])
1058 {
1059 #ifdef NOTIME_FUNC
1060 switch (strlen (ts))
1061 {
1062 case 2:
1063 strcpy (nowstr, " ");
1064 break;
1065 case 3:
1066 strcpy (nowstr, " ");
1067 break;
1068 case 8:
1069 strcpy (nowstr, " ");
1070 break;
1071 default:
1072 strcpy (nowstr, " ");
1073 break;
1074 }
1075 #else
1076 time_t nowbin;
1077 struct tm *nowstruct;
1078 if (time (&nowbin) != (time_t) - 1)
1079 {
1080 nowstruct = localtime (&nowbin);
1081 strftime (nowstr, LINESIZE, ts, nowstruct);
1082 }
1083 #endif
1084 }
1085
1086 /*===============================================
1087 printer utility
1088 */
1089 void
print_llike(MYREAL llike,char * strllike)1090 print_llike (MYREAL llike, char *strllike)
1091 {
1092 if (fabs (llike) > 10e20)
1093 {
1094 sprintf (strllike, "%cInfinity ", llike < 0 ? '-' : ' ');
1095 }
1096 else
1097 sprintf (strllike, "%-10.5f", llike);
1098 }
1099
1100 ///
1101 /// remove traling whitespace from a *string
remove_trailing_blanks(char ** filename)1102 void remove_trailing_blanks(char **filename)
1103 {
1104 int pos;
1105 pos = (long) strlen(*filename);
1106 while (pos>0 && isspace( (*filename)[pos-1]))
1107 pos--;
1108 (*filename)[pos] = '\0';
1109 }
1110
1111 ///
1112 /// read a single word from a string up to the position of
1113 /// the breakchar and return that position
read_word_delim(char * input,char * word,char * delim)1114 long read_word_delim(char *input, char *word, char * delim)
1115 {
1116 long i = 1;
1117 int ch;
1118 long count = 0;
1119 ch = input[0];
1120 if(ch=='\0')
1121 return 0;
1122 while (isspace (ch))
1123 ch = input[i++];
1124 while(strchr(delim, ch)==NULL && ch!=EOF)
1125 {
1126 word[count++] = ch;
1127 ch = input[i++];
1128 }
1129 word[count]='\0';
1130 return i;
1131 }
1132 ///
1133 /// read a single word from a file
1134 /// words are delimited by withespace
read_word(FILE * infile,char * word)1135 int read_word(FILE *infile, char *word)
1136 {
1137 int ch;
1138 long count = 0;
1139 ch = getc(infile);
1140 while (isspace (ch))
1141 ch = getc(infile);
1142 while(strchr(" \t\n\r", ch)==NULL && ch!=EOF)
1143 {
1144 word[count++] = ch;
1145 ch = getc(infile);
1146 }
1147 word[count]='\0';
1148 if(ch==EOF)
1149 return EOF;
1150 else
1151 {
1152 ungetc(ch,infile);
1153 return 0;
1154 }
1155 }
1156
1157 //========================================
1158 // prepends a string to a file
1159 // a blank will be inserted between the string and the
1160 // text in the file
unread_word(FILE * infile,char * word)1161 void unread_word(FILE *infile, char *word)
1162 {
1163 char c = ' ';
1164 long count = (long) strlen(word);
1165 ungetc(c, infile);
1166 while(count > 0)
1167 {
1168 count--;
1169 ungetc(word[count], infile);
1170 }
1171 }
1172
1173 #ifdef STRANGEDEBUG
1174 extern int errno;
1175 #endif
1176
1177 #ifdef ZNZ
1178 /// opens a zip file for reading or writing
1179 void
znzopenfile(znzFile * fp,char * filename,char * mode,int use_compressed)1180 znzopenfile (znzFile * fp, char *filename, char *mode, int use_compressed)
1181 {
1182 //int trials = 0;
1183 znzFile of;
1184 char *file;
1185 char *p;
1186 file = (char *) mycalloc(LINESIZE,sizeof(char));
1187 if ((p = strpbrk (filename, CRLF)) != NULL)
1188 *p = '\0';
1189 remove_trailing_blanks(&filename);
1190 sprintf(file, "%-.*s", (int) LINESIZE, filename);
1191 of = znzopen (file, mode,use_compressed);
1192 if (znz_isnull(of))
1193 error("Could not open zipped bayesallfile");
1194 *fp = of;
1195 myfree(file);
1196 }
1197 #endif
1198
1199
1200 /// opens a file for reading or writing
1201 void
openfile(FILE ** fp,char * filename,char * mode,char * perm)1202 openfile (FILE ** fp, char *filename, char *mode, char *perm)
1203 {
1204 int trials = 0;
1205 FILE *of = NULL;
1206 char *file;
1207 char *p;
1208 #ifdef CAUTIOUS
1209 struct stat sb;
1210 #endif
1211 file = (char *) mycalloc(LINESIZE,sizeof(char));
1212 if ((p = strpbrk (filename, CRLF)) != NULL)
1213 *p = '\0';
1214 remove_trailing_blanks(&filename);
1215 sprintf(file, "%-.*s",(int) LINESIZE, filename);
1216 #ifdef CAUTIOUS
1217 if(cautious)
1218 {
1219 // sb = (stat *) mycalloc(1,sizeof(stat));
1220 if(strchr(mode,'w') && (0 == stat(file, &sb)))
1221 {
1222 file[0] = '\0';
1223 while (file[0] == '\0' && trials++ < 10)
1224 {
1225 printf ("File %s exists, do you want to overwrite it? [YES/No]\n===>",filename);
1226 FGETS (file, LINESIZE, stdin);
1227 if(file[0]=='Y' || file[0] == 'y')
1228 {
1229 file[0] = '\0';
1230 while (file[0] == '\0' && trials++ < 10)
1231 {
1232 printf ("Please enter a new filename>");
1233 FGETS (file, LINESIZE, stdin);
1234 }
1235
1236 }
1237 }
1238 }
1239 // myfree(sb);
1240 }
1241 #endif
1242 while (trials++ < 10)
1243 {
1244 of = fopen (file, mode);
1245 #ifdef STRANGEDEBUG
1246 fprintf(stdout,"original filename -->%s<--\n",filename);
1247 fprintf(stdout,"used filename -->%s<--\n",file);
1248 #endif
1249 if (of!=NULL)
1250 break;
1251 else
1252 {
1253 // extension = strrchr(file,'.');
1254 // if(extension!=NULL && !strncmp(extension,".txt",4))
1255 #ifdef STRANGEDEBUG
1256 fprintf(stdout,"mode -->%s<--\n", mode);
1257 fprintf(stdout,"errorcode -->%i<--\n",errno);
1258 #endif
1259 switch (*mode)
1260 {
1261 case 'r':
1262 #ifdef MPI
1263
1264 printf("%i> ",myID);
1265 #endif
1266
1267 printf ("Cannot read from file \"%s\"\n", file);
1268 file[0] = '\0';
1269 while (file[0] == '\0' && trials++ < 10)
1270 {
1271 #ifdef MPI
1272 printf("%i> ",myID);
1273 #endif
1274
1275 printf ("Please enter a new filename for reading>");
1276 FGETS (file, LINESIZE, stdin);
1277 }
1278 break;
1279 case 'w':
1280 #ifdef MPI
1281
1282 printf("%i> ",myID);
1283 #endif
1284
1285 printf ("Cannot write to file %s\n", file);
1286 file[0] = '\0';
1287 while (file[0] == '\0' && trials++ < 10)
1288 {
1289 #ifdef MPI
1290 printf("%i> ",myID);
1291 #endif
1292
1293 printf ("Please enter a new filename for writing>");
1294 FGETS (file, LINESIZE, stdin);
1295 }
1296 break;
1297 }
1298 }
1299 }
1300 if (trials >= 10)
1301 {
1302 #ifdef MPI
1303 printf("%i> ",myID);
1304 #endif
1305
1306 printf ("You cannot find your file either, so I stop\n\n");
1307 exit (0);
1308 }
1309 *fp = of;
1310 if (perm != NULL)
1311 strcpy (perm, file);
1312 strcpy (filename, file);
1313 myfree(file);
1314 }
1315
get_filename(char ** store,char * value)1316 void get_filename(char **store, char * value)
1317 {
1318 long ls;
1319 long lv;
1320 remove_trailing_blanks(&value);
1321 ls = strlen(*store)+1;
1322 lv = strlen(value)+1;
1323 if (ls < lv)
1324 {
1325 (*store) = (char *) realloc(*store,sizeof(char)*lv);
1326 memset(*store,0,sizeof(char)*lv);
1327 }
1328 strcpy (*store, value);
1329 }
1330
1331
1332 long
read_savesum(world_fmt * world,option_fmt * options,data_fmt * data)1333 read_savesum (world_fmt * world, option_fmt * options, data_fmt * data)
1334 {
1335 long Gmax = 0;
1336 FILE *sumfile = data->sumfile;
1337 timearchive_fmt **ta;
1338 long nrep = 0;
1339 long l, i, j, r;
1340 char input[1024];
1341 long hits;
1342 long tmp;
1343 long z, zz, jj;
1344 MYREAL tmpm;
1345 boolean newstyle = TRUE;
1346 int retval;
1347 #ifdef MPI
1348
1349 long listofthings[6];
1350 #endif
1351
1352 FGETS (input, 1024, sumfile);
1353 if (!strncmp ("# begin genealogy-summary file of migrate", input, 41))
1354 {
1355 FGETS (input, 1024, sumfile);
1356 if (input[0] == '#')
1357 {
1358 newstyle = TRUE;
1359 FGETS (input, 1024, sumfile);
1360 if(input[0]=='M') // added custom migration model 11-25-02
1361 {
1362 options->custm = (char *) myrealloc(options->custm, sizeof(char) * (strlen(input)+1));
1363 strcpy(options->custm,input+6);
1364 FGETS (input, 1024, sumfile);
1365 options->custm2 = (char *) myrealloc(options->custm2, sizeof(char) * (strlen(input)+1));
1366 strcpy(options->custm2,input+6);
1367 FGETS (input, 1024, sumfile);
1368 }
1369 }
1370 else
1371 newstyle = FALSE;
1372 hits =
1373 sscanf (input, "%li %li %li %li %li", &world->loci, &world->numpop,
1374 &world->numpop2, &tmp, &options->replicatenum);
1375 if (hits != 5)
1376 {
1377 nrep = 1;
1378 //xcode hits =
1379 sscanf (input, "%li %li %li", &world->loci, &world->numpop,
1380 &world->numpop2);
1381 options->replicate = FALSE;
1382 options->replicatenum = 0;
1383 }
1384 else
1385 {
1386 nrep = options->replicatenum;
1387 if (nrep == 0)
1388 nrep = 1;
1389 options->replicate = (boolean) tmp;
1390 }
1391 options->muloci = data->loci = world->loci;
1392 data->skiploci =
1393 (boolean *) myrealloc (data->skiploci,
1394 sizeof (boolean) * (data->loci + 1));
1395 memset (data->skiploci, 0, sizeof (boolean) * (data->loci + 1));
1396 read_geofile (data, options, world->numpop);
1397 data->numpop = world->numpop;
1398 options->newpops_numalloc = world->numpop;
1399 options->newpops = (long*) mycalloc(options->newpops_numalloc, sizeof(long));
1400 calculate_newpop_numpop(options,data);
1401 set_plot (options);
1402 world->cold=TRUE;
1403 init_world (world, data, options);
1404 ta = world->atl;
1405 for (l = 0; l < world->loci; l++)
1406 for (r = 0; r < nrep; r++)
1407 {
1408 FGETS (input, 1024, sumfile);
1409 retval = fscanf (sumfile, "%li %li %li\n", &ta[r][l].T, &ta[r][l].numpop,
1410 &ta[r][l].sumtips);
1411 Gmax = (ta[r][l].T > Gmax) ? ta[r][l].T : Gmax;
1412 ta[r][l].allocT = 0;
1413 increase_timearchive (world, l, ta[r][l].T, world->numpop, r);
1414 #ifdef USE_MYREAL_FLOAT
1415 retval = fscanf (sumfile, "%g %g %g\n", &ta[r][l].param_like,
1416 &ta[r][l].thb, &ta[r][l].alpha);
1417 #else
1418 retval = fscanf (sumfile, "%lg %lg %lg\n", &ta[r][l].param_like,
1419 &ta[r][l].thb, &ta[r][l].alpha);
1420 #endif
1421 world->chainlikes[l][r] = ta[r][l].param_like;
1422
1423 for (i = 0; i < world->atl[r][l].T; i++)
1424 {
1425 #ifdef USE_MYREAL_FLOAT
1426 retval = fscanf (sumfile, "%li %f\n", &ta[r][l].tl[i].copies, &ta[r][l].tl[i].lcopies);
1427 #else
1428 retval = fscanf (sumfile, "%li %lf\n", &ta[r][l].tl[i].copies, &ta[r][l].tl[i].lcopies);
1429 #endif
1430 for (j = 0; j < world->numpop; j++)
1431 {
1432 #ifdef USE_MYREAL_FLOAT
1433 retval = fscanf (sumfile, "%f %f %f\n", &ta[r][l].tl[i].km[j],
1434 &ta[r][l].tl[i].kt[j], &ta[r][l].tl[i].p[j]);
1435 #else
1436 retval = fscanf (sumfile, "%lf %lf %lf\n", &ta[r][l].tl[i].km[j],
1437 &ta[r][l].tl[i].kt[j], &ta[r][l].tl[i].p[j]);
1438 #endif
1439 }
1440 if (newstyle)
1441 {
1442 for (j = 0; j < world->numpop2; j++)
1443 {
1444 #ifdef USE_MYREAL_FLOAT
1445 retval = fscanf (sumfile, "%f ", &ta[r][l].tl[i].mindex[j]);
1446 #else
1447 retval = fscanf (sumfile, "%lf ", &ta[r][l].tl[i].mindex[j]);
1448 #endif
1449 }
1450 }
1451 else
1452 {
1453 z = 0;
1454 zz = world->numpop;
1455 for (j = 0; j < world->numpop; j++)
1456 {
1457 for (jj = 0; jj < world->numpop; jj++)
1458 {
1459 #ifdef USE_MYREAL_FLOAT
1460 retval = fscanf (sumfile, "%f ", &tmpm);
1461 #else
1462 retval = fscanf (sumfile, "%lf ", &tmpm);
1463 #endif
1464 if (j == jj)
1465 ta[r][l].tl[i].mindex[z++] = tmpm;
1466 else
1467 ta[r][l].tl[i].mindex[zz++] = tmpm;
1468 }
1469 }
1470 }
1471 }
1472 // ta[r][l].tl[0].lcopies = ta[r][l].tl[0].copies > 1 ? ln_copies(ta[r][l].tl[0].copies)-1. : 0.0;
1473 for (i = 0; i < world->numpop2; i++)
1474 {
1475 #ifdef USE_MYREAL_FLOAT
1476 retval = fscanf (sumfile, "%g %g\n", &ta[r][l].param[i],
1477 &ta[r][l].param0[i]);
1478 #else
1479 retval = fscanf (sumfile, "%lg %lg\n", &ta[r][l].param[i],
1480 &ta[r][l].param0[i]);
1481 #endif
1482 }
1483 log_param0 (ta[r][l].param0, ta[r][l].lparam0, world->numpop2);
1484 //for(i=0;i<ta[l].T;i++)
1485 // {
1486 // fscanf(sumfile,"%lg ",&ta[l].likelihood[i]);
1487 // }
1488 // FGETS(input,1024,sumfile);
1489 #ifdef USE_MYREAL_FLOAT
1490 retval = fscanf (sumfile, "%li %g\n", &ta[r][l].trials, &ta[r][l].normd);
1491 #else
1492 retval = fscanf (sumfile, "%li %lg\n", &ta[r][l].trials, &ta[r][l].normd);
1493 #endif
1494 }
1495 }
1496 else
1497 {
1498 error ("This is not a genealogy-summary file for MIGRATE\n");
1499 }
1500 #ifdef MPI
1501 listofthings[0] = world->loci;
1502 listofthings[1] = world->numpop;
1503 listofthings[2] = world->numpop2;
1504 listofthings[3] = (long) options->replicate;
1505 listofthings[4] = options->replicatenum;
1506 listofthings[5] = Gmax;
1507 //printf("listofthings sent");
1508 fflush(stdout);
1509 MYMPIBCAST (listofthings,5, MPI_LONG, MASTER, comm_world);
1510 #endif
1511
1512 return Gmax;
1513 }
1514
1515 void
write_savesum(world_fmt * world)1516 write_savesum (world_fmt * world)
1517 {
1518 long r, repmax;
1519 FILE *sumfile = world->data->sumfile;
1520 timearchive_fmt **ta = world->atl;
1521
1522 long l, i, j;
1523 if (world->options->replicate)
1524 {
1525 if (world->options->replicatenum == 0)
1526 repmax = world->options->lchains;
1527 else
1528 repmax = world->options->replicatenum;
1529 }
1530 else
1531 repmax = 1;
1532 fprintf (sumfile,
1533 "# begin genealogy-summary file of migrate %s ------\n#\n",
1534 MIGRATEVERSION);
1535 fprintf (sumfile,"Model=%s\n",world->options->custm);
1536 fprintf (sumfile,"Mode2=%s\n",world->options->custm2);
1537 fprintf (sumfile, "%li %li %li %li %li\n", world->loci, world->numpop,
1538 world->numpop2, (long) world->options->replicate, repmax);
1539 for (l = 0; l < world->loci; l++)
1540 {
1541 for (r = 0; r < repmax; r++)
1542 {
1543 fprintf (sumfile,
1544 "%li %li ####### locus %li, replicate %li ################\n",
1545 l, r, l, r);
1546 fprintf(sumfile,"%li %li %li\n",ta[r][l].T, ta[r][l].numpop, ta[r][l].sumtips);
1547
1548 fprintf (sumfile, "%20.20g %20.20g %20.20g\n", ta[r][l].param_like,
1549 ta[r][l].thb, ta[r][l].alpha);
1550 for (i = 0; i < ta[r][l].T; i++)
1551 {
1552 fprintf (sumfile, "%li %f\n", ta[r][l].tl[i].copies,ta[r][l].tl[i].lcopies);
1553 for (j = 0; j < world->numpop; j++)
1554 {
1555 fprintf (sumfile, "%20.20f %20.20f %f\n",
1556 ta[r][l].tl[i].km[j], ta[r][l].tl[i].kt[j],
1557 ta[r][l].tl[i].p[j]);
1558 }
1559 for (j = 0; j < world->numpop2; j++)
1560 {
1561 fprintf (sumfile, "%f ", ta[r][l].tl[i].mindex[j]);
1562 //debug printf (sumfile, "%li ", ta[l].tl[i].l[j]);
1563 }
1564 fprintf (sumfile, "\n");
1565 }
1566 for (i = 0; i < world->numpop2; i++)
1567 {
1568 fprintf (sumfile, "%20.20e %20.20e\n", ta[r][l].param[i],
1569 ta[r][l].param0[i]);
1570 }
1571 // for(i=0;i<ta[l].T;i++)
1572 // {
1573 // printf(sumfile,"%20.20e ",ta[l].likelihood[i]);
1574 // }
1575 // printf(sumfile,"\n");
1576 fprintf (sumfile, "%li %20.20e\n", ta[r][l].trials, ta[r][l].normd);
1577 }
1578 }
1579 fprintf (sumfile, "# end genealogy-summary file of migrate %s ------\n",
1580 MIGRATEVERSION);
1581 fflush (sumfile);
1582 }
1583
1584
1585 /*=======================================================*/
1586
1587
1588
1589
1590 /*--------------------------------
1591 creates the length value in a node
1592 */
1593 MYREAL
lengthof(node * p)1594 lengthof (node * p)
1595 {
1596 if (p->type == 'm')
1597 error ("A migration node was feed into lengthof");
1598 return fabs (p->tyme - crawlback (p)->tyme);
1599 } /* length */
1600
1601
1602 /*------------------------------------------------
1603 Find the next non-migration node starting
1604 with the theNode, returns to backnode which is not
1605 a migration, does NOT return always a top-node!
1606 */
1607 MYINLINE node *
crawlback(const node * theNode)1608 crawlback (const node * theNode)
1609 {
1610 node *tmp = theNode->back;
1611
1612 while (tmp->type == 'm')
1613 {
1614 tmp = tmp->next->back;
1615 }
1616 return tmp;
1617 }
1618
1619 /*--------------------------------------------
1620 returns the last migration node in a branch or
1621 the node if there is no migration node
1622
1623 node *crawl(node * theNode)
1624 {
1625 node *otmp, *tmp = theNode->back;
1626
1627 otmp = theNode;
1628 if (tmp == NULL)
1629 return otmp;
1630 while (tmp->type == 'm') {
1631 otmp = tmp->next;
1632 tmp = tmp->next->back;
1633 if (tmp == NULL)
1634 return otmp;
1635 }
1636 return otmp;
1637 }
1638 */
1639
1640
1641 MYINLINE node *
showtop(node * theNode)1642 showtop (node * theNode)
1643 {
1644 if (theNode == NULL)
1645 return NULL;
1646 else
1647 {
1648 if (theNode->top)
1649 {
1650 return theNode;
1651 }
1652 else
1653 {
1654 if (theNode->next->top)
1655 {
1656 return theNode->next;
1657 }
1658 else
1659 {
1660 return theNode->next->next;
1661 }
1662 }
1663 }
1664
1665 }
1666
1667 /* adjust the time in a node to time */
1668 void
adjust_time(node * theNode,MYREAL tyme)1669 adjust_time (node * theNode, MYREAL tyme)
1670 {
1671 switch (theNode->type)
1672 {
1673 case 'm':
1674 theNode->tyme = theNode->next->tyme = tyme;
1675 set_dirty(theNode);
1676 break;
1677 case 'i':
1678 theNode->tyme = theNode->next->tyme = theNode->next->next->tyme = tyme;
1679 set_dirty(theNode);
1680 break;
1681 case 'r':
1682 break;
1683 case 't':
1684 break;
1685 default:
1686 error ("Wrong node type");
1687 break;
1688 }
1689 }
1690 /* adjust the time in any node to time */
1691 void
adjust_time_all(node * theNode,MYREAL tyme)1692 adjust_time_all (node * theNode, MYREAL tyme)
1693 {
1694 switch (theNode->type)
1695 {
1696 case 'm':
1697 theNode->tyme = theNode->next->tyme = tyme;
1698 set_dirty(theNode);
1699 break;
1700 case 'i':
1701 theNode->tyme = theNode->next->tyme = theNode->next->next->tyme = tyme;
1702 set_dirty(theNode);
1703 break;
1704 case 'r':
1705 theNode->tyme = theNode->next->tyme = theNode->next->next->tyme = tyme;
1706 set_dirty(theNode);
1707 break;
1708 case 't':
1709 theNode->tyme = tyme;
1710 break;
1711 default:
1712 error ("Wrong node type");
1713 break;
1714 }
1715 }
1716
1717 void
insert_migr_node(world_fmt * world,node * up,node * down,migr_table_fmt * migr_table,long * migr_table_counter)1718 insert_migr_node (world_fmt * world, node * up, node * down,
1719 migr_table_fmt * migr_table, long *migr_table_counter)
1720 {
1721 long i;
1722 //, panic;
1723 //node *tmp, *tmp2, *oldNode, *oldNode2, *theNode;
1724 node *theNode;
1725 if (!up->top)
1726 error ("up has to be a top-node");
1727 //xcode theNode = showtop (up)->back;
1728 if (*migr_table_counter > 0 && up->tyme > migr_table[0].time)
1729 {
1730 printf("%i> up->type=%c tyme=%f showtop(up)time=%f mig[0]tabletyme=%f\n", myID, up->type, up->tyme, showtop(up)->tyme, migr_table[0].time);
1731 error
1732 ("insert_migr_node: the first migration node has a wrong time for up");
1733 }
1734 if (migr_table[(*migr_table_counter) - 1].from != down->actualpop)
1735 {
1736 error ("this should never happen -> wrong choice of nodes\n");
1737 (*migr_table_counter)--;
1738 }
1739 if (((*migr_table_counter) > 0)
1740 && (migr_table[(*migr_table_counter) - 1].from != down->actualpop))
1741 error ("problem catched in inser_migr_table");
1742 theNode=up;
1743 for (i = 0; i < (*migr_table_counter); i++)
1744 {
1745 theNode = add_migration(world, theNode,
1746 migr_table[i].from,
1747 migr_table[i].to,
1748 migr_table[i].time - theNode->tyme);
1749 }
1750 if(down->tyme < migr_table[i-1].time)
1751 {
1752 error("Problem with migration addition in coalesce1p()");
1753 }
1754 //printf("%li migration nodelets added\n",(*migr_table_counter) * 2);
1755 down->back = theNode;
1756 theNode->back = down;
1757 /* the stuff below will be excised once I am complete sure that
1758 this recycler above really works -- it seems currently dec 05 07
1759 tmp = (node *) mycalloc (1, sizeof (node));
1760 tmp2 = (node *) mycalloc (1, sizeof (node));
1761 oldNode = up;
1762 theNode = up->back;
1763 panic = 0;
1764 //printf("+");
1765 while (theNode->tyme < migr_table[i].time ****&& panic++ < 10000****)
1766 {
1767 if (theNode->type == 't' )
1768 {
1769 oldNode = theNode;
1770 theNode = theNode->back;
1771 }
1772 else
1773 {
1774 oldNode = theNode->next;
1775 theNode = theNode->next->back;
1776 }
1777 }
1778 tmp->back = oldNode;
1779 oldNode->back = tmp;
1780 tmp->number = -999;
1781 tmp->nayme = NULL;
1782 tmp->tip = FALSE;
1783 tmp->top = FALSE;
1784 tmp->dirty = TRUE;
1785 tmp->id = world->unique_id++;
1786 tmp->tyme = migr_table[i].time;
1787 tmp->type = 'm';
1788 tmp->actualpop = migr_table[i].to;
1789 tmp->pop = migr_table[i].from;
1790 tmp2->tyme = migr_table[i].time;
1791 tmp2->type = 'm';
1792 tmp2->id = world->unique_id++;
1793 tmp2->actualpop = migr_table[i].to;
1794 tmp2->pop = migr_table[i].from;
1795 tmp->next = tmp2;
1796 tmp2->next = tmp;
1797 tmp2->top = 1;
1798 oldNode2 = down;
1799 theNode = down->back;
1800 while (theNode->tyme > migr_table[i].time)
1801 {
1802 oldNode2 = theNode->next;
1803 theNode = theNode->next->back;
1804 }
1805 tmp2->back = oldNode2;
1806 oldNode2->back = tmp2;
1807 }*/
1808 }
1809
1810
1811 void
children(node * mother,node ** brother,node ** sister)1812 children (node * mother, node ** brother, node ** sister)
1813 {
1814 node *m;
1815
1816 m = showtop (mother);
1817
1818 if (m->type == 't')
1819 {
1820 error ("this is a tip, so there are no more child nodes\n");
1821 }
1822 else
1823 {
1824 (*brother) = crawlback (m->next);
1825 (*sister) = crawlback (m->next->next);
1826 }
1827 }
1828 #endif /*not PRIORTEST*/
1829 /* Uses Lanczos-type approximation to ln(gamma) for z > 0. */
1830 /* Reference: */
1831 /* Lanczos, C. 'A precision approximation of the gamma */
1832 /* function', J. SIAM Numer. Anal., B, 1, 86-96, 1964. */
1833 /* Accuracy: About 14 significant digits except for small regions */
1834 /* in the vicinity of 1 and 2. */
1835 /* Programmer: Alan Miller */
1836 /* CSIRO Division of Mathematics & Statistics */
1837 /* Latest revision - 17 April 1988 */
1838 /* translated and modified into C by Peter Beerli 1997 */
1839 MYREAL
mylgamma(MYREAL z)1840 mylgamma (MYREAL z)
1841 {
1842 MYREAL a[9] = { 0.9999999999995183, 676.5203681218835,
1843 -1259.139216722289, 771.3234287757674, -176.6150291498386,
1844 12.50734324009056, -0.1385710331296526, 9.934937113930748e-6,
1845 1.659470187408462e-7
1846 };
1847 MYREAL lnsqrt2pi = 0.9189385332046727;
1848 MYREAL result;
1849 long j;
1850 MYREAL tmp;
1851 if (z <= 0.)
1852 {
1853 return MYREAL_MAX; /*this will kill the receiving calculation */
1854 }
1855 result = 0.;
1856 tmp = z + 7.;
1857 for (j = 9; j >= 2; --j)
1858 {
1859 result += a[j - 1] / tmp;
1860 tmp -= 1.;
1861 }
1862 result += a[0];
1863 result = log (result) + lnsqrt2pi - (z + 6.5) + (z - 0.5) * log (z + 6.5);
1864 return result;
1865 } /* lgamma */
1866
1867 /* ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3
1868 Computation of the Incomplete Gamma Integral
1869 Auxiliary functions required: lgamma() = logarithm of the gamma
1870 function, and alnorm() = algorithm AS66
1871 in Mathematica this is GammaRegularized[a,0,x] === Gamma[a,0,x]/Gamma[a]
1872 */
1873 MYREAL
incompletegamma(MYREAL tx,MYREAL talpha)1874 incompletegamma (MYREAL tx, MYREAL talpha)
1875 {
1876 const MYREAL eps = DBL_EPSILON ;
1877 double gama, d_1, d_2, d_3;
1878 /*static */
1879 double a, b, c, an, rn;
1880 /*static */
1881 double pn1, pn2, pn3, pn4, pn5, pn6, arg;
1882 double x = (double) tx;
1883 double alpha = (double) talpha;
1884
1885 gama = 0.;
1886 /* Check that we have valid values for X and P */
1887 if (alpha <= 0. || x < 0.)
1888 error ("failed in imcompletegamma(): wrong alpha or x\n");
1889 if (fabs (x) < eps)
1890 return (MYREAL) gama;
1891
1892 /* Use a normal approximation if P > PLIMIT */
1893 if (alpha > 1e3)
1894 {
1895 pn1 =
1896 sqrt (alpha) * 3. * (pow (x / alpha, (1. / 3.)) + 1. / (alpha * 9.) -
1897 1.);
1898 gama = alnorm (pn1, FALSE);
1899 return (MYREAL) gama;
1900 }
1901
1902 /* If X is extremely large compared to P then set GAMMAD = 1 */
1903 if (x > 1e8)
1904 {
1905 gama = 1.;
1906 return (MYREAL) gama;
1907 }
1908
1909 if (x <= 1. || x < alpha)
1910 {
1911 /* Use Pearson's series expansion. */
1912 /* (Note that P is not large enough to force overflow in lgamma()). */
1913 arg = alpha * LOG (x) - x - LGAMMA (alpha + 1.);
1914 c = 1.;
1915 gama = 1.;
1916 a = alpha;
1917 while (c > 1e-14)
1918 {
1919 a += 1.;
1920 c = c * x / a;
1921 gama += c;
1922 }
1923 arg += LOG (gama);
1924 gama = 0.;
1925 if (arg >= -88.)
1926 {
1927 gama = EXP (arg);
1928 }
1929
1930 }
1931 else
1932 {
1933 /* Use a continued fraction expansion */
1934 arg = alpha * LOG (x) - x - LGAMMA (alpha);
1935 a = 1. - alpha;
1936 b = a + x + 1.;
1937 c = 0.;
1938 pn1 = 1.;
1939 pn2 = x;
1940 pn3 = x + 1.;
1941 pn4 = x * b;
1942 gama = pn3 / pn4;
1943 for (;;)
1944 {
1945 a += 1.;
1946 b += 2.;
1947 c += 1.;
1948 an = a * c;
1949 pn5 = b * pn3 - an * pn1;
1950 pn6 = b * pn4 - an * pn2;
1951 if (fabs (pn6) > 0.)
1952 {
1953 rn = pn5 / pn6;
1954 /* Computing MIN */
1955 d_2 = 1e-14, d_3 = rn * 1e-14;
1956 if ((d_1 = gama - rn, fabs (d_1)) <= MIN (d_2, d_3))
1957 {
1958 arg += LOG (gama);
1959 gama = 1.;
1960 if (arg >= -88.)
1961 {
1962 gama = 1. - EXP (arg);
1963 }
1964 return (MYREAL) gama;
1965 }
1966 gama = rn;
1967 }
1968 pn1 = pn3;
1969 pn2 = pn4;
1970 pn3 = pn5;
1971 pn4 = pn6;
1972 if (fabs (pn5) >= 1e37)
1973 {
1974 /* Re-scale terms in continued fraction if terms are large */
1975 pn1 /= 1e37;
1976 pn2 /= 1e37;
1977 pn3 /= 1e37;
1978 pn4 /= 1e37;
1979 }
1980 }
1981 }
1982 return (MYREAL) gama;
1983 } /* incompletegamma() */
1984
1985
1986 /* calculation is replaced by the correct function in
1987 polygamma.c (which is a translation of a fortran program by amos
1988
1989 driver for the polygamma calculation */
1990 MYREAL
polygamma(long n,MYREAL z)1991 polygamma (long n, MYREAL z)
1992 {
1993 MYREAL ans;
1994 long nz, ierr;
1995 dpsifn (&z, &n, 1, 1, &ans, &nz, &ierr);
1996 if (n == 0)
1997 return -ans;
1998 else
1999 return ans;
2000 }
2001
2002 /*-------------------------------------------------------*/
2003 /* nrcheck subroutine (used in damped newton raphson proc */
2004 /* syntax: nrcheck(matrix,inversematrix,ncols=nrows,returnval1,returnval2) */
2005 /* mai 95 PB */
2006 boolean
nrcheck(MYREAL ** m,MYREAL ** tm,MYREAL * v,long nrows,MYREAL * r1,MYREAL * r2,boolean do_newton)2007 nrcheck (MYREAL **m, MYREAL **tm, MYREAL *v, long nrows, MYREAL *r1,
2008 MYREAL *r2, boolean do_newton)
2009 {
2010 long i, j, k;
2011 MYREAL *tmp, *tmp2, tmp3 = 0.0, tmp4 = 0.0;
2012 tmp = (MYREAL *) mycalloc (1, sizeof (MYREAL) * nrows);
2013 tmp2 = (MYREAL *) mycalloc (1, sizeof (MYREAL) * nrows);
2014 /*first evaluate r1 */
2015 (*r1) = (*r2) = 0.0;
2016 for (i = 0; i < nrows; i++)
2017 {
2018 (*r1) += v[i] * v[i];
2019 }
2020 /* T */
2021 for (j = 0; j < nrows; j++)
2022 { /* g . G */
2023 for (k = 0; k < nrows; k++)
2024 {
2025 tmp[j] += v[k] * m[j][k];
2026 tmp2[j] += v[k] * tm[j][k];
2027 }
2028 }
2029 /* T */
2030 for (i = 0; i < nrows; i++)
2031 { /* g . G . g */
2032 (*r2) += tmp[i] * v[i];
2033 tmp3 += tmp2[i] * v[i];
2034 }
2035 tmp4 = LOG (fabs ((*r1)));
2036 tmp4 = tmp4 + tmp4 - LOG (fabs ((*r2)));
2037 tmp4 = ((*r2) < 0 ? -1 : 1) * EXP (tmp4);
2038 myfree(tmp);
2039 if (do_newton && (tmp3 > (tmp4 > 0 ? tmp4 : 0)))
2040 {
2041 memcpy (v, tmp2, sizeof (MYREAL) * nrows);
2042 myfree(tmp2);
2043 return TRUE;
2044 }
2045 myfree(tmp2);
2046 return FALSE;
2047 }
2048
2049
2050 /*-------------------------------------------------------*/
2051 /* Matrix inversion subroutine */
2052 /* The passed matrix will be replaced by its inverse!!!!! */
2053 /* Gauss-Jordan reduction -- invert matrix a in place, */
2054 /* overwriting previous contents of a. On exit, matrix a */
2055 /* contains the inverse. */
2056 void
invert_matrix(MYREAL ** a,long nsize)2057 invert_matrix (MYREAL **a, long nsize)
2058 {
2059 long i, j;
2060 long *indeks;
2061 MYREAL *column, **result;
2062 indeks = (long *) mymalloc (sizeof (long) * nsize);
2063 column = (MYREAL *) mymalloc (sizeof (MYREAL) * nsize);
2064 result = (MYREAL **) mymalloc (sizeof (MYREAL *) * nsize);
2065 for (i = 0; i < nsize; i++)
2066 {
2067 result[i] = (MYREAL *) mymalloc (sizeof (MYREAL) * nsize);
2068 }
2069 lu_decomp (a, indeks, nsize);
2070 for (j = 0; j < nsize; j++)
2071 {
2072 memset (column, 0, sizeof (MYREAL) * nsize);
2073 column[j] = 1.0;
2074 lu_substitution (a, indeks, column, nsize);
2075 for (i = 0; i < nsize; i++)
2076 result[i][j] = column[i];
2077 }
2078 for (i = 0; i < nsize; i++)
2079 {
2080 memcpy (a[i], result[i], sizeof (MYREAL) * nsize);
2081 myfree(result[i]);
2082 }
2083 myfree(result);
2084 myfree(column);
2085 myfree(indeks);
2086 }
2087
2088 /*=======================================================*/
2089
2090 /*-------------------------------------------------------*/
2091 /* LU decomposition */
2092 /* after Dahlquist et al. 1974 and Press et al. 1988 */
2093 /* the method's uses Crout's procedure and the pivoting */
2094 /* described in Press et al. */
2095 /* Syntax: lu_decomp(matrix, indeks, nrows) */
2096 /* matrix will be destroyed and filled with the two */
2097 /* triangular matrices, indeks is the index vector for the */
2098 /* pivoting and the row change in case of 0 pivot values */
2099 /* nrows is the number of rows and columns in matrix */
2100 /* april 95 PB */
2101 void
lu_decomp(MYREAL ** m,long * indeks,long nrows)2102 lu_decomp (MYREAL **m, long *indeks, long nrows)
2103 {
2104 long i, j, k, p, kmax = -1;
2105 MYREAL *max_row_vals, big, summ, pivot, bigt;
2106 max_row_vals = (MYREAL *) mycalloc (1, sizeof (MYREAL) * nrows);
2107 for (i = 0; i < nrows; i++)
2108 {
2109 big = 0.0;
2110 for (j = 0; j < nrows; j++)
2111 {
2112 if ((bigt = fabs (m[i][j])) > big)
2113 big = bigt;
2114 }
2115 max_row_vals[i] = 1.0 / big;
2116 if (big == 0.0)
2117 {
2118 error ("Singular matrix detected in lu_decomp\n");
2119 }
2120 }
2121 for (i = 0; i < nrows; i++)
2122 {
2123 for (k = 0; k < i; k++)
2124 { /* upper half of matrix */
2125 summ = m[k][i];
2126 for (p = 0; p < k; p++)
2127 summ -= m[k][p] * m[p][i];
2128 m[k][i] = summ;
2129 }
2130 big = 0.0;
2131 for (k = i; k < nrows; k++)
2132 { /* lower half of matrix */
2133 summ = m[k][i];
2134 for (p = 0; p < i; p++)
2135 summ -= m[k][p] * m[p][i];
2136 m[k][i] = summ;
2137 pivot = fabs (summ) /**max_row_vals[k]*/ ;
2138 /* printf(stdout,"i=%li,pivot=%f,big=%f\n",i,pivot,big); */
2139 if (pivot >= big)
2140 {
2141 big = pivot;
2142 kmax = k;
2143 }
2144 }
2145 if (i != kmax)
2146 {
2147 for (p = 0; p < nrows; p++)
2148 {
2149 pivot = m[kmax][p];
2150 m[kmax][p] = m[i][p];
2151 m[i][p] = pivot;
2152 }
2153 max_row_vals[kmax] = max_row_vals[i];
2154 }
2155 indeks[i] = kmax;
2156 if (m[i][i] == 0.0)
2157 m[i][i] = SMALL_VALUE;
2158 if (i != nrows - 1)
2159 {
2160 pivot = 1. / m[i][i];
2161 for (k = i + 1; k < nrows; k++)
2162 m[k][i] *= pivot;
2163 }
2164 }
2165 myfree(max_row_vals);
2166 } /* end of lu_decomp */
2167
2168 /*-------------------------------------------------------*/
2169 /* LU substitution */
2170 /* after Dahlquist et al. 1974 and Press et al. 1988 */
2171 /* needs first the evaluation LU decomposition */
2172 /* Syntax: lu_substition(matrix, indeks, vector, nrows) */
2173 /* matrix = LU decomposed matrix, indeks = order of matrix */
2174 /* vector = value vector, nrows = number of rows/columns */
2175 /* april 95 PB */
2176 void
lu_substitution(MYREAL ** m,long * indeks,MYREAL * v,long nrows)2177 lu_substitution (MYREAL **m, long *indeks, MYREAL *v, long nrows)
2178 {
2179 long i, j;
2180 MYREAL summ;
2181 for (i = 0; i < nrows; i++)
2182 {
2183 summ = v[indeks[i]];
2184 v[indeks[i]] = v[i];
2185 for (j = 0; j < i; j++)
2186 summ -= m[i][j] * v[j];
2187 v[i] = summ;
2188 }
2189 for (i = nrows - 1; i >= 0; i--)
2190 {
2191 summ = v[i];
2192 for (j = i + 1; j < nrows; j++)
2193 summ -= m[i][j] * v[j];
2194 v[i] = summ / m[i][i];
2195 }
2196 }
2197
2198
2199 /* Algorithm AS66 Applied Statistics (1973) vol22 no.3
2200 Evaluates the tail area of the standardised normal curve
2201 from x to infinity if upper is .true. or
2202 from minus infinity to x if upper is .false. */
2203 MYREAL
alnorm(MYREAL x,int up)2204 alnorm (MYREAL x, int up)
2205 {
2206 /* Initialized data */
2207 /* *** machine dependent constants ????????????? */
2208 /*static */ MYREAL zero = 0.;
2209 /*static */
2210 MYREAL a1 = 5.75885480458;
2211 /*static */
2212 MYREAL a2 = 2.62433121679;
2213 /*static */
2214 MYREAL a3 = 5.92885724438;
2215 /*static */
2216 MYREAL b1 = -29.8213557807;
2217 /*static */
2218 MYREAL b2 = 48.6959930692;
2219 /*static */
2220 MYREAL c1 = -3.8052e-8;
2221 /*static */
2222 MYREAL c2 = 3.98064794e-4;
2223 /*static */
2224 MYREAL c3 = -.151679116635;
2225 /*static */
2226 MYREAL c4 = 4.8385912808;
2227 /*static */
2228 MYREAL c5 = .742380924027;
2229 /*static */
2230 MYREAL one = 1.;
2231 /*static */
2232 MYREAL c6 = 3.99019417011;
2233 /*static */
2234 MYREAL d1 = 1.00000615302;
2235 /*static */
2236 MYREAL d2 = 1.98615381364;
2237 /*static */
2238 MYREAL d3 = 5.29330324926;
2239 /*static */
2240 MYREAL d4 = -15.1508972451;
2241 /*static */
2242 MYREAL d5 = 30.789933034;
2243 /*static */
2244 MYREAL half = .5;
2245 /*static */
2246 MYREAL ltone = 7.;
2247 /*static */
2248 MYREAL utzero = 18.66;
2249 /*static */
2250 MYREAL con = 1.28;
2251 /*static */
2252 MYREAL p = .398942280444;
2253 /*static */
2254 MYREAL q = .39990348504;
2255 /*static */
2256 MYREAL r = .398942280385;
2257
2258 /*static */
2259 MYREAL y, result;
2260
2261 if (x < zero)
2262 {
2263 up = !up;
2264 x = -x;
2265 }
2266 if (x <= ltone || (up && x <= utzero))
2267 {
2268 y = half * x * x;
2269 if (x > con)
2270 {
2271 result =
2272 r * EXP (-y) / (x + c1 +
2273 d1 / (x + c2 +
2274 d2 / (x + c3 +
2275 d3 / (x + c4 +
2276 d4 / (x + c5 +
2277 d5 / (x + c6))))));
2278 return ((!up) ? one - result : result);
2279 }
2280 result =
2281 half - x * (p - q * y / (y + a1 + b1 / (y + a2 + b2 / (y + a3))));
2282 return ((!up) ? one - result : result);
2283 }
2284 else
2285 {
2286 return ((!up) ? 1.0 : 0.);
2287 }
2288 /*fake */ return -99;
2289 } /* alnorm */
2290
2291 /* dpsifn.c -- translated by f2c (version 19950808).
2292 and hand-patched by Peter Beerli Seattle, 1996
2293 SUBROUTINE DPSIFN (X, N, KODE, M, ANS, NZ, IERR)
2294
2295 C***BEGIN PROLOGUE DPSIFN
2296 C***PURPOSE Compute derivatives of the Psi function.
2297 C***LIBRARY SLATEC
2298 C***CATEGORY C7C
2299 C***TYPE MYREAL PRECISION (PSIFN-S, DPSIFN-D)
2300 C***KEYWORDS DERIVATIVES OF THE GAMMA FUNCTION, POLYGAMMA FUNCTION,
2301 C PSI FUNCTION
2302 C***AUTHOR Amos, D. E., (SNLA)
2303 C***DESCRIPTION
2304 C
2305 C The following definitions are used in DPSIFN:
2306 C
2307 C Definition 1
2308 C PSI(X) = d/dx (ln(GAMMA(X)), the first derivative of
2309 C the log GAMMA function.
2310 C Definition 2
2311 C K K
2312 C PSI(K,X) = d /dx (PSI(X)), the K-th derivative of PSI(X).
2313 C ___________________________________________________________________
2314 C DPSIFN computes a sequence of SCALED derivatives of
2315 C the PSI function; i.e. for fixed X and M it computes
2316 C the M-member sequence
2317 C
2318 C ((-1)**(K+1)/GAMMA(K+1))*PSI(K,X)
2319 C for K = N,...,N+M-1
2320 C
2321 C where PSI(K,X) is as defined above. For KODE=1, DPSIFN returns
2322 C the scaled derivatives as described. KODE=2 is operative only
2323 C when K=0 and in that case DPSIFN returns -PSI(X) + LN(X). That
2324 C is, the logarithmic behavior for large X is removed when KODE=2
2325 C and K=0. When sums or differences of PSI functions are computed
2326 C the logarithmic terms can be combined analytically and computed
2327 C separately to help retain significant digits.
2328 C
2329 C Note that CALL DPSIFN(X,0,1,1,ANS) results in
2330 C ANS = -PSI(X)
2331 C
2332 C Input X is MYREAL PRECISION
2333 C X - Argument, X .gt. 0.0D0
2334 C N - First member of the sequence, 0 .le. N .le. 100
2335 C N=0 gives ANS(1) = -PSI(X) for KODE=1
2336 C -PSI(X)+LN(X) for KODE=2
2337 C KODE - Selection parameter
2338 C KODE=1 returns scaled derivatives of the PSI
2339 C function.
2340 C KODE=2 returns scaled derivatives of the PSI
2341 C function EXCEPT when N=0. In this case,
2342 C ANS(1) = -PSI(X) + LN(X) is returned.
2343 C M - Number of members of the sequence, M.ge.1
2344 C
2345 C Output ANS is MYREAL PRECISION
2346 C ANS - A vector of length at least M whose first M
2347 C components contain the sequence of derivatives
2348 C scaled according to KODE.
2349 C NZ - Underflow flag
2350 C NZ.eq.0, A normal return
2351 C NZ.ne.0, Underflow, last NZ components of ANS are
2352 C set to zero, ANS(M-K+1)=0.0, K=1,...,NZ
2353 C IERR - Error flag
2354 C IERR=0, A normal return, computation completed
2355 C IERR=1, Input error, no computation
2356 C IERR=2, Overflow, X too small or N+M-1 too
2357 C large or both
2358 C IERR=3, Error, N too large. Dimensioned
2359 C array TRMR(NMAX) is not large enough for N
2360 C
2361 C The nominal computational accuracy is the maximum of unit
2362 C roundoff (=D1MACH(4)) and 1.0D-18 since critical constants
2363 C are given to only 18 digits.
2364 C
2365 C PSIFN is the single precision version of DPSIFN.
2366 C
2367 C *Long Description:
2368 C
2369 C The basic method of evaluation is the asymptotic expansion
2370 C for large X.ge.XMIN followed by backward recursion on a two
2371 C term recursion relation
2372 C
2373 C W(X+1) + X**(-N-1) = W(X).
2374 C
2375 C This is supplemented by a series
2376 C
2377 C SUM( (X+K)**(-N-1) , K=0,1,2,... )
2378 C
2379 C which converges rapidly for large N. Both XMIN and the
2380 C number of terms of the series are calculated from the unit
2381 C roundoff of the machine environment.
2382 C
2383 C***REFERENCES Handbook of Mathematical Functions, National Bureau
2384 C of Standards Applied Mathematics Series 55, edited
2385 C by M. Abramowitz and I. A. Stegun, equations 6.3.5,
2386 C 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964.
2387 C D. E. Amos, A portable Fortran subroutine for
2388 C derivatives of the Psi function, Algorithm 610, ACM
2389 C Transactions on Mathematical Software 9, 4 (1983),
2390 C pp. 494-502.
2391 C***ROUTINES CALLED D1MACH, I1MACH
2392 C***REVISION HISTORY (YYMMDD)
2393 C 820601 DATE WRITTEN
2394 C 890531 Changed all specific intrinsics to generic. (WRB)
2395 C 890911 Removed unnecessary intrinsics. (WRB)
2396 C 891006 Cosmetic changes to prologue. (WRB)
2397 C 891006 REVISION DATE from Version 3.2
2398 C 891214 Prologue converted to Version 4.0 format. (BAB)
2399 C 920501 Reformatted the REFERENCES section. (WRB)
2400 C***END PROLOGUE DPSIFN
2401
2402
2403 */
2404
2405 /*static*/ long fifteen = 15;
2406 /*static*/
2407 long sixteen = 16;
2408 /*static*/
2409 long five = 5;
2410 /*static*/
2411 long four = 4;
2412 /*static*/
2413 long fourteen = 14;
2414
2415 MYREAL
d1mach(long i)2416 d1mach (long i)
2417 {
2418 //#ifdef MYREAL == float
2419 //const MYREAL eps = FLT_EPSILON ;
2420 //const MYREAL numbermax = FLT_MAX;
2421 //const MYREAL numbermin = FLT_MIN;
2422 //#else
2423 const MYREAL eps = DBL_EPSILON ;
2424 const MYREAL numbermax = MYREAL_MAX;
2425 const MYREAL numbermin = DBL_MIN;
2426 //#endif
2427
2428 switch (i)
2429 {
2430 case 1:
2431 return numbermin;
2432 case 2:
2433 return numbermax;
2434 case 3:
2435 return eps / FLT_RADIX;
2436 case 4:
2437 return eps;
2438 case 5:
2439 return log10 ((MYREAL) FLT_RADIX);
2440 }
2441 usererror ("invalid argument: d1mach(%ld)\n", i);
2442 return 0; /* for compilers that complain of missing return values */
2443 }
2444
2445 long
i1mach(long i)2446 i1mach (long i)
2447 {
2448 switch (i)
2449 {
2450 case 1:
2451 return 5; /* standard input */
2452 case 2:
2453 return 6; /* standard output */
2454 case 3:
2455 return 7; /* standard punch */
2456 case 4:
2457 return 0; /* standard error */
2458 case 5:
2459 return 32; /* bits per integer */
2460 case 6:
2461 return 1; /* Fortran 77 value */
2462 case 7:
2463 return 2; /* base for integers */
2464 case 8:
2465 return 31; /* digits of integer base */
2466 case 9:
2467 return LONG_MAX;
2468 case 10:
2469 return FLT_RADIX;
2470 case 11:
2471 return FLT_MANT_DIG;
2472 case 12:
2473 return FLT_MIN_EXP;
2474 case 13:
2475 return FLT_MAX_EXP;
2476 case 14:
2477 return DBL_MANT_DIG;
2478 case 15:
2479 return DBL_MIN_EXP;
2480 case 16:
2481 return DBL_MAX_EXP;
2482 }
2483 usererror ("invalid argument: i1mach(%ld)\n", i);
2484 return 0; /* for compilers that complain of missing return values */
2485 }
2486
2487 int
dpsifn(MYREAL * x,long * n,long kode,long m,MYREAL * ans,long * nz,long * ierr)2488 dpsifn (MYREAL *x, long *n, long kode, long m, MYREAL *ans, long *nz,
2489 long *ierr)
2490 {
2491 /* Initialized data */
2492
2493 /*static */ long nmax = 100;
2494 /*static */
2495 MYREAL b[22] = { 1., -.5, .166666666666666667,
2496 -.0333333333333333333, .0238095238095238095, -.0333333333333333333,
2497 .0757575757575757576, -.253113553113553114, 1.16666666666666667,
2498 -7.09215686274509804, 54.9711779448621554, -529.124242424242424,
2499 6192.1231884057971, -86580.2531135531136, 1425517.16666666667,
2500 -27298231.067816092, 601580873.900642368, -15116315767.0921569,
2501 429614643061.166667, -13711655205088.3328, 488332318973593.167,
2502 -19296579341940068.1
2503 };
2504
2505 /* System generated locals */
2506 long i1, i2;
2507 MYREAL d1, d2;
2508
2509
2510 /* Local variables */
2511 /*static */
2512 MYREAL elim, xinc, xmin, tols, xdmy, yint, trmr[100], rxsq;
2513 /*static */
2514 long i__, j, k;
2515 /*static */
2516 MYREAL s, t, slope, xdmln, wdtol;
2517 /*static */
2518 MYREAL t1, t2;
2519 /*static */
2520 long fn;
2521 /*static */
2522 MYREAL ta;
2523 /*static */
2524 long mm, nn, np;
2525 /*static */
2526 MYREAL fx, tk;
2527 /*static */
2528 long mx, nx;
2529 /*static */
2530 MYREAL xm, tt, xq, den, arg, fln, r1m4, r1m5, eps, rln, tol,
2531 xln, trm[22], tss, tst;
2532 int i;
2533 for(i=0;i<22;i++)
2534 trm[i]=0.0;
2535 for(i=0;i<100;i++)
2536 trmr[i]=0.0;
2537
2538 /* Parameter adjustments */
2539 --ans;
2540
2541 /* Function Body */
2542 /* ----------------------------------------------------------------------- */
2543 /* BERNOULLI NUMBERS */
2544 /* ----------------------------------------------------------------------- */
2545
2546 /* ***FIRST EXECUTABLE STATEMENT DPSIFN */
2547 *ierr = 0;
2548 *nz = 0;
2549 if (*x <= 0.)
2550 {
2551 *ierr = 1;
2552 }
2553 if (*n < 0)
2554 {
2555 *ierr = 1;
2556 }
2557 if (kode < 1 || kode > 2)
2558 {
2559 *ierr = 1;
2560 }
2561 if (m < 1)
2562 {
2563 *ierr = 1;
2564 }
2565 if (*ierr != 0)
2566 {
2567 return 0;
2568 }
2569 mm = m;
2570 /* Computing MIN */
2571 //xcode i1 = -fifteen;
2572 //xcode i2 = sixteen;
2573 nx = MIN (-i1mach (fifteen), i1mach (sixteen));
2574 r1m5 = d1mach (five);
2575 r1m4 = d1mach (four) * .5;
2576 wdtol = MAX (r1m4, 5e-19);
2577 /* ----------------------------------------------------------------------- */
2578 /* ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT */
2579 /* ----------------------------------------------------------------------- */
2580 elim = (nx * r1m5 - 3.) * 2.302;
2581 xln = LOG (*x);
2582 L41:
2583 nn = *n + mm - 1;
2584 fn = nn;
2585 t = (fn + 1) * xln;
2586 /* ----------------------------------------------------------------------- */
2587 /* OVERFLOW AND UNDERFLOW TEST FOR SMALL AND LARGE X */
2588 /* ----------------------------------------------------------------------- */
2589 if (fabs (t) > elim)
2590 {
2591 goto L290;
2592 }
2593 if (*x < wdtol)
2594 {
2595 goto L260;
2596 }
2597 /* ----------------------------------------------------------------------- */
2598 /* COMPUTE XMIN AND THE NUMBER OF TERMS OF THE SERIES, FLN+1 */
2599 /* ----------------------------------------------------------------------- */
2600 rln = r1m5 * i1mach (fourteen);
2601 rln = MIN (rln, 18.06);
2602 fln = MAX (rln, 3.) - 3.;
2603 yint = fln * .4 + 3.5;
2604 slope = fln * (fln * 6.038e-4 + .008677) + .21;
2605 xm = yint + slope * fn;
2606 mx = (long) xm + 1;
2607 xmin = (MYREAL) mx;
2608 if (*n == 0)
2609 {
2610 goto L50;
2611 }
2612 xm = rln * -2.302 - MIN (0., xln);
2613 arg = xm / *n;
2614 arg = MIN (0., arg);
2615 eps = EXP (arg);
2616 xm = 1. - eps;
2617 if (fabs (arg) < .001)
2618 {
2619 xm = -arg;
2620 }
2621 fln = *x * xm / eps;
2622 xm = xmin - *x;
2623 if (xm > 7. && fln < 15.)
2624 {
2625 goto L200;
2626 }
2627 L50:
2628 xdmy = *x;
2629 xdmln = xln;
2630 xinc = 0.;
2631 if (*x >= xmin)
2632 {
2633 goto L60;
2634 }
2635 nx = (long) (*x);
2636 xinc = xmin - nx;
2637 xdmy = *x + xinc;
2638 xdmln = LOG (xdmy);
2639 L60:
2640 /* ----------------------------------------------------------------------- */
2641 /* GENERATE W(N+MM-1,X) BY THE ASYMPTOTIC EXPANSION */
2642 /* ----------------------------------------------------------------------- */
2643 t = fn * xdmln;
2644 t1 = xdmln + xdmln;
2645 t2 = t + xdmln;
2646 /* Computing MAX */
2647 d1 = fabs (t), d2 = fabs (t1), d1 = MAX (d1, d2), d2 = fabs (t2);
2648 tk = MAX (d1, d2);
2649 if (tk > elim)
2650 {
2651 goto L380;
2652 }
2653 tss = EXP (-t);
2654 tt = .5 / xdmy;
2655 t1 = tt;
2656 tst = wdtol * tt;
2657 if (nn != 0)
2658 {
2659 t1 = tt + 1. / fn;
2660 }
2661 rxsq = 1. / (xdmy * xdmy);
2662 ta = rxsq * .5;
2663 t = (fn + 1) * ta;
2664 s = t * b[2];
2665 if (fabs (s) < tst)
2666 {
2667 goto L80;
2668 }
2669 tk = 2.;
2670 for (k = 4; k <= 22; ++k)
2671 {
2672 t = t * ((tk + fn + 1) / (tk + 1.)) * ((tk + fn) / (tk + 2.)) * rxsq;
2673 trm[k - 1] = t * b[k - 1];
2674 if ((d1 = trm[k - 1], fabs (d1)) < tst)
2675 {
2676 goto L80;
2677 }
2678 s += trm[k - 1];
2679 tk += 2.;
2680 /* L70: */
2681 }
2682 L80:
2683 s = (s + t1) * tss;
2684 if (xinc == 0.)
2685 {
2686 goto L100;
2687 }
2688 /* ----------------------------------------------------------------------- */
2689 /* BACKWARD RECUR FROM XDMY TO X */
2690 /* ----------------------------------------------------------------------- */
2691 nx = (long) xinc;
2692 np = nn + 1;
2693 if (nx > nmax)
2694 {
2695 goto L390;
2696 }
2697 if (nn == 0)
2698 {
2699 goto L160;
2700 }
2701 xm = xinc - 1.;
2702 fx = *x + xm;
2703 /* ----------------------------------------------------------------------- */
2704 /* THIS LOOP SHOULD NOT BE CHANGED. FX IS ACCURATE WHEN X IS SMALL */
2705 /* ----------------------------------------------------------------------- */
2706 i1 = nx;
2707 for (i__ = 1; i__ <= i1; ++i__)
2708 {
2709 i2 = -np;
2710 trmr[i__ - 1] = pow (fx, (MYREAL) i2);
2711 s += trmr[i__ - 1];
2712 xm += -1.;
2713 fx = *x + xm;
2714 /* L90: */
2715 }
2716 L100:
2717 ans[mm] = s;
2718 if (fn == 0)
2719 {
2720 goto L180;
2721 }
2722 /* ----------------------------------------------------------------------- */
2723 /* GENERATE LOWER DERIVATIVES, J.LT.N+MM-1 */
2724 /* ----------------------------------------------------------------------- */
2725 if (mm == 1)
2726 {
2727 return 0;
2728 }
2729 i1 = mm;
2730 for (j = 2; j <= i1; ++j)
2731 {
2732 --fn;
2733 tss *= xdmy;
2734 t1 = tt;
2735 if (fn != 0)
2736 {
2737 t1 = tt + 1. / fn;
2738 }
2739 t = (fn + 1) * ta;
2740 s = t * b[2];
2741 if (fabs (s) < tst)
2742 {
2743 goto L120;
2744 }
2745 tk = (MYREAL) (fn + 4);
2746 for (k = 4; k <= 22; ++k)
2747 {
2748 trm[k - 1] = trm[k - 1] * (fn + 1) / tk;
2749 if ((d1 = trm[k - 1], fabs (d1)) < tst)
2750 {
2751 goto L120;
2752 }
2753 s += trm[k - 1];
2754 tk += 2.;
2755 /* L110: */
2756 }
2757 L120:
2758 s = (s + t1) * tss;
2759 if (xinc == 0.)
2760 {
2761 goto L140;
2762 }
2763 if (fn == 0)
2764 {
2765 goto L160;
2766 }
2767 xm = xinc - 1.;
2768 fx = *x + xm;
2769 i2 = nx;
2770 for (i__ = 1; i__ <= i2; ++i__)
2771 {
2772 trmr[i__ - 1] *= fx;
2773 s += trmr[i__ - 1];
2774 xm += -1.;
2775 fx = *x + xm;
2776 /* L130: */
2777 }
2778 L140:
2779 mx = mm - j + 1;
2780 ans[mx] = s;
2781 if (fn == 0)
2782 {
2783 goto L180;
2784 }
2785 /* L150: */
2786 }
2787 return 0;
2788 /* ----------------------------------------------------------------------- */
2789 /* RECURSION FOR N = 0 */
2790 /* ----------------------------------------------------------------------- */
2791 L160:
2792 i1 = nx;
2793 for (i__ = 1; i__ <= i1; ++i__)
2794 {
2795 s += 1. / (*x + nx - i__);
2796 /* L170: */
2797 }
2798 L180:
2799 if (kode == 2)
2800 {
2801 goto L190;
2802 }
2803 ans[1] = s - xdmln;
2804 return 0;
2805 L190:
2806 if (xdmy == *x)
2807 {
2808 return 0;
2809 }
2810 xq = xdmy / *x;
2811 ans[1] = s - LOG (xq);
2812 return 0;
2813 /* ----------------------------------------------------------------------- */
2814 /* COMPUTE BY SERIES (X+K)**(-(N+1)) , K=0,1,2,... */
2815 /* ----------------------------------------------------------------------- */
2816 L200:
2817 nn = (long) fln + 1;
2818 np = *n + 1;
2819 t1 = (*n + 1) * xln;
2820 t = EXP (-t1);
2821 s = t;
2822 den = *x;
2823 i1 = nn;
2824 for (i__ = 1; i__ <= i1; ++i__)
2825 {
2826 den += 1.;
2827 i2 = -np;
2828 trm[i__ - 1] = pow (den, (MYREAL) i2);
2829 s += trm[i__ - 1];
2830 /* L210: */
2831 }
2832 ans[1] = s;
2833 if (*n != 0)
2834 {
2835 goto L220;
2836 }
2837 if (kode == 2)
2838 {
2839 ans[1] = s + xln;
2840 }
2841 L220:
2842 if (mm == 1)
2843 {
2844 return 0;
2845 }
2846 /* ----------------------------------------------------------------------- */
2847 /* GENERATE HIGHER DERIVATIVES, J.GT.N */
2848 /* ----------------------------------------------------------------------- */
2849 tol = wdtol / 5.;
2850 i1 = mm;
2851 for (j = 2; j <= i1; ++j)
2852 {
2853 t /= *x;
2854 s = t;
2855 tols = t * tol;
2856 den = *x;
2857 i2 = nn;
2858 for (i__ = 1; i__ <= i2; ++i__)
2859 {
2860 den += 1.;
2861 trm[i__ - 1] /= den;
2862 s += trm[i__ - 1];
2863 if (trm[i__ - 1] < tols)
2864 {
2865 goto L240;
2866 }
2867 /* L230: */
2868 }
2869 L240:
2870 ans[j] = s;
2871 /* L250: */
2872 }
2873 return 0;
2874 /* ----------------------------------------------------------------------- */
2875 /* SMALL X.LT.UNIT ROUND OFF */
2876 /* ----------------------------------------------------------------------- */
2877 L260:
2878 i1 = -(*n) - 1;
2879 ans[1] = pow (*x, (MYREAL) i1);
2880 if (mm == 1)
2881 {
2882 goto L280;
2883 }
2884 k = 1;
2885 i1 = mm;
2886 for (i__ = 2; i__ <= i1; ++i__)
2887 {
2888 ans[k + 1] = ans[k] / *x;
2889 ++k;
2890 /* L270: */
2891 }
2892 L280:
2893 if (*n != 0)
2894 {
2895 return 0;
2896 }
2897 if (kode == 2)
2898 {
2899 ans[1] += xln;
2900 }
2901 return 0;
2902 L290:
2903 if (t > 0.)
2904 {
2905 goto L380;
2906 }
2907 *nz = 0;
2908 *ierr = 2;
2909 return 0;
2910 L380:
2911 ++(*nz);
2912 ans[mm] = 0.;
2913 --mm;
2914 if (mm == 0)
2915 {
2916 return 0;
2917 }
2918 goto L41;
2919 L390:
2920 *nz = 0;
2921 *ierr = 3;
2922 return 0;
2923 } /* dpsifn_ */
2924
2925
2926
2927
2928 MYREAL
rannor(MYREAL mean,MYREAL sd)2929 rannor (MYREAL mean, MYREAL sd)
2930 {
2931 MYREAL r1, r2;
2932 r1 = RANDUM ();
2933 r2 = RANDUM ();
2934 return sd * sqrt (-2. * LOG (r1)) * cos (TWOPI * r2) + mean;
2935 }
2936
2937
2938 char
lowercase(int c)2939 lowercase (int c)
2940 {
2941 return (char) tolower (c);
2942 }
2943
2944 char
uppercase(int c)2945 uppercase (int c)
2946 {
2947 return (char) toupper (c);
2948 }
2949
upper(char * from,char ** to)2950 void upper(char *from, char **to)
2951 {
2952 long i=0;
2953 while(from[i] != '\0')
2954 {
2955 (*to)[i] = from[i];
2956 ++i;
2957 }
2958 (*to)[++i] = '\0';
2959 }
2960
2961 MYREAL
find_chi(long df,MYREAL prob)2962 find_chi (long df, MYREAL prob)
2963 {
2964 double a, b, m;
2965 double xb = 200.0;
2966 double xa = 0.0;
2967 double xm = 5.;
2968 double dprob = (double) prob;
2969 a = probchi (df, xa);
2970 m = probchi (df, xm);
2971 b = probchi (df, xb);
2972 while (fabs (m - dprob) > EPSILON)
2973 {
2974 if (m < dprob)
2975 {
2976 b = m;
2977 xb = xm;
2978 }
2979 else
2980 {
2981 a = m;
2982 xa = xm;
2983 }
2984 xm = (-(b * xa) + prob * xa + a * xb - dprob * xb) / (a - b); //(xa + xb)/2.;
2985
2986 m = probchi (df, xm);
2987 }
2988 return (MYREAL) xm;
2989 }
2990
2991
2992 MYREAL
probchi(long df,double chi)2993 probchi (long df, double chi)
2994 {
2995 const MYREAL eps = DBL_EPSILON ;
2996 double prob;
2997 double v = ((MYREAL) df) / 2.;
2998
2999 if (chi > eps && v > eps)
3000 {
3001 //lg = EXP (LGAMMA (v));
3002 prob = 1. - incompletegamma (chi / 2., v);
3003 }
3004 else
3005 prob = 1.0;
3006 // printf("prob=%f v=%f chi=%f lg(v/2)=%f ig(chi/2,v/2)=%f\n",
3007 // prob,v,chi,lg, incompletegamma(chi/2.,v/2.));
3008
3009 return prob;
3010 }
3011
3012 MYREAL
probchiboundary(MYREAL chi,long zeros,long all)3013 probchiboundary (MYREAL chi, long zeros, long all)
3014 {
3015 long nonzeros = all - zeros;
3016
3017 //xcode MYREAL a, b, m;
3018 MYREAL m;
3019 MYREAL xb = 1.0; //1.0-EPSILON/1000.;
3020 MYREAL xa = 0.0; //EPSILON/1000.;
3021 MYREAL xm = 0.51;
3022 if(all==0)
3023 return 1.;
3024 if (zeros == 0)
3025 {
3026 return probchi (all, chi);
3027 }
3028 //xcode a = chiboundary (zeros, nonzeros, xa);
3029 m = chiboundary (zeros, nonzeros, xm);
3030 //xcode b = chiboundary (zeros, nonzeros, xb);
3031 while (fabs (m - chi) > EPSILON
3032 && (fabs (xa - xm) > EPSILON && fabs (xb - xm) > EPSILON))
3033 {
3034 if (m < chi)
3035 {
3036 //xcode b = m;
3037 xb = xm;
3038 }
3039 else
3040 {
3041 //xcode a = m;
3042 xa = xm;
3043 }
3044 xm = /*(-(b * xa) + chi * xa + a * xb - chi * xb) / (a - b); */
3045 (xa + xb) / 2.;
3046
3047 m = chiboundary (zeros, nonzeros, xm);
3048 }
3049 return xm;
3050 }
3051
3052 MYREAL
chiboundary(long zeros,long nonzeros,MYREAL alpha)3053 chiboundary (long zeros, long nonzeros, MYREAL alpha)
3054 {
3055 // MYREAL prob;
3056 MYREAL sum = 0.;
3057 long i;
3058 long k = zeros;
3059 MYREAL freq;
3060 MYREAL summ;
3061 //printf("z=%li a=%4.2f ", k,alpha);
3062 for (i = 0; i <= zeros; i++)
3063 {
3064 // sum = zerovec[i] * (i+nonzeros) == 0 ? 0. : find_chi(i+nonzeros,alpha);
3065 freq = EXP (logfac (k) - logfac (k - i) - logfac (i) - LOG2 * k);
3066 summ = (i + nonzeros) == 0 ? 0. : find_chi (i + nonzeros, alpha);
3067 // printf(" %.2f(%.4f)",freq*pow(2.,k),summ);
3068 sum += freq * summ;
3069 }
3070 //printf("\n");
3071 return sum;
3072 }
3073
3074
3075 MYREAL
chisquare(long df,MYREAL alpha)3076 chisquare (long df, MYREAL alpha)
3077 {
3078 const MYREAL table05[] =
3079 {
3080 3.84146, 5.99147, 7.81473, 9.48773, 11.0705, 12.5916
3081 };
3082 const MYREAL table01[] =
3083 {
3084 6.63490, 9.21034, 11.3449, 13.2767, 15.0863, 16.8119
3085 };
3086 if( df<=6)
3087 {
3088 if (alpha == 0.05)
3089 return table05[df - 1];
3090 if (alpha == 0.01)
3091 return table01[df - 1];
3092 }
3093 return -100000000;
3094 // return generalized_gamma(0.5,df/2.0,1.0,alpha);
3095 }
3096
3097 MYREAL
calc_sum(MYREAL * vector,long n)3098 calc_sum (MYREAL *vector, long n)
3099 {
3100 long i;
3101 MYREAL summ = 0.0;
3102 for (i = 0; i < n; i++)
3103 summ += vector[i];
3104 return summ;
3105 }
3106
3107 //==========================================
3108 // searching and finding
3109
3110 boolean
find(long i,long * list,long listlen)3111 find (long i, long *list, long listlen)
3112 {
3113 long j;
3114 for (j = 0; j < listlen; j++)
3115 {
3116 if (i == list[j])
3117 return TRUE;
3118 }
3119 return FALSE;
3120 }
3121
3122 //====================================================
3123 // conversion between the different parameter schemes
3124 // returns the begining of mig_.i
3125 long
mstart(long pop,long numpop)3126 mstart (long pop, long numpop)
3127 {
3128 return numpop + pop * numpop - pop;
3129 }
3130
3131 // returns the end of mig_.i
3132 long
mend(long pop,long numpop)3133 mend (long pop, long numpop)
3134 {
3135 return numpop + pop * numpop - pop + numpop - 1;
3136 }
3137
3138 ///
3139 /// return the first element of population pop
3140 long
mmstart(long pop,long numpop)3141 mmstart (long pop, long numpop)
3142 {
3143 return pop * (numpop);
3144 }
3145
3146
3147 ///
3148 /// return the element past the last element in of population pop
3149 long
mmend(long pop,long numpop)3150 mmend (long pop, long numpop)
3151 {
3152 return pop * numpop + numpop;
3153 }
3154
3155 ///
3156 /// Returns the location in a full matrix given the abbreviated matrix
3157 /// given i,j coordinates in an abbreviated matrix {d,d,...,d,a,a,a,...,a,b,b,...,b, ...}
3158 /// it returns the position j* in a linearized full matrix with diagonal element "-"
3159 /// the linearized matrix would look like this {-, a,a,a,...,a, b, - , b, .... ,b, c, c, -, c, ....},
3160 /// the position of i is the same for the abbreviated and the full matrix, and is not returned because
3161 /// the calling function knows this already.
3162 long
mm2m(long frompop,long topop,long numpop)3163 mm2m (long frompop, long topop, long numpop)
3164 {
3165 if (frompop == topop)
3166 return (frompop);
3167 if (frompop < topop)
3168 return numpop + topop * (numpop - 1) + frompop;
3169 else
3170 return numpop + topop * (numpop - 1) + (frompop - 1);
3171 }
3172
3173 ///
3174 /// Calulates the j and i from a linear abbreviated matrix
3175 /// position z in {d,d,d,...,d,a,a,...,a,b,b,...,b, c, ...} we know how many populations c (columns or rows)
3176 /// exist and wnat to calculate the i,j coordinates in the c x c matrix {{d,a,a,...,a},{b,d,b,...},...}
3177 /// in the function z above is im and i = frompop, j is topop
3178 /// frompop and topop contain the result
3179 void
m2mm(long i,long numpop,long * frompop,long * topop)3180 m2mm (long i, long numpop, long *frompop, long *topop)
3181 {
3182 if (i < numpop)
3183 {
3184 *frompop = i;
3185 *topop = i;
3186 return;
3187 }
3188 else
3189 {
3190 (*topop) = (long) (i - numpop) / (numpop - 1);
3191 (*frompop) = i - numpop - (*topop) * (numpop - 1);
3192 if (*frompop >= *topop)
3193 *frompop += 1;
3194 }
3195 }
3196
3197
3198 ///
3199 /// position i in linear abbreviated array
3200 long
m2mml(long i,long numpop)3201 m2mml (long i, long numpop)
3202 {
3203 long topop, frompop;
3204
3205 if (i < numpop)
3206 {
3207 return i * numpop + i;
3208 }
3209 else
3210 {
3211 topop = (long) (i - numpop) / (numpop - 1);
3212 frompop = i - numpop - (topop) * (numpop - 1);
3213 if (frompop >= topop)
3214 frompop += 1;
3215 return numpop * topop + frompop;
3216 }
3217 }
3218
3219
3220 long
mml2m(long pos,long numpop)3221 mml2m (long pos, long numpop)
3222 {
3223 long topop = 0, frompop = 0, i = 1;
3224 if (pos == 0)
3225 return 0;
3226 while (pos > numpop * (i++))
3227 topop++;
3228 frompop = pos - topop * numpop;
3229 return mm2m (frompop, topop, numpop);
3230 }
3231
3232
3233 long
m2mml2(long i,long topop,long numpop)3234 m2mml2 (long i, long topop, long numpop)
3235 {
3236 long frompop;
3237
3238 if (i < numpop)
3239 {
3240 return i * numpop + i;
3241 }
3242 else
3243 {
3244 frompop = i - numpop - (topop) * (numpop - 1);
3245 if (frompop >= topop)
3246 frompop += 1;
3247 return numpop * topop + frompop;
3248 }
3249 }
3250
3251 #ifndef PRIORTEST
set_paramstr(char * paramstr,long j,long numpop,boolean usem)3252 void set_paramstr(char *paramstr, long j, long numpop, boolean usem)
3253 {
3254 long frompop;
3255 long topop;
3256 const long numpop2 = numpop * numpop;
3257 if (j < numpop)
3258 sprintf(paramstr,"Theta_%-3li",j+1);
3259 else
3260 {
3261 if (j < numpop2)
3262 {
3263 m2mm (j, numpop, &frompop, &topop);
3264 if(usem)
3265 {
3266 sprintf(paramstr, "M_%li->%li", frompop+1, topop+1);
3267 }
3268 else
3269 {
3270 sprintf(paramstr, "xN_%lim_%li->%li", topop+1, frompop+1, topop+1);
3271 }
3272 }
3273 else
3274 {
3275 sprintf(paramstr, "Rate");
3276 }
3277 }
3278 }
3279
3280 void
gamma_rates(MYREAL * rate,MYREAL * probcat,long categs,char * input)3281 gamma_rates (MYREAL *rate, MYREAL *probcat, long categs, char *input)
3282 {
3283 long i;
3284 MYREAL alpha = MYREAL_MAX;
3285 MYREAL value;
3286 while (!isdigit (*input) && *input != '\0')
3287 input++;
3288 if ((value = strtod (input, (char **) NULL)) > 0)
3289 alpha = value;
3290 initgammacat (categs, alpha, 1., rate, probcat);
3291 for (i = 0; i < categs; i++)
3292 {
3293 probcat[i] = EXP (probcat[i]);
3294 }
3295 // calc_gamma (alpha, rate, categs);
3296 }
3297
3298 /* calculation of rate values following a gamma distribution for
3299 given probability values */
3300 void
calc_gamma(MYREAL alpha,MYREAL * gama,long categs)3301 calc_gamma (MYREAL alpha, MYREAL *gama, long categs)
3302 {
3303 long i, panic;
3304 MYREAL low, mid, high, xlow, xhigh, tmp, freq = 0, x = 10, elements =
3305 (MYREAL) categs;
3306 freq = -(0.5 / elements); /*so we have midpoints instead of endpoints */
3307 for (i = 0; i < categs; i++)
3308 {
3309 low = 0;
3310 mid = incompletegamma (10., alpha);
3311 high = 1.;
3312 freq += 1. / (elements);
3313 if (freq < mid)
3314 {
3315 high = mid;
3316 xlow = 0;
3317 xhigh = 10.;
3318 x = 5.;
3319 }
3320 else
3321 {
3322 low = mid;
3323 xhigh = 1e10;
3324 xlow = 10.;
3325 x = 1e5;
3326 }
3327 panic = 0;
3328 while (panic++ < 1000 && fabs (low - high) > 0.0001 && x > 0.000000001)
3329 {
3330 mid = incompletegamma (x, alpha);
3331 if (freq < mid)
3332 {
3333 high = mid;
3334 tmp = x;
3335 x = (x + xlow) / 2.;
3336 xhigh = tmp;
3337 }
3338 else
3339 {
3340 low = mid;
3341 tmp = x;
3342 x = (x + xhigh) / 2.;
3343 xlow = tmp;
3344 }
3345 }
3346 gama[i] = x / alpha;
3347 //Debug
3348 // printf (stderr, " %li> %f\n", i, gama[i]);
3349
3350 if (x >= 10e10)
3351 {
3352 error ("calc_gamma(): x is too big");
3353 }
3354 }
3355 }
3356
3357
3358 void
fprintf2(FILE * file,long filesize,const char * fmt,...)3359 fprintf2(FILE *file, long filesize, const char *fmt, ...)
3360 {
3361 char *p;
3362 va_list ap;
3363 long bufsize;
3364 long pallocsize = filesize+strlen(fmt)+1;
3365 p = (char *) mycalloc(pallocsize,sizeof(char));
3366 va_start(ap, fmt);
3367 bufsize = vsprintf(p, fmt, ap);
3368 if(bufsize>=pallocsize)
3369 error("failed in printf2()");
3370 fprintf(file,"%s", p);
3371 va_end(ap);
3372 myfree(p);
3373 }
3374
3375
3376 void
print_line(FILE * outfile,char c,long nn,long flag)3377 print_line (FILE * outfile, char c, long nn, long flag)
3378 {
3379 long i, start = 0;
3380 switch (flag)
3381 {
3382 case START:
3383 start = 2;
3384 FPRINTF (outfile, "=--");
3385 break;
3386 case STOP:
3387 start = 2;
3388 FPRINTF (outfile, "==-");
3389 break;
3390 default:
3391 start = 0;
3392 }
3393 for (i = start; i < nn; i++)
3394 {
3395 FPRINTF (outfile, "%c",c);
3396 }
3397 FPRINTF(outfile, "\n");
3398 }
3399
3400
3401 void
sprint_line(char * buffer,char c,long nn,long flag)3402 sprint_line
3403 (char *buffer, char c, long nn, long flag)
3404 {
3405 char ch[2];
3406 long i, start = 0;
3407 char fp[LINESIZE];
3408 buffer[0] = '\0';
3409 ch[0] = c;
3410 ch[1] = '\0';
3411 switch (flag)
3412 {
3413 case START:
3414 start = 2;
3415 sprintf (fp, "=%c%c", c, c);
3416 strcat (buffer, fp);
3417 break;
3418 case STOP:
3419 start = 2;
3420 sprintf (fp, "==%c", c);
3421 strcat (buffer, fp);
3422 break;
3423 default:
3424 start = 0;
3425 }
3426 for (i = start; i < nn; i++)
3427 {
3428 strcat (buffer, ch);
3429 }
3430 strcat (buffer, "\n");
3431 }
3432
3433 char
sgetc(char ** buffer)3434 sgetc (char **buffer)
3435 {
3436 char ch;
3437 ch = **buffer;
3438 (*buffer)++;
3439 return ch;
3440 }
3441
3442 /// line-end transparent string gets command
3443 char *
sgets(char * s,int size,char ** stream)3444 sgets (char *s, int size, char **stream)
3445 {
3446 long ch = '\0';
3447 long counter = 0;
3448 while (counter < size - 1)
3449 {
3450 ch = **stream;
3451 (*stream)++;
3452 switch (ch)
3453 {
3454 case '\0':
3455 case '\r':
3456 case '\n':
3457 s[counter] = '\0';
3458 return s;
3459 default:
3460 s[counter] = (char) ch;
3461 break;
3462 }
3463 counter++;
3464 }
3465 return s;
3466 }
3467 /// line-end transparent string gets command,
3468 /// this needs and end on the stream but allocates memory for s
3469 char *
sgets_safe(char ** s,long * size,char ** stream)3470 sgets_safe (char **s, long *size, char **stream)
3471 {
3472 boolean notdone=TRUE;
3473 long ch = '\0';
3474 long counter = 0;
3475 while (notdone)
3476 {
3477 ch = **stream;
3478 (*stream)++;
3479 switch (ch)
3480 {
3481 case '\0':
3482 case '\r':
3483 case '\n':
3484 (*s)[counter] = '\0';
3485 return *s;
3486 default:
3487 (*s)[counter] = (char) ch;
3488 break;
3489 }
3490 if(counter >= *size)
3491 {
3492 *size *= 2;
3493 *s = (char*) realloc(*s, *size * sizeof(char));
3494 }
3495 counter++;
3496 }
3497 return *s;
3498 }
3499
3500
3501 /// sticks a text into the printing buffer
add_to_buffer(char * fp,long * bufsize,char ** buffer,long * allocbufsize)3502 void add_to_buffer(char *fp, long *bufsize, char **buffer, long *allocbufsize)
3503 {
3504 long fpsize = (long) strlen (fp) + 1;
3505 if(*allocbufsize <= (*bufsize + fpsize))
3506 {
3507 *allocbufsize += 100 * fpsize;
3508 (*buffer) = (char *) myrealloc (*buffer, (*allocbufsize) * sizeof (char));
3509 }
3510 (*bufsize) += sprintf((*buffer) + (*bufsize),"%s",fp);
3511 }
3512
3513
3514 /// sticks a text into the printing buffer
print_to_buffer(char ** buffer,long * maxbufsize,char * tempbuffer,long * pos,const char * fmt,...)3515 long print_to_buffer(char **buffer, long *maxbufsize, char *tempbuffer, long *pos, const char *fmt, ...)
3516 {
3517 long mypos=0;
3518 char *p = tempbuffer;
3519 va_list ap;
3520 //p = (char *) mycalloc(1024,sizeof(char));
3521 va_start(ap, fmt);
3522 mypos = vsprintf(p, fmt, ap);
3523 va_end(ap);
3524 if((*pos + mypos) < (*maxbufsize))
3525 {
3526 (*pos) += sprintf((*buffer) + (*pos), "%s",p);
3527 // if(*pos + mypos >= *maxbufsize )
3528 // {
3529 // printf("%i> pos=%li + mypos=%li > maxbufsize=%li \n",myID, *pos, mypos, *maxbufsize);
3530 // }
3531 }
3532 else
3533 {
3534 *maxbufsize = *pos + 4 * mypos; // add some extra space
3535 (*buffer) = (char *) myrealloc ((*buffer), (*maxbufsize) * sizeof (char));
3536 (*pos) += sprintf((*buffer) + (*pos), "%s",p);
3537 }
3538 // myfree(p);
3539 return (*pos);
3540 }
3541
3542 /// sticks a text into the warning buffer that is printed at the end of the PDF and the end of the TEXT file
record_warnings(world_fmt * world,const char * fmt,...)3543 void record_warnings(world_fmt * world, const char *fmt, ...)
3544 {
3545 long mypos=0;
3546 char *p = (char *) calloc(LINESIZE,sizeof(char));
3547 va_list ap;
3548 if(myID==MASTER && world->cold)
3549 {
3550 va_start(ap, fmt);
3551 mypos = vsprintf(p, fmt, ap);
3552 va_end(ap);
3553 if((world->warningsize + mypos) < world->warningallocsize)
3554 {
3555 world->warningsize += sprintf(world->warning + world->warningsize, "%s\n",p);
3556 }
3557 else
3558 {
3559 world->warningallocsize = world->warningsize + 4 * mypos; // add some extra space
3560 if(world->warning != NULL)
3561 world->warning = (char *) myrealloc (world->warning, world->warningallocsize * sizeof (char));
3562 else
3563 world->warning = (char *) mycalloc (world->warningallocsize, sizeof (char));
3564 world->warningsize += sprintf(world->warning + world->warningsize, "%s\n",p);
3565 }
3566 }
3567 myfree(p);
3568 }
3569
3570
3571
3572
print_warning2(FILE * file,world_fmt * world)3573 void print_warning2(FILE *file, world_fmt *world)
3574 {
3575 char paragraph[] = "This section reports potential problems with your run, but such reporting is often not very accurate. Whith many parameters in a multilocus analysi\
3576 s, it is very common that some parameters for some loci will not be very informative, triggering suggestions (for example to increase the prior ran\
3577 ge) that are not sensible. This suggestion tool will improve with time, therefore do not blindly follow its suggestions. If some parameters are fla\
3578 gged, inspect the tables carefully and judge wether an action is required. For example, if you run a Bayesian inference with sequence data, for mac\
3579 roscopic species there is rarely the need to increase the prior for Theta beyond 0.1; but if you use microsatellites it is rather common that your \
3580 prior distribution for Theta should have a range from 0.0 to 100 or more. With many populations (>3) it is also very common that some migration rou\
3581 tes are estimated poorly because the data contains little or no information for that route. Increasing the range will not help in such situations, \
3582 reducing number of parameters may help in such situations.\0";
3583 long i=0;
3584 long z=0;
3585 char *section;
3586 section = (char*) mycalloc(LINESIZE,sizeof(char));
3587 fprintf(file,"\nPOTENTIAL PROBLEMS\n");
3588 print_line(file,'-',90,10);
3589 while (paragraph[i] != '\0')
3590 {
3591 section[z]=paragraph[i];
3592 section[z+1] = '\0';
3593 z++;
3594 i++;
3595 if(z >= 90)
3596 {
3597 while(section[z-1]!=' ')
3598 {
3599 z--;
3600 i--;
3601 }
3602 section[z]='\0';
3603 fprintf(file,"%s\n",section);
3604 section[0]='\0';
3605 z=0;
3606 }
3607 }
3608 fprintf(file,"%s\n",section);
3609 print_line(file,'-',90,10);
3610 if (world->warning[0]=='\0')
3611 fprintf(file,"No warning was recorded during the run\n\n");
3612 else
3613 fprintf(file,"%s", world->warning);
3614 print_line(file,'-',90,10);
3615 myfree(section);
3616 }
3617
3618
print_stored_warnings(world_fmt * world)3619 void print_stored_warnings(world_fmt *world)
3620 {
3621 if(world->warningsize>0)
3622 {
3623 print_warning2(world->outfile,world);
3624 if(world->options->progress)
3625 {
3626 print_warning2(stdout,world);
3627 }
3628 }
3629 }
3630
3631
3632 // true should continue/shortcut loop
3633 // false should do the rest of the loop, and reset j0 to j
3634 // selects variable to work on, used on Bayesian context
shortcut(long j0,bayes_fmt * bayes,long * j)3635 boolean shortcut(long j0,bayes_fmt *bayes, long *j)
3636 {
3637 if(bayes->map[j0][1] == INVALID)
3638 return TRUE;
3639 else
3640 {
3641 *j = bayes->map[j0][1];
3642 }
3643 if(*j < j0)
3644 {
3645 return TRUE;
3646 }
3647 else
3648 {
3649 return FALSE;
3650 }
3651 }
3652 #endif
3653
3654 // calculates a rough approximation to log(val)
3655 // max error = 0.0007
3656 // Submitted by Laurent de Soras, posted on 30 March 2001
3657 // http://www.flipcode.com/cgi-bin/msg.cgi?showThread=Tip-Fastlogfunction&forum=totd&id=-1
3658 // Fast log() Function, by Laurent de Soras:
3659 // Here is a code snippet to replace the slow log() function...
3660 // It just performs an approximation, but the maximum error is below 0.007.
3661 // Speed gain is about x5, and probably could be increased by tweaking the assembly code.
3662 // The function is based on floating point coding.
3663 // It's easy to get floor (log2(N)) by isolating exponent part.
3664 // We can refine the approximation by using the mantissa. This function returns log2(N)
3665 /*
3666 MYINLINE float fast_log2 (float vval)
3667 {
3668 float vv;
3669 int * const exp_ptr = (int *) (&vval);
3670 int x = *exp_ptr;
3671 const int log_2 = ((x >> 23) & 255) - 128;
3672 x &= ~(255 << 23);
3673 x += 127 << 23;
3674 *exp_ptr = x;
3675
3676 vval = ((-1.0/3) * vval + 2) * vval - 2.0/3; // (1)
3677 vv = vval + log_2;
3678 return (vv);
3679 }
3680
3681 // The line (1) computes 1+log2(m), m ranging from 1 to 2. The proposed
3682 // formula is a 3rd degree polynomial keeping first derivate
3683 // continuity. Higher degree could be used for more accuracy. For faster
3684 // results, one can remove this line, if accuracy is not the matter (it
3685 // gives some linear interpolation between powers of 2).
3686 //Now we got log2(N), we have to multiply it by ln(2) to get the natural log :
3687
3688
3689 MYINLINE float fast_log (float vval)
3690 {
3691 float v = fast_log2 (vval) * 0.69314718f;
3692 //fprintf(stdout,"val= %f fast_log=%f log=%f\n", (float) vval, (float) v, log(vval));
3693 return v;
3694 }
3695
3696 //#define myEXPA (1048576 / 0.693147180559945309417232121458)
3697 #define myEXPA 1512775.39519518569383584038231
3698 #define myEXPC 60801
3699
3700 MYINLINE double fast_exp(double y)
3701 {
3702 union
3703 {
3704 double d;
3705 // weird the implementation says ifdef but that does no work at all
3706 #ifndef LITTLE_ENDIAN
3707 struct { int j, i; } n;
3708 #else
3709 struct { int i, j; } n;
3710 #endif
3711 } eco;
3712 eco.n.i = (int)(myEXPA*(y)) + (1072693248 - myEXPC);
3713 eco.n.j = 0;
3714 return eco.d;
3715 }
3716
3717 */
3718