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