1 /*
2
3 HyPhy - Hypothesis Testing Using Phylogenies.
4
5 Copyright (C) 1997-now
6 Core Developers:
7 Sergei L Kosakovsky Pond (sergeilkp@icloud.com)
8 Art FY Poon (apoon42@uwo.ca)
9 Steven Weaver (sweaver@temple.edu)
10
11 Module Developers:
12 Lance Hepler (nlhepler@gmail.com)
13 Martin Smith (martin.audacis@gmail.com)
14
15 Significant contributions from:
16 Spencer V Muse (muse@stat.ncsu.edu)
17 Simon DW Frost (sdf22@cam.ac.uk)
18
19 Permission is hereby granted, free of charge, to any person obtaining a
20 copy of this software and associated documentation files (the
21 "Software"), to deal in the Software without restriction, including
22 without limitation the rights to use, copy, modify, merge, publish,
23 distribute, sublicense, and/or sell copies of the Software, and to
24 permit persons to whom the Software is furnished to do so, subject to
25 the following conditions:
26
27 The above copyright notice and this permission notice shall be included
28 in all copies or substantial portions of the Software.
29
30 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
31 OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
32 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
33 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
34 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
35 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
36 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
37
38 */
39
40 #include <math.h>
41 #include <ctype.h>
42 #include <float.h>
43
44 #include "calcnode.h"
45 #include "function_templates.h"
46
47 #include "category.h"
48 #include "batchlan.h"
49 #include "tree.h"
50
51 #include "global_things.h"
52 #include "hbl_env.h"
53
54 using namespace hy_global;
55 using namespace hy_env;
56
57 //#define _UBER_VERBOSE_MX_UPDATE_DUMP
58 //#define _UBER_VERBOSE_MX_UPDATE_DUMP_EVAL 1
59
60 //_______________________________________________________________________________________________
61
_CalcNode()62 _CalcNode::_CalcNode () {
63 theProbs = nil; // default constructor, doesn't do much
64 compExp = nil;
65 matrixCache = nil;
66 flags = 0;
67 cBase = 0;
68 }
69
70 //_______________________________________________________________________________________________
_CalcNode(_String name,_String parms,int codeBase,_VariableContainer * theP,_AVLListXL * aCache)71 _CalcNode::_CalcNode (_String name, _String parms, int codeBase, _VariableContainer* theP, _AVLListXL* aCache):_VariableContainer (name, "", theP) {
72 // construct a node from a string of the form
73 // matrix name, <optional comma separated variable declarations, inititalizations>
74 // also should be passed the pointer to a container tree
75 InitializeCN (parms, codeBase, theP, aCache);
76 }
77
78 //_______________________________________________________________________________________________
_CalcNode(_CalcNode * sourceNode,_VariableContainer * theP)79 _CalcNode::_CalcNode (_CalcNode* sourceNode, _VariableContainer* theP):_VariableContainer (sourceNode->ContextFreeName(), "", theP) {
80 _String model = *sourceNode->GetModelName();
81 InitializeCN (model, 0, theP);
82 if (model.nonempty()) { // copy model parameter values
83 CopyMatrixParameters(sourceNode, true);
84 }
85 }
86
87 //_______________________________________________________________________________________________
InitializeCN(_String const & parms,int,_VariableContainer * theP,_AVLListXL * aCache)88 void _CalcNode::InitializeCN ( _String const& parms, int, _VariableContainer* theP, _AVLListXL* aCache) {
89 if (theIndex < 0) return;
90
91 cBase = 0;
92 flags = 0;
93 theProbs = nil;
94 compExp = nil;
95 matrixCache = nil;
96
97 _List parameters;
98 _ElementaryCommand::ExtractConditions(parms, 0, parameters, ',');
99
100 InitializeVarCont (kEmptyString, *(_String*)parameters.GetItem(0), theP, aCache);
101
102 // parse and instantiate all model expressions
103 parameters.ForEach([=] (BaseRef expression, unsigned long) -> void {
104 _Formula fg (*(_String*)expression, this);
105 }, (GetModelIndex() == HY_NO_MODEL && parms.nonempty()) ? 0L : 1L);
106
107 // attach the variables to the lists inside the node
108 ScanAndAttachVariables();
109
110 // check for category variables
111
112 //_SimpleList local_cat_vars;
113
114 ForEachLocalVariable(iVariables, [&] (long var_idx, long ref_idx, unsigned long idx) -> void {
115 if (ref_idx >= 0) {
116 /* TODO:
117 this used to do with local category variables
118 NOT TESTED and is LIKELY BROKEN
119 replaced with an error message;
120 */
121 _Variable *test_var = (_Variable*)LocateVar(ref_idx);
122 if (test_var->IsCategory()) {
123 HandleApplicationError(_String ("Local instance to category variables are not supported, ") & test_var->GetName()->Enquote());
124 return;
125 }
126 }
127 });
128
129 if (gVariables) {
130
131 _SimpleList moved_to_cat;
132
133 gVariables->Each ([&] (long var_idx, unsigned long idx) -> void {
134 _Variable *test_v = LocateVar(var_idx);
135 if (test_v->IsCategory()) {
136 categoryVariables << var_idx;
137 categoryIndexVars << -1;
138 moved_to_cat << idx;
139 }
140 });
141
142 if (moved_to_cat.nonempty()) {
143 gVariables->DeleteList(moved_to_cat);
144 if (gVariables->empty()) {
145 DeleteAndZeroObject(gVariables);
146 }
147 }
148 }
149
150 variablePtrs.Replace(theIndex, this, true);
151
152 /** TODO check equivalence **
153
154 BaseRef temp = variablePtrs(theIndex);
155 variablePtrs[theIndex]=this->makeDynamic();
156 DeleteObject(temp);
157 */
158
159 }
160
161 //__________________________________________________________________________________
SetModel(long modelID,_AVLListXL * varCache)162 void _CalcNode::SetModel (long modelID, _AVLListXL* varCache) {
163 _VariableContainer::SetModel (modelID, varCache);
164 }
165
166 //_______________________________________________________________________________________________
167
SetDependance(long var_index)168 long _CalcNode::SetDependance (long var_index) {
169 var_index = _VariableContainer::SetDependance (var_index);
170 if (var_index >= 0L) {
171 /** if the new constraint includes a category variable,
172 it needs to be marked as such in this CalcNode
173 */
174
175
176 PopulateAndSort ([&] (_AVLList& avl) -> void {
177 LocateVar (var_index)->ScanForVariables (avl,true);
178 }).Each ( [&] (long var_index, unsigned long idx) -> void {
179 if (LocateVar(var_index)->IsCategory() &&(categoryVariables >> var_index)) {
180 categoryIndexVars<<-1;
181 }
182 });
183
184 // also clear out previously computed matrix exponentials
185 if (compExp) {
186 DeleteAndZeroObject(compExp);
187 } /*else {
188 if (matrixCache) {
189 }
190 }*/
191 }
192 return var_index;
193 }
194
195 //_______________________________________________________________________________________________
196
SetCodeBase(int codeBase)197 void _CalcNode::SetCodeBase (int codeBase) {
198 if (codeBase>0) {
199 if (codeBase != cBase || !theProbs) {
200 if (theProbs) {
201 delete [] theProbs;
202 }
203 theProbs = new hyFloat [codeBase];
204 cBase = codeBase;
205 }
206
207 theProbs[0] = 1.0;
208 }
209 }
210
211 //_______________________________________________________________________________________________
SetCompMatrix(long categID)212 void _CalcNode::SetCompMatrix (long categID) {
213 compExp = GetCompExp (categID);
214 }
215
216 //_______________________________________________________________________________________________
217
ProcessTreeBranchLength(_String const & branch_length)218 hyFloat _CalcNode::ProcessTreeBranchLength (_String const& branch_length) {
219 hyFloat res = -1.;
220
221 if (branch_length.nonempty()) {
222 if (branch_length.char_at(0UL)==':') {
223 res = branch_length.Cut(1L,-1L).to_float();
224 } else {
225 res = branch_length.to_float ();
226 }
227 res = MAX (res, 1e-10);
228 }
229
230 return res;
231 }
232
233 //_______________________________________________________________________________________________
234
Clear(void)235 void _CalcNode::Clear (void) {
236 if (compExp) {
237 DeleteAndZeroObject (compExp);
238 }
239 if (theProbs) {
240 delete [] theProbs;
241 theProbs = nil;
242 }
243 _VariableContainer::Clear();
244 }
245
246 //_______________________________________________________________________________________________
247
~_CalcNode(void)248 _CalcNode::~_CalcNode (void) {
249 Clear();
250 }
251
252 //_______________________________________________________________________________________________
253
FreeUpMemory(long)254 long _CalcNode::FreeUpMemory (long) {
255 long res = 0L;
256 if (compExp) {
257 res = compExp->GetMySize();
258 DeleteAndZeroObject (compExp);
259 }
260 return res;
261 }
262
263 //__________________________________________________________________________________
264
RemoveModel(void)265 void _CalcNode::RemoveModel (void) {
266
267 categoryVariables.Clear();
268 categoryIndexVars.Clear();
269 remapMyCategories.Clear();
270
271 Clear();
272 }
273
274 //__________________________________________________________________________________
GetBranchSpec(void)275 _String* _CalcNode::GetBranchSpec (void) {
276
277 _StringBuffer * res = new _StringBuffer (32UL);
278 *res << GetModelName();
279
280 ForEachLocalVariable(iVariables, [&] (long var_index, long ref_index, unsigned long idx) -> void {
281 if (idx == 0UL) {
282 (*res) << (res->nonempty() ? ',' : '{');
283 } else {
284 (*res) << ',';
285 }
286 _Variable * av = LocateVar (var_index);
287 if (ref_index >= 0L) {
288 res->AppendAnAssignmentToBuffer(LocateVar (ref_index)->GetName(),
289 new _String (av->Value()));
290 } else {
291 res->AppendAnAssignmentToBuffer(av->GetName(),
292 new _String (av->Value()));
293 }
294 });
295
296 ForEachLocalVariable(dVariables, [&] (long var_index, long ref_index, unsigned long) -> void {
297 if (ref_index < 0L) {
298 (*res) << (res->nonempty() ? ',' : '{');
299
300 _Variable * av = LocateVar (var_index);
301 res->AppendAnAssignmentToBuffer(av->GetName(),
302 av->GetFormulaString(kFormulaStringConversionNormal),
303 kAppendAnAssignmentToBufferFree | kAppendAnAssignmentToBufferAssignment);
304 }
305 });
306
307 if (res->nonempty()) {
308 (*res) << '}';
309 }
310
311 res->TrimSpace();
312 return res;
313 }
314
315
316 //__________________________________________________________________________________
317
ReplaceModel(_String & modelName,_VariableContainer * parent_tree)318 void _CalcNode::ReplaceModel (_String& modelName, _VariableContainer* parent_tree) {
319 // TODO: SLKP 20171203 this is FUGLY
320 RemoveModel ();
321
322 _TheTree * parent_tree_object = (_TheTree*)parent_tree;
323
324 long index_in_parent = parent_tree_object->flatCLeaves._SimpleList::Find ((long)this);
325 _List * container_object = nil;
326
327 if (index_in_parent >= 0) {
328 parent_tree_object->flatCLeaves.Replace (index_in_parent, nil, false);
329 container_object = &parent_tree_object->flatCLeaves;
330 } else {
331 index_in_parent = parent_tree_object->flatTree._SimpleList::Find ((long)this);
332 if (index_in_parent >= 0) {
333 parent_tree_object->flatTree.Replace (index_in_parent, nil, false);
334 container_object = &parent_tree_object->flatTree;
335 }
336 }
337
338 long my_var_index = theIndex;
339
340 DeleteVariable (theIndex, false); // this will clean up all the node.xx variables
341 InitializeCN (modelName, 0 , parent_tree);
342
343 if (container_object) {
344 _Variable * new_node = LocateVar(my_var_index);
345 container_object->Replace (index_in_parent, new_node, false);
346 new_node->AddAReference();
347 }
348 }
349
350 //_______________________________________________________________________________________________
MatchSubtree(_CalcNode * mNode)351 bool _CalcNode::MatchSubtree (_CalcNode* mNode) {
352 // TODO: SLKP 20171203 possible deprecate
353 node <long>* myNode = LocateMeInTree (),
354 * matchNode = mNode->LocateMeInTree ();
355 if (myNode&&matchNode) {
356 return myNode->compare_subtree(matchNode);
357 }
358 return false;
359 }
360
361 //_______________________________________________________________________________________________
362
ComputeBranchLength(void)363 hyFloat _CalcNode::ComputeBranchLength (void) {
364
365 static const _String kLargeMatrixBranchLengthDimension ("LARGE_MATRIX_BRANCH_LENGTH_MODIFIER_DIMENSION"),
366 kLargeMatrixBranchLengthModifier ("LARGE_MATRIX_BRANCH_LENGTH_MODIFIER");
367
368
369 if (GetModelIndex() < 0) {
370 return Value();
371 }
372
373 HBLObjectRef stencil = (HBLObjectRef)hy_env::EnvVariableGet(hy_env::branch_length_stencil, STRING | ASSOCIATIVE_LIST);
374
375 if (stencil) {
376 if (stencil->ObjectClass () == STRING) {
377 if (((_FString*)stencil)->get_str() == kStringSuppliedLengths) {
378 return Value();
379 }
380 } else {
381 _AssociativeList *lookup = (_AssociativeList*)stencil;
382 _String lookup_name = ContextFreeName();
383 _Constant * value = (_Constant*)lookup->GetByKey (lookup_name, NUMBER);
384 if (value) {
385 return value->Value();
386 }
387 }
388 }
389
390
391 _Matrix *freqMx,
392 *theMx;
393
394 bool mbf;
395
396 RetrieveModelComponents (theModel, theMx, freqMx, mbf);
397
398 if ( ! (freqMx && theMx)) {
399 return Value();
400 }
401
402 freqMx = (_Matrix*)freqMx->ComputeNumeric();
403
404 hyFloat result = 0.0;
405
406 IntergrateOverAssignments (categoryVariables, true, [&] (long current_cat, hyFloat weight) -> void {
407 result += fabs(ComputeModelMatrix()->ExpNumberOfSubs (freqMx, mbf))*weight;
408 });
409
410 if (freqMx->GetSize() > hy_env::EnvVariableGetNumber(kLargeMatrixBranchLengthDimension, 20.)) {
411 result /= hy_env::EnvVariableGetNumber(kLargeMatrixBranchLengthModifier, 3.);
412 }
413
414 return result;
415 }
416
417 //_______________________________________________________________________________________________
418
operator [](unsigned long i)419 hyFloat& _CalcNode::operator[] (unsigned long i) {
420 // TODO SLKP 20171203, possibly DEPRECATE?
421 return theProbs [i];
422 }
423
424 //_______________________________________________________________________________________________
425
toStr(unsigned long)426 BaseRef _CalcNode::toStr (unsigned long) {
427 _StringBuffer * res = new _StringBuffer (64L);
428 (*res) << theName << '(' << _String(CountIndependents()) << ',' << _String(CountDependents()) << ')';
429 res->TrimSpace();
430 return res;
431 }
432
433 //__________________________________________________________________________________
434
Duplicate(BaseRefConst theO)435 void _CalcNode::Duplicate (BaseRefConst theO) {
436 _VariableContainer::Duplicate (theO);
437 cBase = 0;
438 compExp = nil;
439 matrixCache = nil;
440 theProbs = nil;
441 }
442
443 //_______________________________________________________________________________________________
444
HasChanged(bool)445 bool _CalcNode::HasChanged(bool) {
446
447 return _VariableContainer::HasChanged() || categoryVariables.Any([&] (long cat_idx, unsigned long) -> bool {
448 return LocateVar (cat_idx)->HasChanged();
449 });
450 }
451
452
453 //_______________________________________________________________________________________________
454
NeedNewCategoryExponential(long catID) const455 bool _CalcNode::NeedNewCategoryExponential(long catID) const {
456
457 if (_VariableContainer::NeedToExponentiate(catID>=0)) {
458 return true;
459 }
460
461 if (catID==-1) {
462 return !compExp || categoryVariables.Any([&] (long cat_idx, unsigned long) -> bool {
463 return LocateVar (cat_idx)->HasChanged();
464 });
465 } else {
466 return !GetCompExp(catID) || categoryVariables.Any([&] (long cat_idx, unsigned long i) -> bool {
467 return ((_CategoryVariable*)LocateVar (cat_idx))->HaveParametersChanged(remapMyCategories.list_data[catID*(categoryVariables.countitems()+1)+i+1]);
468 });
469 }
470 return false;
471 }
472
473 //_______________________________________________________________________________________________
RecomputeMatrix(long categID,long totalCategs,_Matrix * storeRateMatrix,_List * queue,_SimpleList * tags,_List * bufferedOps)474 bool _CalcNode::RecomputeMatrix (long categID, long totalCategs, _Matrix* storeRateMatrix, _List* queue, _SimpleList* tags, _List* bufferedOps)
475 {
476 // assumed that NeedToExponentiate was called prior to this function
477
478 //_Variable* curVar, *locVar;
479
480 _SimpleList * var_lists [2] = {iVariables, dVariables};
481 for (_SimpleList* iterable : var_lists) {
482 ForEachLocalVariable(iterable,CopyModelParameterValue);
483 }
484
485 #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
486 // if (1|| likeFuncEvalCallCount == _UBER_VERBOSE_MX_UPDATE_DUMP_LF_EVAL && gVariables) {
487 for (unsigned long i=0; i<gVariables->lLength; i++) {
488 _Variable* curVar = LocateVar(gVariables->GetElement(i));
489 fprintf (stderr, "[_CalcNode::RecomputeMatrix] Node %s, var %s, value = %15.12g\n", GetName()->get_str(), curVar->GetName()->get_str(), curVar->Compute()->Value());
490 }
491 //}
492 #endif
493
494 /*
495 for (unsigned long i=0; i<categoryVariables.lLength; i++) {
496 if (categoryIndexVars.list_data[i]<0) {
497 continue;
498 }
499 curVar = LocateVar (categoryIndexVars.list_data[i]);
500 locVar = LocateVar (categoryVariables.list_data[i]);
501 curVar->SetValue(locVar->Compute());
502 } */
503
504 if (!storeRateMatrix) {
505 if (totalCategs>1) {
506 #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
507 fprintf (stderr, "[_CalcNode::RecomputeMatrix] Deleting category %ld for node %s at %p\n", categID, GetName()->get_str(), GetCompExp(categID));
508 #endif
509 if (clear_exponentials()) {
510 DeleteObject(GetCompExp(categID, true));
511 }
512 } else {
513 if (compExp) {
514 if (clear_exponentials ()) {
515 DeleteAndZeroObject(compExp);
516 }
517 }
518 }
519 }
520
521 bool isExplicitForm = HasExplicitFormModel ();
522
523 if (isExplicitForm && bufferedOps) {
524 _Matrix * bufferedExp = (_Matrix*)GetExplicitFormModel()->Compute (0,nil, bufferedOps);
525 #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
526 fprintf (stderr, "[_CalcNode::RecomputeMatrix] Setting (buffered) category %ld/%ld for node %s\n", categID, totalCategs, GetName()->get_str());
527 #endif
528 SetCompExp ((_Matrix*)bufferedExp->makeDynamic(), totalCategs>1?categID:-1);
529 return false;
530 }
531
532 unsigned long previous_length = queue && tags ? queue->lLength: 0;
533 _Matrix * myModelMatrix = GetModelMatrix(queue,tags);
534
535 if (isExplicitForm && !myModelMatrix) { // matrix exponentiations got cached
536 if (queue && queue->lLength > previous_length) {
537 return true;
538 } else {
539 HandleApplicationError ("Internal error");
540 return false;
541 }
542 } else {
543
544 if (myModelMatrix->MatrixType()!=_POLYNOMIAL_TYPE) {
545 _Matrix *temp = nil;
546 if (isExplicitForm) {
547 temp = (_Matrix*)myModelMatrix->makeDynamic();
548 } else {
549 temp = (_Matrix*)myModelMatrix->MultByFreqs(theModel, true);
550 }
551
552 // copy updated model (local) constrained parameters to their external references
553 ForEachLocalVariable(dVariables, [&] (long var_idx, long ref_idx, unsigned long) -> void {
554 if (ref_idx >= 0) {
555 _Variable * model_var = LocateVar (ref_idx);
556 if (!model_var -> IsIndependent()) {
557 _Variable * param_var = LocateVar (var_idx);
558 if (param_var->IsIndependent()) {
559 param_var->SetValue(param_var->Compute(),true,true,NULL);
560 }
561 }
562 }
563 });
564
565
566 if (storeRateMatrix) {
567 storeRateMatrix->Duplicate(temp);
568 return isExplicitForm;
569 }
570
571
572 if (queue) {
573 (*queue) << temp;
574 if (tags) {
575 (*tags) << (isExplicitForm);
576 }
577 return isExplicitForm;
578 }
579
580 #ifdef _UBER_VERBOSE_MX_UPDATE_DUMP
581 fprintf (stderr, "[_CalcNode::RecomputeMatrix] Setting category %ld/%ld for node %s\n", categID, totalCategs, GetName()->get_str());
582 #endif
583 SetCompExp ((_Matrix*)(isExplicitForm?temp:temp->Exponentiate(1., true)), totalCategs>1?categID:-1);
584
585 } else {
586 compExp = (_Matrix*)myModelMatrix->Evaluate(false);
587 }
588 }
589 return false;
590 }
591
592 //_______________________________________________________________________________________________
SetCompExp(_Matrix * m,long catID,bool do_exponentiation)593 void _CalcNode::SetCompExp (_Matrix* m, long catID, bool do_exponentiation) {
594
595 _Matrix ** store_exp_here;
596
597 if (catID >= 0 && matrixCache) {
598 if (remapMyCategories.lLength) {
599 catID = remapMyCategories.list_data[catID*(categoryVariables.lLength+1)];
600 }
601 store_exp_here = &(matrixCache[catID]);
602 } else {
603 store_exp_here = &compExp;
604 }
605
606 if (do_exponentiation) {
607 compExp = m->Exponentiate(1., true, *store_exp_here);
608 reuse_exponentials ();
609 } else {
610 compExp = m;
611 }
612 *store_exp_here = compExp;
613 }
614 //_______________________________________________________________________________________________
ComputeModelMatrix(bool)615 _Matrix* _CalcNode::ComputeModelMatrix (bool) {
616 // assumed that NeedToExponentiate was called prior to this function
617 _SimpleList * var_lists [2] = {iVariables, dVariables};
618 for (_SimpleList* iterable : var_lists) {
619 ForEachLocalVariable(iterable,CopyModelParameterValue);
620 }
621
622 _Matrix * modelMx = GetModelMatrix();
623 if (modelMx && modelMx->ObjectClass()==MATRIX && modelMx->MatrixType()!=_POLYNOMIAL_TYPE) {
624 return (_Matrix*)modelMx->ComputeNumeric();
625 }
626
627 return nil;
628 }
629
630 //_______________________________________________________________________________________________
631
GetCompExp(long catID,bool doClear) const632 _Matrix* _CalcNode::GetCompExp (long catID, bool doClear) const {
633 if (catID==-1) {
634 return compExp;
635 } else {
636
637 if (remapMyCategories.lLength) {
638 catID = remapMyCategories.list_data[catID * (categoryVariables.lLength+1)];
639 }
640
641 _Matrix* ret = matrixCache?matrixCache[catID]:compExp;
642
643 if (doClear && matrixCache) {
644 if (matrixCache[catID] && matrixCache[catID]->CanFreeMe()) {
645 matrixCache[catID] = nil;
646 }
647 }
648 return ret;
649 }
650 }
651
652 //_______________________________________________________________________________________________
653
makeDynamic(void) const654 BaseRef _CalcNode::makeDynamic(void) const {
655 _CalcNode* res = new (_CalcNode);
656 res->_VariableContainer::Duplicate (this);
657 res->categoryVariables.Duplicate (&categoryVariables);
658 res->categoryIndexVars.Duplicate (&categoryIndexVars);
659 res->theValue = theValue;
660 res->cBase = cBase;
661 res->flags = flags;
662 if (cBase) {
663 res->theProbs = new hyFloat [cBase];
664 CopyArray(res->theProbs, theProbs, cBase);
665 } else {
666 res->theProbs = nil;
667 }
668 res->compExp = compExp;
669 if (compExp) {
670 compExp->AddAReference();
671 }
672 return res;
673 }
674
675 //_______________________________________________________________________________________________
676
SetupCategoryMap(_List & containerVariables,_SimpleList & classCounter,_SimpleList & multipliers)677 void _CalcNode::SetupCategoryMap (_List& containerVariables, _SimpleList& classCounter, _SimpleList& multipliers) {
678 // TODO 20171203, SLKP: this needs review
679
680 long totalCategories = classCounter.Element(-1),
681 globalCatCount = containerVariables.countitems(),
682 localCategories = 1L,
683 catCount = categoryVariables.countitems(),
684 entriesPerCat = 1L+catCount;
685
686 //for (long k = 0; k<categoryVariables.lLength;k++)
687 // printf ("%ld\n", categoryVariables(k));//, ((_Variable*)categoryVariables(k))->GetName()->sData);
688
689 if (catCount == 0L) {
690 remapMyCategories.Clear();
691 } else {
692 remapMyCategories.Populate (totalCategories*entriesPerCat,0,0);
693
694 _SimpleList remappedIDs,
695 rateMultiplers (catCount,1,0),
696 categoryValues (globalCatCount,0,0);
697
698 /* find where in the global list the categories for this node are */
699 for (long myCatID = 0L; myCatID < catCount; myCatID++) {
700 long index = containerVariables.FindPointer(LocateVar(categoryVariables.get(myCatID)));
701 if (index < 0) {
702 HandleApplicationError ("Internal error in SetupCategoryMap. Please report to sergeilkp@icloud.com");
703 }
704 localCategories *= classCounter.get(index);
705 remappedIDs << index;
706 }
707
708 /* generate multiplicative offsets for each category; a move by one category value
709 in variable myCatID changes the compositve index by "offset"
710 */
711 for (long myCatID = catCount-2L; myCatID >= 0L; myCatID--) {
712 rateMultiplers.list_data[myCatID] = rateMultiplers.list_data[myCatID+1]*classCounter.list_data[remappedIDs.list_data[myCatID+1]];
713 }
714
715 for (long currentRateCombo = 0L; currentRateCombo < totalCategories; currentRateCombo++) {
716 long copyRateCombo = currentRateCombo;
717 for (long variableID = 0L; variableID < globalCatCount; variableID++) {
718 categoryValues.list_data[variableID] = copyRateCombo / multipliers.list_data[variableID];
719 copyRateCombo = copyRateCombo%multipliers.list_data[variableID];
720 //printf ("%d %d %d %d\n", currentRateCombo, variableID, multipliers.list_data[variableID], categoryValues.list_data[variableID]);
721 }
722
723 long localCatID = 0L;
724
725 for (long localVariableID = 0; localVariableID<catCount; localVariableID++) {
726 localCatID += rateMultiplers.list_data[localVariableID] * categoryValues.list_data[remappedIDs.list_data[localVariableID]];
727 }
728
729 long offset = currentRateCombo * entriesPerCat;
730 remapMyCategories.list_data[offset] = localCatID;
731 //printf ("[%ld] = %ld (%ld)\n", offset, localCatID, );
732
733 offset++;
734 for (long localVariableID = 0; localVariableID<catCount; localVariableID++) {
735 remapMyCategories[offset++] = categoryValues.list_data[remappedIDs.list_data[localVariableID]];
736 }
737
738 }
739 }
740
741 //printf ("Node remap at %s yielded %s\n", GetName()->sData, _String((_String*)remapMyCategories.toStr()).sData);
742
743 }
744
745 //_______________________________________________________________________________________________
746
LocateMeInTree(void) const747 node<long>* _CalcNode::LocateMeInTree (void) const {
748
749 _String parentName = ParentObjectName (),
750 myName = ContextFreeName();
751
752 return ((_TreeTopology*)FetchVar(LocateVarByName(parentName)))->FindNodeByName(&myName);
753
754 }
755
756
757 //_______________________________________________________________________________________________
758
ConvertToSimpleMatrix(void) const759 void _CalcNode::ConvertToSimpleMatrix (void) const {
760 _Formula * mf = GetExplicitFormModel();
761 if (mf) {
762 mf->ConvertMatrixArgumentsToSimpleOrComplexForm (false);
763 } else {
764 _Matrix * mm [2] = {GetModelMatrix(), GetFreqMatrix()};
765 for (_Matrix * m : mm) {
766 if (m) {
767 m->MakeMeSimple();
768 }
769 }
770 }
771 }
772
773 //_______________________________________________________________________________________________
774
ConvertFromSimpleMatrix(void)775 void _CalcNode::ConvertFromSimpleMatrix (void) {
776 _Formula * mf = GetExplicitFormModel();
777 if (mf) {
778 mf->ConvertMatrixArgumentsToSimpleOrComplexForm (true);
779 } else {
780 _Matrix * mm [2] = {GetModelMatrix(), GetFreqMatrix()};
781 for (_Matrix * m : mm) {
782 if (m) {
783 m->MakeMeGeneral();
784 }
785 }
786 }
787 }
788
789 //_______________________________________________________________________________________________
790
RecurseMC(long varToConstrain,node<long> * whereAmI,bool first,char rooted)791 _Formula* _CalcNode::RecurseMC (long varToConstrain, node<long>* whereAmI, bool first, char rooted) {
792 // TODO 20171203, SLKP: this needs review
793
794 long descendants = whereAmI->get_num_nodes(),
795 f = iVariables?iVariables->FindStepping(varToConstrain,2,1):-1,
796 start = 0;
797
798 if (f<0 && !first) {
799 HandleApplicationError (_String ("Molecular clock constraint has failed, since variable ")
800 &LocateVar(varToConstrain)->GetName()->Enquote()
801 &" is not an independent member of the node "
802 &GetName()->Enquote()
803 );
804 return nil;
805 }
806
807
808 if (descendants == 0) {
809 if (first) {
810 return nil;
811 } else {
812 return new _Formula (LocateVar(iVariables->get(f-1)),true);
813 }
814 }
815
816 if (first && (!whereAmI->get_parent()) && (rooted == ROOTED_LEFT)) {
817 descendants --;
818 }
819 if (first && (!whereAmI->get_parent()) && (rooted == ROOTED_RIGHT)) {
820 start++;
821 }
822
823 // internal node - must do some work
824
825 _Formula** nodeConditions = new _Formula * [descendants-start];
826
827 for (long k=start+1; k<=descendants; k++) {
828 node<long>* downWeGo = whereAmI->go_down(k);
829 if (!(nodeConditions[k-1-start] = map_node_to_calcnode(downWeGo)->RecurseMC (varToConstrain, downWeGo))) {
830 for (long f2 = 0; f2 < k-start-1; f2++) {
831 delete nodeConditions[f2];
832 }
833
834 delete[] nodeConditions;
835 return nil;
836 }
837 }
838
839 // all the conditions have been written. now check how we should resolve them
840
841 long k;
842
843 for (k=0; k<descendants-start; k++)
844 if ((nodeConditions[k])->Length()>1) {
845 break;
846 }
847
848 if (k==descendants-start) { // all underlying branches are "simple"
849 for (long n=1; n<descendants-start; n++) {
850 //printf ("Setting simple constraint at %s %s\n", LocateVar(nodeConditions[n]->GetIthTerm(0)->GetAVariable())->GetName()->get_str(),
851 // _String((_String*)nodeConditions[0]->toStr(kFormulaStringConversionNormal, nil, true)).get_str());
852 LocateVar (nodeConditions[n]->GetIthTerm(0)->GetAVariable())->SetFormula (*nodeConditions[0]);
853 delete (nodeConditions[n]);
854 nodeConditions[n] = nil;
855 }
856 k = 0;
857 } else {
858 long l;
859 for ( l=k+1; l<descendants-start; l++)
860 if (nodeConditions[l]->Length()>1) {
861 break;
862 }
863
864 if (l==descendants-start) // all but one underlying branches are "simple"
865 for (long n=0; n<descendants-start; n++) {
866 if (n==k) {
867 continue;
868 }
869 //printf ("Setting semi-simple constraint at %s : %s\n", LocateVar(nodeConditions[n]->GetIthTerm(0)->GetAVariable())->GetName()->get_str(),
870 // _String((_String*)nodeConditions[k]->toStr(kFormulaStringConversionNormal, nil, true)).get_str());
871
872 LocateVar (nodeConditions[n]->GetIthTerm(0)->GetAVariable())->SetFormula (*nodeConditions[k]);
873 delete (nodeConditions[n]);
874 nodeConditions[n] = nil;
875 }
876 // really bad bongos! must solve for non-additive constraint
877 else
878 for (long l=0; l<descendants-start; l++) {
879 if (l==k) {
880 continue;
881 }
882 if (nodeConditions[l]->Length()==1) {
883 //printf ("Setting simple non-additive at %s %s\n", LocateVar(nodeConditions[l]->GetIthTerm(0)->GetAVariable())->GetName()->get_str(),
884 // _String((_String*)nodeConditions[k]->toStr(kFormulaStringConversionNormal, nil, true)).get_str());
885
886 LocateVar (nodeConditions[l]->GetIthTerm(0)->GetAVariable())->SetFormula (*nodeConditions[k]);
887 } else { // solve for a non-additive constraint
888 _Variable* nonAdd = LocateVar (nodeConditions[l]->GetIthTerm(0)->GetAVariable());
889 nodeConditions[l]->GetList().Delete(0);
890 _Formula newConstraint;
891 newConstraint.Duplicate(nodeConditions[k]);
892 for (long m=0; m<nodeConditions[l]->GetList().lLength; m++) {
893 _Operation* curOp = (_Operation*)(*nodeConditions[l]).GetList()(m);
894 if (curOp->GetNoTerms()) {
895 newConstraint.GetList().AppendNewInstance(new _Operation (HY_OP_CODE_SUB, 2L));
896 } else {
897 newConstraint.GetList()<<curOp;
898 }
899 }
900 delete (nodeConditions[l]);
901 nodeConditions[l] = nil;
902 nonAdd->SetFormula(newConstraint);
903 //printf ("Setting complex non-additive at %s %s\n", nonAdd->GetName()->get_str(),
904 // _String((_String*)newConstraint.toStr(kFormulaStringConversionNormal, nil, true)).get_str());
905 }
906 }
907 }
908
909
910 if(!first) {
911 _Formula *result = nodeConditions[k];
912 _Operation *newVar = new _Operation;
913 newVar->SetAVariable (iVariables->list_data[f-1]);
914 result->GetList().AppendNewInstance( newVar);
915 result->GetList().AppendNewInstance(new _Operation (HY_OP_CODE_ADD, 2L));
916
917 delete [] nodeConditions;
918 return result;
919 }
920
921 for (long k=0; k<descendants-start; k++)
922 if (nodeConditions[k]) {
923 delete nodeConditions[k];
924 }
925
926 delete [] nodeConditions;
927 return nil;
928 }
929
930 //_______________________________________________________________________________________________
931
ParentTree(void)932 _VariableContainer* _CalcNode::ParentTree(void) {
933 _String parentTree = ParentObjectName();
934 return (_VariableContainer* )FetchObjectFromVariableByType(&parentTree, TREE);
935 }
936
937
938
939
940
941