1 /*------------------------------------------------------
2  Maximum likelihood estimation
3  of migration rate  and effectice population size
4  using a Metropolis-Hastings Monte Carlo algorithm
5  -------------------------------------------------------
6  Unique event polymorphism   R O U T I N E S
7 
8  Peter Beerli 2000, Seattle
9  beerli@fsu.edu
10 
11 Copyright 2002 Peter Beerli and Joseph Felsenstein
12 
13  Permission is hereby granted, free of charge, to any person obtaining
14  a copy of this software and associated documentation files (the
15  "Software"), to deal in the Software without restriction, including
16  without limitation the rights to use, copy, modify, merge, publish,
17  distribute, sublicense, and/or sell copies of the Software, and to
18  permit persons to whom the Software is furnished to do so, subject
19  to the following conditions:
20 
21  The above copyright notice and this permission notice shall be
22  included in all copies or substantial portions of the Software.
23 
24  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
25  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
26  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
27  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
28  ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
29  CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
30  WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
31 
32 $Id: uep.c 2158 2013-04-29 01:56:20Z beerli $
33 */
34 /* \file uep.c
35 Unique event polymorphism manipulation
36 */
37 #ifdef UEP
38 #include "migration.h"
39 #include "tools.h"
40 #include "random.h"
41 #include "uep.h"
42 #include "sighandler.h"
43 #include "migrate_mpi.h"
44 
45 #ifdef DMALLOC_DEBUG
46 #include <dmalloc.h>
47 #endif
48 
49 #define AVERAGE    0  /* used for time of uep */
50 #define POPULATION 1  /* used for time of uep given that it raised in */
51 /*   pop x */
52 
53 node *first_uep (node * p, node * nop, long uepsites);
54 boolean uep_check (node * root, int *uep, long uepsites);
55 node *first_uep2 (node * p, node * root, long uepsites);
56 MYREAL sum_ueplike (proposal_fmt * proposal);
57 
58 void print_uep (FILE * out, world_fmt * world, long copies);
59 MYREAL *calc_uep_popprob (world_fmt * world, long copies, long rootstate);
60 MYREAL *calc_uep_timeprob (world_fmt * world, long copies, int type,
61                            long rootstate);
62 long normalize_uep (world_fmt * world);
63 void analyze_uep (world_fmt * world);
64 void print_ancestor (FILE * out, world_fmt * world, long copies);
65 void update_uepanc (world_fmt * world);
66 void check_uep_root (node * p, world_fmt * world);
67 
68 void
setup_uep(world_fmt * world,option_fmt * options)69 setup_uep (world_fmt * world, option_fmt * options)
70 {
71     if (options->uep)
72     {
73         alloc_ueplike (&(world->ueplike), world->data->uepsites,
74                        &(world->ueplikestore), world->options->lsteps,
75                        &(world->ueprootstore), world->numpop, world->loci);
76         alloc_ueptime (&(world->ueptime), world->data->uepsites,
77                        &(world->ueptimestore), world->options->lsteps,
78                        world->numpop, world->loci);
79         world->oldrootuep =
80             (long *) mycalloc (world->data->uepsites, sizeof (long));
81 
82         alloc_uepanc (&(world->uepanc), world->data->uepsites, world->numpop,
83                       world->loci);
84     }
85 }
86 
87 void
destroy_uep(world_fmt * world,long storesize)88 destroy_uep (world_fmt * world, long storesize)
89 {
90     long loci = world->loci;
91     long locus;
92     long i;
93     long j;
94   if (world->options->uep)
95     {
96       myfree(world->uepanc[0]);
97       myfree(world->uepanc);
98       myfree(world->ueplike[0]);
99       myfree(world->ueplike);
100       for (locus = 0; locus < loci; locus++)
101 	{
102 	  for (i = 0; i < storesize; ++i)
103 	    {
104 	      for (j = 0; j < world->options->lsteps; ++j)
105 		{
106 		  myfree(world->ueplikestore[locus][i][j]);
107 		  myfree(world->ueptime[j].populations);
108 		  myfree(world->ueptime[j].ueptime);
109 		}
110 	      myfree(world->ueplikestore[locus][i]);
111 	      myfree(world->ueptimestore[locus][i]);
112 	      myfree(world->ueprootstore[locus][i]);
113 	    }
114 	  myfree(world->ueplikestore[locus]);
115 	  myfree(world->ueptimestore[locus]);
116 	  myfree(world->ueprootstore[locus]);
117 	}
118       myfree(world->ueplikestore);
119       myfree(world->ueptimestore);
120       myfree(world->ueprootstore);
121       for (j = 0; j < world->options->lsteps; ++j)
122 	{
123 	  myfree(world->ueptime[j].populations);
124 	  myfree(world->ueptime[j].ueptime);
125 	}
126       myfree(world->ueptime);
127       myfree(world->oldrootuep);
128     }
129 }
130 
131 void
pseudonu_twostate(proposal_fmt * proposal,ueparray_fmt * xx1,MYREAL * lx1,MYREAL v1,ueparray_fmt * xx2,MYREAL * lx2,MYREAL v2)132 pseudonu_twostate (proposal_fmt * proposal,
133                    ueparray_fmt * xx1, MYREAL *lx1, MYREAL v1,
134                    ueparray_fmt * xx2, MYREAL *lx2, MYREAL v2)
135 {
136     //  static long count=0;
137     long j;
138     MYREAL pija11, pija12, pija21, pija22;
139     //  MYREAL x3m = -MYREAL_MAX;
140     MYREAL r = proposal->world->options->uepmu;
141     MYREAL s = proposal->world->options->uepnu;
142     MYREAL rs = r + s;
143     pair *x1;
144     pair *x2;
145     v1 = EXP (-v1 * rs);
146     v2 = EXP (-v2 * rs);
147     r = r / rs;
148     s = s / rs;
149     for (j = 0; j < proposal->world->data->uepsites; j++)
150     {
151         x1 = &(xx1->s[j]);
152         x2 = &(xx2->s[j]);
153         pija11 = (v1 * r + s) * (*x1)[0] + (s - v1 * s) * (*x1)[1];
154         pija12 = (v2 * r + s) * (*x2)[0] + (s - v2 * s) * (*x2)[1];
155         pija21 = (s - v1 * s) * (*x1)[0] + (v1 * s + r) * (*x1)[1];
156         pija22 = (s - v2 * s) * (*x2)[0] + (v2 * s + r) * (*x2)[1];
157         (*x1)[0] = pija11 * pija12;
158         (*x1)[1] = pija21 * pija22;
159     }
160     //*lx1 += lx2;
161     //count++;
162     //if (count == SCALEINTERVAL)
163     // {
164     // count = 0;
165     //if ((*xx1)[0] > x3m)
166     //      x3m = (*xx1)[0];
167     //      if ((*xx1)[1] > x3m)
168     //      x3m = (*xx1)[1];
169     //      (*xx1)[0] /= x3m;
170     //     (*xx1)[1] /= x3m;
171     //      *lx1 += log (x3m);
172     //    }
173 }
174 
175 
176 void
twostate_nuview(node * mother,world_fmt * world,const long locus)177 twostate_nuview (node * mother, world_fmt * world, const long locus)
178 {
179     //  static long count=0;
180     long j;
181     node *d1 = NULL, *d2 = NULL;
182     MYREAL v1, v2;
183     MYREAL pija11, pija12, pija21, pija22;
184     pair *xx1, *xx2;
185     //  MYREAL x3m = -MYREAL_MAX;
186     pair *xx3;
187     MYREAL r = world->options->uepmu;
188     MYREAL s = world->options->uepnu;
189     MYREAL rs = r + s;
190     children (mother, &d1, &d2);
191     v1 = EXP (-d1->v * rs);
192     v2 = EXP (-d2->v * rs);
193     r = r / rs;
194     s = s / rs;
195     for (j = 0; j < world->data->uepsites; j++)
196     {
197         xx1 = &d1->ux.s[j];
198         xx2 = &d2->ux.s[j];
199         xx3 = &mother->ux.s[j];
200 
201         pija11 = (v1 * r + s) * (*xx1)[0] + (s - v1 * s) * (*xx1)[1];
202         pija12 = (v2 * r + s) * (*xx2)[0] + (s - v2 * s) * (*xx2)[1];
203         (*xx3)[0] = pija11 * pija12;
204         pija21 = (s - v1 * s) * (*xx1)[0] + (v1 * s + r) * (*xx1)[1];
205         pija22 = (s - v2 * s) * (*xx2)[0] + (v2 * s + r) * (*xx2)[1];
206         (*xx3)[0] = pija21 * pija22;
207     }
208     //  count++;
209     //  mother->scale[0] = d1->scale[0] + d2->scale[0];
210     //  if (count == SCALEINTERVAL)
211     //    {
212     //      count = 0;
213     //      if (xx3[0] > x3m)
214     //      x3m = xx3[0];
215     //      if (xx3[1] > x3m)
216     //      x3m = xx3[1];
217     //      xx3[0] /= x3m;
218     //      xx3[1] /= x3m;
219     //      mother->scale[0] += log (x3m);
220     //    }
221 }
222 
223 void
alloc_uepanc(long *** uepanc,long size,long numpop,long loci)224 alloc_uepanc (long ***uepanc, long size, long numpop, long loci)
225 {
226     long locus;
227     (*uepanc) = (long **) mycalloc (loci, sizeof (long *));
228     (*uepanc)[0] = (long *) mycalloc (loci * 2 * size * numpop, sizeof (long));
229     for (locus = 1; locus < loci; locus++)
230       (*uepanc)[locus] = (*uepanc)[0] + 2 * size * numpop;
231 }
232 
233 void
alloc_ueplike(MYREAL *** ueplike,long size,MYREAL ***** ueplikestore,long storesize,long **** ueprootstore,long numpop,long loci)234 alloc_ueplike (MYREAL ***ueplike, long size,
235                MYREAL *****ueplikestore, long storesize,
236                long ****ueprootstore, long numpop, long loci)
237 {
238     long i, j, locus;
239     (*ueplike) = (MYREAL **) mycalloc (size, sizeof (MYREAL *));
240     (*ueplike)[0] = (MYREAL *) mycalloc (size * numpop, sizeof (MYREAL));
241     for (i = 1; i < size; ++i)
242         (*ueplike)[i] = (*ueplike)[0] + i * numpop;
243 
244     (*ueprootstore) = (long ***) mycalloc (loci, sizeof (long **));
245     (*ueplikestore) = (MYREAL ****) mycalloc (loci, sizeof (MYREAL ***));
246 
247     for (locus = 0; locus < loci; locus++)
248     {
249         (*ueprootstore)[locus] = (long **) mycalloc (storesize, sizeof (long *));
250         (*ueplikestore)[locus] =
251             (MYREAL ***) mycalloc (storesize, sizeof (MYREAL **));
252         for (i = 0; i < storesize; ++i)
253         {
254             (*ueprootstore)[locus][i] = (long *) mycalloc (size, sizeof (long));
255             (*ueplikestore)[locus][i] =
256                 (MYREAL **) mycalloc (size, sizeof (MYREAL *));
257             for (j = 0; j < size; ++j)
258                 (*ueplikestore)[locus][i][j] =
259                     (MYREAL *) mycalloc (numpop, sizeof (MYREAL));
260         }
261     }
262 }
263 
264 void
alloc_ueptime(ueptime_fmt ** ueptime,long size,ueptime_fmt **** ueptimestore,long storesize,long numpop,long loci)265 alloc_ueptime (ueptime_fmt ** ueptime, long size,
266                ueptime_fmt **** ueptimestore,
267                long storesize, long numpop, long loci)
268 {
269     long i, j, locus;
270     (*ueptime) = (ueptime_fmt *) mycalloc (size, sizeof (ueptime_fmt));
271     for (j = 0; j < size; ++j)
272     {
273         (*ueptime)[j].size = 3;
274         (*ueptime)[j].populations = (long *) mycalloc (3, sizeof (long));
275         (*ueptime)[j].ueptime = (MYREAL *) mycalloc (3, sizeof (MYREAL));
276     }
277     (*ueptimestore) = (ueptime_fmt ***) mycalloc (loci, sizeof (ueptime_fmt **));
278     for (locus = 0; locus < loci; locus++)
279     {
280         (*ueptimestore)[locus] = (ueptime_fmt **) mycalloc (storesize,
281                                  sizeof (ueptime_fmt
282                                          *));
283         (*ueptimestore)[locus][0] =
284             (ueptime_fmt *) mycalloc (size * storesize, sizeof (ueptime_fmt));
285         for (i = 1; i < storesize; ++i)
286         {
287             (*ueptimestore)[locus][i] = (*ueptimestore)[locus][0] + i * size;
288         }
289         for (i = 0; i < storesize; ++i)
290         {
291             for (j = 0; j < size; ++j)
292             {
293                 (*ueptimestore)[locus][i][j].size = 3;
294                 (*ueptimestore)[locus][i][j].populations =
295                     (long *) mycalloc (3, sizeof (long));
296                 (*ueptimestore)[locus][i][j].ueptime =
297                     (MYREAL *) mycalloc (3, sizeof (MYREAL));
298             }
299         }
300     }
301 }
302 
303 void
allocate_uep(node * p,world_fmt * world,char datatype,boolean withtips)304 allocate_uep (node * p, world_fmt * world, char datatype, boolean withtips)
305 {
306     if (p->type != 't')
307     {
308         if (p->next->back != NULL)
309             allocate_uep (crawlback (p->next), world, datatype, withtips);
310         if (p->next->next->back != NULL)
311             allocate_uep (crawlback (p->next->next), world, datatype, withtips);
312         //      p->ux.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair));
313     }
314     else
315     {
316         if (withtips)
317         {
318             //   p->ux.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair));
319         }
320     }
321 }
322 
323 void
constrain_distance_uep(int ** uep,long uepsites,MYREAL ** distm,long tips)324 constrain_distance_uep (int **uep, long uepsites, MYREAL **distm, long tips)
325 {
326     long i, j;
327     //  long ueppatt = howmany_uep_pattern(uep,uepsites, tips);
328     //  long ueppatt = 2;
329     for (i = 0; i < tips; ++i)
330     {
331         for (j = 0; j < i; ++j)
332         {
333             if (uep[i][0] == uep[j][0])
334             {
335                 distm[i][j] = 0.;
336                 distm[j][i] = 0.;
337             }
338             else
339             {
340                 distm[i][j] = 1.;
341                 distm[j][i] = 1.;
342             }
343         }
344     }
345 }
346 
347 void
set_uepvalues(data_fmt * data,node ** treenode,long tii,long ii)348 set_uepvalues (data_fmt * data, node ** treenode, long tii, long ii)
349 {
350     long j;
351     for (j = 0; j < data->uepsites; j++)
352     {
353         switch (data->uep[ii][j])
354         {
355         case '0':
356             treenode[tii]->ux.s[j][0] = 1.;
357             treenode[tii]->ux.s[j][1] = 0.;
358             break;
359         case '1':
360             treenode[tii]->ux.s[j][0] = 0.;
361             treenode[tii]->ux.s[j][1] = 1.;
362             break;
363         case '?':
364         default:
365             treenode[tii]->ux.s[j][0] = treenode[tii]->ux.s[j][1] = 1.;
366             break;
367         }
368     }
369 }
370 
371 void
make_uep_values(world_fmt * world,data_fmt * data,long locus)372 make_uep_values (world_fmt * world, data_fmt * data, long locus)
373 {
374     long pop, ind, ii = 0, poptips;
375     long halfuep = world->sumtips / 2;
376     long tmp = 0;
377     node **treenode = world->nodep;
378     if (strchr (SEQUENCETYPES, world->options->datatype))
379     {
380         for (ii = 0; ii < world->sumtips; ii++)
381         {
382             set_uepvalues (data, treenode, ii, ii);
383         }
384     }
385     else
386     {
387         for (pop = 0; pop < world->numpop; pop++)
388         {
389             poptips = data->numind[pop][0] / 2;
390             for (ind = 0; ind < poptips; ind++)
391             {
392                 set_uepvalues (data, treenode, tmp + ind, ii);
393                 set_uepvalues (data, treenode, tmp + ind + poptips,
394                                halfuep + ii);
395                 ii++;
396             }
397             tmp += poptips;
398         }
399     }
400 }
401 
402 MYREAL
pseudo_tl_uep(ueparray_fmt * xx1,ueparray_fmt * xx2,MYREAL v1,MYREAL v2,proposal_fmt * proposal,world_fmt * world)403 pseudo_tl_uep (ueparray_fmt * xx1, ueparray_fmt * xx2, MYREAL v1, MYREAL v2,
404                proposal_fmt * proposal, world_fmt * world)
405 {
406     MYREAL summ;
407     long i;
408     MYREAL tterm;
409     pair *x1;
410     summ = 0.0;
411     for (i = 0; i < world->data->uepsites; i++)
412     {
413         x1 = &(xx1->s[i]);
414         tterm =
415             world->options->uepfreq0 * (*x1)[0] + world->options->uepfreq1 * (*x1)[1];
416         summ += (LOG (tterm) + proposal->mf[i]);
417     }
418     return summ;
419 }
420 
421 
422 MYREAL
treelike_uep(world_fmt * world,long locus)423 treelike_uep (world_fmt * world, long locus)
424 {
425     MYREAL tterm;
426     MYREAL summ;
427     long i;
428     //  MYREAL scale;
429     node *p;
430     pair *x1;
431     p = crawlback (world->root->next);
432     summ = 0.0;
433     for (i = 0; i < world->data->uepsites; i++)
434     {
435         x1 = &(p->ux.s[i]);
436         //        scale = p->uepscale[i];
437         tterm =
438             world->options->uepfreq0 * (*x1)[0] + world->options->uepfreq1 * (*x1)[1];
439         summ += (LOG (tterm)); //+ scale);
440     }
441     return summ;
442 }    /* treelikelihood */
443 
444 
445 // uses the conditional likelihoods to set
446 // the actual state at a node
447 void
update_uep(node * p,world_fmt * world)448 update_uep (node * p, world_fmt * world)
449 {
450     long i;
451     if (p->type != 't')
452     {
453         if (p->next->back != NULL)
454             update_uep (crawlback (p->next), world);
455         if (p->next->next->back != NULL)
456             update_uep (crawlback (p->next->next), world);
457 
458         if (p->type != 'r')
459         {
460             for (i = 0; i < world->data->uepsites; ++i)
461                 p->uep[i] = (p->ux.s[i][0] < p->ux.s[i][1]) ? '1' : '0';
462         }
463     }
464 }
465 
466 void
check_uep_root(node * p,world_fmt * world)467 check_uep_root (node * p, world_fmt * world)
468 {
469     long i;
470     for (i = 0; i < world->data->uepsites; ++i)
471     {
472         if (p->uep[i] > 1)
473         {
474             p->uep[i] = (char) RANDINT (0, 1);
475         }
476         world->oldrootuep[i] = p->uep[i];
477     }
478     //  memcpy(world->oldrootuep,p->uep,sizeof(long)*world->data->uepsites);
479 }
480 
481 boolean
is_success_pseudo_uep(proposal_fmt * proposal)482 is_success_pseudo_uep (proposal_fmt * proposal)
483 {
484     boolean o, ob, t, tb;
485     long uepsites = proposal->world->data->uepsites;
486     //if(!proposal->world->in_last_chain && proposal->world->options->uep_last))
487     //  return TRUE;
488     node *root = proposal->world->root->next->back;
489 
490     if (proposal->world->options->uep)
491     {
492         o = uep_check (root, proposal->origin->uep, uepsites);
493         ob = uep_check (root, proposal->oback->uep, uepsites);
494         t = uep_check (root, proposal->target->uep, uepsites);
495         tb =
496             uep_check (root, showtop (crawlback (proposal->target))->uep,
497                        uepsites);
498         if ((o && t && ob && tb) || (!o && !t && !ob && !tb))
499             return TRUE;
500         if (o && ob && t && !tb)
501             return TRUE;
502         if (o && !ob && t && tb)
503             return TRUE;
504         if (o && !ob && !t && !tb)
505             return TRUE;
506         if (!o && !ob && t && !tb)
507             return TRUE;
508     }
509     return FALSE;
510 }
511 
512 boolean
is_success_pseudo_uepOLD(proposal_fmt * proposal)513 is_success_pseudo_uepOLD (proposal_fmt * proposal)
514 {
515     boolean o, ob, t;
516     //if(!proposal->world->in_last_chain && proposal->world->options->uep_last))
517     //  return TRUE;
518     node *root = proposal->world->root->next->back;
519     if (proposal->world->options->uep)
520     {
521         o =
522             uep_check (root, proposal->origin->uep,
523                        proposal->world->data->uepsites);
524         ob =
525             uep_check (root, proposal->oback->uep,
526                        proposal->world->data->uepsites);
527         t =
528             uep_check (root, proposal->target->uep,
529                        proposal->world->data->uepsites);
530         if (o && t)
531             return TRUE;
532         else
533         {
534             if (o && !t && !ob)
535             {
536                 return TRUE;
537             }
538             else
539             {
540                 if (!o && !t && !ob)
541                     return TRUE;
542             }
543         }
544     }
545     return FALSE;
546 }
547 
548 boolean
uep_check(node * root,int * uep,long uepsites)549 uep_check (node * root, int *uep, long uepsites)
550 {
551     long i;
552     for (i = 0; i < uepsites; ++i)
553         if (uep[i] == root->uep[i])
554             return FALSE;
555     return TRUE;
556 }
557 
558 boolean
uep_compare(int * target,int * origin,long uepsites)559 uep_compare (int *target, int *origin, long uepsites)
560 {
561     long i;
562 
563     for (i = 0; i < uepsites; ++i)
564         if (target[i] != origin[i])
565             return FALSE;
566     return TRUE;
567 }
568 
569 
570 void
store_uep(world_fmt * world)571 store_uep (world_fmt * world)
572 {
573     long i;
574     if (world->options->uep)
575     {
576         if (world->in_last_chain)
577         {
578             for (i = 0; i < world->data->uepsites; ++i)
579             {
580                 memcpy (world->ueplikestore[world->locus][world->G][i],
581                         world->ueplike[i], sizeof (MYREAL) * world->numpop);
582                 memcpy (world->ueptimestore[world->locus][world->G][i].
583                         populations, world->ueptime[i].populations,
584                         sizeof (long) * world->ueptime[i].size);
585                 memcpy (world->ueptimestore[world->locus][world->G][i].ueptime,
586                         world->ueptime[i].ueptime,
587                         sizeof (MYREAL) * world->ueptime[i].size);
588                 world->ueptimestore[world->locus][world->G][i].size =
589                     world->ueptime[i].size;
590             }
591         }
592     }
593 }
594 
595 void
print_uep(FILE * out,world_fmt * world,long copies)596 print_uep (FILE * out, world_fmt * world, long copies)
597 {
598     MYREAL *prob, *tyme, *timep;
599     long pop;
600     FPRINTF (out, "\n\n\nUnique event polymorphism\n");
601     FPRINTF (out, "=========================\n\n");
602     FPRINTF (out, "Probabilities and time of UEP\n");
603     FPRINTF (out, "Ancestral State of UEP is either 0 or 1\n");
604     FPRINTF (out, "-------------------------------------------------\n");
605     FPRINTF (out, "Populations  Probability   Time(*)       Time(**)\n");
606     FPRINTF (out, "-------------------------------------------------\n");
607     prob = calc_uep_popprob (world, copies, -1); // don't care if anc uep=1/0
608     tyme = calc_uep_timeprob (world, copies, AVERAGE, -1);
609     timep = calc_uep_timeprob (world, copies, POPULATION, -1);
610     for (pop = 0; pop < world->numpop; ++pop)
611     {
612         FPRINTF (out, "  %3li        %2.6f      %2.6f      %2.6f\n",
613                  pop + 1, prob[pop], tyme[pop], timep[pop]);
614     }
615     FPRINTF (out, "-------------------------------------------------\n");
616     FPRINTF (out, "(*)  Assuming that mutation arose on branch to UEP\n");
617     FPRINTF (out, "(**) Mutation arose in specific population.\n");
618 
619     FPRINTF (out, "\nAncestral State of UEP is 0\n");
620     FPRINTF (out, "-------------------------------------------------\n");
621     FPRINTF (out, "Populations  Probability   Time(*)       Time(**)\n");
622     FPRINTF (out, "-------------------------------------------------\n");
623     prob = calc_uep_popprob (world, copies, 1); // don't care if anc uep=1
624     tyme = calc_uep_timeprob (world, copies, AVERAGE, 1);
625     timep = calc_uep_timeprob (world, copies, POPULATION, 1);
626     for (pop = 0; pop < world->numpop; ++pop)
627     {
628         FPRINTF (out, "  %3li        %2.6f      %2.6f      %2.6f\n",
629                  pop + 1, prob[pop], tyme[pop], timep[pop]);
630     }
631     FPRINTF (out, "-------------------------------------------------\n");
632     FPRINTF (out, "Probability=Prob(is in population | uep_anc=0)\n");
633     FPRINTF (out, "(*), (**) see above\n");
634 
635     FPRINTF (out, "\nAncestral State of UEP is 1\n");
636     FPRINTF (out, "-------------------------------------------------\n");
637     FPRINTF (out, "Populations  Probability   Time(*)       Time(**)\n");
638     FPRINTF (out, "-------------------------------------------------\n");
639     prob = calc_uep_popprob (world, copies, 0); // don't care if anc uep=0
640     tyme = calc_uep_timeprob (world, copies, AVERAGE, 0);
641     timep = calc_uep_timeprob (world, copies, POPULATION, 0);
642     for (pop = 0; pop < world->numpop; ++pop)
643     {
644         FPRINTF (out, "  %3li        %2.6f      %2.6f      %2.6f\n",
645                  pop + 1, prob[pop], tyme[pop], timep[pop]);
646     }
647     FPRINTF (out, "-------------------------------------------------\n");
648     FPRINTF (out, "Probability=Prob(is in population | uep_anc=1)\n");
649     FPRINTF (out, "(*), (**) see above\n");
650     myfree(prob);
651     myfree(tyme);
652     myfree(timep);
653 }
654 
655 MYREAL *
calc_uep_popprob(world_fmt * world,long copies,long rootstate)656 calc_uep_popprob (world_fmt * world, long copies, long rootstate)
657 {
658     long tree, pop, j, locus;
659     MYREAL *result;
660     result = (MYREAL *) mycalloc (world->numpop, sizeof (MYREAL));
661     for (locus = 0; locus < world->loci; locus++)
662     {
663         for (tree = 0; tree < world->atl[world->repstop - 1][locus].T - 1;
664                 tree++)
665         {
666             for (pop = 0; pop < world->numpop; pop++)
667             {
668                 for (j = 0; j < world->data->uepsites; ++j)
669                 {
670                     if (rootstate != world->ueprootstore[locus][tree][j])
671                         result[pop] += world->ueplikestore[locus][tree][j][pop];
672                 }
673             }
674         }
675     }
676     for (pop = 0; pop < world->numpop; pop++)
677     {
678         result[pop] /= copies;
679     }
680     return result;
681 }
682 
683 MYREAL
ueptime_report(ueptime_fmt * ueptime,int type,long pop)684 ueptime_report (ueptime_fmt * ueptime, int type, long pop)
685 {
686     long i;
687     MYREAL summ = 0.0;
688     long elem = 0;
689 
690     switch (type)
691     {
692     case AVERAGE:
693         return (ueptime->ueptime[ueptime->size - 1] + ueptime->ueptime[0]) / 2.;
694     case POPULATION:
695         summ = 0.0;
696         elem = 0;
697         for (i = 0; i < ueptime->size - 1; ++i)
698         {
699             if (ueptime->populations[i] == pop)
700             {
701                 summ += ueptime->ueptime[i] + ueptime->ueptime[i + 1];
702                 elem += 2;
703             }
704         }
705         if (elem > 0)
706             return summ / (MYREAL) elem;
707         else
708             return 0.0;
709         break;
710     }
711     return 0.0;
712 }
713 
714 MYREAL *
calc_uep_timeprob(world_fmt * world,long copies,int type,long rootstate)715 calc_uep_timeprob (world_fmt * world, long copies, int type, long rootstate)
716 {
717     long pop, tree, j, locus;
718     MYREAL *result;
719     result = (MYREAL *) mycalloc (world->numpop, sizeof (MYREAL));
720     for (locus = 0; locus < world->loci; locus++)
721     {
722         for (tree = 0; tree < world->atl[world->repstop - 1][locus].T - 1;
723                 tree++)
724         {
725             for (pop = 0; pop < world->numpop; pop++)
726             {
727                 for (j = 0; j < world->data->uepsites; ++j)
728                 {
729                     if (rootstate != world->ueprootstore[locus][tree][j])
730                         result[pop] += world->ueplikestore[locus][tree][j][pop] *
731                                        ueptime_report (&world->ueptimestore[locus][tree][j],
732                                                        type, pop);
733                 }
734             }
735         }
736     }
737     for (pop = 0; pop < world->numpop; pop++)
738     {
739         result[pop] /= copies;
740     }
741     return result;
742 }
743 
744 long
normalize_uep(world_fmt * world)745 normalize_uep (world_fmt * world)
746 {
747     long tree, pop, j, locus;
748     MYREAL summ;
749     long copies = 0;
750     for (locus = 0; locus < world->loci; locus++)
751     {
752         for (tree = 0; tree < world->atl[world->repstop - 1][locus].T - 1;
753                 tree++)
754         {
755             summ = 0.0;
756             copies += world->atl[world->repstop - 1][locus].tl[tree].copies;
757             for (pop = 0; pop < world->numpop; pop++)
758             {
759                 for (j = 0; j < world->data->uepsites; ++j)
760                     summ += world->ueplikestore[locus][tree][j][pop];
761             }
762             for (pop = 0; pop < world->numpop; pop++)
763             {
764                 for (j = 0; j < world->data->uepsites; ++j)
765                     world->ueplikestore[locus][tree][j][pop] *=
766                         ((MYREAL) world->atl[world->repstop - 1][locus].tl[tree].
767                          copies) / summ;
768             }
769         }
770     }
771     return copies;
772 }
773 
774 void
analyze_uep(world_fmt * world)775 analyze_uep (world_fmt * world)
776 {
777     long copies;
778     copies = normalize_uep (world);
779     print_uep (world->outfile, world, copies);
780     print_ancestor (world->outfile, world, copies);
781 }
782 
783 void
show_uep_store(world_fmt * world)784 show_uep_store (world_fmt * world)
785 {
786     long i, j, z, steps, locus;
787     FILE *file;
788     char filename[20] = "uepout";
789     openfile (&file, filename, "w+",  NULL);
790     for (locus = 0; locus < world->loci; locus++)
791     {
792         for (steps = 0; steps < world->G; ++steps)
793         {
794             for (i = 0; i < world->numpop; ++i)
795             {
796                 for (j = 0; j < world->data->uepsites; ++j)
797                     FPRINTF (file, "%f ",
798                              world->ueplikestore[locus][steps][j][i]);
799             }
800             FPRINTF (file, "\n");
801         }
802         for (steps = 0; steps < world->G; ++steps)
803         {
804             for (j = 0; j < world->data->uepsites; ++j)
805             {
806                 FPRINTF (file, "%li ",
807                          world->ueptimestore[locus][steps][j].size);
808                 for (z = 0; z < world->ueptimestore[locus][steps][j].size; ++z)
809                 {
810                     FPRINTF (file, "%li ",
811                              world->ueptimestore[locus][steps][j].
812                              populations[z]);
813                 }
814                 FPRINTF (file, "\n");
815                 FPRINTF (file, "%li ",
816                          world->ueptimestore[locus][steps][j].size);
817                 for (z = 0; z < world->ueptimestore[locus][steps][j].size; ++z)
818                 {
819                     FPRINTF (file, "%f ",
820                              world->ueptimestore[locus][steps][j].ueptime[z]);
821                 }
822                 FPRINTF (file, "\n");
823             }
824         }
825     }
826 }
827 
828 void
print_ancestor(FILE * out,world_fmt * world,long copies)829 print_ancestor (FILE * out, world_fmt * world, long copies)
830 {
831     char uepsite[10];
832     long us, pop, half;
833     MYREAL subtotal = 0.0;
834     MYREAL total = 0;
835     long locus;
836     MYREAL ancsum = 0.;
837     for (locus = 0; locus < world->loci; locus++)
838     {
839         for (us = 0; us < 2 * world->numpop * world->data->uepsites; us++)
840             total += world->uepanc[locus][us];
841     }
842     FPRINTF (out, "\n\nUEP probabilities at the MRCA\n");
843     FPRINTF (out, "--------------------------------------\n");
844     FPRINTF (out, "UEP allele    Population   Probability\n");
845     FPRINTF (out, "--------------------------------------\n");
846     half = world->numpop * world->data->uepsites;
847     for (us = 0; us < world->data->uepsites; us++)
848     {
849         if (world->data->uepsites != 1)
850             sprintf (uepsite, "%3li:", us + 1);
851         else
852             sprintf (uepsite, "%4s", "     ");
853         for (pop = 0; pop < world->numpop; pop++)
854         {
855             ancsum = 0;
856             for (locus = 0; locus < world->loci; locus++)
857             {
858                 subtotal += world->uepanc[locus][pop * us + pop];
859                 ancsum += world->uepanc[locus][pop * us + pop];
860             }
861             FPRINTF (out, "%4s 0            %3li     %f\n", uepsite, pop + 1,
862                      ancsum / total);
863         }
864         FPRINTF (out, "%4s 0            All     %f\n", uepsite,
865                  subtotal / total);
866         FPRINTF (out, "--------------------------------------\n");
867         subtotal = 0.;
868         for (pop = 0; pop < world->numpop; pop++)
869         {
870             ancsum = 0;
871             for (locus = 0; locus < world->loci; locus++)
872             {
873                 subtotal += world->uepanc[locus][half + pop * us + pop];
874                 ancsum += world->uepanc[locus][half + pop * us + pop];
875             }
876             FPRINTF (out, "%4s 1            %3li     %f\n", uepsite, pop + 1,
877                      ancsum / total);
878         }
879         FPRINTF (out, "%4s 1            All     %f\n", uepsite,
880                  subtotal / total);
881         FPRINTF (out, "--------------------------------------\n");
882     }
883 }
884 
885 void
calc_ueplike(node * p,world_fmt * world,MYREAL ** ueplike)886 calc_ueplike (node * p, world_fmt * world, MYREAL **ueplike)
887 {
888     boolean done;
889     int temp;
890     long i, save_i = 0;
891     node *d1, *d2;
892     if (p->type != 't')
893     {
894         if (p->next->back != NULL)
895             calc_ueplike (crawlback (p->next), world, ueplike);
896         if (p->next->next->back != NULL)
897             calc_ueplike (crawlback (p->next->next), world, ueplike);
898 
899         if (p->type != 'r')
900         {
901             d1 = crawlback (p->next);
902             d2 = crawlback (p->next->next);
903             done = FALSE;
904             for (i = 0; i < world->data->uepsites; ++i)
905             {
906                 temp = (d1->uep[i]=='0'  ? 0 : 1 )+ (d2->uep[i]=='0' ? 0 : 1);
907                 if (temp == 1)
908                 {
909                     if (!done) //shortcut as long we have only one UEP
910                     {
911 		      if(world->options->verbose)
912                         printf ("CALCUEP: at <%li> with time %f (%f,%f)\n",p->id,p->tyme,d1->tyme,d2->tyme);
913                         collect_ueplike (p, d1, d2, i, world, ueplike[i]);
914                         save_i = i;
915                     }
916                     else
917                         memcpy (ueplike[i], ueplike[save_i],
918                                 sizeof (MYREAL) * world->numpop);
919                     done = TRUE;
920                 }
921             }
922         }
923     }
924 }
925 
926 void
fill_ueptime(ueptime_fmt * ueptime,node * p,node * last)927 fill_ueptime (ueptime_fmt * ueptime, node * p, node * last)
928 {
929     long i = 1;
930     long count = 0;
931     node *nod;
932 
933     ueptime->populations[0] = last->actualpop;
934     ueptime->ueptime[0] = last->tyme;
935     while ((nod = showtop (last->back)) != p)
936     {
937         last = nod;
938         count++;
939     }
940     count += 2;
941     if (count > ueptime->size)
942     {
943         ueptime->populations = (long *) myrealloc (ueptime->populations,
944                                sizeof (long) * count);
945         ueptime->ueptime = (MYREAL *) myrealloc (ueptime->ueptime,
946                                                sizeof (MYREAL) * count);
947         ueptime->size = count;
948     }
949     while ((nod = showtop (last->back)) != p)
950     {
951         ueptime->populations[i] = nod->actualpop;
952         ueptime->ueptime[i] = nod->tyme;
953         last = nod;
954         i++;
955     }
956     ueptime->populations[i] = nod->actualpop;
957     ueptime->ueptime[i] = nod->tyme;
958     ueptime->size = i + 1;
959 }
960 
961 void
collect_ueplike(node * p,node * d1,node * d2,long uepsite,world_fmt * world,MYREAL * ueplike)962 collect_ueplike (node * p, node * d1, node * d2, long uepsite,
963                  world_fmt * world, MYREAL *ueplike)
964 {
965     node *nod, *last;
966     //  long pop;
967     //  MYREAL interval = world->root->next->back->tyme - p->tyme;
968     //  printf("collect_ueplike: %f %f %f\n",world->root->next->back->tyme,p->tyme,interval);
969     memset (ueplike, 0, sizeof (MYREAL) * world->numpop);
970     //  if (d1->uep[uepsite] == 1)
971     if (d1->uep[uepsite] != world->root->next->back->uep[uepsite])
972         last = d1;
973     else
974         last = d2;
975     fill_ueptime (&world->ueptime[uepsite], p, last);
976     // = (p->tyme + last->tyme) / 2.;
977     while ((nod = showtop (last->back)) != p)
978     {
979         //      printf("real: %li> %f %f\n",nod->id,nod->tyme, nod->tyme - last->tyme);
980         ueplike[nod->actualpop] += nod->tyme - last->tyme;
981         last = nod;
982     }
983     ueplike[nod->actualpop] += nod->tyme - last->tyme;
984     //for(pop=0;pop<world->numpop;++pop)
985     //{
986     //  if(ueplike[pop]!=0.)
987     //    ueplike[pop] *= EXP (world->treelen);
988     //}
989     //      FPRINTF(stdout,"%f %f %f %li - ",ueplike[pop], ueplike[pop]-interval,
990     //      nod->tyme, nod->id);
991     //  }
992     // FPRINTF(stdout,"real\n");
993 }
994 
995 
996 node *
first_uep(node * p,node * nop,long uepsites)997 first_uep (node * p, node * nop, long uepsites)
998 {
999     node *bt = showtop (crawlback (p));
1000     while (uep_compare (p->uep, bt->uep, uepsites) && bt != nop)
1001     {
1002         p = bt;
1003         bt = showtop (crawlback (p));
1004     }
1005     return p;
1006 }
1007 
1008 node *
first_uep2(node * p,node * root,long uepsites)1009 first_uep2 (node * p, node * root, long uepsites)
1010 {
1011     node *pn, *pnn;
1012     if (p->type != 'r' && p->type != 't')
1013     {
1014         pn = p->next;
1015         pnn = p->next->next;
1016         p = first_uep2 (crawlback (pn), root, uepsites);
1017         if (uep_check (root, p->uep, uepsites))
1018             return p;
1019         p = first_uep2 (crawlback (pnn), root, uepsites);
1020         if (uep_check (root, p->uep, uepsites))
1021             return p;
1022     }
1023     return p;
1024 }
1025 
1026 MYREAL
sum_ueplike(proposal_fmt * proposal)1027 sum_ueplike (proposal_fmt * proposal)
1028 {
1029     long u, i;
1030     MYREAL like = 1.;
1031     MYREAL poplike;
1032 
1033     for (u = 0; u < proposal->world->data->uepsites; ++u)
1034     {
1035         poplike = 0.;
1036         for (i = 0; i < proposal->world->numpop; ++i)
1037         {
1038             if (proposal->ueplike[u][i] != 0.0)
1039                 poplike += proposal->ueplike[u][i];
1040         }
1041         like *= poplike;
1042     }
1043     if (proposal->world->options->ueprate > 0.0)
1044     {
1045         return LOG (like) -
1046                (proposal->world->options->ueprate * proposal->treelen);
1047     }
1048     else
1049         return LOG (like);
1050 }
1051 
1052 
1053 MYREAL
adjust_uep_base(proposal_fmt * proposal,MYREAL interval,MYREAL oldinterval)1054 adjust_uep_base (proposal_fmt * proposal, MYREAL interval, MYREAL oldinterval)
1055 {
1056     long u, pop;
1057     // MYREAL expT = EXP (-(proposal->treelen-proposal->world->treelen));
1058     for (u = 0; u < proposal->world->data->uepsites; ++u)
1059     {
1060         memcpy (proposal->ueplike[u], proposal->world->ueplike[u],
1061                 sizeof (MYREAL) * proposal->world->numpop);
1062         //                  printf("adjust_base   ");
1063         for (pop = 0; pop < proposal->world->numpop; ++pop)
1064         {
1065             if (proposal->world->ueplike[u][pop] != 0.0)
1066                 proposal->ueplike[u][pop] -= oldinterval - interval;
1067             //      proposal->ueplike[u][pop] *= expT;
1068             //                      printf("%f ", proposal->ueplike[u][pop]);
1069         }
1070         //            printf("- %f %f\n",interval, oldinterval);
1071     }
1072     return sum_ueplike (proposal);
1073 }
1074 
1075 MYREAL
adjust_uep_target(node * first,node * firstb,proposal_fmt * proposal,MYREAL interval,MYREAL oldinterval)1076 adjust_uep_target (node * first, node * firstb, proposal_fmt * proposal,
1077                    MYREAL interval, MYREAL oldinterval)
1078 {
1079     long u, pop;
1080     node *last, *p;
1081     MYREAL lasttime, endtime;
1082     //MYREAL expT = EXP (-(proposal->treelen-proposal->world->treelen));
1083     boolean o =
1084         uep_check (proposal->world->root->next->back, proposal->origin->uep,
1085                    proposal->world->data->uepsites);
1086     for (u = 0; u < proposal->world->data->uepsites; ++u)
1087     {
1088         memset (proposal->ueplike[u], 0,
1089                 sizeof (MYREAL) * proposal->world->numpop);
1090         last = first;
1091         p = showtop (first->back);
1092         if (o)
1093         {
1094             lasttime = proposal->time;
1095             endtime = firstb->tyme;
1096         }
1097         else
1098         {
1099             lasttime = last->tyme;
1100             endtime = proposal->time;
1101         }
1102         while (p->tyme < lasttime)
1103         {
1104             last = p;
1105             p = showtop (last->back);
1106         }
1107         while (p->tyme <= endtime)
1108         {
1109             proposal->ueplike[u][p->actualpop] += p->tyme - lasttime;
1110             last = p;
1111             lasttime = p->tyme;
1112             p = showtop (last->back);
1113         }
1114         //                  printf("adjust_target ");
1115         for (pop = 0; pop < proposal->world->numpop; ++pop)
1116         {
1117             if (proposal->world->ueplike[u][pop] != 0.0)
1118                 proposal->ueplike[u][pop] -= oldinterval - interval;
1119             //  proposal->ueplike[u][pop] *= expT;
1120             //               printf("%f ", proposal->ueplike[u][pop]);
1121         }
1122         //                 printf("- %f %f\n",interval, oldinterval);
1123     }
1124     return sum_ueplike (proposal);
1125 }
1126 
1127 MYREAL
adjust_uep_origin(node * first,node * firstb,proposal_fmt * proposal,MYREAL interval,MYREAL oldinterval)1128 adjust_uep_origin (node * first, node * firstb, proposal_fmt * proposal,
1129                    MYREAL interval, MYREAL oldinterval)
1130 {
1131     long u, i, pop;
1132     MYREAL lasttime;
1133     //  node *last, *p;
1134     migr_table_fmt *mt = proposal->migr_table;
1135     long mtc = proposal->migr_table_counter;
1136     //MYREAL expT = EXP (-(proposal->treelen-proposal->world->treelen));
1137     //  printf("@@@@@@@ adjust branch below origin @@@@@@@@@\n");
1138     for (u = 0; u < proposal->world->data->uepsites; ++u)
1139     {
1140         // memcpy(proposal->ueplike[u], proposal->world->ueplike[u],
1141         //             sizeof(MYREAL)*proposal->world->numpop);
1142         memset (proposal->ueplike[u], 0,
1143                 sizeof (MYREAL) * proposal->world->numpop);
1144         lasttime = first->tyme;
1145         for (i = 0; i < mtc; ++i)
1146         {
1147             proposal->ueplike[u][mt[i].to] += mt[i].time - lasttime;
1148             lasttime = mt[i].time;
1149         }
1150         proposal->ueplike[u][proposal->target->actualpop] +=
1151             proposal->time - lasttime;
1152         //                  printf("adjust_target ");
1153         for (pop = 0; pop < proposal->world->numpop; ++pop)
1154         {
1155             if (proposal->world->ueplike[u][pop] != 0.0)
1156                 proposal->ueplike[u][pop] -= oldinterval - interval;
1157             //  proposal->ueplike[u][pop] *= expT;
1158             //               printf("%f ", proposal->ueplike[u][pop]);
1159         }
1160         //        printf("- %f %f\n",interval, oldinterval);
1161     }
1162     return sum_ueplike (proposal);
1163 }
1164 
1165 MYREAL
adjust_uep_oback(node * first,node * firstb,proposal_fmt * proposal,MYREAL interval,MYREAL oldinterval)1166 adjust_uep_oback (node * first, node * firstb, proposal_fmt * proposal,
1167                   MYREAL interval, MYREAL oldinterval)
1168 {
1169     long u, pop;
1170     MYREAL lasttime = 0.;
1171     node *last, *p;
1172     //MYREAL expT = EXP (-(proposal->treelen-proposal->world->treelen));
1173     if (first->tyme < proposal->time) // target == osister
1174         first = proposal->oback;
1175     for (u = 0; u < proposal->world->data->uepsites; ++u)
1176     {
1177         //memcpy(proposal->ueplike[u], proposal->world->ueplike[u],
1178         //             sizeof(MYREAL)*proposal->world->numpop);
1179         memset (proposal->ueplike[u], 0,
1180                 sizeof (MYREAL) * proposal->world->numpop);
1181         last = first;
1182         p = showtop (first->back);
1183         if (last != proposal->oback)
1184             lasttime = last->tyme;
1185         else
1186             lasttime = proposal->time;
1187         while (p != firstb)
1188         {
1189             proposal->ueplike[u][p->actualpop] += p->tyme - lasttime;
1190             last = p;
1191             p = showtop (last->back);
1192             if (p == proposal->oback)
1193                 p = showtop (p->back);
1194             lasttime = last->tyme;
1195         }
1196         proposal->ueplike[u][p->actualpop] += p->tyme - lasttime;
1197         //       printf(" adjust_oback ");
1198         for (pop = 0; pop < proposal->world->numpop; ++pop)
1199         {
1200             if (proposal->world->ueplike[u][pop] != 0.0)
1201                 proposal->ueplike[u][pop] -= oldinterval - interval;
1202             //  proposal->ueplike[u][pop] *= expT;
1203             //                printf("%f ", proposal->ueplike[u][pop]);
1204         }
1205         //                  printf("- %f %f\n",interval, oldinterval);
1206     }
1207     return sum_ueplike (proposal);
1208 }
1209 
1210 MYREAL
adjust_uep_osister(node * first,node * firstb,proposal_fmt * proposal,MYREAL interval,MYREAL oldinterval)1211 adjust_uep_osister (node * first, node * firstb, proposal_fmt * proposal,
1212                     MYREAL interval, MYREAL oldinterval)
1213 {
1214     long u, pop;
1215     MYREAL lasttime = 0.;
1216     node *last, *p;
1217     //MYREAL expT = EXP (-(proposal->treelen-proposal->world->treelen));
1218     for (u = 0; u < proposal->world->data->uepsites; ++u)
1219     {
1220         memset (proposal->ueplike[u], 0,
1221                 sizeof (MYREAL) * proposal->world->numpop);
1222         last = first;
1223         p = showtop (first->back);
1224         lasttime = last->tyme;
1225         while (p != firstb)
1226         {
1227             proposal->ueplike[u][p->actualpop] += p->tyme - lasttime;
1228             last = p;
1229             p = showtop (last->back);
1230             lasttime = last->tyme;
1231         }
1232         proposal->ueplike[u][p->actualpop] += p->tyme - lasttime;
1233         //            printf("adjust_osister ");
1234         for (pop = 0; pop < proposal->world->numpop; ++pop)
1235         {
1236             if (proposal->world->ueplike[u][pop] != 0.0)
1237                 proposal->ueplike[u][pop] -= oldinterval - interval;
1238             //  proposal->ueplike[u][pop] *= expT;
1239             //        printf("%f ", proposal->ueplike[u][pop]);
1240         }
1241         // printf("- %f %f\n",interval, oldinterval);
1242     }
1243     return sum_ueplike (proposal);
1244 }
1245 
1246 MYREAL
adjust_uep_osistertarget(node * first,node * firstb,proposal_fmt * proposal,MYREAL interval,MYREAL oldinterval)1247 adjust_uep_osistertarget (node * first, node * firstb,
1248                           proposal_fmt * proposal, MYREAL interval,
1249                           MYREAL oldinterval)
1250 {
1251     long u, pop;
1252     MYREAL lasttime = 0.;
1253     node *last, *p;
1254     //MYREAL expT = EXP (-(proposal->treelen-proposal->world->treelen));
1255     for (u = 0; u < proposal->world->data->uepsites; ++u)
1256     {
1257         memset (proposal->ueplike[u], 0,
1258                 sizeof (MYREAL) * proposal->world->numpop);
1259         last = first;
1260         p = showtop (first->back);
1261         lasttime = last->tyme;
1262         while (p != firstb)
1263         {
1264             proposal->ueplike[u][p->actualpop] += p->tyme - lasttime;
1265             last = p;
1266             p = showtop (last->back);
1267             if (last != proposal->oback)
1268                 lasttime = last->tyme;
1269             else
1270                 lasttime = proposal->time;
1271         }
1272         proposal->ueplike[u][p->actualpop] += p->tyme - lasttime;
1273         printf ("adjust_osistertarget ");
1274         for (pop = 0; pop < proposal->world->numpop; ++pop)
1275         {
1276             if (proposal->world->ueplike[u][pop] != 0.0)
1277                 proposal->ueplike[u][pop] -= oldinterval - interval;
1278             //  proposal->ueplike[u][pop] *= expT;
1279             printf ("%f ", proposal->ueplike[u][pop]);
1280         }
1281         printf ("- %f %f\n", interval, oldinterval);
1282     }
1283     return sum_ueplike (proposal);
1284 }
1285 
1286 MYREAL
adjust_uep_osistertarget2(node * first,node * firstb,proposal_fmt * proposal,MYREAL interval,MYREAL oldinterval)1287 adjust_uep_osistertarget2 (node * first, node * firstb,
1288                            proposal_fmt * proposal, MYREAL interval,
1289                            MYREAL oldinterval)
1290 {
1291     long u, pop, actualpop;
1292     MYREAL lasttime = 0., endtime;
1293     node *last, *p;
1294     // MYREAL expT = EXP (-(proposal->treelen-proposal->world->treelen));
1295     for (u = 0; u < proposal->world->data->uepsites; ++u)
1296     {
1297         memset (proposal->ueplike[u], 0,
1298                 sizeof (MYREAL) * proposal->world->numpop);
1299         last = first;
1300         p = showtop (first->back);
1301         actualpop = p->actualpop;
1302         lasttime = last->tyme;
1303         endtime = proposal->time;
1304         while (p != proposal->realtarget)
1305         {
1306             proposal->ueplike[u][p->actualpop] += p->tyme - last->tyme;
1307             last = p;
1308             p = showtop (last->back);
1309         }
1310         proposal->ueplike[u][p->actualpop] += p->tyme - last->tyme;
1311         proposal->ueplike[u][p->pop] += endtime - p->tyme;
1312         //                  printf("adjust_osistertarget ");
1313         for (pop = 0; pop < proposal->world->numpop; ++pop)
1314         {
1315             if (proposal->world->ueplike[u][pop] != 0.0)
1316                 proposal->ueplike[u][pop] -= oldinterval - interval;
1317             //  proposal->ueplike[u][pop] *= expT;
1318             //                printf("%f ", proposal->ueplike[u][pop]);
1319         }
1320         //      printf("- %f %f\n",interval, oldinterval);
1321     }
1322     return sum_ueplike (proposal);
1323 }
1324 
1325 MYREAL
pseudo_ueplikelihood(world_fmt * world,proposal_fmt * proposal)1326 pseudo_ueplikelihood (world_fmt * world, proposal_fmt * proposal)
1327 {
1328     MYREAL interval = 0., oldinterval = 0.; //  node *first, *firstb;
1329     node *root = proposal->root->next->back;
1330     long uepsites = proposal->world->data->uepsites;
1331 
1332     boolean o, t, ob, tb;
1333 
1334     node *firstb;
1335     node *first;
1336     //memcpy(root->uep, world->oldrootuep, sizeof(long)*uepsites);
1337     o = uep_check (root, proposal->origin->uep, uepsites);
1338     t = uep_check (root, proposal->target->uep, uepsites);
1339     ob = uep_check (root, proposal->oback->uep, uepsites);
1340     tb = uep_check (root, showtop (crawlback (proposal->target))->uep,
1341                     uepsites);
1342     first = first_uep2 (root, root, uepsites);
1343     firstb =
1344         showtop (crawlback
1345                  (first = first_uep (first, proposal->world->root, uepsites)));
1346     proposal->firstuep = first;
1347     if (proposal->world->options->ueprate > 0.0)
1348         proposal->treelen =
1349             calc_pseudotreelength (proposal, proposal->world->treelen);
1350     /*  if(root==proposal->target)
1351        {
1352        oldinterval = root->tyme - firstb->tyme;
1353        interval = proposal->time - firstb->tyme;
1354        }
1355        else
1356        {
1357        if(root==proposal->oback)
1358        {
1359        newroot = crawlback(proposal->oback->next);
1360        if(newroot==proposal->origin)
1361        newroot = crawlback(proposal->oback->next->next);
1362        oldinterval = root->tyme - firstb->tyme;
1363        if(proposal->time > newroot->tyme)
1364        interval = proposal->time - firstb->tyme;
1365        else
1366        interval = newroot->tyme - firstb->tyme;
1367        }
1368        else
1369        {
1370        oldinterval= root->tyme - firstb->tyme;
1371        interval = root->tyme - firstb->tyme;
1372        }
1373        } */
1374     if (o && t && ob && tb)
1375     {
1376         if (proposal->oback == first)
1377         {
1378 
1379             return adjust_uep_oback (proposal->osister, firstb, proposal,
1380                                      interval, oldinterval);
1381         }
1382         else
1383         {
1384             //      FPRINTF(stdout,"+");
1385             return adjust_uep_base (proposal, interval, oldinterval);
1386         }
1387     }
1388     else
1389     {
1390         if (!o && !t && !ob && !tb)
1391         {
1392             if (proposal->osister == first)
1393             {
1394                 //      FPRINTF(stdout,"s");
1395                 if (proposal->target == first)
1396                     return adjust_uep_osistertarget (proposal->osister, firstb,
1397                                                      proposal, interval,
1398                                                      oldinterval);
1399                 else
1400                 {
1401                     if (proposal->target == proposal->oback)
1402                         return adjust_uep_osistertarget2 (proposal->osister,
1403                                                           firstb, proposal,
1404                                                           interval, oldinterval);
1405                     else
1406                         return adjust_uep_osister (proposal->osister,
1407                                                    showtop (crawlback (firstb)),
1408                                                    proposal, interval,
1409                                                    oldinterval);
1410                 }
1411             }
1412             else
1413             {
1414                 //      FPRINTF(stdout,"-");
1415                 return adjust_uep_base (proposal, interval, oldinterval);
1416             }
1417         }
1418         else
1419         {
1420             if (proposal->target == first)
1421             {
1422                 //    if(proposal->oback == firstb)
1423                 // adjust_uep_target_oback(first,firstb,proposal,
1424                 //                      interval,oldinterval);
1425                 //else
1426                 return adjust_uep_target (first, firstb, proposal, interval,
1427                                           oldinterval);
1428             }
1429             else
1430             {
1431                 if (proposal->oback == firstb && proposal->origin != first)
1432                     return adjust_uep_target (first, showtop (crawlback (firstb)),
1433                                               proposal, interval, oldinterval);
1434                 else
1435                 {
1436                     if (proposal->origin == first)
1437                         return adjust_uep_origin (proposal->origin,
1438                                                   proposal->oback, proposal,
1439                                                   interval, oldinterval);
1440                     //  else
1441                     //    error ("do not know what to do with oback!=firstb");
1442                 }
1443             }
1444         }
1445     }
1446     return -MYREAL_MAX;  //never come here
1447 }
1448 
1449 MYREAL
ueplikelihood(world_fmt * world)1450 ueplikelihood (world_fmt * world)
1451 {
1452     long u, i;
1453     MYREAL like = 1.;
1454     MYREAL poplike = 0;
1455     calc_ueplike (world->root->next->back, world, world->ueplike);
1456     //  printf("update_likeli ");
1457     for (u = 0; u < world->data->uepsites; ++u)
1458     {
1459         poplike = 0.;
1460         for (i = 0; i < world->numpop; ++i)
1461         {
1462             poplike += world->ueplike[u][i];
1463             //      printf("%f ", world->ueplike[u][i]);
1464         }
1465         like *= poplike;
1466         //      printf("\n");
1467     }
1468     if (world->options->ueprate > 0.0)
1469     {
1470         return LOG (like) - (world->options->ueprate * world->treelen);
1471     }
1472     else
1473         return LOG (like);
1474 }
1475 
1476 void
show_uep_time(node * p,world_fmt * world)1477 show_uep_time (node * p, world_fmt * world)
1478 {
1479     if (p->type != 'r')
1480     {
1481         if (p->type == 'm')
1482         {
1483             printf ("%li> - %f %f\n", p->id, p->tyme,
1484                     world->root->next->back->tyme - p->tyme);
1485             show_uep_time (showtop (p->back), world);
1486         }
1487         else
1488         {
1489             printf ("%li> %c %f %f\n", p->id, p->uep[0], p->tyme,
1490                     world->root->next->back->tyme - p->tyme);
1491             show_uep_time (showtop (p->back), world);
1492         }
1493     }
1494 }
1495 void
show_uep_time2(node * p,world_fmt * world)1496 show_uep_time2 (node * p, world_fmt * world)
1497 {
1498     if (p->type == 't')
1499         return;
1500     if (p->type != 'r')
1501     {
1502         if (p->type == 'm')
1503         {
1504             show_uep_time2 (showtop (p->next->back), world);
1505             printf ("%li> - %f %f\n", p->id, p->tyme,
1506                     world->root->next->back->tyme - p->tyme);
1507         }
1508         else
1509         {
1510             show_uep_time2 (showtop (p->next->back), world);
1511             show_uep_time2 (showtop (p->next->next->back), world);
1512             printf ("%li> %c %f %f\n", p->id, p->uep[0], p->tyme,
1513                     world->root->next->back->tyme - p->tyme);
1514 
1515         }
1516     }
1517 }
1518 
1519 void
update_uepanc(world_fmt * world)1520 update_uepanc (world_fmt * world)
1521 {
1522     long i, half, pop, popi;
1523     half = world->data->uepsites * world->numpop;
1524     for (i = 0; i < world->data->uepsites; ++i)
1525     {
1526         pop = world->root->next->back->actualpop;
1527         popi = pop * i + pop;
1528         world->ueprootstore[world->locus][world->G][i] =
1529             world->root->next->back->uep[i];
1530         if (world->root->next->back->uep[i] == 0)
1531             world->uepanc[world->locus][popi] += 1;
1532         else
1533             world->uepanc[world->locus][half + popi] += 1;
1534     }
1535 }
1536 
1537 void
copy_uepx(proposal_fmt * proposal,ueparray_fmt xx1,ueparray_fmt xx2)1538 copy_uepx (proposal_fmt * proposal, ueparray_fmt xx1, ueparray_fmt xx2)
1539 {
1540     memcpy (xx1.s, xx2.s, sizeof (pair)*proposal->world->data->uepsites);
1541 }
1542 
1543 #endif
1544 
1545 
1546 
1547 
1548 
1549 
1550