1 /*
2 (c) Peter Beerli 2010-12 Tallhassee
3 beerli@fsu.edu
4
5 Permission is hereby granted, free of charge, to any person obtaining a copy
6 of this software and associated documentation files (the "Software"), to deal
7 in the Software without restriction, including without limitation the rights
8 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
9 of the Software, and to permit persons to whom the Software is furnished to do
10 so, subject to the following conditions:
11
12 The above copyright notice and this permission notice shall be included in all copies
13 or substantial portions of the Software.
14
15 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
16 INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
17 PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
18 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
19 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
20 OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21
22 */
23 #ifdef BEAGLE
24 #include "beagle.h"
25 #include "migration.h"
26 #include "tree.h"
27 #include "tools.h"
28 #include "sighandler.h"
29
30 #define BEAGLE_PARTIALS 7
31
32 extern long myID;
33
34 void debug_beagle(beagle_fmt *beagle);
35 long new_id(long id, long sumtips);
36 /*-----------------------------------------------------------------------------
37 | This function sets up the beagle library and initializes all data members.
38 */
init_beagle(world_fmt * world,long locus)39 void init_beagle(world_fmt *world, long locus)
40 {
41
42 beagle_fmt *beagle = world->beagle;
43 beagle->instance_handle = (int *) mycalloc(world->nummutationmodels[locus],sizeof(int));
44 beagle->numallocpartials = 4 * world->sumtips*world->mutationmodels[0].numstates*world->mutationmodels[0].numsites;
45 beagle->partials = (double *) mycalloc(beagle->numallocpartials, sizeof(double));
46
47 beagle->numallocbranches = 200;
48 beagle->branch_lengths = (double *) mycalloc(beagle->numallocbranches, sizeof(double));
49 beagle->numallocoperations = 200;
50 beagle->operations = (int *) mycalloc(beagle->numallocoperations * BEAGLE_PARTIALS,sizeof(int));//parent,parentIDscalingfactor, lchild,lchildtrans,rchild,rtrans
51 beagle->branch_indices = (int *) mycalloc(beagle->numallocoperations * 2, sizeof(int));
52
53 beagle->scalingfactorsindices= (int *) mycalloc(2*(2 * world->sumtips-1), sizeof(int));
54 beagle->scalingfactorscount = 0; //actual count of which scale factors are updated? this was world->sumtips-1;
55
56 beagle->numallocallyweights = world->mutationmodels[0].numsites;
57 beagle->allyweights = (int *) mycalloc(beagle->numallocallyweights,sizeof(int));
58 // this is already set and definitely should be set here, left as reminder
59 // that the current setting in makevalues() is a hack
60 // world->mutationmodels[world->sublocistarts[locus]].basefreqs = (double *) calloc(world->mutationmodels[world->sublocistarts[locus]].numstates,sizeof(double));
61 }
62
set_beagle_instances(world_fmt * world,boolean instance)63 void set_beagle_instances(world_fmt *world, boolean instance)
64 {
65 int resource=1;//cpu=0, gpu=1
66 long i;
67 long ii;
68 long locus = world->locus;
69 beagle_fmt *beagle = world->beagle;
70 for(i=0;i<world->nummutationmodels[locus];i++)
71 {
72 ii = world->sublocistarts[locus] + i;
73 beagle->instance_handle[i] = beagleCreateInstance(world->sumtips, // number of tips
74 2 * (2 * world->sumtips-1), // total buffers (= total nodes in tree)
75 // twice to accomodate rejections.
76 0, // number of compact state representation [funny tips]
77 world->mutationmodels[ii].numstates, // number of states (nucleotides etc)
78 world->mutationmodels[ii].numsites, // number of site patterns
79 1, // number of rate matrices eigen-decomp buffers
80 2*(2 * world->sumtips - 2),// number of rate matrix buffers (=number of branches)
81 world->mutationmodels[ii].numsiterates,// categoryCount
82 0, //2*(2 * world->sumtips - 1), number of scaling buffers
83 &resource, // List of resources (?)
84 1, // Number of resources
85 0L, // preferenceFlags (see BeagleFlags)
86 0L // requirementFlags (see BeagleFlags)
87 );
88
89 if (beagle->instance_handle[i] < 0)
90 usererror("createInstance returned a negative instance handle (and that's not good)");
91 // initialize the instance
92 int code = beagleInitializeInstance(beagle->instance_handle[i], NULL);
93
94 if (code < 0)
95 {
96 usererror("initializeInstance returned a negative error code (and that is bad)\n\n");
97 }
98 }
99 }
100
101
102 long
set_branch_index(node * p,long * bid)103 set_branch_index (node * p, long *bid)
104 {
105 long bbid = *bid;
106 if (p->type != 't')
107 {
108 bbid = set_branch_index (crawlback (p->next), &bbid);
109 bbid = set_branch_index (crawlback (p->next->next), &bbid);
110 }
111 p->bid = bbid;
112 //printf("%li -- %li\n",p->id, bbid);
113 return bbid+1;
114 } /* set_branch_index */
115
116
adjust_beagle(beagle_fmt * beagle)117 void adjust_beagle(beagle_fmt *beagle)
118 {
119 beagle->numallocoperations = beagle->numoperations + 10;
120 beagle->operations = (int *) myrealloc(beagle->operations, sizeof(int) * BEAGLE_PARTIALS * beagle->numallocoperations);
121 beagle->branch_lengths = (double *) myrealloc(beagle->branch_lengths, sizeof(double) * 2 * beagle->numallocoperations);
122 beagle->branch_indices = (int *) myrealloc(beagle->branch_indices, beagle->numallocoperations * 2 * sizeof(int));
123 beagle->numallocbranches = beagle->numallocoperations * 2;
124 }
125
126
prepare_beagle_instances(node * theNode,node * left,node * right,beagle_fmt * beagle)127 void prepare_beagle_instances(node *theNode, node * left, node *right, beagle_fmt *beagle)
128 {
129 long parentid = theNode->id;
130 long rightid = right->id;
131 long leftid = left->id;
132
133 long rightbid = right->bid;
134 long leftbid = left->bid;
135
136 double leftbranch = left->v;
137 double rightbranch = right->v;
138 long ii = beagle->numoperations;
139 long i = ii * BEAGLE_PARTIALS ;
140 long j = ii * 2;
141 if(beagle->numoperations >= beagle->numallocoperations)
142 {
143 adjust_beagle(beagle);
144 }
145 beagle->branch_indices[j] = leftbid;
146 beagle->branch_indices[j+1] = rightbid;
147
148 beagle->operations[i] = parentid;
149 beagle->operations[i+1] = parentid; //BEAGLE_OP_NONE; //parentid;
150 beagle->operations[i+2] = parentid; //BEAGLE_OP_NONE; //parentid;
151 beagle->operations[i+3] = leftid;
152 beagle->operations[i+4] = leftbid;
153 beagle->operations[i+5] = rightid;
154 beagle->operations[i+6] = rightbid;
155
156 beagle->branch_lengths[j] = leftbranch;
157 beagle->branch_lengths[j+1] = rightbranch;
158
159 beagle->scalingfactorsindices[ii] = parentid;
160 beagle->scalingfactorscount += 1;
161 beagle->numbranches += 2;
162 beagle->numoperations++;
163 //#ifdef BEAGLEDEBUG
164 printf("%li> TREE: %i {%i %i %i %i %i %i %i} {%i %i} {%f %f}\n",myID, beagle->numoperations,
165 beagle->operations[i],
166 beagle->operations[i+1] ,
167 beagle->operations[i+2] ,
168 beagle->operations[i+3] ,
169 beagle->operations[i+4] ,
170 beagle->operations[i+5] ,
171 beagle->operations[i+6] ,
172 beagle->branch_indices[j],
173 beagle->branch_indices[j+1],
174 beagle->branch_lengths[j],
175 beagle->branch_lengths[j+1]);
176 //#endif
177 // printf("c");
178 }
179
prepare_beagle_instances_proposal(proposal_fmt * proposal,long trueparentid,long leftid,long leftbid,double leftbranch,long rightid,long rightbid,double rightbranch,beagle_fmt * beagle)180 void prepare_beagle_instances_proposal(proposal_fmt *proposal, long trueparentid, long leftid, long leftbid, double leftbranch,
181 long rightid, long rightbid, double rightbranch, beagle_fmt *beagle)
182 {
183 long i = beagle->numoperations * BEAGLE_PARTIALS;
184 long j = beagle->numoperations * 2;
185 long nodep_boundary = proposal->world->sumtips * 2 - 1;
186 long parentid = trueparentid > nodep_boundary ? trueparentid - nodep_boundary : trueparentid + nodep_boundary;
187 if(beagle->numoperations >= beagle->numallocoperations)
188 {
189 adjust_beagle(beagle);
190 }
191 printf("PINT: (%li, %li) l:%li,%f, r:%li,%f\n",trueparentid,parentid,leftid,leftbranch,rightid,rightbranch);
192 beagle->branch_indices[j] = leftbid;
193 beagle->branch_indices[j+1] = rightbid;
194 beagle->operations[i] = parentid;
195 beagle->operations[i+1] = parentid; //BEAGLE_OP_NONE;//parentid;
196 beagle->operations[i+2] = parentid; //BEAGLE_OP_NONE; //parentid;
197 beagle->operations[i+3] = leftid;
198 beagle->operations[i+4] = leftbid;
199 beagle->operations[i+5] = rightid;
200 beagle->operations[i+6] = rightbid;
201
202 beagle->branch_lengths[j] = leftbranch;
203 beagle->branch_lengths[j+1] = rightbranch;
204 beagle->numbranches += 2;
205 beagle->numoperations++;
206 //#ifdef BEAGLEDEBUG
207 printf("%li> PROP: %i {%i %i %i %i %i %i %i} {%i %i} {%f %f}\n",myID,beagle->numoperations,
208 beagle->operations[i],
209 beagle->operations[i+1] ,
210 beagle->operations[i+2] ,
211 beagle->operations[i+3] ,
212 beagle->operations[i+4] ,
213 beagle->operations[i+5] ,
214 beagle->operations[i+6] ,
215 beagle->branch_indices[j],
216 beagle->branch_indices[j+1],
217 beagle->branch_lengths[j],
218 beagle->branch_lengths[j+1]);
219 //#endif
220 //printf("p");
221 }
222
223 void
evaluate_beagle_instances_proposal(proposal_fmt * proposal,node * mother,node * newdaughter,long newdaughter_id,long newdaughter_bid,MYREAL v)224 evaluate_beagle_instances_proposal (proposal_fmt * proposal,
225 node * mother,
226 node * newdaughter, long newdaughter_id, long newdaughter_bid,
227 MYREAL v)
228 {
229 node *nn = NULL, *d1 = NULL, *d2 = NULL, *oldnn = NULL;
230
231 if (mother->type != 'r')
232 {
233 children (mother, &d1, &d2);
234 if (d1 == newdaughter)
235 {
236 d1 = d2;
237 }
238 long id0, id1, id2;
239 long bid1, bid2;
240 double v1, v2;
241 id0 = mother->id;
242 id1 = new_id(newdaughter_id, proposal->world->sumtips);
243 bid1 = newdaughter_bid;
244 v1 = v;
245 id2 = d1->id;
246 bid2 = d1->bid;
247 v2 = d1->v;
248 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
249
250 oldnn = mother;
251 nn = showtop (crawlback (mother));
252 while (nn->type != 'r')
253 {
254 children (nn, &d1, &d2);
255 if (d1 == oldnn)
256 {
257 d1 = d2;
258 d2 = oldnn;
259 }
260 if(d2->id != oldnn->id)
261 {
262 warning("One of the children is not what it should be\n");
263 }
264 id0 = nn->id;
265 id1 = new_id(d2->id,proposal->world->sumtips);;
266 bid1 = d2->bid;
267 v1 = d2->v;
268 id2 = d1->id;
269 bid2 = d1->bid;
270 v2 = d1->v;
271 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
272
273 oldnn = nn;
274 nn = showtop (crawlback (nn));
275 }
276 }
277 }
278
279
280
281
fill_beagle_instances(world_fmt * world)282 void fill_beagle_instances(world_fmt *world)
283 {
284 beagle_fmt *beagle = world->beagle;
285 unsigned long i;
286 unsigned long ii;
287 unsigned long j;
288 //unsigned long zz;
289 unsigned long site;
290 unsigned int numsites;
291 unsigned int numstates;
292 unsigned int numweights=0;
293 int code;
294 long locus = world->locus;
295 // set partials for all tipnodes
296 // use z to advance through the partial array
297 unsigned long z = 0L;
298
299 for(i=0; i<world->nummutationmodels[locus]; i++)
300 {
301 ii = world->sublocistarts[locus] + i;
302 numstates = world->mutationmodels[ii].numstates;
303 numsites = world->mutationmodels[ii].numsites;
304 // will currently fail with more complex data
305 if(numsites + numweights < beagle->numallocallyweights)
306 {
307 beagle->numallocallyweights += numsites * 2;
308 beagle->allyweights = (int *) myrealloc(beagle->allyweights,sizeof(int) * beagle->numallocallyweights);
309 }
310 memcpy(beagle->allyweights+numweights, world->data->seq[0]->aliasweight, sizeof(int) * numsites);
311 numweights += numsites;
312 for (j = 0; j < world->sumtips; ++j)
313 {
314 for(site=0;site<numsites;site++)
315 {
316 if(z > beagle->numallocpartials)
317 {
318 beagle->numallocpartials += numsites * numstates;
319 beagle->partials = (double *) myrealloc(beagle->partials,sizeof(double) * beagle->numallocpartials);
320 }
321 memcpy(&beagle->partials[z],
322 &(world->nodep[j]->x.s[site][0][0]),
323 sizeof(double) * numstates);
324 beagle->numpartials++;
325 #if 0
326 fprintf(stdout,"Site %li Partial: %li (%f ", site, j,
327 beagle->partials[z]);
328 for(zz=1;zz<world->mutationmodels[ii].numstates;zz++)
329 {
330 fprintf(stdout,"%f ",beagle->partials[z+zz]);
331 }
332 fprintf(stdout,")\n")
333 #endif
334 // advance z
335 z += numstates;
336 }
337 code = beagleSetTipPartials(
338 beagle->instance_handle[i], // instance
339 j, // indicator for tips
340 &beagle->partials[z - numsites*numstates]);// inPartials
341 if (code != 0)
342 usererror("setTipPartials encountered a problem");
343 }
344 }
345
346
347 for(i=0; i<world->nummutationmodels[locus]; i++)
348 {
349 ii = world->sublocistarts[locus] + i;
350
351 code = beagleSetCategoryRates(beagle->instance_handle[i], world->mutationmodels[ii].siterates);
352 if (code != 0)
353 usererror("setCategoryRates encountered a problem");
354
355 code = beagleSetEigenDecomposition(beagle->instance_handle[i], // instance
356 0, // eigenIndex,
357 (const double *)world->mutationmodels[ii].eigenvectormatrix, // inEigenVectors,
358 (const double *)world->mutationmodels[ii].inverseeigenvectormatrix,// inInverseEigenVectors,
359 world->mutationmodels[i].eigenvalues); // inEigenValues
360
361 if (code != 0)
362 usererror("setEigenDecomposition encountered a problem");
363 }
364 }
365 /*-----------------------------------------------------------------------------
366 | Calculates the log likelihood by calling the beagle functions
367 | updateTransitionMatrices, updatePartials and calculateEdgeLogLikelihoods.
368 */
calcLnL(world_fmt * world,boolean instance)369 double calcLnL(world_fmt *world, boolean instance)
370 {
371 beagle_fmt *beagle = world->beagle;
372 double logL = 0.0;
373 unsigned long i;
374 unsigned long j;
375 //unsigned long z;
376 long locus = world->locus;
377 long ii;
378 double *patternloglike = (double *) calloc(world->maxnumpattern[locus],sizeof(double));
379 double *outlike = (double *) calloc(10*world->maxnumpattern[locus],sizeof(double));
380 int rootIndex = beagle->operations[BEAGLE_PARTIALS * (beagle->numoperations-1)];//world->root->next->back->id;
381 for(i=0; i<world->nummutationmodels[locus]; i++)
382 {
383 ii = world->sublocistarts[locus] + i;
384 int code = beagleUpdateTransitionMatrices(beagle->instance_handle[i], // instance,
385 0, // eigenIndex,
386 (const int *) beagle->branch_indices, // indicators transitionrates for each branch,
387 NULL, // firstDerivativeIndices,
388 NULL, // secondDervativeIndices,
389 beagle->branch_lengths, // edgeLengths,
390 beagle->numbranches); // number branches to update, count
391
392 if (code != 0)
393 usererror("updateTransitionMatrices encountered a problem");
394
395 int cumulativeScalingFactorIndex = 0; //BEAGLE_OP_NONE; //this would be the index of the root scaling location
396
397 beagleResetScaleFactors(beagle->instance_handle[i],
398 cumulativeScalingFactorIndex);
399
400 beagleAccumulateScaleFactors(beagle->instance_handle[i],
401 beagle->scalingfactorsindices,
402 beagle->scalingfactorscount,
403 cumulativeScalingFactorIndex);
404
405
406 code = beagleUpdatePartials((const int *) &beagle->instance_handle[i], // instance
407 1, // instanceCount
408 beagle->operations, // operations
409 beagle->numoperations, // operationCount
410 cumulativeScalingFactorIndex);//BEAGLE_OP_NONE); // connected to accumulate....
411 #ifdef BEAGLEDEBUG
412 for(j=0;j<2*(world->sumtips * 2 - 1); j++)
413 {
414 beagleGetPartials(beagle->instance_handle[i],j,BEAGLE_OP_NONE, outlike);
415 if(j==world->sumtips * 2 - 1)
416 printf("-----------------------\n");
417 printf("%li: {%f, %f, %f, %f}\n",j, outlike[0],outlike[1],outlike[2],outlike[3]);
418 }
419 #endif
420 if (code != 0)
421 usererror("updatePartials encountered a problem");
422
423
424
425 if(beagle->weights==NULL)
426 beagle->weights = (double *) mycalloc(1,sizeof(double));
427 // else
428 // beagle->weights = (double *) myrealloc(beagle->weights,beagle->numoperations * sizeof(double));
429
430 beagle->weights[0]= 1.0;
431
432 // calculate the site likelihoods at the root node
433 code = beagleCalculateRootLogLikelihoods(beagle->instance_handle[i], // instance
434 (const int *) &rootIndex, // bufferIndices
435 (const double *) world->mutationmodels[ii].siteprobs, // weights
436 (const double *) world->mutationmodels[ii].basefreqs,// stateFrequencies
437 &cumulativeScalingFactorIndex, //scalingfactors index,
438 1, // count is this correct
439 patternloglike); // outLogLikelihoods
440 //trash from function above &beagle->scalingfactorscount,//size of the scaling factor index
441 if (code != 0)
442 usererror("calculateRootLogLikelihoods encountered a problem");
443
444 for (j = 0; j < world->mutationmodels[ii].numsites; j++)
445 {
446 logL += beagle->allyweights[j] * patternloglike[j];
447 // printf("%.1f ",patternloglike[j]);
448 }
449 }
450 printf("Log LnL=%f (instance=%li)\n",logL,(long) instance);
451 #ifdef BEAGLEDEBUG
452 debug_beagle(beagle);
453 #endif
454 myfree(patternloglike);
455 myfree(outlike);
456 return logL;
457 }
458
force_beagle_recalculate(world_fmt * world,long locus)459 double force_beagle_recalculate(world_fmt *world, long locus)
460 {
461 reset_beagle(world->beagle);
462 set_all_dirty(world->root->next, crawlback (world->root->next), world, locus);
463 smooth (world->root->next, crawlback (world->root->next), world, locus);
464 return treelikelihood(world);
465 }
466
467 ///
468 /// start out at the root node and traverse up and match the ids with the new ids
change_beagle(node * theNode,beagle_fmt * beagle,long sumtips)469 void change_beagle(node *theNode, beagle_fmt *beagle, long sumtips)
470 {
471 // move node ids;
472 // this works only just after the pseudolikelihood step + the acceptlike with TRUE outcome
473 long i;
474 long parentid;
475 long oldparentid;
476 long nodep_boundary = sumtips * 2 - 1;
477 if (theNode->type != 't')
478 {
479 change_beagle(crawlback(theNode->next), beagle,sumtips);
480 change_beagle(crawlback(theNode->next->next),beagle, sumtips);
481 }
482 for(i=0;i<beagle->numoperations;i++)
483 {
484 parentid = beagle->operations[i*BEAGLE_PARTIALS];
485 oldparentid = parentid < nodep_boundary ? parentid + nodep_boundary : parentid - nodep_boundary;
486 if(theNode->id == oldparentid)
487 {
488 theNode->id = parentid;
489 }
490 }
491 }
492
new_id(long id,long sumtips)493 long new_id(long id, long sumtips)
494 {
495 const long nodep_boundary = 2 * sumtips - 1;
496 return (id < nodep_boundary ? id + nodep_boundary : id - nodep_boundary);
497 }
498
reset_beagle(beagle_fmt * beagle)499 void reset_beagle(beagle_fmt *beagle)
500 {
501 beagle->numoperations = 0;
502 beagle->numbranches = 0;
503 beagle->scalingfactorscount = 0;
504 }
505
beagle_stop(world_fmt ** universe,long usize)506 void beagle_stop(world_fmt **universe, long usize)
507 {
508 long u;
509 for(u=0;u< usize; u++)
510 {
511 world_fmt *world = universe[u];
512 beagle_fmt *beagle = world->beagle;
513 unsigned long i;
514 long locus = world->locus;
515 for(i=0; i<world->nummutationmodels[locus]; i++)
516 {
517 beagleFinalizeInstance(beagle->instance_handle[i]);
518 }
519 }
520 }
521
522
set_beagle_dirty(node * origin,node * target,node * mrca)523 void set_beagle_dirty(node *origin, node *target, node *mrca)
524 {
525 node *nn = origin;
526 while((nn=showtop(crawlback(nn)))!=mrca)
527 {
528 set_dirty(nn);
529 }
530 nn = target;
531 while((nn=showtop(crawlback(nn)))!=mrca)
532 {
533 set_dirty(nn);
534 }
535 }
536
debug_beagle(beagle_fmt * beagle)537 void debug_beagle(beagle_fmt *beagle)
538 {
539 printf("----beagle content---------------------\n");
540 printf("operations : %p\n",beagle->operations);
541 printf(" alloc : %i\n",beagle->numallocoperations);
542 printf(" num : %i\n",beagle->numoperations);
543 printf("branch indices: %p\n",beagle->branch_indices);
544 printf(" alloc : %i\n",beagle->numallocbranches);
545 printf(" num : %i\n",beagle->numbranches);
546 printf("----beagle content end-----------------\n\n");
547 }
548
549
550
551 #endif /*BEAGLE*/
552