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