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
45 #include "global_things.h"
46 #include "global_object_lists.h"
47 #include "tree.h"
48 #include "tree_iterator.h"
49 #include "hbl_env.h"
50 #include "category.h"
51 #include "likefunc.h"
52
53 const _String kTreeErrorMessageEmptyTree ("Cannot construct empty trees");
54
55
56 using namespace hy_global;
57 using namespace hy_env;
58 using namespace hyphy_global_objects;
59
60 #define TREE_V_SHIFT 8.0
61 #define TREE_H_SHIFT 10.0
62
63
64 hyFloat _TheTree::_timesCharWidths[256]= { // Hardcoded relative widths of all 255 characters in the Times font, for the use of PSTreeString
65 0,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0,0.25098,0.721569,0.721569,0.721569,0,0.721569,0.721569,
66 0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0.721569,0,0.721569,0.721569,
67 0.25098,0.333333,0.407843,0.501961,0.501961,0.831373,0.776471,0.180392,0.333333,0.333333,0.501961,0.564706,0.25098,0.333333,0.25098,0.278431,
68 0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.278431,0.278431,0.564706,0.564706,0.564706,0.443137,
69 0.921569,0.721569,0.666667,0.666667,0.721569,0.611765,0.556863,0.721569,0.721569,0.333333,0.388235,0.721569,0.611765,0.890196,0.721569,0.721569,
70 0.556863,0.721569,0.666667,0.556863,0.611765,0.721569,0.721569,0.945098,0.721569,0.721569,0.611765,0.333333,0.278431,0.333333,0.470588,0.501961,
71 0.333333,0.443137,0.501961,0.443137,0.501961,0.443137,0.333333,0.501961,0.501961,0.278431,0.278431,0.501961,0.278431,0.776471,0.501961,0.501961,
72 0.501961,0.501961,0.333333,0.388235,0.278431,0.501961,0.501961,0.721569,0.501961,0.501961,0.443137,0.478431,0.2,0.478431,0.541176,0.721569,
73 0.721569,0.721569,0.666667,0.611765,0.721569,0.721569,0.721569,0.443137,0.443137,0.443137,0.443137,0.443137,0.443137,0.443137,0.443137,0.443137,
74 0.443137,0.443137,0.278431,0.278431,0.278431,0.278431,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,0.501961,
75 0.501961,0.4,0.501961,0.501961,0.501961,0.34902,0.454902,0.501961,0.760784,0.760784,0.980392,0.333333,0.333333,0.54902,0.890196,0.721569,
76 0.713725,0.54902,0.54902,0.54902,0.501961,0.576471,0.494118,0.713725,0.823529,0.54902,0.27451,0.27451,0.309804,0.768627,0.666667,0.501961,
77 0.443137,0.333333,0.564706,0.54902,0.501961,0.54902,0.611765,0.501961,0.501961,1,0.25098,0.721569,0.721569,0.721569,0.890196,0.721569,
78 0.501961,1,0.443137,0.443137,0.333333,0.333333,0.54902,0.494118,0.501961,0.721569,0.168627,0.745098,0.333333,0.333333,0.556863,0.556863,
79 0.501961,0.25098,0.333333,0.443137,1,0.721569,0.611765,0.721569,0.611765,0.611765,0.333333,0.333333,0.333333,0.333333,0.721569,0.721569,
80 0.788235,0.721569,0.721569,0.721569,0.721569,0.278431,0.333333,0.333333,0.333333,0.333333,0.333333,0.333333,0.333333,0.333333,0.333333,0.333333
81 },
82 _TheTree::_maxTimesCharWidth = 0.980392;
83
84 _String const _TheTree::kTreeOutputLabel ( "TREE_OUTPUT_BRANCH_LABEL"),
85 _TheTree::kTreeOutputTLabel ( "TREE_OUTPUT_BRANCH_TLABEL"),
86 _TheTree::kTreeOutputFSPlaceH ( "__FONT_SIZE__");
87
88
89 #define DEGREES_PER_RADIAN 57.29577951308232286465
90
91
92 //#define _UBER_VERBOSE_MX_UPDATE_DUMP
93 //#define _UBER_VERBOSE_MX_UPDATE_DUMP_LF_EVAL 1
94
95 extern long likeFuncEvalCallCount,
96 matrix_exp_count;
97
98
99 hyFloat _lfScalerPower = 64,
100 _lfScalerUpwards = pow(2.,_lfScalerPower),
101 _lfScalingFactorThreshold = 1./_lfScalerUpwards,
102 _logLFScaler = _lfScalerPower *log(2.),
103 _lfMaxScaler = sqrt(DBL_MAX * 1.e-10),
104 _lfMinScaler = 1./_lfMaxScaler;
105
106
107
108 _Vector _scalerMultipliers,
109 _scalerDividers;
110
111 #ifdef _SLKP_DEBUG_
112
113 /*----------------------------------------------------------------------------------------------------------*/
echoNodeList(_SimpleList & theNodes,_SimpleList & leaves,_SimpleList & iNodes,_SimpleList & flatParents)114 void echoNodeList (_SimpleList& theNodes, _SimpleList& leaves, _SimpleList& iNodes, _SimpleList& flatParents) {
115 for (long n = 0; n < theNodes.lLength; n++) {
116 bool is_leaf = theNodes(n)<leaves.lLength;
117 node<long>* nd = (node<long>*)(is_leaf?leaves(theNodes(n)):iNodes(theNodes(n)-leaves.lLength));
118 long children = 0;
119 if (!is_leaf) {
120 long myIndex = theNodes(n)-leaves.lLength;
121 flatParents.Each ([&](long v, unsigned long)->void {if (v == myIndex) children++;});
122 }
123 printf ("%ld %ld %s (%d)\n", n, theNodes(n), LocateVar(nd->in_object)->GetName()->get_str(), children);
124 }
125 }
126
127 #endif
128
129 /*----------------------------------------------------------------------------------------------------------*/
130
_computeReductionScaler(hyFloat currentScaler,hyFloat sum,long & didScale)131 __attribute__((__always_inline__)) hyFloat _computeReductionScaler (hyFloat currentScaler, hyFloat sum, long& didScale) {
132 if (currentScaler > _lfMinScaler) {
133 didScale = -1;
134 sum *= _lfScalingFactorThreshold;
135
136 hyFloat tryScale2 = currentScaler * _lfScalingFactorThreshold,
137 scaler = _lfScalingFactorThreshold;
138
139 while (sum > _lfScalerUpwards && tryScale2 > _lfMinScaler) {
140 sum *= _lfScalingFactorThreshold;
141 tryScale2 *= _lfScalingFactorThreshold;
142 scaler *= _lfScalingFactorThreshold;
143 didScale --;
144 }
145
146 return scaler;
147 }
148 return NAN;
149 }
150
_computeBoostScaler(hyFloat currentScaler,hyFloat sum,long & didScale)151 __attribute__((__always_inline__)) hyFloat _computeBoostScaler (hyFloat currentScaler, hyFloat sum, long& didScale) {
152 if (currentScaler < _lfMaxScaler) {
153 didScale = 1;
154 sum *= _lfScalerUpwards;
155
156 hyFloat tryScale2 = currentScaler * _lfScalerUpwards,
157 scaler = _lfScalerUpwards;
158
159 while (sum < _lfScalingFactorThreshold && tryScale2 < _lfMaxScaler) {
160 sum *= _lfScalerUpwards;
161 tryScale2 *= _lfScalerUpwards;
162 scaler *= _lfScalerUpwards;
163 didScale ++;
164 }
165 return scaler;
166 } /*else {
167 #pragma omp critical
168 printf ("Overflow in _computeBoostScaler %15.12lg\n", currentScaler);
169 }*/
170 return NAN;
171 }
172 /*----------------------------------------------------------------------------------------------------------*/
173
acquireScalerMultiplier(long s)174 hyFloat acquireScalerMultiplier (long s) {
175 if (s>0) {
176 if (s >= _scalerMultipliers.get_used())
177 for (long k = _scalerMultipliers.get_used(); k <= s; k++) {
178 _scalerMultipliers.Store (exp (-_logLFScaler * k));
179 }
180 return _scalerMultipliers.theData[s];
181 }
182 s = -s;
183 if (s >= _scalerDividers.get_used())
184 for (long k = _scalerDividers.get_used(); k <= s; k++) {
185 _scalerDividers.Store (exp (_logLFScaler * k));
186 }
187 return _scalerDividers.theData[s];
188 }
189
190
191 //_______________________________________________________________________________________________
192
_TheTree()193 _TheTree::_TheTree () {
194 categoryCount = 1L;
195 aCache = nil;
196 } // default constructor - doesn't do much
197
198
199
200 //_______________________________________________________________________________________________
201
~_TheTree(void)202 _TheTree::~_TheTree (void) {
203 DeleteObject (aCache);
204
205 }
206
207 //_______________________________________________________________________________________________
208
PurgeTree(void)209 void _TheTree::PurgeTree (void) {
210 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
211
212 while (_CalcNode * iterator = ti.Next()) {
213 DeleteVariable (*iterator->GetName());
214 }
215
216 theRoot->delete_tree (true);
217 }
218
219 //_______________________________________________________________________________________________
_TheTree(_String const & name,_String const & parms,bool make_a_copy)220 _TheTree::_TheTree (_String const & name, _String const & parms, bool make_a_copy):_TreeTopology (&name) {
221 PreTreeConstructor (make_a_copy);
222 _TreeTopologyParseSettings settings = CollectParseSettings();
223 settings.AllocateCache();
224
225 if (_AssociativeList * meta = MainTreeConstructor (parms, settings)) {
226 PostTreeConstructor (make_a_copy, meta);
227 }
228 }
229
230
231 //_______________________________________________________________________________________________
_TheTree(_String const & name,_TreeTopology * top)232 _TheTree::_TheTree (_String const & name, _TreeTopology* top):_TreeTopology (&name) {
233 PreTreeConstructor (false);
234 if (top->theRoot) {
235 isDefiningATree = kTreeIsBeingParsed;
236 theRoot = top->theRoot->duplicate_tree ();
237
238 _TreeTopologyParseSettings parse_settings = _TreeTopology::CollectParseSettings();
239 parse_settings.AllocateCache();
240
241 ConditionalTraverser (
242 [&] (node<long>* iterator, node_iterator<long> const& ni) -> bool {
243 hyFloat stored_branch_length = top->GetBranchLength (iterator);
244
245 _String nodeVS ,
246 nodeName (top->GetNodeName (iterator, false)),
247 *nodeSpec = map_node_to_calcnode(iterator)->GetBranchSpec();
248
249 if (!nodeName.IsValidIdentifier(fIDAllowFirstNumeric)) {
250 nodeName = nodeName.ConvertToAnIdent (fIDAllowFirstNumeric);
251 }
252
253 if (stored_branch_length != 0.0) {
254 nodeVS = stored_branch_length;
255 }
256
257 FinalizeNode (iterator, 0, nodeName, *nodeSpec, nodeVS, parse_settings);
258 DeleteObject (nodeSpec);
259 return false;
260 },
261 true // SLKP 20180309 : TODO check to see if it is necessary to traverse the root
262 );
263
264 isDefiningATree = kTreeNotBeingDefined;
265 PostTreeConstructor (false, nil);
266 } else {
267 HandleApplicationError ("Can't create an empty tree");
268 return;
269 }
270 }
271
272 //_______________________________________________________________________________________________
273
PreTreeConstructor(bool)274 void _TheTree::PreTreeConstructor (bool) {
275 rooted = UNROOTED;
276 categoryCount = 1;
277 aCache = new _AVLListXL (new _SimpleList);
278 }
279
280 //_______________________________________________________________________________________________
281
delete_associated_calcnode(node<long> * n) const282 void _TheTree::delete_associated_calcnode (node<long> * n) const {
283 DeleteVariable(*map_node_to_calcnode(n)->GetName());
284 }
285
286
287 //_______________________________________________________________________________________________
288
PostTreeConstructor(bool make_copy,_AssociativeList * meta)289 void _TheTree::PostTreeConstructor (bool make_copy, _AssociativeList* meta) {
290 // TODO: SLKP 20180311, extensive duplication with the same function from _TreeToology, consider
291 // a possible refactor
292
293 if (aCache) {
294 DeleteObject (aCache->dataList);
295 DeleteAndZeroObject(aCache);
296 }
297
298 auto variable_handler = [&] (void) -> void {
299 /** TODO SLKP 20171211, make sure the semantics are unchanged */
300 // existing variable is a CalcNode
301 variablePtrs.Replace (get_index(), make_copy ? this->makeDynamic() : this, false);
302 setParameter (WrapInNamespace (_TreeTopology::kMeta, GetName()), meta ? meta : new _MathObject, nil, false);
303 //printf ("makecopy = %d [%ld]\n", make_copy, this->SingleReference());
304 };
305
306 bool accept_rooted = EnvVariableTrue(accept_rooted_trees);
307 try {
308 if (theRoot->get_num_nodes() <= 2) { // rooted tree - check
309 if (accept_rooted == false) {
310
311 //long node_index = theRoot->get_data();
312 bool recurse = false;
313
314 if (theRoot->get_num_nodes() == 2) {
315 for (int i = 1; i<=2; i++) {
316 node<long> *node_temp = theRoot->go_down(i);
317 if (node_temp->get_num_nodes()) { // an internal node - make it a root
318 delete_associated_calcnode(theRoot);
319 node_temp->detach_parent();
320 node_temp->add_node(*theRoot->go_down(3-i));
321 delete theRoot;
322 theRoot = node_temp;
323 //delete_associated_calcnode (theRoot);
324 rooted = i == 1 ? ROOTED_LEFT : ROOTED_RIGHT;
325 ReportWarning (_String("Rooted topology. Removing one branch - the ") & (i==1 ? "left" : "right") & " root child has been promoted to be the new root");
326 break;
327 }
328 }
329
330 if (rooted==UNROOTED) {
331 ReportWarning ("One branch tree supplied - hopefully this IS what you meant to do.");
332 node<long> *delete_me = theRoot->go_down(2);
333 //delete_associated_calcnode(theRoot);
334 delete_me->detach_parent();
335 theRoot->detach_child(2);
336 delete_associated_calcnode(delete_me);
337 delete delete_me;
338 rooted = ROOTED_LEFT;
339 ReportWarning ("RIGHT child collapsed to the root");
340 //delete_associated_calcnode(theRoot);
341 }
342 } else {
343 if (theRoot->get_num_nodes() == 0) {
344 //delete this;
345 throw (kTreeErrorMessageEmptyTree);
346 }
347 node<long> *node_temp = theRoot->go_down(1);
348 node_temp->detach_parent();
349 delete_associated_calcnode(theRoot);
350 delete theRoot;
351 theRoot = node_temp;
352 ReportWarning ("The root has a single child, which is be promoted to the root");
353 recurse = true;
354 }
355
356 if (recurse) {
357 PostTreeConstructor (make_copy, meta);
358 return;
359 }
360 }
361 }
362 if (!theRoot) {
363 //delete this;
364 throw _String ("Invalid tree/topology string specification.");
365 }
366 } catch (const _String& e) {
367 HandleApplicationError (e);
368 }
369
370 variable_handler ();
371 }
372
373 //_______________________________________________________________________________________________
_TheTree(_String const & name,_TheTree * otherTree)374 _TheTree::_TheTree (_String const &name, _TheTree* otherTree):_TreeTopology (&name) {
375 PreTreeConstructor (false);
376 if (otherTree->theRoot) {
377 isDefiningATree = kTreeIsBeingParsed;
378 theRoot = otherTree->theRoot->duplicate_tree ();
379
380 node_iterator<long> tree_iterator (theRoot, _HY_TREE_TRAVERSAL_POSTORDER);
381 while (node<long>*topTraverser = tree_iterator.Next()) {
382 _CalcNode *sourceNode = (_CalcNode*)LocateVar(topTraverser->in_object),
383 copiedNode (sourceNode, this);
384 topTraverser->init (copiedNode.theIndex);
385 }
386
387 isDefiningATree = kTreeNotBeingDefined;
388 PostTreeConstructor (false, nil);
389 } else {
390 HandleApplicationError(kTreeErrorMessageEmptyTree);
391 return;
392 }
393 }
394
395 //_______________________________________________________________________________________________
396
makeDynamic(void) const397 BaseRef _TheTree::makeDynamic (void) const {
398 _TheTree* res = new _TheTree;
399 res->_CalcNode::Duplicate (this);
400 res->rooted = rooted;
401 res->categoryCount = 1;
402 res->theRoot = CopyTreeStructure (theRoot,true);
403 return res;
404 }
405
406
407 //_______________________________________________________________________________________________
408
makeDynamicCopy(_String const * replacementName) const409 BaseRef _TheTree::makeDynamicCopy (_String const* replacementName) const {
410 _TheTree* res = new _TheTree;
411
412 res->rooted = rooted;
413 if (theRoot) {
414 _String rn = *replacementName&'.';
415 res->theRoot = DuplicateTreeStructure (theRoot, &rn, true);
416 } else {
417 res->theRoot = nil;
418 }
419
420 res->SetIndex(variableNames.GetXtra(LocateVarByName (*replacementName)));
421 res->theName = new _String (*replacementName);
422 res->theName->AddAReference();
423 return res;
424 }
425
426 //_______________________________________________________________________________________________
427
DuplicateTreeStructure(node<long> * source_node_object,_String const * new_name,bool) const428 node<long>* _TheTree::DuplicateTreeStructure (node <long>* source_node_object, _String const* new_name, bool) const {
429
430 // TODO SLKP 20180311: this needs to be reviewed; lots of old ugly code
431
432 node<long>* locNode = new node<long>;
433
434 for (unsigned long i=1UL; i <= source_node_object->get_num_nodes(); i++) {
435 locNode->add_node(*DuplicateTreeStructure (source_node_object->go_down(i), new_name , false));
436 }
437
438 _String replace_me = *GetName()&'.',
439 *temp;
440
441 _CalcNode* sourceNode = (_CalcNode*)map_node_to_calcnode(source_node_object)->makeDynamic();
442 _String new_node_name = sourceNode->GetName()->Replace(replace_me, *new_name, false);
443 _Variable node_var (new_node_name);
444 InsertVar (&node_var);
445 locNode->init (node_var.get_index());
446
447 sourceNode->ForEachLocalVariable(sourceNode->iVariables, [&] (long var_idx, long ref_idx, long array_index) {
448 _Variable new_node_var (LocateVar (var_idx)->GetName()->Replace(replace_me, *new_name, false));
449 (*sourceNode->iVariables)[array_index] = new_node_var.get_index();
450 });
451
452 sourceNode->ForEachLocalVariable(sourceNode->dVariables, [&] (long var_idx, long ref_idx, long array_index) {
453 _Variable new_node_var (LocateVar (var_idx)->GetName()->Replace(replace_me, *new_name, false));
454 (*sourceNode->dVariables)[array_index] = new_node_var.get_index();
455 _Variable * v = LocateVar (new_node_var.get_index());
456 //_Formula const *existing_constraint = LocateVar (var_idx)->get_constraint();
457 // TODO : this is not a robust way to rename variables
458 _String * existing_formula = LocateVar (var_idx)->GetFormulaString(kFormulaStringConversionNormal);
459 *existing_formula = existing_formula->Replace (replace_me, *new_name,true);
460 _Formula new_constraint (*existing_formula);
461 v->SetFormula (new_constraint);
462 DeleteObject (existing_formula);
463 });
464
465 return locNode;
466 }
467
468 //_______________________________________________________________________________________________
469
_MainTreeConstructor_error(const _String & error,const _String & tree_string,unsigned long index)470 bool _MainTreeConstructor_error (const _String& error, const _String& tree_string, unsigned long index) {
471 isDefiningATree = kTreeNotBeingDefined;
472 HandleApplicationError ( error & ", in the following string context " &
473 tree_string.Cut(index>31L?index-32L:0L,index)&
474 " <ERROR HERE> "&
475 tree_string.Cut(index+1L,tree_string.length()-index>32L?index+32L:-1L)
476 );
477
478 return false;
479 }
480
481
482 //_______________________________________________________________________________________________
483
484
FinalizeNode(node<long> * nodie,long number,_String nodeName,_String const & nodeParameters,_String & nodeValue,_TreeTopologyParseSettings const & settings)485 _String _TheTree::FinalizeNode (node<long>* nodie, long number , _String nodeName, _String const& nodeParameters, _String& nodeValue, _TreeTopologyParseSettings const& settings)
486 {
487
488 static const _String kCommentSuffix ("_comment");
489
490 bool isAutoGenerated = nodeName.empty() || (settings.ingore_user_inode_names && nodie->get_num_nodes()>0);
491
492 if (isAutoGenerated) {
493 nodeName = settings.inode_prefix & number;
494 } else {
495 if (!nodeName.IsValidIdentifier(fIDAllowFirstNumeric)) {
496 _String new_name = nodeName.ConvertToAnIdent(fIDAllowFirstNumeric);
497 ReportWarning (_String ("Automatically renamed ") & nodeName.Enquote() & " to " & new_name.Enquote() & " in order to create a valid HyPhy identifier");
498 nodeName = new_name;
499 }
500 }
501
502 _String node_parameters = nodeParameters;
503
504 if (nodie == theRoot) {
505 node_parameters = kEmptyString;
506 nodeValue = kEmptyString;
507 } else {
508 if (node_parameters.empty() && lastMatrixDeclared !=-1L ) {
509 node_parameters=* hyphy_global_objects::GetObjectNameByType (HY_BL_MODEL, lastMatrixDeclared, false);
510 }
511
512 if (node_parameters.nonempty()) {
513 ReportWarning (_StringBuffer (128UL) << "Model " << '"' << node_parameters << '"' << " assigned to \"" << nodeName << '"');
514 } else {
515 ReportWarning (_String("No nodel was assigned to ")& nodeName.Enquote());
516 }
517 }
518
519 hyTreeDefinitionPhase saveIDT = isDefiningATree;
520 isDefiningATree = kTreeNodeBeingCreated;
521 _CalcNode cNt (nodeName,node_parameters, 4, this, aCache);
522 isDefiningATree = saveIDT;
523 nodie->init (cNt.get_index());
524
525 _Constant val (ProcessTreeBranchLength(nodeValue));
526
527 if (nodeValue.nonempty() && settings.accept_user_lengths) {
528 if (cNt.CountIndependents () == 1UL) { // can assign default values
529 bool use_direct_value = true;
530 if (settings.auto_convert_lengths) {
531 long nodeModelID = cNt.GetModelIndex();
532 if (nodeModelID != HY_NO_MODEL && !IsModelOfExplicitForm(nodeModelID)) {
533 _Formula * expressionToSolveFor = nil;
534 long already_converted = settings.parser_cache->FindLong (nodeModelID);
535
536 if (already_converted < 0) {
537 _Variable * tV, * tV2;
538 bool mByF;
539 RetrieveModelComponents (nodeModelID, tV, tV2, mByF);
540 _String * result = ((_Matrix*)tV->GetValue())->BranchLengthExpression((_Matrix*)tV2->GetValue(),mByF);
541 if (result->nonempty()) {
542 expressionToSolveFor = new _Formula (*result);
543 for (unsigned long cc = 0; cc < cNt.categoryVariables.countitems(); cc++) {
544 _CategoryVariable * thisCC = cNt.get_ith_category (cc);
545 thisCC -> SetValue (new _Constant(thisCC->Mean()), false, true, NULL);
546 }
547 settings.parser_cache->Insert ((BaseRef)nodeModelID, (long)expressionToSolveFor, false, false);
548 }
549 DeleteObject (result);
550 } else {
551 expressionToSolveFor = (_Formula*)settings.parser_cache->GetXtra(already_converted);
552 }
553
554 if (expressionToSolveFor != nil) {
555 _Variable * solveForMe = LocateVar (cNt.iVariables->list_data[1]);
556 hyFloat modelP = expressionToSolveFor->Brent (solveForMe,solveForMe->GetLowerBound(), solveForMe->GetUpperBound(), 1e-6, nil, val.Value());
557 ReportWarning (_String("Branch parameter of ") & nodeName.Enquote() &" set to " & modelP);
558 cNt.GetIthIndependent(0) ->SetValue(new _Constant (modelP), false,true, NULL);
559 use_direct_value = false;
560 }
561 }
562 }
563
564 if (use_direct_value) {
565 cNt.GetIthIndependent(0)->SetValue (&val,true,true,NULL);
566 ReportWarning (_String("Branch parameter of ") & nodeName&" set to " & nodeValue);
567 }
568 } else {
569 ReportWarning (_StringBuffer (128) << nodeName << " has " << _String ((long)(cNt.iVariables?cNt.iVariables->lLength/2:0)) << " parameters - branch length not assigned");
570 }
571 }
572
573 _CalcNode *nodeVar = (_CalcNode*)LocateVar(cNt.get_index());
574
575 if (nodeVar == NULL) return kEmptyString;
576
577 nodeVar->SetValue (&val,true,true,NULL);
578
579 node_parameters.Clear();
580 nodeValue.Clear();
581
582 nodeVar->categoryVariables.TrimMemory();
583 nodeVar->categoryIndexVars.TrimMemory();
584 nodeVar->_VariableContainer::TrimMemory();
585
586 return nodeName;
587 }
588
589 //_______________________________________________________________________________________________
590
GetNodeModel(node<long> * n) const591 _String const* _TheTree::GetNodeModel (node<long>* n) const {
592 return map_node_to_calcnode(n)->GetModelName ();
593 }
594
595 //_______________________________________________________________________________________________
596
GetNodeName(node<long> * n,bool fullName) const597 const _String _TheTree::GetNodeName (node<long>* n, bool fullName) const {
598 return GetNodeName(map_node_to_calcnode(n), fullName);
599 }
600
601 //_______________________________________________________________________________________________
602
GetNodeName(_CalcNode const * cn,bool fullName) const603 const _String _TheTree::GetNodeName (_CalcNode const* cn, bool fullName) const {
604 if (fullName) {
605 return *cn->GetName();
606 }
607 return cn->GetName()->Cut (GetName()->length()+1L,kStringEnd);
608 }
609
610
611 //_______________________________________________________________________________________________
612
toStr(unsigned long)613 BaseRef _TheTree::toStr (unsigned long) {
614 return _TreeTopology::toStr();
615 }
616
617
618 //__________________________________________________________________________________
619
CompileListOfModels(_SimpleList & l)620 void _TheTree::CompileListOfModels (_SimpleList& l) {
621 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
622 while (_CalcNode* iterator = ti.Next()) {
623 long modelID = iterator->GetModelIndex();
624 if (modelID != HY_NO_MODEL && l.Find(modelID) == -1 ) {
625 l << modelID;
626 }
627 }
628 }
629
630 //_______________________________________________________________________________________________
SetCompMatrices(long catID) const631 void _TheTree::SetCompMatrices (long catID) const {
632 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot);
633 while (_CalcNode* iterator = ti.Next()) {
634 iterator->SetCompMatrix (catID);
635 }
636 }
637
638 //__________________________________________________________________________________
639
SetUp(void)640 void _TheTree::SetUp (void) {
641
642 flatTree.Clear();
643 flatNodes.Clear();
644 flatLeaves.Clear();
645 flatCLeaves.Clear();
646 flatParents.Clear();
647
648 _SimpleList flatINodeParents;
649
650 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
651
652 while (_CalcNode* iterator = ti.Next()) {
653 if (ti.IsAtLeaf()) {
654 flatCLeaves << iterator;
655 flatLeaves << (long)ti.GetNode();
656 flatParents << (long)ti.GetNode()->get_parent();
657 } else {
658 flatTree<<iterator;
659 flatNodes<< (long)ti.GetNode();
660 flatINodeParents << (long)ti.GetNode()->get_parent();
661 }
662
663 }
664
665 flatParents << flatINodeParents;
666
667 _SimpleList parentlist (flatNodes),
668 indexer (flatNodes.lLength,0,1);
669
670 SortLists (&parentlist,&indexer);
671 flatParents.Each ([&] (long value, unsigned long index) -> void {
672 if (value) {
673 flatParents[index] = indexer.get(parentlist.BinaryFind(value));
674 } else {
675 flatParents[index] = -1L;
676 }
677 });
678
679 /*for (unsigned long k=0UL; k<flatParents.countitems(); k++) {
680 if (flatParents.get(k)) { // not root
681 flatParents[k] = indexer.get(parentlist.BinaryFind(flatParents.get(k)));
682 } else {
683 flatParents[k] = -1L;
684 }
685 }*/
686
687 }
688
689 //__________________________________________________________________________________
690
AllBranchesHaveModels(long matchSize) const691 bool _TheTree::AllBranchesHaveModels (long matchSize) const {
692 // TODO SLKP 20180313: possible deprecation
693
694 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot);
695
696 while (_CalcNode* iterator = ti.Next()) {
697 _Matrix * model = iterator->GetModelMatrix();
698 if (! (model && (matchSize < 0L || model->GetHDim()==matchSize))) {
699 return false;
700 }
701 }
702
703 return true;
704 }
705
706
707 //__________________________________________________________________________________
708
ExecuteSingleOp(long opCode,_List * arguments,_hyExecutionContext * context,HBLObjectRef cache)709 HBLObjectRef _TheTree::ExecuteSingleOp (long opCode, _List* arguments, _hyExecutionContext* context, HBLObjectRef cache) {
710
711 switch (opCode) {
712 case HY_OP_CODE_TYPE: // Type
713 return Type(cache);
714 break;
715 }
716
717 _MathObject * arg0 = _extract_argument (arguments, 0UL, false);
718
719 if (arg0) {
720 switch (opCode) {
721 case HY_OP_CODE_TEXTREESTRING: // TEXTreeString
722 return TEXTreeString(arg0, cache);
723 }
724
725 _MathObject * arg1 = _extract_argument (arguments, 1UL, false);
726
727 if (arg1) {
728 switch (opCode) {
729 case HY_OP_CODE_PSTREESTRING: //PlainTreeString
730 return PlainTreeString(arg0,arg1, cache);
731 }
732 }
733 }
734
735 switch (opCode) {
736 case HY_OP_CODE_TEXTREESTRING:
737 case HY_OP_CODE_PSTREESTRING:
738 WarnWrongNumberOfArguments (this, opCode,context, arguments);
739 return new _MathObject;
740 }
741
742 return _TreeTopology::ExecuteSingleOp (opCode,arguments,context);
743
744 }
745
746
747
748
749 //__________________________________________________________________________________
GetBranchLengthString(node<long> * n,bool get_expression) const750 _String const _TheTree::GetBranchLengthString (node<long> * n, bool get_expression) const {
751
752 _CalcNode* tree_node = map_node_to_calcnode(n);
753
754 if (get_expression) {
755
756 bool mbf;
757
758 _Variable *mm,
759 *fv;
760
761 RetrieveModelComponents(tree_node->GetModelIndex(), mm, fv, mbf);
762
763 if (mm && fv && mm->ObjectClass() == MATRIX && fv->ObjectClass() == MATRIX) {
764 return _String (((_Matrix*)mm->GetValue())->BranchLengthExpression((_Matrix*)fv->GetValue(),mbf));
765 } else {
766 return kEmptyString;
767 }
768
769 } else {
770 return tree_node->ComputeBranchLength();
771 }
772 }
773
774 //__________________________________________________________________________________
GetBranchLength(node<long> * n) const775 hyFloat _TheTree::GetBranchLength (node<long> * n) const {
776 return map_node_to_calcnode(n)->ComputeBranchLength();
777 }
778
779
780 //__________________________________________________________________________________
GetBranchValue(node<long> * n) const781 _String const _TheTree::GetBranchValue (node<long> * n) const {
782 hyFloat t = map_node_to_calcnode(n)->Value();
783 if (!CheckEqual(t, -1)) {
784 return t;
785 } else {
786 return kEmptyString;
787 }
788 }
789
790
791
792 //__________________________________________________________________________________
GetBranchVarValue(node<long> * n,long idx) const793 _String const _TheTree::GetBranchVarValue (node<long> * n, long idx) const {
794 _CalcNode * tree_node = map_node_to_calcnode(n);
795
796 long iVarValue = tree_node->iVariables->FindStepping(idx,2,1);
797 if (iVarValue>0) {
798 // model parameter
799 return _String(LocateVar (tree_node->iVariables->get(iVarValue-1))->Value());
800 } else {
801 // user parameter
802 _String query = _String('.') & *LocateVar(idx)->GetName();
803
804 long local_idx = tree_node->FindLocalVariable (tree_node->iVariables, [&] (long local, long template_var, unsigned long i) -> bool {
805 return LocateVar (local)->GetName()->EndsWith(query);
806 });
807
808 if (local_idx != kNotFound) {
809 return tree_node->GetIthIndependent(local_idx)->GetName();
810 }
811 }
812 return kEmptyString;
813 }
814
815
816 //__________________________________________________________________________________
817
AlignNodes(node<nodeCoord> * theNode) const818 void _TheTree::AlignNodes (node<nodeCoord>* theNode) const {
819 long k = theNode->get_num_nodes();
820 if (k) {
821 theNode->in_object.v = (theNode->go_down(1)->in_object.v+theNode->go_down(k)->in_object.v)/2.0;
822 theNode->in_object.h = 0;
823 for (; k; k--) {
824 hyFloat t = theNode->go_down(k)->in_object.h;
825 if (t<theNode->in_object.h) {
826 theNode->in_object.h = t;
827 }
828 }
829 theNode->in_object.h -= TREE_H_SHIFT;
830 } else {
831 theNode->in_object.v = 0;
832 theNode->in_object.h = 0;
833 }
834 }
835 //__________________________________________________________________________________
836
AlignedTipsMapping(node<long> * iterator,hyFloat & current_offset,bool first,bool respectRoot) const837 node<nodeCoord>* _TheTree::AlignedTipsMapping (node<long>* iterator, hyFloat& current_offset, bool first, bool respectRoot) const{
838 // 20180615 TODO: SLKP this needs review and possibly deprecation
839 if (first) {
840 current_offset = 0.0;
841 long descendants = theRoot->get_num_nodes();
842 node<nodeCoord>* aRoot = new node<nodeCoord>; // reslove rootedness here
843 aRoot->in_object.varRef = -1;
844 if (rooted == UNROOTED || !respectRoot) {
845 for (long k=1L; k<=descendants; k++) {
846 aRoot->add_node(*AlignedTipsMapping(theRoot->go_down(k), current_offset));
847 }
848 AlignNodes (aRoot);
849 return aRoot;
850 } else {
851 node<nodeCoord>* aChild = new node<nodeCoord>;
852 aChild->in_object.varRef = -1;
853 if (rooted == ROOTED_LEFT) {
854 aRoot->add_node (*aChild);
855 for (long k=1L; k<descendants; k++) {
856 aChild->add_node(*AlignedTipsMapping(theRoot->go_down(k), current_offset));
857 }
858 aRoot->add_node(*AlignedTipsMapping(theRoot->go_down(descendants), current_offset));
859 } else {
860 aRoot->add_node(*AlignedTipsMapping(theRoot->go_down(1), current_offset));
861 for (long k=2L; k<=descendants; k++) {
862 aChild->add_node(*AlignedTipsMapping(theRoot->go_down(k), current_offset));
863 }
864 aRoot->add_node (*aChild);
865 }
866 AlignNodes (aChild);
867 AlignNodes (aRoot);
868 return aRoot;
869 }
870 } else {
871 node<nodeCoord>*aNode = new node<nodeCoord>;
872 long descendants = iterator->get_num_nodes();
873 for (long k=1L; k<=descendants; k++) {
874 aNode->add_node(*AlignedTipsMapping(iterator->go_down(k), current_offset));
875 }
876 if (!descendants) { // terminal node
877 aNode->in_object.v = current_offset;
878 aNode->in_object.h = 0;
879 current_offset+=TREE_V_SHIFT;
880 } else {
881 AlignNodes (aNode);
882 }
883 aNode->in_object.varRef = iterator->in_object;
884 return aNode;
885 }
886 }
887
888
889 //__________________________________________________________________________________
890
ScaledBranchReMapping(node<nodeCoord> * theNode,hyFloat tw) const891 void _TheTree::ScaledBranchReMapping (node<nodeCoord>* theNode, hyFloat tw) const
892 {
893 theNode->in_object.h -= tw;
894 for (long k=1; k<=theNode->get_num_nodes(); k++) {
895 ScaledBranchReMapping (theNode->go_down(k), tw);
896 }
897 }
898
899 //__________________________________________________________________________________
900
DetermineBranchLengthGivenScalingParameter(long varRef,_String & matchString,hyTopologyBranchLengthMode mapMode) const901 hyFloat _TheTree::DetermineBranchLengthGivenScalingParameter (long varRef, _String& matchString, hyTopologyBranchLengthMode mapMode) const {
902 if (mapMode == kTopologyBranchLengthNone) { // was 3
903 return 1.;
904 }
905
906 _CalcNode * tree_node = (_CalcNode*)LocateVar(varRef);
907
908 hyFloat branchLength = HY_REPLACE_BAD_BRANCH_LENGTH_WITH_THIS;
909
910 if (mapMode==kTopologyBranchLengthExpectedSubs) { // was 1
911 return tree_node->ComputeBranchLength();
912 } else if (mapMode==kTopologyBranchLengthUserLengths) { // was 2
913 branchLength = tree_node->Value();
914 if (branchLength<=0.0) {
915 branchLength = HY_REPLACE_BAD_BRANCH_LENGTH_WITH_THIS;
916 }
917 } else {
918 auto match_pattern = [&] (long local_idx, long template_index, unsigned long index) -> void {
919 _Variable * local_parameter = LocateVar(local_idx);
920 if (local_parameter->GetName()->EndsWith(matchString)) {
921 branchLength = local_parameter->Compute()->Value();
922 if (branchLength <= 0.) {
923 branchLength = HY_REPLACE_BAD_BRANCH_LENGTH_WITH_THIS;
924 }
925 throw (0);
926 }
927 };
928
929 try {
930 tree_node->ForEachLocalVariable(tree_node->iVariables, match_pattern);
931 tree_node->ForEachLocalVariable(tree_node->dVariables, match_pattern);
932 } catch (int) {
933 return branchLength;
934 }
935 }
936 return branchLength;
937 }
938
939
940 //__________________________________________________________________________________
941
ScaledBranchMapping(node<nodeCoord> * theParent,_String * scalingParameter,long locDepth,long & depth,hyTopologyBranchLengthMode mapMode) const942 node<nodeCoord>* _TheTree::ScaledBranchMapping (node<nodeCoord>* theParent, _String* scalingParameter, long locDepth, long& depth, hyTopologyBranchLengthMode mapMode) const{
943 // 20180615 TODO: SLKP this needs further review and possibly deprecation
944
945 // run a pass of aligned tip mapping then perform one more pass from the root to the children
946 // pre-order to remap the length of branches.
947
948 static hyFloat treeWidth;
949 bool wasRoot = !theParent;
950
951 if (!theParent) { // init
952 hyFloat vertical_offset;
953 theParent = AlignedTipsMapping (theRoot, vertical_offset, true, true);
954 theParent->in_object.h = 0.0;
955 treeWidth = 0.;
956 }
957
958 node<nodeCoord>* currentN;
959 long descendants = theParent->get_num_nodes(),
960 k = 1,
961 j,
962 b = -1;
963
964 hyFloat branchLength = HY_REPLACE_BAD_BRANCH_LENGTH_WITH_THIS;
965
966 for (; k<=descendants; k++) {
967 currentN = theParent->go_down(k);
968 j = currentN->in_object.varRef;
969
970 if (j>=0) {
971
972 branchLength = currentN->in_object.bL = DetermineBranchLengthGivenScalingParameter(j,*scalingParameter,mapMode);
973 branchLength += theParent->in_object.h;
974
975 StoreIfGreater(treeWidth, branchLength);
976
977 theParent->go_down (k)->in_object.h = branchLength;
978 ScaledBranchMapping (theParent->go_down(k), scalingParameter, locDepth+1, depth, mapMode);
979
980 } else {
981 theParent->go_down (k)->in_object.h = 0;
982 b = k;
983 }
984
985 }
986
987 if (k==descendants+1) {
988 if (locDepth>=depth) {
989 depth = locDepth+1;
990 }
991 }
992
993 if (wasRoot) {
994 if (b>0 && descendants==2) {
995 j = (b==1)? 2 : 1;
996
997 ScaledBranchReMapping (theParent->go_down(j),0);
998 theParent->go_down(b)->in_object.h = 0;
999 ScaledBranchMapping (theParent->go_down(b), scalingParameter, locDepth, depth, mapMode);
1000 }
1001
1002 ScaledBranchReMapping (theParent, treeWidth);
1003 return theParent;
1004 }
1005 return nil;
1006 }
1007
1008 //__________________________________________________________________________________
1009
RadialBranchMapping(node<long> * referenceNode,node<nodeCoord> * parentNode,_String * scalingParameter,hyFloat anglePerTip,long & currentTipID,hyFloat & maxRadius,hyTopologyBranchLengthMode mapMode)1010 node<nodeCoord>* _TheTree::RadialBranchMapping (node<long>* referenceNode, node<nodeCoord>* parentNode, _String* scalingParameter, hyFloat anglePerTip, long& currentTipID, hyFloat& maxRadius, hyTopologyBranchLengthMode mapMode)
1011 {
1012 // 20180618 TODO: SLKP this needs review and possibly deprecation
1013
1014 // label 1 stores current radial distance from the root
1015 // label 2 stores the angle of the line to this node
1016 // h and v store the Cartesian coordinates
1017
1018 node <nodeCoord>* current_node = new node <nodeCoord>;
1019
1020 hyFloat branchL = 0.,
1021 referenceL = 0.;
1022
1023 if (parentNode == nil) {
1024 current_node->in_object.label1 = 0.0;
1025 current_node->in_object.label2 = 0.0;
1026 } else {
1027 referenceL = parentNode->in_object.label1;
1028 branchL = DetermineBranchLengthGivenScalingParameter(referenceNode->in_object,*scalingParameter,mapMode);
1029 }
1030
1031 long children = referenceNode->get_num_nodes();
1032
1033 current_node->in_object.label1 = referenceL + branchL;
1034 if (children == 0) {
1035 current_node->in_object.label2 = anglePerTip * currentTipID++;
1036 } else {
1037 hyFloat angleSum = 0.;
1038 for (long n = 1; n <= children; n++) {
1039 node<nodeCoord>* newChild = RadialBranchMapping (referenceNode->go_down(n), current_node, scalingParameter, anglePerTip, currentTipID, maxRadius, mapMode);
1040 current_node->add_node(*newChild);
1041 angleSum += newChild->in_object.label2;
1042 }
1043 current_node->in_object.label2 = angleSum / children;
1044 }
1045
1046 current_node->in_object.h = current_node->in_object.label1*cos(current_node->in_object.label2);
1047 current_node->in_object.v = current_node->in_object.label1*sin(current_node->in_object.label2);
1048 maxRadius = MAX(maxRadius, current_node->in_object.label1);
1049 current_node->in_object.varRef = referenceNode->in_object;
1050 current_node->in_object.bL = branchL;
1051
1052 return current_node;
1053 }
1054
1055
1056
1057 //__________________________________________________________________________________
1058
AssignLabelsToBranches(node<nodeCoord> * theParent,_String * scalingParameter,bool below)1059 void _TheTree::AssignLabelsToBranches (node<nodeCoord>* theParent, _String* scalingParameter, bool below) {
1060 // 20180618 TODO: SLKP this needs review and possibly deprecation
1061
1062 bool wasRoot = !(theParent->parent);
1063 long descendants = theParent->get_num_nodes(),b=-1;
1064
1065 hyFloat branchLength = HY_REPLACE_BAD_BRANCH_LENGTH_WITH_THIS;
1066
1067 hyTopologyBranchLengthMode mapMode;
1068 _String matchString = DetermineBranchLengthMappingMode(scalingParameter, mapMode);
1069
1070 for (long k = 1; k<=descendants; k++) {
1071 node<nodeCoord>* currentN = theParent->go_down(k);
1072 long j = currentN->in_object.varRef;
1073 if (j>=0) {
1074 branchLength = DetermineBranchLengthGivenScalingParameter(j, matchString, mapMode);
1075 if (below) {
1076 currentN->in_object.label2 = branchLength;
1077 } else {
1078 currentN->in_object.label1 = branchLength;
1079 }
1080 AssignLabelsToBranches (theParent->go_down(k), scalingParameter, below);
1081 } else {
1082 if (below) {
1083 currentN->in_object.label2 = 0.;
1084 } else {
1085 currentN->in_object.label1 = 0.;
1086 }
1087 b = k;
1088 AssignLabelsToBranches (theParent->go_down(k), scalingParameter, below);
1089 }
1090
1091 }
1092
1093 if (wasRoot) {
1094 if ( b>0 && descendants==2 ) {
1095
1096 b = b == 1 ? 2 : 1;
1097
1098 if (below) {
1099 theParent->in_object.label2 = theParent->go_down(b)->in_object.label2/2.;
1100 theParent->go_down(b)->in_object.label2/=2.;
1101 } else {
1102 theParent->in_object.label1 = theParent->go_down(b)->in_object.label1/2.;
1103 theParent->go_down(b)->in_object.label1/=2.;
1104 }
1105 }
1106 }
1107 }
1108
1109 //__________________________________________________________________________________
1110
PSStringWidth(_String const & s)1111 hyFloat _TheTree::PSStringWidth (_String const& s) {
1112
1113 hyFloat nnWidth = 0.;
1114
1115 s.Each([&] (unsigned char cc, unsigned long) -> void {
1116 nnWidth += _timesCharWidths [cc];
1117 });
1118
1119 return nnWidth;
1120 }
1121
1122 //__________________________________________________________________________________
PlainTreeString(HBLObjectRef p,HBLObjectRef p2,HBLObjectRef cache)1123 HBLObjectRef _TheTree::PlainTreeString (HBLObjectRef p, HBLObjectRef p2, HBLObjectRef cache) {
1124 // 20180618 TODO: SLKP this needs review and possibly deprecation
1125
1126
1127 static const _String kTreeOutputEmbed ("TREE_OUTPUT_EMBED"),
1128 kTreeOutputOptions ("TREE_OUTPUT_OPTIONS"),
1129 kTreeOutputBackground ( "TREE_OUTPUT_BACKGROUND"),
1130 kTreeOutputRightMargin ( "TREE_OUTPUT_RIGHT_MARGIN"),
1131 kTreeOutputXtraMargin ( "TREE_OUTPUT_XTRA_MARGIN"),
1132 kTreeOutputSymbols ( "TREE_OUTPUT_SYMBOLS"),
1133 kTreeOutputSymbolSize ( "TREE_OUTPUT_SYMBOL_SIZE"),
1134 kTreeOutputExtraPS ( "TREE_OUTPUT_EXTRA_POSTSCRIPT"),
1135 kTreeOutputPrefixPS ( "TREE_OUTPUT_PREFIX_POSTSCRIPT"),
1136 kTreeOutputLayout ( "TREE_OUTPUT_LAYOUT");
1137
1138
1139 auto computeChordLength = [] (hyFloat l, hyFloat angle, hyFloat* maxCoord = nil) -> hyFloat {
1140 hyFloat sinV = sin(angle),
1141 cosV = cos(angle);
1142
1143 if (maxCoord) {
1144 maxCoord[0] = MAX (maxCoord[0], cosV*l);
1145 maxCoord[1] = MIN (maxCoord[1], cosV*l);
1146 maxCoord[2] = MAX (maxCoord[2], sinV*l);
1147 maxCoord[3] = MIN (maxCoord[3], sinV*l);
1148 }
1149
1150 return l/MAX(fabs(sinV),fabs(cosV));
1151 };
1152
1153 _TreeTopologyParseSettings parse_settings = _TreeTopology::CollectParseSettings();
1154
1155 _StringBuffer * res = new _StringBuffer (512UL);
1156 long treeLayout = 0;
1157
1158 if (p&&(p->ObjectClass()==STRING)) {
1159 if (p2&&(p2->ObjectClass()==MATRIX)) {
1160 node<nodeCoord>* newRoot,
1161 *currentNd;
1162
1163 bool doEmbed = false;
1164 bool doSymbol = false;
1165 _FString *extraPS = nil,
1166 *prefixPS = nil;
1167 long symbolsize = 3;
1168
1169
1170 _AssociativeList * toptions = (_AssociativeList*)FetchObjectFromVariableByType (&kTreeOutputOptions,ASSOCIATIVE_LIST);
1171
1172 if (toptions) {
1173 HBLObjectRef lc = toptions->GetByKey (kTreeOutputLayout, NUMBER);
1174 if (lc) {
1175 treeLayout = lc->Value();
1176 }
1177 lc = toptions->GetByKey (kTreeOutputEmbed, NUMBER);
1178 if (lc) {
1179 doEmbed = lc->Value();
1180 }
1181 lc = toptions->GetByKey(kTreeOutputSymbols, NUMBER);
1182 if ( lc ) {
1183 doSymbol = lc->Value();
1184 }
1185
1186 lc = toptions->GetByKey(kTreeOutputExtraPS, STRING);
1187 if ( lc ) {
1188 extraPS = (_FString*)lc->Compute();
1189 }
1190
1191 lc = toptions->GetByKey(kTreeOutputPrefixPS, STRING);
1192 if ( lc ) {
1193 prefixPS = (_FString*)lc->Compute();
1194 }
1195
1196 lc = toptions->GetByKey(kTreeOutputSymbolSize, NUMBER);
1197 if ( lc ) {
1198 symbolsize = lc->Value();
1199 }
1200 }
1201
1202 _String* theParam = (_String*) p->toStr(),
1203 t;
1204
1205 bool scaling = theParam->length(),
1206 doLabelWidth = true;
1207
1208 long tipCount = 0L,
1209 fontSize = -1L;
1210
1211 _Matrix* dimMatrix = ((_Matrix*)p2->Compute());
1212
1213
1214 hyFloat hScale = 1.0,
1215 vScale = 1.0,
1216 labelWidth = 0.,
1217
1218 treeHeight = (*dimMatrix)(0,1),
1219 treeWidth = (*dimMatrix)(0,0),
1220
1221 treeRotation = (dimMatrix->GetVDim()>2)?(*dimMatrix)(0,2):0.0,
1222 treeRadius ,
1223 treeArcAngle ,
1224 totalTreeL = 0.,
1225
1226 mappedTreeHeight = 0.0,
1227 mappedTreeHeight2 = 1e+100;
1228
1229 // letter size in points
1230 if (treeLayout == 1) {
1231 treeRadius = treeWidth;
1232 treeArcAngle = treeHeight;
1233
1234 if (treeRadius <= 0.0) {
1235 treeRadius = 300.;
1236 }
1237 if (treeArcAngle <= 0.0 || treeArcAngle > 360.) {
1238 treeArcAngle = 360.0;
1239 }
1240
1241 treeWidth = treeHeight = treeRadius * 2.;
1242
1243 treeRotation = (treeRotation - floor(treeRotation/360.)*360.)/DEGREES_PER_RADIAN;
1244 } else {
1245 if (treeHeight <= 0.0) {
1246 treeHeight = 792.;
1247 }
1248 if (treeWidth <= 0.0) {
1249 treeWidth = 576.;
1250 }
1251 }
1252
1253
1254
1255 (*res)<< (_String("%!\n%% PS file for the tree '")&*GetName()&"'.\n%% Generated by "& GetVersionString() & " on " & GetTimeStamp ());
1256 if (treeLayout == 1) {
1257 (*res) << "% Radial layout\n"
1258 << "/righttext {dup newpath 0 0 moveto false charpath closepath pathbbox pop exch pop exch sub 4 -1 roll exch sub 3 -1 roll newpath moveto show} def\n";
1259 }
1260 if (!doEmbed) {
1261 (*res)<<"<< /PageSize [" << _String(treeWidth+15) /*wayne added 15 to make trees sit inside the page */
1262 <<' '
1263 <<_String(treeHeight+15)
1264 <<"] >> setpagedevice\n";
1265 }
1266
1267
1268 long xtraChars = 0;
1269
1270 if (toptions) {
1271 _Matrix* rgbColor = (_Matrix*)(toptions)->GetByKey (kTreeOutputBackground, MATRIX);
1272 if (rgbColor) {
1273 t = _String((*rgbColor)(0,0)) & " " & _String((*rgbColor)(0,1)) & " " & _String((*rgbColor)(0,2)) & " setrgbcolor\nnewpath\n";
1274 (*res) << t
1275 << "0 0 moveto\n"
1276 << "0 "
1277 << _String(treeHeight)
1278 << " lineto\n"
1279 << _String(treeWidth)
1280 << ' '
1281 << _String(treeHeight)
1282 << " lineto\n"
1283 << _String(treeWidth)
1284 << " 0 lineto\n"
1285 << "closepath\nfill\n0 0 0 setrgbcolor\n";
1286 }
1287
1288 if ( doSymbol ) {
1289 /*add some symbol drawing postscript functions */
1290 (*res) << "/size "
1291 << _String(symbolsize)
1292 << " def\n"
1293 << "/box { 0 -0.5 size mul rmoveto\n"
1294 << "1 size mul 0 size mul rlineto\n"
1295 << "0 size mul 1 size mul rlineto\n"
1296 << "-1 size mul 0 size mul rlineto\n"
1297 << "closepath\n"
1298 << "} def\n"
1299 << "/triangle { size size 0.5 mul rlineto 0 size -1 mul rlineto closepath } def\n"
1300 << "/circle {currentpoint exch 0.5 size mul add exch 0.5 size mul 180 540 arc\n"
1301 << "closepath\n"
1302 << "} def\n"
1303 << "/diamond { 0 -0.5 size mul rmoveto 0.5 size mul 0 rmoveto 45 rotate 0.707107 size mul 0 rlineto 0 size 0.707107 mul rlineto -0.707107 size mul 0 rlineto -45 rotate closepath} def\n";
1304 }
1305
1306
1307 _Constant* fontSizeIn = (_Constant*)(toptions)->GetByKey (kTreeOutputFSPlaceH, NUMBER);
1308 if (fontSizeIn) {
1309 fontSize = fontSizeIn->Value();
1310 }
1311
1312 fontSizeIn = (_Constant*)(toptions)->GetByKey (kTreeOutputRightMargin, NUMBER);
1313 if (fontSizeIn) {
1314 treeWidth = MAX(treeWidth/4+1,treeWidth-fontSizeIn->Value());
1315 doLabelWidth = false;
1316 }
1317
1318 if ((fontSizeIn = (_Constant*)(toptions)->GetByKey (kTreeOutputXtraMargin, NUMBER))) {
1319 xtraChars = fontSizeIn->Value();
1320 }
1321 }
1322
1323 hyTopologyBranchLengthMode mapMode;
1324 _String scalingStringMatch;
1325
1326 if (scaling) {
1327 scalingStringMatch = DetermineBranchLengthMappingMode (theParam, mapMode);
1328 }
1329
1330 if (treeLayout == 1) {
1331
1332 long tipID = 0;
1333 EdgeCount (tipCount, tipID);
1334 tipID = 0;
1335 hScale = 0.;
1336 newRoot = RadialBranchMapping (theRoot,nil,&scalingStringMatch,treeArcAngle*pi_const/(180.*tipCount),tipID,hScale,mapMode);
1337 totalTreeL = hScale;
1338 } else {
1339 if (scaling) {
1340 newRoot = ScaledBranchMapping (nil, &scalingStringMatch, 0, tipCount,mapMode);
1341 } else {
1342 hyFloat offset;
1343 newRoot = AlignedTipsMapping (theRoot, offset, true);
1344 }
1345
1346 hScale = -treeWidth/newRoot->in_object.h;
1347 }
1348
1349 currentNd = NodeTraverser (newRoot);
1350
1351 while (currentNd) {
1352 if (currentNd->in_object.v > mappedTreeHeight) {
1353 mappedTreeHeight = currentNd->in_object.v;
1354 }
1355 if (currentNd->in_object.v < mappedTreeHeight2) {
1356 mappedTreeHeight2 = currentNd->in_object.v;
1357 }
1358
1359 if (!currentNd->get_num_nodes()) {
1360 tipCount++;
1361 }
1362
1363 currentNd = NodeTraverser ((node<nodeCoord>*)nil);
1364 }
1365
1366 // compute the size of the font, based on the spacing between branches.
1367 // 36 is max and 2 is min;
1368
1369 if (fontSize<0) {
1370 fontSize = treeHeight/(tipCount+1+(scaling>0)*1.5);
1371
1372 if (fontSize>36) {
1373 fontSize = 36;
1374 } else if (fontSize<2) {
1375 fontSize = 2;
1376 }
1377 }
1378
1379 // now recompute the width of the of the tree, including the labels
1380 // do not allow names to take up more that 1/2 of the width
1381 // try to keep full length but if the names are too wide,
1382 // shrink the size of the font
1383
1384 currentNd = NodeTraverser (newRoot);
1385
1386 hyFloat plotBounds[4];
1387
1388 if (treeLayout == 1) {
1389 plotBounds [0] = plotBounds[1] = plotBounds [2] = plotBounds[3] = 0.;
1390 }
1391
1392 while (currentNd) {
1393 if (currentNd->in_object.varRef>=0) {
1394 _String nodeName (*LocateVar (currentNd->in_object.varRef)->GetName());
1395 nodeName.Trim (nodeName.FindBackwards ('.',0,-1)+1,-1);
1396
1397 _AssociativeList * nodeOptions = nil;
1398 if (toptions) {
1399 nodeOptions = (_AssociativeList *) toptions->GetByKey (nodeName, ASSOCIATIVE_LIST);
1400 }
1401
1402
1403 HBLObjectRef nodeLabel = nodeOptions?nodeOptions->GetByKey (_TheTree::kTreeOutputLabel,STRING):nil,
1404 nodeTLabel = nodeOptions?nodeOptions->GetByKey (_TheTree::kTreeOutputTLabel,STRING):nil;
1405
1406 hyFloat nnWidth = 0.0;
1407
1408 if (nodeLabel) {
1409 nnWidth = 0.0;
1410 } else if (nodeTLabel) {
1411 nnWidth = 1.+PSStringWidth (((_FString*)nodeTLabel->Compute())->get_str());
1412 //printf ("%g\n", nnWidth);
1413 } else if (currentNd->get_num_nodes() == 0) {
1414 nnWidth = 1.+PSStringWidth (nodeName);
1415 }
1416
1417 nnWidth += _maxTimesCharWidth * xtraChars;
1418 nnWidth *= fontSize;
1419
1420 if (treeLayout == 1) {
1421 currentNd->in_object.label2 += treeRotation;
1422 hyFloat chordLength = computeChordLength (treeRadius, currentNd->in_object.label2,plotBounds),
1423 overflow = MAX(0., treeRadius +
1424 (nnWidth - chordLength) *
1425 hScale / MAX(HY_REPLACE_BAD_BRANCH_LENGTH_WITH_THIS,currentNd->in_object.label1));
1426
1427 //nnWidth + currentNd->in_object.label1 * vScale-chordLength);
1428
1429
1430
1431 if (overflow > treeRadius*.5) { // shift too big
1432 chordLength -= currentNd->in_object.label1/(2.*hScale) * treeRadius;
1433
1434 chordLength = chordLength / nnWidth;
1435 fontSize = fontSize * chordLength;
1436
1437 if (fontSize<2) {
1438 fontSize = 2;
1439 }
1440
1441 nnWidth = treeRadius*.5;
1442 } else {
1443 nnWidth = overflow;
1444 }
1445 } else {
1446 if (nnWidth>treeWidth*.5) {
1447 fontSize = fontSize / (2.*nnWidth/treeWidth);
1448 if (fontSize<2) {
1449 fontSize = 2;
1450 }
1451 nnWidth = treeWidth*.5;
1452 }
1453 }
1454 if (nnWidth > labelWidth) {
1455 labelWidth = nnWidth;
1456 }
1457 }
1458 currentNd = NodeTraverser ((node<nodeCoord>*)nil);
1459 }
1460
1461 if (!doLabelWidth) {
1462 labelWidth = 0;
1463 }
1464
1465 if (scaling) {
1466 treeHeight -= 3.0*fontSize;
1467 }
1468
1469 if (treeLayout == 1) {
1470 hScale = (treeRadius-labelWidth)/hScale;
1471 vScale = 0.;
1472 } else {
1473 hScale = -(treeWidth-labelWidth)/newRoot->in_object.h;
1474 vScale = (treeHeight-2*fontSize)/(mappedTreeHeight - mappedTreeHeight2);
1475 }
1476 t = _String ("/Times-Roman findfont\n") & fontSize & " scalefont\nsetfont\n";
1477 (*res)<<&t;
1478
1479 nodeCoord dummy;
1480
1481 long lw = fontSize/6+1;
1482
1483 (*res) << _String(lw) << " setlinewidth\n1 setlinecap\n";
1484
1485 if (prefixPS) {
1486 (*res) << prefixPS->get_str().Replace (kTreeOutputFSPlaceH, _String(fontSize), true);
1487 }
1488
1489 if (treeLayout == 1) {
1490 //newRoot->in_object.h = -plotBounds[1];
1491 //newRoot->in_object.v = -plotBounds[3];
1492 //hScale *= 2.*treeRadius / MAX(plotBounds[0]-plotBounds[1],plotBounds[2]-plotBounds[3]);
1493 newRoot->in_object.h = treeRadius;
1494 newRoot->in_object.v = treeRadius;
1495 vScale = 0.;
1496
1497 TreePSRecurse (newRoot, (*res), hScale, vScale, ceil(treeWidth), ceil(treeHeight), fontSize/2, labelWidth-fontSize/2, parse_settings, toptions, 1, &treeRadius);
1498 } else {
1499 TreePSRecurse (newRoot, (*res), hScale, vScale, ceil(treeWidth), ceil(treeHeight), fontSize/2, labelWidth-fontSize/2, parse_settings, toptions);
1500 }
1501
1502 if (scaling) { /* ruler */
1503 if (fontSize < 8) {
1504 fontSize = 8;
1505 (*res) << (_String ("/Times-Roman findfont\n") & fontSize & " scalefont\nsetfont\n");
1506 }
1507
1508 if (treeWidth < 4*fontSize) { // enforce minimal ruler width
1509 treeWidth = 4*fontSize;
1510 }
1511
1512 vScale = exp (log(10.)*floor (log10(treeLayout==1?totalTreeL:treeWidth/hScale))); // compute the scaling factor
1513 if (vScale*hScale < 10.0) {
1514 vScale *= 10.;
1515 }
1516
1517 _String rulerLabel (vScale, "%5.2g");
1518
1519 while (vScale*hScale > (treeLayout==1?treeRadius/3:treeWidth-newRoot->in_object.h)) {
1520 vScale *= 0.5;
1521 rulerLabel = vScale;
1522 }
1523
1524 while (PSStringWidth(rulerLabel)*fontSize>vScale*hScale-3) {
1525 vScale *= 2.0;
1526 rulerLabel = vScale;
1527 }
1528
1529 long lm = newRoot->in_object.h,
1530 rm = lw+vScale*hScale;
1531
1532 if (treeLayout == 1) {
1533 treeHeight = 2*treeRadius - 2*fontSize;
1534 lm = treeWidth - fontSize - rm;
1535 rm = treeWidth - fontSize - lw;
1536 }
1537
1538 (*res) << "newpath\n"
1539 << (_String (lm) & ' ' & _String ((long)treeHeight+2*lw) & " moveto\n")
1540 << (_String (rm) & ' ' & (long)(treeHeight+2*lw) & " lineto\nstroke\nnewpath\n") // main horizontal bar
1541 << (_String (lm) & ' ' & _String ((long)treeHeight) & " moveto\n")
1542 << (_String (lm) & ' ' & (long)(treeHeight+4*lw) & " lineto\nstroke\nnewpath\n")// left notch
1543 << ( _String (rm) & ' ' & (long)(treeHeight) & " moveto\n")
1544 << (_String (rm) & ' ' & (long)(treeHeight+4*lw) & " lineto\nstroke\nnewpath\n") // right notch
1545 << (_String (lm+lw) & ' ' & _String ((long)treeHeight+3*lw) & " moveto\n") // main text
1546 << (_String ('(') & _String (rulerLabel) & ") show\n");
1547 }
1548
1549 newRoot->delete_tree ();
1550 delete newRoot;
1551
1552 if (extraPS) {
1553 (*res) << extraPS->get_str().Replace (kTreeOutputFSPlaceH, _String(fontSize), true);
1554 }
1555
1556 if (!doEmbed) {
1557 (*res)<< "showpage";
1558 }
1559 DeleteObject (theParam);
1560 } else {
1561 ReportWarning ("An invalid 3rd parameter was passed to PSTreeString.");
1562 }
1563 } else {
1564 ReportWarning ("An invalid 2nd parameter was passed to PSTreeString.");
1565 }
1566 res->TrimSpace();
1567 return _returnStringOrUseCache(res, cache);
1568 }
1569
1570
1571
1572
1573 //_______________________________________________________________________________________________
1574
CountTreeCategories(void)1575 long _TheTree::CountTreeCategories (void) {
1576 // 20180618 TODO: SLKP this is an expensive call (rescans all category variables recursively);
1577 // could use an optimization pass
1578
1579 categoryVariables.Clear();
1580 _AVLList cVA (&categoryVariables);
1581 ScanForCVariables (cVA);
1582 cVA.ReorderList ();
1583 categoryCount = 1L;
1584
1585 categoryVariables.Each ([&] (long v, unsigned long) -> void {
1586 categoryCount *= ((_CategoryVariable*)LocateVar(v))->GetNumberOfIntervals();
1587 });
1588 return categoryCount;
1589 }
1590
1591
1592 //__________________________________________________________________________________
1593
TreePSRecurse(node<nodeCoord> * iterator,_StringBuffer & res,hyFloat hScale,hyFloat vScale,long hSize,long vSize,long halfFontSize,long shift,const _TreeTopologyParseSettings & settings,_AssociativeList * outOptions,char layout,hyFloat * xtra) const1594 void _TheTree::TreePSRecurse (node<nodeCoord>* iterator, _StringBuffer &res, hyFloat hScale, hyFloat vScale,
1595 long hSize, long vSize, long halfFontSize, long shift, const _TreeTopologyParseSettings& settings, _AssociativeList* outOptions,
1596 char layout, hyFloat * xtra) const {
1597
1598 // TODO: SLKP 20180803, this needs review and possible deprecation
1599
1600
1601 static _String const kTreeOutputThickness ( "TREE_OUTPUT_BRANCH_THICKNESS"),
1602 kTreeOutputLinecap ( "TREE_OUTPUT_BRANCH_LINECAP"),
1603 kTreeOutputSplit ( "TREE_OUTPUT_BRANCH_SPLIT"),
1604 kTreeOutputNotchesColor ( "TREE_OUTPUT_BRANCH_NOTCHES_COLOR"),
1605 kTreeOutputNotches ( "TREE_OUTPUT_BRANCH_NOTCHES"),
1606 kTreeOutputColor ( "TREE_OUTPUT_BRANCH_COLOR"),
1607 kTreeOutputDash ( "TREE_OUTPUT_BRANCH_DASH"),
1608 kTreeOutputOLabel ( "TREE_OUTPUT_OVER_BRANCH"),
1609 kTreeOutputNNPlaceH ( "__NODE_NAME__");
1610
1611
1612
1613
1614 // 20180618 TODO: SLKP this needs review and possibly deprecation
1615
1616 unsigned long descendants = iterator->get_num_nodes();
1617 bool is_leaf = descendants == 0UL;
1618 //lineW = halfFontSize/3+1;
1619
1620 hyFloat vc,
1621 hc,
1622 vcl,
1623 hcl,
1624 hc1,
1625 hc2;
1626
1627 _String t,
1628 varName,
1629 colorString ("0 0 0 setrgbcolor\n");
1630
1631 if (!iterator->is_root()) {
1632 if (iterator->in_object.varRef >=0) {
1633 varName = *map_node_to_calcnode(iterator)->GetName();
1634 }
1635 } else {
1636 hc = GetRoot().in_object;
1637 if (hc >= 0.) {
1638 varName = (*LocateVar(hc)->GetName());
1639 }
1640 if (layout == 1) {
1641 res << (_String (iterator->in_object.h) & ' ' & _String (iterator->in_object.v) & " translate\n");
1642 }
1643 }
1644
1645 varName.Trim(varName.FindBackwards ('.',0L,-1L)+1L,-1L);
1646
1647 _AssociativeList * nodeOptions = outOptions ? (_AssociativeList *) outOptions->GetByKey (varName, ASSOCIATIVE_LIST) : nil;
1648
1649 HBLObjectRef nodeLabel = nodeOptions?nodeOptions->GetByKey (_TheTree::kTreeOutputLabel,STRING):nil,
1650 nodeTLabel = nodeOptions?nodeOptions->GetByKey (_TheTree::kTreeOutputTLabel,STRING):nil;
1651
1652 if (layout == 1) {
1653 hcl = iterator->in_object.h*hScale;
1654 vcl = iterator->in_object.v*hScale;
1655 } else {
1656 vcl = vSize-iterator->in_object.v*vScale,
1657 hcl = hSize+iterator->in_object.h*hScale-shift;
1658 }
1659
1660 if (is_leaf || nodeLabel)
1661 // terminal node or default label
1662 {
1663 t = kEmptyString;
1664 hyFloat myAngle = layout==1?iterator->in_object.label2*DEGREES_PER_RADIAN:0.0;
1665 if (layout == 1) {
1666 res << (_String(myAngle) & " rotate\n");
1667 vc = 0;
1668 hcl = (vScale + iterator->in_object.bL)*hScale;
1669 hc = hcl + halfFontSize;
1670 } else {
1671 vc = vcl-halfFontSize;
1672 hc = hcl+halfFontSize;
1673 }
1674
1675 if (!nodeTLabel) {
1676 res<<"newpath\n";
1677
1678 if (nodeLabel) {
1679 res << (_String(hcl) & ' ' & _String (vc) & " moveto\n");
1680 t = ((_FString*)nodeLabel->Compute())->get_str().Replace (kTreeOutputNNPlaceH, varName, true).Replace (kTreeOutputFSPlaceH, _String(halfFontSize*2), true) & '\n';
1681 }
1682
1683 if (is_leaf && t.empty()) {
1684 // generate the default label
1685 if (layout == 1 && myAngle > 90. && myAngle < 270.) {
1686 hyFloat xt = hc-halfFontSize/2,
1687 yt = vc-2*halfFontSize/3;
1688 t = _String (xt) & _String (" 0 translate 180 rotate 0 ") & _String (yt) & _String ('(') & varName & ") righttext 180 rotate -" & xt & " 0 translate\n";
1689 } else {
1690 res<< (_String(hc-halfFontSize/2) & ' ' & _String (vc-2*halfFontSize/3) & " moveto\n");
1691 t = _String ('(') & varName & ") show\n";
1692 }
1693 }
1694 res<<&t;
1695 }
1696
1697 if (is_leaf) {
1698 iterator->in_object.h = hc-halfFontSize;
1699 }
1700
1701 if (layout == 1) {
1702 res << (_String(-iterator->in_object.label2*DEGREES_PER_RADIAN) & " rotate\n");
1703 }
1704 }
1705
1706 long minChildHC = 0x0fffffff,
1707 newV = 0;
1708
1709 if (!is_leaf) {
1710 vc = vSize - iterator->in_object.v*vScale;
1711 hc = hSize + iterator->in_object.h*hScale-shift;
1712
1713 nodeCoord childCoord;
1714 for (long k = 1; k<=descendants; k++) {
1715 node<nodeCoord>* child = iterator->go_down(k);
1716 TreePSRecurse (child, res, hScale, (layout==1)?vScale+iterator->in_object.bL:vScale, hSize, vSize,halfFontSize,shift,settings, outOptions,layout,xtra);
1717 if (k==1) {
1718 hc1 = layout==1?child->in_object.label2:child->in_object.v;
1719 }
1720 if (k==descendants) {
1721 hc2 = layout==1?child->in_object.label2:child->in_object.v;
1722 }
1723 }
1724
1725 char doVLines = 3;
1726
1727 for (long k = 1; k<=descendants; k++) {
1728 node<nodeCoord>* child = iterator->go_down(k);
1729
1730 if (child->in_object.varRef>=0) {
1731 t = map_node_to_calcnode(child)->ContextFreeName();
1732 } else {
1733 t = kEmptyString;
1734 }
1735
1736 newV += child->in_object.v;
1737
1738 _AssociativeList * childOptions = nil;
1739
1740 hyFloat splitBranch = -1.,
1741 lineWP = 0.0;
1742
1743 _String childColor,
1744 notchColor,
1745 childDash,
1746 linewidth1,
1747 linewidth2,
1748 linecap1,
1749 linecap2;
1750
1751 _String const * blabelString = nil;
1752
1753
1754 _Matrix * notches = nil,
1755 * multiColor = nil;
1756
1757 if (layout == 1) {
1758 res << (_String(child->in_object.label2*DEGREES_PER_RADIAN) & " rotate\n");
1759 }
1760
1761 if (outOptions) {
1762 childOptions = (_AssociativeList *) outOptions->GetByKey (t, ASSOCIATIVE_LIST);
1763 if (childOptions) {
1764 HBLObjectRef keyVal = childOptions->GetByKey (kTreeOutputThickness,NUMBER);
1765 if (keyVal) {
1766 lineWP = keyVal->Compute()->Value();
1767 linewidth1 = _String("currentlinewidth ") & lineWP & " setlinewidth\n";
1768 linewidth2 = "setlinewidth\n";
1769 }
1770 keyVal = childOptions->GetByKey (kTreeOutputLinecap,NUMBER);
1771 if (keyVal) {
1772 linecap1 = _String("currentlinecap ") & (long)keyVal->Compute()->Value() & " setlinecap\n";
1773 linecap2 = "setlinecap\n";
1774 }
1775 keyVal = childOptions->GetByKey (kTreeOutputSplit,NUMBER);
1776 if (keyVal) {
1777 splitBranch = keyVal->Compute()->Value();
1778 }
1779
1780 keyVal = childOptions->GetByKey (kTreeOutputNotches,MATRIX);
1781 if (keyVal) {
1782 notches = (_Matrix*)(((_Matrix*)keyVal->Compute())->ComputeNumeric())->makeDynamic();
1783 }
1784
1785 keyVal = childOptions->GetByKey (kTreeOutputColor,MATRIX);
1786 if (keyVal) {
1787 _Matrix* rgbColor = (_Matrix*)keyVal->Compute();
1788 if (rgbColor->GetHDim() > 1 && rgbColor->GetVDim () == 4 && layout != 1) {
1789 multiColor = (_Matrix*)rgbColor->makeDynamic();
1790 } else {
1791 childColor = _String((*rgbColor)(0,0)) & " " & _String((*rgbColor)(0,1)) & " " & _String((*rgbColor)(0,2)) & " setrgbcolor\n";
1792 }
1793 }
1794
1795 keyVal = childOptions->GetByKey (kTreeOutputNotchesColor,MATRIX);
1796 if (keyVal) {
1797 _Matrix* rgbColor = (_Matrix*)keyVal->Compute();
1798 notchColor = _String((*rgbColor)(0,0)) & " " & _String((*rgbColor)(0,1)) & " " & _String((*rgbColor)(0,2)) & " setrgbcolor\n";
1799 }
1800
1801 keyVal = childOptions->GetByKey (kTreeOutputDash,MATRIX);
1802 if (keyVal) {
1803 _Matrix* dash = (_Matrix*)keyVal->Compute();
1804 childDash = _String('[') & _String((*dash)(0,0)) & " " & _String((*dash)(0,1)) & "] " & _String((*dash)(0,2)) & " setdash\n";
1805 }
1806 keyVal = childOptions->GetByKey (kTreeOutputOLabel,STRING);
1807 if (keyVal) {
1808 blabelString = &((_FString*)keyVal->Compute())->get_str();
1809 }
1810 }
1811 }
1812
1813
1814 if (blabelString) {
1815 if (layout == 1) {
1816 t = _String(iterator->in_object.label1*hScale) & " 0 moveto\n";
1817 } else {
1818 t = _String(hc) & ' ' & _String (child->in_object.v) & " moveto\n";
1819 }
1820 res<<&t
1821 << (blabelString->Replace (kTreeOutputNNPlaceH, varName, true).Replace (_TheTree::kTreeOutputFSPlaceH, _String(halfFontSize*2), true)) << '\n'
1822 <<'\n';
1823 }
1824
1825 res << childDash << childColor << linewidth1 << linecap1;
1826
1827 if (layout == 1) {
1828 if (child->get_num_nodes()) { // internal node
1829 res << "newpath\n" << ((_String("0 0 ") & child->in_object.label1*hScale & ' ' & child->in_object.v & ' ' & (child->in_object.auxD) & " arc \n")) << "stroke\n";
1830 }
1831 } else {
1832 if (childOptions && descendants == 2) {
1833 res << "newpath\n"
1834 << ((_String(hc) & ' ' & _String (0.5*(hc1+hc2)) & " moveto\n"))
1835 << ((_String(hc) & ' ' & _String (k==1?hc1:hc2) & " lineto\n"))
1836 << "stroke\n";
1837
1838 doVLines -= k;
1839 }
1840
1841 }
1842
1843 res << "newpath\n";
1844 if (layout == 1) {
1845 res << (_String(child->in_object.label1*hScale) & " 0 moveto\n");
1846 t = _String(-child->in_object.bL*hScale) & " 0 rlineto\n";
1847 } else {
1848
1849 hyFloat lineWidthInset = 0.0;
1850
1851 //if (lineWP > 0.0)
1852 // lineWidthInset = (lineWP-lineW)*.5;
1853
1854 if (multiColor) {
1855 hyFloat span = child->in_object.h - hc + 2*lineWidthInset,
1856 currentX = hc - lineWidthInset;
1857
1858 res << (_String(currentX) & ' ' & _String (child->in_object.v) & " moveto\n");
1859 for (long seg = 0; seg < multiColor->GetHDim(); seg++) {
1860 res << (_String((*multiColor)(seg,0)) & " " & _String((*multiColor)(seg,1)) & " " & _String((*multiColor)(seg,2)) & " setrgbcolor\n");
1861 hyFloat mySpan = span*(*multiColor)(seg,3);
1862 res << (_String(mySpan) & " 0 rlineto\n");
1863 if (seg < multiColor->GetHDim()-1) {
1864 currentX += mySpan;
1865 res << "stroke\nnewpath\n"
1866 << ((_String(currentX) & ' ' & _String (child->in_object.v) & " moveto\n"));
1867 }
1868 }
1869 DeleteObject (multiColor);
1870 t = kEmptyString;
1871 } else {
1872 res<< (_String(hc - lineWidthInset) & ' ' & _String (child->in_object.v) & " moveto\n");
1873 t = _String(child->in_object.h) & ' ' & _String (child->in_object.v + lineWidthInset) & " lineto\n";
1874 }
1875 if (layout == 1) {
1876 minChildHC = MIN(child->in_object.label1,minChildHC);
1877 } else {
1878 minChildHC = MIN(child->in_object.h,minChildHC);
1879 }
1880 }
1881 res <<t
1882 << "stroke\n"
1883 << linecap2
1884 << linewidth2;
1885
1886 if (childDash.nonempty()) {
1887 res << "[] 0 setdash\n";
1888 }
1889
1890 if (splitBranch >= 0.0 && splitBranch <= 1.0) {
1891 res << "newpath\n";
1892 hyFloat x,
1893 y;
1894
1895 if (layout == 1) {
1896 x = (child->in_object.label1 - child->in_object.bL*splitBranch)*hScale;
1897 y = 0.;
1898 } else {
1899 x = hc+(child->in_object.h-hc)*(1.-splitBranch);
1900 y = child->in_object.v;
1901 }
1902
1903 res<< (_String(x) & ' ' & _String (y) & " " & halfFontSize & " 0 360 arc\n") << "fill\n";
1904 }
1905
1906 if (notches) {
1907 notches->CheckIfSparseEnough(true);
1908 res << notchColor;
1909 for (long l = 0; l < notches->GetSize(); l++) {
1910 hyFloat aNotch = (*notches)[l];
1911 if (aNotch >= 0. && aNotch <= 1.) {
1912 res << "newpath\n";
1913 hyFloat x,
1914 y;
1915
1916 if (layout == 1) {
1917 x = (child->in_object.label1 - child->in_object.bL*aNotch)*hScale;
1918 y = 0.;
1919 } else {
1920 x = hc+(child->in_object.h-hc)*(1.-aNotch);
1921 y = child->in_object.v;
1922 }
1923
1924 res << (_String(x-0.5*halfFontSize) & ' ' & _String (y-0.5*halfFontSize) & " moveto ")
1925 << (_String(x+0.5*halfFontSize) & ' ' & _String (y+0.5*halfFontSize) & " lineto\n")
1926 << (_String(x-0.5*halfFontSize) & ' ' & _String (y+0.5*halfFontSize) & " moveto ")
1927 << (_String(x+0.5*halfFontSize) & ' ' & _String (y-0.5*halfFontSize) & " lineto\n")
1928 << "stroke\n";
1929 }
1930 }
1931 DeleteObject (notches);
1932 }
1933
1934 res << &colorString;
1935 if (layout == 1) {
1936 res << (_String(-child->in_object.label2*DEGREES_PER_RADIAN) & " rotate\n");
1937 }
1938 }
1939
1940 if (layout == 0 && doVLines) {
1941 _String linewidth1,
1942 linewidth2,
1943 linecap1,
1944 linecap2;
1945
1946 if (nodeOptions) {
1947 HBLObjectRef keyVal = nodeOptions->GetByKey (kTreeOutputThickness,NUMBER);
1948 if (keyVal) {
1949 hyFloat lineWP = keyVal->Compute()->Value();
1950 linewidth1 = _String("currentlinewidth ") & lineWP & " setlinewidth\n";
1951 linewidth2 = "setlinewidth\n";
1952 }
1953 keyVal = nodeOptions->GetByKey (kTreeOutputLinecap,NUMBER);
1954 if (keyVal) {
1955 linecap1 = _String("currentlinecap ") & (long)keyVal->Compute()->Value() & " setlinecap\n";
1956 linecap2 = "setlinecap\n";
1957 }
1958 }
1959
1960 res << linecap1
1961 << linewidth1
1962 << "newpath\n";
1963
1964 if (doVLines == 3) {
1965 t = _String(hc) & ' ' & _String (hc2) & " moveto\n";
1966 } else {
1967 t = _String(hc) & ' ' & _String (0.5*(hc1+hc2)) & " moveto\n";
1968 }
1969 res<< t;
1970 if (doVLines == 3) {
1971 t = _String(hc) & ' ' & _String (hc1) & " lineto\n";
1972 } else {
1973 t = _String(hc) & ' ' & _String (doVLines==1?hc1:hc2) & " lineto\n";
1974 }
1975 res<< t
1976 << "stroke\n"
1977 << &linewidth2
1978 << &linecap2;
1979 }
1980
1981 if (layout == 0) {
1982 iterator->in_object.h = hc;
1983 iterator->in_object.v = newV/descendants;
1984 } else {
1985 iterator->in_object.auxD = (hc2-iterator->in_object.label2)*DEGREES_PER_RADIAN;
1986 if (!(iterator->is_root())) {
1987 iterator->in_object.v = (hc1-iterator->in_object.label2)*DEGREES_PER_RADIAN;
1988 }
1989 }
1990 } else {
1991 iterator->in_object.v = vc;
1992 }
1993
1994 if (nodeTLabel) {
1995 t = ((_FString*)nodeTLabel->Compute())->get_str();
1996 if (t.nonempty()) {
1997 hyFloat scF = 2.*halfFontSize;
1998
1999 if (layout == 1) {
2000 res << (_String(iterator->in_object.label2*DEGREES_PER_RADIAN) & " rotate\n");
2001 } else {
2002 if (!is_leaf) {
2003 vcl = iterator->in_object.v;
2004 }
2005 }
2006
2007
2008 if (scF != 2.*halfFontSize) {
2009 res << (_String("/Times-Roman findfont\n")& scF & " scalefont\nsetfont\n");
2010 }
2011
2012 res << "newpath\n";
2013 if (layout == 1) {
2014 res << (_String(iterator->in_object.label1 * hScale + halfFontSize) & ' ' & _String (-2*halfFontSize/3) & " moveto\n");
2015 } else {
2016 if (!is_leaf) {
2017 hc = hcl + scF*0.5;
2018 vc = vcl - scF*0.33;
2019
2020 } else {
2021 hc -= scF*0.25;
2022 vc -= scF*0.33;
2023 }
2024 res << (_String(hc) & ' ' & _String (vc) & " moveto\n");
2025 }
2026
2027 res << '(' << t << ") show \n";
2028
2029 if (scF != 2.*halfFontSize) {
2030 res << (_String("/Times-Roman findfont\n")& halfFontSize*2 & " scalefont\nsetfont\n");
2031 }
2032
2033 if (layout == 1) {
2034 res << (_String(-iterator->in_object.label2*DEGREES_PER_RADIAN) & " rotate\n");
2035 }
2036 }
2037 }
2038
2039 if (colorString.nonempty()) {
2040 res << "0 0 0 setrgbcolor\n";
2041 }
2042
2043 if (iterator->is_root() == false && layout == 1) {
2044 res << (_String (-iterator->in_object.h) & ' ' & _String (-iterator->in_object.v) & " translate\n");
2045 }
2046 }
2047
2048 //__________________________________________________________________________________
2049
SetUpMatrices(long categCount)2050 void _TheTree::SetUpMatrices (long categCount) {
2051 //fprintf (stderr, "[_TheTree::SetUpMatrices] %ld\n", categCount);
2052
2053 categoryCount = Maximum (categCount,1L);
2054
2055 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
2056
2057 while (_CalcNode* iterator = ti.Next()) {
2058 if (iterator->IsConstant()) {
2059 iterator->varFlags |= HY_VC_NO_CHECK;
2060 }
2061 iterator->ConvertToSimpleMatrix();
2062
2063 if (categoryCount==1L) {
2064 iterator->matrixCache = nil;
2065 } else {
2066 iterator->matrixCache = new _Matrix* [categoryCount];
2067 InitializeArray (iterator->matrixCache, categCount, (_Matrix*)nil);
2068 }
2069 }
2070 }
2071
2072
2073 //__________________________________________________________________________________
2074
CleanUpMatrices(void)2075 void _TheTree::CleanUpMatrices (void) {
2076 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
2077
2078 if (categoryCount == 1L) {
2079 while (_CalcNode* iterator = ti.Next()) {
2080
2081 // mod 05/03/2003 - uncomment next 5 lines
2082 // this breaks after ReplicateConstraint or MolecularClock is called
2083 // WTF?
2084
2085 iterator->ConvertFromSimpleMatrix();
2086
2087 if (iterator->compExp) {
2088 DeleteObject (iterator->compExp);
2089 iterator->compExp = nil;
2090 }
2091
2092 iterator->varFlags &= HY_VC_CLR_NO_CHECK;
2093 }
2094 } else {
2095 while (_CalcNode* iterator = ti.Next()) {
2096 iterator->ConvertFromSimpleMatrix();
2097
2098 for (long i=0; i<categoryCount; i++) {
2099 DeleteAndZeroObject(iterator->matrixCache[i]);
2100 }
2101
2102 delete [] iterator->matrixCache;
2103 iterator->matrixCache = nil;
2104 iterator->compExp = nil;
2105 iterator->varFlags &= HY_VC_CLR_NO_CHECK;
2106 }
2107 categoryCount = 1;
2108 }
2109 }
2110
2111 //__________________________________________________________________________________
2112
RemoveModel(void)2113 void _TheTree::RemoveModel (void) {
2114 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
2115 while (_CalcNode* iterator = ti.Next()) {
2116 iterator->RemoveModel();
2117 }
2118 categoryCount = 1;
2119 }
2120
2121
2122
2123 //__________________________________________________________________________________
2124
ScanContainerForVariables(_AVLList & l,_AVLList & l2,_AVLListX * tagger,long weight,_AVLListX * map_variables_to_nodes) const2125 void _TheTree::ScanContainerForVariables (_AVLList& l,_AVLList& l2, _AVLListX * tagger, long weight, _AVLListX * map_variables_to_nodes) const {
2126 /**
2127 map_variables_to_nodes, if supplied, will be used to map variable IDs to "post-order" node index of nodes that they affect, IF they only affect one node; otherwise they will be tagged with '-1'
2128 */
2129
2130 unsigned long traversal_order = 0UL,
2131 leaf_index = 0UL,
2132 int_index = 0UL;
2133
2134 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot);
2135 while (_CalcNode* iterator = ti.Next()) {
2136 bool is_leaf = ti.IsAtLeaf();
2137
2138
2139 iterator->ScanContainerForVariables(l,l2, tagger, weight + flatNodes.lLength + flatLeaves.lLength - traversal_order, map_variables_to_nodes, is_leaf ? leaf_index : int_index + flatLeaves.lLength);
2140
2141 if (is_leaf) {
2142 leaf_index ++;
2143 } else {
2144 int_index ++;
2145 }
2146 traversal_order ++;
2147 }
2148 }
2149
2150 //__________________________________________________________________________________
2151
ScanAndAttachVariables(void) const2152 void _TheTree::ScanAndAttachVariables (void) const {
2153 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot);
2154 while (_CalcNode* iterator = ti.Next()) {
2155 iterator->ScanAndAttachVariables();
2156 }
2157 }
2158
2159 //__________________________________________________________________________________
2160
ScanForDVariables(_AVLList & l,_AVLList & l2) const2161 void _TheTree::ScanForDVariables (_AVLList& l,_AVLList& l2) const {
2162 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot);
2163 while (_CalcNode* iterator = ti.Next()) {
2164 iterator->ScanForDVariables(l,l2);
2165 }
2166 }
2167
2168 //__________________________________________________________________________________
2169
ScanForGVariables(_AVLList & li,_AVLList & ld,_AVLListX * tagger,long weight) const2170 void _TheTree::ScanForGVariables (_AVLList& li, _AVLList& ld, _AVLListX * tagger, long weight) const {
2171 _SimpleList cL;
2172 _AVLList cLL (&cL);
2173 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot );
2174
2175 while (_CalcNode* iterator = ti.Next()) {
2176
2177 _Formula *explicitFormMExp = iterator->GetExplicitFormModel ();
2178 _Matrix *modelM = explicitFormMExp?nil:iterator->GetModelMatrix();
2179
2180 if ((explicitFormMExp && cLL.Find ((BaseRef)explicitFormMExp) < 0) || (modelM && cLL.Find(modelM) < 0)) {
2181 _SimpleList temp;
2182 {
2183 _AVLList tempA (&temp);
2184 if (modelM) {
2185 modelM->ScanForVariables(tempA, true);
2186 } else {
2187 explicitFormMExp->ScanFForVariables(tempA, true, false, true, true);
2188 }
2189 tempA.ReorderList();
2190 }
2191 for (unsigned long i=0; i<temp.lLength; i++) {
2192 long p = temp.list_data[i];
2193 _Variable* v = LocateVar (p);
2194 if (v&&v->IsGlobal()) {
2195 if(v->IsIndependent()) {
2196 li.Insert ((BaseRef)p);
2197 if (tagger) {
2198 tagger->UpdateValue((BaseRef)p, weight, 0);
2199 }
2200 } else {
2201 ld.Insert ((BaseRef)p);
2202 }
2203 }
2204 }
2205 cLL.Insert (modelM?(BaseRef)modelM:(BaseRef)explicitFormMExp);
2206 }
2207 iterator -> ScanForGVariables(li,ld);
2208 }
2209 }
2210
2211 //__________________________________________________________________________________
2212
ScanForCVariables(_AVLList & lcat) const2213 void _TheTree::ScanForCVariables (_AVLList& lcat) const {
2214 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot);
2215 while (_CalcNode* iterator = ti.Next()) {
2216 for (unsigned long i = 0UL; i < iterator->categoryVariables.lLength; i++) {
2217 lcat.Insert ((BaseRef)iterator->categoryVariables.Element(i));
2218 }
2219 }
2220 }
2221
2222 //__________________________________________________________________________________
2223
HasChanged(bool)2224 bool _TheTree::HasChanged (bool) {
2225 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot);
2226 while (_CalcNode* iterator = ti.Next()) {
2227 if (iterator->HasChanged()) {
2228 return true;
2229 }
2230 }
2231 return false;
2232 }
2233
2234 //__________________________________________________________________________________
2235
HasChanged2(void)2236 bool _TheTree::HasChanged2 (void) {
2237
2238 for (unsigned long k = 0; k < categoryVariables.lLength; k++) {
2239 if (((_CategoryVariable*)LocateVar(categoryVariables.Element(k)))->HaveParametersChanged()) {
2240 return true;
2241 }
2242 }
2243
2244 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER | fTreeIteratorTraversalSkipRoot);
2245 while (_CalcNode* iterator = ti.Next()) {
2246 if (iterator->_VariableContainer::HasChanged()) {
2247 return true;
2248 }
2249 }
2250 return false;
2251 }
2252
2253
2254 //_______________________________________________________________________________________________
InitializeTreeFrequencies(_Matrix * mx,bool setDim)2255 void _TheTree::InitializeTreeFrequencies (_Matrix *mx, bool setDim) {
2256 // this will take the matrix of frequencies and
2257 // 1) use its dimensions to initialize tree freq holders
2258 // 2) place global frequencies into the tree holder for later use by the pruning algo
2259 // must be called before any tree pruning computations are started
2260
2261 unsigned long vecDim = mx->GetHDim()*mx->GetVDim();
2262 // theModel = mx;
2263
2264 if (setDim) {
2265 SetTreeCodeBase (vecDim);
2266 } else {
2267 CopyArray (theProbs, mx->theData, vecDim);
2268 }
2269 }
2270
2271 //_______________________________________________________________________________________________
SetTreeCodeBase(long b)2272 void _TheTree::SetTreeCodeBase (long b) {
2273 // this will take the matrix of frequencies and
2274 // 1) use its dimensions to initialize tree freq holders
2275 // 2) place global frequencies into the tree holder for later use by the pruning algo
2276 // must be called before any tree pruning computations are started
2277 SetCodeBase (b);
2278 if (cBase>0) {
2279 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
2280 while (_CalcNode* iterator = ti.Next()) {
2281 iterator -> SetCodeBase (b);
2282 }
2283 }
2284 }
2285
2286 //_______________________________________________________________________________________________
2287
IsLinkedToALF(long & pid) const2288 long _TheTree::IsLinkedToALF (long& pid) const {
2289 return likeFuncList.FindOnCondition([&] (BaseRefConst value, unsigned long index) -> bool {
2290 if (likeFuncNamesList.GetItem(index)) {
2291 pid = ((_LikelihoodFunction const*)value)->DependOnTree (*GetName());
2292 if (pid >= 0) {
2293 return true;
2294 }
2295 }
2296 return false;
2297 });
2298 }
2299
2300
2301
2302 //_______________________________________________________________________________________________
2303
IntPopulateLeaves(_DataSetFilter const * dsf,long site_index) const2304 bool _TheTree::IntPopulateLeaves (_DataSetFilter const* dsf, long site_index) const {
2305 // assign proper values to leaf conditional probability vectors
2306 bool site_has_all_gaps = true;
2307
2308 _String* buffer = dsf->MakeSiteBuffer();
2309
2310 for (long leaf_index = 0; leaf_index<flatLeaves.lLength; leaf_index++) {
2311 _CalcNode * iterator = (_CalcNode*)flatCLeaves.GetItem(leaf_index);
2312 dsf->RetrieveState(site_index, leaf_index, *buffer, false);
2313
2314 site_has_all_gaps &= (dsf->Translate2Frequencies (*buffer, iterator->theProbs, true)<0); // ambig
2315 site_has_all_gaps &= (!ArrayAny (iterator->theProbs, cBase, [](hyFloat x, unsigned long) {return x == 0.0; })); //completely unresolved
2316 map_node_to_calcnode (((node <long>*)flatLeaves.GetElement(leaf_index))->parent)->cBase = -1;
2317 }
2318
2319 DeleteObject (buffer);
2320 return site_has_all_gaps;
2321 }
2322
2323
2324 //_______________________________________________________________________________________________
2325
RecoverNodeSupportStates(_DataSetFilter const * dsf,long site_index,_Matrix & resultMatrix)2326 void _TheTree::RecoverNodeSupportStates (_DataSetFilter const* dsf, long site_index, _Matrix& resultMatrix) {
2327 // TODO SLKP 20180803 this needs to be moved to using standard log likelihood calculations
2328
2329 // assume current values of all parameters
2330 // return 2 sets of vectors for each branch
2331 // - top-down conditionals
2332 // - bottom-up conditionals
2333 // resultMatrix is assumed to contain
2334 // uniqueSites X (flatLeaves.lLength+flatTree.lLength)*cBase*2 X categoryCount
2335
2336 long globalShifter = (flatLeaves.countitems()+flatTree.countitems())*cBase,
2337 catShifer = dsf->GetPatternCount() * 2 * globalShifter;
2338
2339 IntPopulateLeaves (dsf, site_index);
2340
2341 /* pass 1; populate top-down vectors */
2342 /* ugly top-bottom algorithm for debuggability and compactness */
2343
2344 // SLKP 20180803 -- this is a very ugly hack to avoid reitroducing CalcNode->nodeIndex
2345 // this should be taken care of when RecoverNodeSupportStates is re-engineered
2346 _SimpleList _avl_storage;
2347 _AVLListX node_to_index (&_avl_storage);
2348 long leaf_count = 0L, int_node_count = flatLeaves.countitems();
2349
2350 node_iterator<long> tree_iterator (theRoot, _HY_TREE_TRAVERSAL_POSTORDER);
2351 while (node<long>*topTraverser = tree_iterator.Next()) {
2352 _CalcNode *sourceNode = map_node_to_calcnode (topTraverser);
2353 if (topTraverser->is_leaf()) {
2354 node_to_index.Insert (sourceNode, leaf_count++, false, false);
2355 } else {
2356 node_to_index.Insert (sourceNode, int_node_count++, false, false);
2357 }
2358 }
2359
2360
2361 for (long catCount = 0L; catCount < categoryCount; catCount ++) {
2362 hyFloat * currentStateVector = resultMatrix.theData + 2*globalShifter*site_index + catShifer*catCount,
2363 * vecPointer = currentStateVector;
2364
2365 for (long nodeCount = 0L; nodeCount<flatCLeaves.lLength; nodeCount++) {
2366 hyFloat *leafVec = ((_CalcNode*)(((BaseRef*)flatCLeaves.list_data)[nodeCount]))->theProbs;
2367 CopyArray(vecPointer, leafVec, cBase);
2368 vecPointer += cBase;
2369 }
2370
2371 // TODO SLKP 20180703: ugly fix for underflow which WON'T work if category count > 1
2372
2373 for (long iNodeCount = 0L; iNodeCount < flatTree.lLength - 1; iNodeCount++) {
2374 node<long>* thisINode = (node<long>*)flatNodes.list_data[iNodeCount];
2375
2376 hyFloat sum = 0.;
2377
2378 for (long cc = 0L; cc < cBase; cc++) {
2379 hyFloat tmp = 1.0;
2380
2381 for (long nc = 0; nc < thisINode->get_num_nodes(); nc++) {
2382 hyFloat tmp2 = 0.0;
2383 _CalcNode * child = map_node_to_calcnode(thisINode->go_down(nc+1));
2384
2385 hyFloat * childSupport = currentStateVector + node_to_index.GetDataByKey(child) *cBase,
2386 * transMatrix = child->GetCompExp(categoryCount>1?catCount:(-1))->theData + cc*cBase;
2387
2388 for (long cc2 = 0; cc2 < cBase; cc2++) {
2389 tmp2 += transMatrix[cc2] * childSupport[cc2];
2390 }
2391
2392 tmp *= tmp2;
2393 }
2394 vecPointer[cc] = tmp;
2395 sum += tmp;
2396 }
2397
2398 if (sum < _lfScalingFactorThreshold && sum > 0.0) {
2399 for (long cc = 0; cc < cBase; cc++) {
2400 vecPointer[cc] *= _lfScalerUpwards;
2401 }
2402 }
2403
2404 vecPointer += cBase;
2405
2406 }
2407 RecoverNodeSupportStates2 (&GetRoot(),currentStateVector+globalShifter,currentStateVector,categoryCount>1?catCount:(-1), node_to_index);
2408 }
2409 /* pass 2; populate bottom-up vectors */
2410 /* for this we need to traverse the tree pre-order */
2411 /* because speed is not much of a concern, use a recursive call for compactness */
2412
2413 }
2414
2415 //_______________________________________________________________________________________________
2416
RecoverNodeSupportStates2(node<long> * thisNode,hyFloat * resultVector,hyFloat * forwardVector,long catID,_AVLListX & lookup)2417 void _TheTree::RecoverNodeSupportStates2 (node<long>* thisNode, hyFloat * resultVector, hyFloat* forwardVector, long catID, _AVLListX& lookup) {
2418
2419 _CalcNode * thisNodeC = map_node_to_calcnode (thisNode);
2420 hyFloat * vecPointer = resultVector + lookup.GetDataByKey(thisNodeC) * cBase;
2421
2422 if (thisNode->parent) {
2423 if (thisNode->parent->parent) {
2424 hyFloat sum = 0.;
2425 for (long cc = 0; cc < cBase; cc++) {
2426 hyFloat tmp = 1.0;
2427 for (long nc = 0; nc < thisNode->parent->get_num_nodes(); nc++) {
2428 hyFloat tmp2 = 0.0;
2429 _CalcNode * child = map_node_to_calcnode(thisNode->parent->go_down (nc+1));
2430 bool invert = (child == thisNodeC);;
2431 if (invert) {
2432 child = map_node_to_calcnode (thisNode->parent);
2433 }
2434
2435 hyFloat * childSupport = invert?resultVector + cBase*lookup.GetDataByKey(child)
2436 :forwardVector + lookup.GetDataByKey(child)*cBase,
2437 * transMatrix = child->GetCompExp(catID)->theData + cc*cBase;
2438
2439 for (long cc2 = 0; cc2 < cBase; cc2++) {
2440 tmp2 += transMatrix[cc2] * childSupport[cc2];
2441 }
2442
2443 tmp *= tmp2;
2444 }
2445 vecPointer[cc] = tmp;
2446 sum += tmp;
2447 }
2448 if (sum < _lfScalingFactorThreshold && sum > 0.0) {
2449 for (long cc = 0; cc < cBase; cc++) {
2450 vecPointer[cc] *= _lfScalerUpwards;
2451 }
2452 }
2453
2454 } else {
2455 for (long cc = 0; cc < cBase; cc++,vecPointer++) {
2456 hyFloat tmp = 1.0;
2457 for (long nc = 0; nc < thisNode->parent->get_num_nodes(); nc++) {
2458 hyFloat tmp2 = 0.0;
2459 _CalcNode * child = ((_CalcNode*)((BaseRef*)variablePtrs.list_data)[thisNode->parent->get_node(nc+1)->in_object]);
2460 if (child != thisNodeC) {
2461 hyFloat * childSupport = forwardVector + lookup.GetDataByKey(child)*cBase,
2462 * transMatrix = child->GetCompExp(catID)->theData + cc*cBase;
2463
2464 for (long cc2 = 0; cc2 < cBase; cc2++) {
2465 tmp2 += transMatrix[cc2] * childSupport[cc2];
2466 }
2467
2468 tmp *= tmp2;
2469 }
2470 }
2471 *vecPointer = tmp;
2472 }
2473 }
2474 } else {
2475 InitializeArray (vecPointer, cBase, 1.0);
2476 }
2477
2478 for (long nc = 0; nc < thisNode->get_num_nodes(); nc++) {
2479 RecoverNodeSupportStates2 (thisNode->get_node(nc+1),resultVector,forwardVector,catID,lookup);
2480 }
2481 }
2482
2483
2484 //_______________________________________________________________________________________________
ConstructNodeToIndexMap(bool doINodes) const2485 _AVLListX* _TheTree::ConstructNodeToIndexMap (bool doINodes) const {
2486 _SimpleList * nodes = new _SimpleList;
2487 const _SimpleList * whichL = doINodes?&flatNodes:&flatLeaves;
2488 _AVLListX * result = new _AVLListX (nodes);
2489
2490 for (unsigned long pistolero = 0UL; pistolero < whichL->lLength; pistolero++) {
2491 result->Insert ((BaseRef)whichL->list_data[pistolero], pistolero, false);
2492 }
2493
2494 return result;
2495
2496 }
2497
2498 //_______________________________________________________________________________________________
MapPostOrderToInOrderTraversal(_SimpleList & storeHere,bool doINodes) const2499 void _TheTree::MapPostOrderToInOrderTraversal (_SimpleList& storeHere, bool doINodes) const {
2500 _AVLListX* nodeMapper = ConstructNodeToIndexMap (doINodes);
2501
2502 _TreeIterator ti (this, doINodes ? _HY_TREE_TRAVERSAL_PREORDER : _HY_TREE_TRAVERSAL_POSTORDER);
2503
2504 unsigned long allNodeCount = 0UL;
2505
2506 storeHere.Populate (doINodes?flatTree.lLength:flatLeaves.lLength, 0, 0);
2507
2508 while (_CalcNode* iterator = ti.Next()) {
2509 bool isTip = ti.IsAtLeaf();
2510 if ( isTip && !doINodes || !isTip && doINodes) {
2511 storeHere.list_data[nodeMapper->GetXtra (nodeMapper->Find((BaseRef)(ti.GetNode())))] = allNodeCount++;
2512 }
2513 }
2514
2515 nodeMapper->DeleteAll(false);
2516 DeleteObject (nodeMapper);
2517 }
2518
2519 //_______________________________________________________________________________________________
MarkDone(void)2520 void _TheTree::MarkDone (void) {
2521 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
2522 while (_CalcNode* iterator = ti.Next()) {
2523 iterator -> MarkDone();
2524 }
2525 }
2526
2527 //_______________________________________________________________________________________________
2528
ComputeReleafingCostChar(_DataSetFilter const * dsf,long firstIndex,long secondIndex,const _SimpleList * child_count) const2529 long _TheTree::ComputeReleafingCostChar (_DataSetFilter const* dsf, long firstIndex, long secondIndex, const _SimpleList* child_count) const {
2530
2531 const char *pastState = dsf->GetColumn(firstIndex),
2532 *thisState = dsf->GetColumn(secondIndex);
2533
2534 long rootIndex = flatTree.lLength - 1,
2535 theCost = child_count->list_data[rootIndex],
2536 offset = flatLeaves.countitems();
2537
2538 bool * marked_nodes = (bool*) calloc ( flatTree.lLength, sizeof(bool));
2539 //InitializeArray(marked_nodes, flatTree.lLength, false);
2540 //memset (marked_nodes, 0, sizeof(bool) * flatTree.lLength);
2541
2542 for (long node_index = 0L; node_index < flatLeaves.lLength; node_index++) {
2543 long seq_index = dsf->theNodeMap.list_data[node_index];
2544 if (thisState[seq_index] != pastState[seq_index]) {
2545 /*long parentIndex = flatParents.list_data[node_index];
2546 while (marked_nodes[parentIndex] == false && parentIndex < rootIndex) {
2547 marked_nodes[parentIndex] = true;
2548 theCost += child_count->list_data [parentIndex];
2549 parentIndex = flatParents.list_data[parentIndex+offset];
2550 }*/
2551 marked_nodes [flatParents.list_data[node_index]] = true;
2552 }
2553 }
2554
2555 // the root will always be changed, so don't need to check it
2556 for (long i=0; i<flatTree.lLength - 1; i++) {
2557 if (marked_nodes[i]) {
2558 marked_nodes [flatParents.list_data [i + offset]] = true;
2559 theCost += child_count->list_data [i];
2560 }
2561 }
2562 free (marked_nodes);
2563
2564 return theCost;
2565
2566 }
2567
2568 //_______________________________________________________________________________________________
2569
ClearConstraints(void)2570 void _TheTree::ClearConstraints (void) {
2571 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
2572 while (_CalcNode* iterator = ti.Next()) {
2573 iterator->ClearConstraints();
2574 }
2575 }
2576
2577
2578 //_______________________________________________________________________________________________
2579
ComputeReleafingCost(_DataSetFilter const * dsf,long firstIndex,long secondIndex,_SimpleList * traversalTags,long orderIndex,_SimpleList const * childCount) const2580 long _TheTree::ComputeReleafingCost (_DataSetFilter const* dsf, long firstIndex, long secondIndex, _SimpleList* traversalTags, long orderIndex, _SimpleList const* childCount) const {
2581
2582 long filterL = dsf->GetPatternCount();
2583
2584
2585
2586 bool * marked_nodes = (bool*) calloc (flatTree.lLength, sizeof(bool));
2587 //memset (marked_nodes, 0, sizeof(bool) * flatTree.lLength);
2588 //InitializeArray(marked_nodes, flatTree.lLength, false);
2589
2590
2591 dsf->CompareTwoSitesCallback (firstIndex, secondIndex, [marked_nodes, this] (long idx, unsigned long di)->void{
2592 marked_nodes [this->flatParents.list_data[di]] = true;
2593 });
2594
2595 long theCost = 0,
2596 rootIndex = flatTree.lLength - 1;
2597
2598 // don't check the root; it is always "dirty"
2599
2600 if (childCount) {
2601 if (traversalTags && orderIndex) {
2602 for (long i=0; i<rootIndex; i++) {
2603 if (marked_nodes[i]) {
2604 marked_nodes[flatParents.list_data[flatLeaves.lLength + i]] = true;
2605 theCost += childCount->list_data[i];
2606 } else {
2607 long theIndex = filterL * i + orderIndex;
2608 traversalTags->list_data[theIndex/_HY_BITMASK_WIDTH_] |= bitMaskArray.masks[theIndex%_HY_BITMASK_WIDTH_];
2609 }
2610 }
2611 }
2612
2613 for (long i=0; i<rootIndex; i++) {
2614 if (marked_nodes[i]) {
2615 marked_nodes[flatParents.list_data[flatLeaves.lLength + i]] = true;
2616 theCost += childCount->list_data[i];
2617 }
2618 }
2619 } else {
2620 for (long i=0; i<rootIndex; i++) {
2621 if (marked_nodes[i]) {
2622 marked_nodes[flatParents.list_data[flatLeaves.lLength + i]] = true;
2623 theCost += ((node <long>*)(flatNodes.list_data[i]))->get_num_nodes();
2624 } else if (traversalTags && orderIndex) {
2625 long theIndex = filterL * i + orderIndex;
2626 traversalTags->list_data[theIndex/_HY_BITMASK_WIDTH_] |= bitMaskArray.masks[theIndex%_HY_BITMASK_WIDTH_];
2627 }
2628 }
2629 }
2630
2631 theCost += ((node <long>*)(flatNodes.list_data[rootIndex]))->get_num_nodes();
2632
2633 free (marked_nodes);
2634
2635 return theCost;
2636
2637 }
2638
2639
2640 //_______________________________________________________________________________________________
2641
MolecularClock(_String const & baseNode,_List & varsToConstrain) const2642 void _TheTree::MolecularClock (_String const& baseNode, _List& varsToConstrain) const {
2643 // TODO SLKP 20180803 : this needs a review
2644
2645 node<long>* topNode = nil;
2646
2647 if (baseNode.empty()) { // called Molecular Clock on the entire tree
2648 topNode = &GetRoot();
2649 _String* childNameP;
2650 if (rooted == ROOTED_LEFT) { // run separate constraint on the right child of the root
2651 MolecularClock (map_node_to_calcnode(theRoot->go_down (theRoot->get_num_nodes()))->ContextFreeName(), varsToConstrain);
2652
2653 } else if (rooted == ROOTED_RIGHT) {
2654 MolecularClock (map_node_to_calcnode(theRoot->go_down (1))->ContextFreeName(), varsToConstrain);
2655 }
2656 } else {
2657 topNode = FindNodeByName (&baseNode);
2658 }
2659
2660 if (!topNode) {
2661 HandleApplicationError (_String ("Molecular clock constraint has failed, since node '")
2662 &baseNode
2663 &"' is not a part of tree '"
2664 &*GetName() & "'");
2665 } else
2666 for (unsigned long k=1UL; k<varsToConstrain.lLength; k++) {
2667 long varIndex = LocateVarByName (*(_String*)varsToConstrain (k));
2668
2669 if (varIndex<0) {
2670 HandleApplicationError (_String ("Molecular clock constraint has failed, since variable' ") &*(_String*)varsToConstrain (k) &"' is undefined.");
2671 return ;
2672 }
2673 map_node_to_calcnode(topNode)->RecurseMC (variableNames.GetXtra(varIndex), topNode, true, rooted);
2674 }
2675 }
2676
2677
2678 //_______________________________________________________________________________________________
2679
AddNodeNamesToDS(_DataSet * ds,bool doTips,bool doInternals,char dOrS) const2680 void _TheTree::AddNodeNamesToDS (_DataSet* ds, bool doTips, bool doInternals, char dOrS) const {
2681 // TODO SLKP 20180803 needs review
2682 if (dOrS == 2 && doTips && doInternals) {
2683 AddNodeNamesToDS (ds, false, true, 0);
2684 AddNodeNamesToDS (ds, true, false, 0);
2685 return;
2686 }
2687
2688 RetrieveNodeNames (doTips, doInternals,dOrS ? _HY_TREE_TRAVERSAL_POSTORDER : _HY_TREE_TRAVERSAL_PREORDER).ForEach ([&ds] (BaseRef name, unsigned long )->void {
2689 ds->AddName (*(_String*)name);
2690 });
2691 }
2692
2693 // LF COMPUTE FUNCTIONS
2694 // TODO SLKP 20180803 these all could use a review
2695
2696 /*----------------------------------------------------------------------------------------------------------*/
ExponentiateMatrices(_List & expNodes,long tc,long catID)2697 void _TheTree::ExponentiateMatrices (_List& expNodes, long tc, long catID) {
2698 _List matrixQueue, nodesToDo;
2699
2700 _SimpleList isExplicitForm ((unsigned long)expNodes.countitems());
2701 bool hasExpForm = false;
2702
2703 for (unsigned long nodeID = 0; nodeID < expNodes.lLength; nodeID++) {
2704 long didIncrease = matrixQueue.lLength;
2705 _CalcNode* thisNode = (_CalcNode*) expNodes(nodeID);
2706 if (thisNode->RecomputeMatrix (catID, categoryCount, nil, &matrixQueue,&isExplicitForm)) {
2707 hasExpForm = true;
2708 }
2709 #ifdef _UBER_VERBOSE_DUMP
2710 if (likeFuncEvalCallCount == _UBER_VERBOSE_DUMP)
2711 printf ("NodeID %ld (%s). Old length %ld, new length %ld (%ld)\n", nodeID, thisNode->GetName()->sData, didIncrease,matrixQueue.lLength, isExplicitForm.lLength);
2712 #endif
2713 if ((didIncrease = (matrixQueue.lLength - didIncrease))) {
2714 for (long copies = 0; copies < didIncrease; copies++) {
2715 nodesToDo << thisNode;
2716 }
2717 }
2718 }
2719
2720 //printf ("%ld %d\n", nodesToDo.lLength, hasExpForm);
2721
2722 unsigned long matrixID;
2723
2724 _List * computedExponentials = hasExpForm? new _List (matrixQueue.lLength) : nil;
2725
2726 #ifdef _OPENMP
2727 unsigned long nt = cBase<20?1:(MIN(tc, matrixQueue.lLength / 3 + 1));
2728 hy_global::matrix_exp_count += matrixQueue.lLength;
2729 #endif
2730
2731 #ifdef _OPENMP
2732 #if _OPENMP>=201511
2733 #pragma omp parallel for default(shared) private (matrixID) schedule(monotonic:guided) proc_bind(spread) if (nt>1) num_threads (nt)
2734 #else
2735 #if _OPENMP>=200803
2736 #pragma omp parallel for default(shared) private (matrixID) schedule(guided) proc_bind(spread) if (nt>1) num_threads (nt)
2737 #endif
2738 #endif
2739 #endif
2740 for (matrixID = 0; matrixID < matrixQueue.lLength; matrixID++) {
2741 if (isExplicitForm.list_data[matrixID] == 0 || !hasExpForm) { // normal matrix to exponentiate
2742 ((_CalcNode*) nodesToDo(matrixID))->SetCompExp ((_Matrix*)matrixQueue(matrixID), catID, true);
2743 } else {
2744 (*computedExponentials) [matrixID] = ((_Matrix*)matrixQueue(matrixID))->Exponentiate(1., true);
2745 }
2746 }
2747
2748 if (computedExponentials) {
2749 _CalcNode * current_node = nil;
2750 _List buffered_exponentials;
2751
2752 for (unsigned long mx_index = 0; mx_index < nodesToDo.lLength; mx_index++) {
2753 if (isExplicitForm.list_data[mx_index]) {
2754 _CalcNode *next_node = (_CalcNode*) nodesToDo (mx_index);
2755 //printf ("%x %x\n", current_node, next_node);
2756 if (next_node != current_node) {
2757 if (current_node) {
2758 current_node->RecomputeMatrix (catID, categoryCount, nil, nil, nil, &buffered_exponentials);
2759 }
2760 current_node = next_node;
2761 buffered_exponentials.Clear(true);
2762 buffered_exponentials.AppendNewInstance((*computedExponentials)(mx_index));
2763 }
2764 else {
2765 buffered_exponentials.AppendNewInstance((*computedExponentials)(mx_index));
2766 }
2767 } else {
2768 if (current_node) {
2769 current_node->RecomputeMatrix (catID, categoryCount, nil, nil, nil, &buffered_exponentials);
2770 }
2771 current_node = nil;
2772 }
2773 }
2774 if (current_node) {
2775 current_node->RecomputeMatrix (catID, categoryCount, nil, nil, nil, &buffered_exponentials);
2776 }
2777 DeleteObject(computedExponentials);
2778 #ifdef _UBER_VERBOSE_DUMP_MATRICES
2779 if (likeFuncEvalCallCount == _UBER_VERBOSE_DUMP) {
2780 fprintf (stderr, "\n T_MATRIX = {");
2781 for (unsigned long nodeID = 0; nodeID < flatLeaves.lLength + flatTree.lLength - 1; nodeID++) {
2782 bool isLeaf = nodeID < flatLeaves.lLength;
2783
2784 _CalcNode * current_node = isLeaf? (((_CalcNode**) flatCLeaves.list_data)[nodeID]):
2785 (((_CalcNode**) flatTree.list_data) [nodeID - flatLeaves.lLength]);
2786 if (nodeID) {
2787 fprintf (stderr, ",");
2788 }
2789 fprintf (stderr, "\n\"%s\":%s", current_node->GetName()->get_str(), _String((_String*)current_node->GetCompExp()->toStr()).get_str());
2790
2791 }
2792 fprintf (stderr, "\n};\n");
2793 }
2794 #endif
2795
2796
2797 }
2798 }
2799
2800 /*----------------------------------------------------------------------------------------------------------*/
2801
DetermineNodesForUpdate(_SimpleList & updateNodes,_List * expNodes,long catID,long addOne,bool canClear,_AVLListX * var_to_node_map,_AVLList * changed_variables)2802 long _TheTree::DetermineNodesForUpdate (_SimpleList& updateNodes, _List* expNodes, long catID, long addOne, bool canClear,
2803 _AVLListX * var_to_node_map, _AVLList * changed_variables) {
2804 _CalcNode *currentTreeNode;
2805 long lastNodeID = -1;
2806 unsigned long tagged_node_count = 0UL;
2807
2808 // look for nodes with model changes and mark the path up to the root as needing an update
2809
2810 #define DIRECT_INDEX(N) (flatParents.list_data[N]+flatLeaves.lLength)
2811
2812 auto _handle_node = [&] (long node_id, bool do_list)->void {
2813 bool isLeaf = node_id < flatLeaves.lLength;
2814
2815 if (isLeaf) {
2816 currentTreeNode = (((_CalcNode**) flatCLeaves.list_data)[node_id]);
2817 } else {
2818 currentTreeNode = (((_CalcNode**) flatTree.list_data) [node_id - flatLeaves.lLength]);
2819 }
2820
2821 if (currentTreeNode->NeedNewCategoryExponential (catID)) {
2822 if (expNodes) {
2823 (*expNodes) << currentTreeNode;
2824 lastNodeID = node_id;
2825 } else {
2826 currentTreeNode->RecomputeMatrix (catID, categoryCount, nil);
2827 }
2828
2829 if (do_list) {
2830 nodesToUpdate.list_data[node_id] = 2;
2831 tagged_node_count++;
2832 }
2833 }
2834
2835 if (do_list && nodesToUpdate.list_data[node_id]) {
2836 long parent_index = DIRECT_INDEX(node_id);
2837 //if (nodesToUpdate.list_data[parent_index] == 0) {
2838 nodesToUpdate.list_data[parent_index] = 2;
2839 tagged_node_count++;
2840 //}
2841 }
2842 };
2843
2844 bool already_done = false;
2845 if (changed_variables && var_to_node_map && changed_variables->countitems() < 8 && forceRecalculationOnTheseBranches.empty()) {
2846 already_done = true;
2847 //_SimpleList nodes;
2848 updateNodes.RequestSpace (changed_variables->countitems());
2849 _AVLList uniques (&updateNodes);
2850
2851 for (long i = 0L; i < changed_variables->dataList->lLength; i++) {
2852 long var_index = changed_variables->dataList->list_data[i];
2853 long node_index = var_to_node_map->FindAndGetXtra ((BaseRefConst)var_index,-1L);
2854 if (node_index < 0) {
2855 already_done = false;
2856 break;
2857 }
2858 //updateNodes << node_index;
2859 uniques.InsertNumber(node_index);
2860 }
2861 //already_done = false;
2862 //if (likeFuncEvalCallCount == 7) {
2863 // ObjectToConsole(&updateNodes); NLToConsole();
2864 //}
2865 if (already_done) {
2866 if (addOne >= 0) {
2867 uniques.InsertNumber(addOne);
2868 }
2869 long i = 0L;
2870 for (; i < updateNodes.lLength; i++) {
2871 _handle_node (updateNodes.get (i), false);
2872 }
2873 //if (likeFuncEvalCallCount == 7) {
2874 // ObjectToConsole(&updateNodes); NLToConsole();
2875 //}
2876
2877
2878 for (long i2 = 0; i2 < i; i2++) {
2879 long my_index = updateNodes.get (i2),
2880 parent_index = flatParents.list_data[my_index];
2881
2882 while (parent_index >= 0) {
2883 long dir_index = parent_index + flatLeaves.lLength;
2884 //if (uniques.InsertNumber(dir_index) >= 0) { // performed insertion
2885 // node<long>* parent_node = (node<long>*)flatNodes.get (parent_index);
2886 //}
2887 uniques.InsertNumber(dir_index);
2888 parent_index = flatParents.list_data[dir_index];
2889 }
2890 }
2891
2892 _SimpleList children;
2893
2894 /*nodesToUpdate.Populate (flatTree.lLength, 0, 0);
2895
2896 for (long i = 0; i < updateNodes.lLength; i++) {
2897 long tagged_node = updateNodes.get (i);
2898 if (tagged_node >= flatLeaves.lLength) {
2899 nodesToUpdate.list_data [tagged_node - flatLeaves.lLength] = true;
2900 }
2901 }
2902
2903 for (long i = 0; i < flatParents.lLength - 1; i++) {
2904 long my_parent = flatParents.list_data [i];
2905 if (nodesToUpdate.list_data[my_parent]) {
2906 children << i;
2907 }
2908 }*/
2909
2910 for (long i = 0; i < flatParents.lLength - 1; i++) {
2911 long my_parent = flatParents.list_data [i] + flatLeaves.lLength;
2912 if (uniques.FindLong (my_parent) >= 0) {
2913 children << i;
2914 }
2915 }
2916
2917 for (long i = 0; i < children.lLength; i++) {
2918 uniques.InsertNumber(children.get (i));
2919 }
2920
2921 if (uniques.countitems() <= 32) {
2922 updateNodes.BubbleSort();
2923 } else {
2924 uniques.ReorderList();
2925 }
2926
2927 updateNodes.Pop();
2928
2929 } else {
2930 updateNodes.Clear(false);
2931 }
2932 }
2933
2934
2935 if (!already_done) {
2936 nodesToUpdate.Populate (flatLeaves.lLength + flatTree.lLength, 0, 0);
2937
2938 if (addOne >= 0) {
2939 nodesToUpdate.list_data[addOne] = 2;
2940 tagged_node_count++;
2941 }
2942
2943 if (forceRecalculationOnTheseBranches.nonempty()) {
2944 forceRecalculationOnTheseBranches.Each ([this] (long value, unsigned long) -> void {
2945 this->nodesToUpdate.list_data [value] = 2L;
2946 });
2947 tagged_node_count += forceRecalculationOnTheseBranches.countitems();
2948
2949 if (canClear) {
2950 forceRecalculationOnTheseBranches.Clear();
2951 }
2952 }
2953
2954 for (unsigned long nodeID = 0UL; nodeID < nodesToUpdate.lLength - 1UL; nodeID++) {
2955 _handle_node (nodeID, true);
2956 }
2957
2958
2959 // one more pass to pick up all DIRECT descendants of changed internal nodes
2960
2961 for (unsigned long nodeID = 0UL; nodeID < nodesToUpdate.lLength - 1UL; nodeID++)
2962 if (nodesToUpdate.list_data[nodeID] == 0 && nodesToUpdate.list_data[DIRECT_INDEX(nodeID)] == 2) {
2963 nodesToUpdate.list_data[nodeID] = 1;
2964 tagged_node_count++;
2965 }
2966
2967
2968 // write out all changed nodes
2969 updateNodes.RequestSpace(tagged_node_count);
2970 // 20200610: for larger trees, this helps with reducing memory allocations
2971 for (unsigned long nodeID = 0UL; nodeID < nodesToUpdate.lLength - 1UL; nodeID++) {
2972 if (nodesToUpdate.list_data[nodeID]) {
2973 updateNodes << nodeID;
2974 }
2975 }
2976 }
2977
2978 /*printf ("%ld ", likeFuncEvalCallCount);
2979 ObjectToConsole(&updateNodes);
2980 printf ("\n");*/
2981
2982 if (expNodes && expNodes->countitems() == 1) {
2983 return lastNodeID;
2984 }
2985
2986 return -1;
2987 }
2988
2989 /*----------------------------------------------------------------------------------------------------------*/
2990
FillInConditionals(_DataSetFilter const * theFilter,hyFloat * iNodeCache,_SimpleList * tcc)2991 void _TheTree::FillInConditionals (_DataSetFilter const* theFilter, hyFloat* iNodeCache, _SimpleList* tcc)
2992 // this utility function will simply fill in all the conditional probability vectors for internal nodes,
2993 // including those that were skipped due to column sorting optimization
2994 // this is useful to avoid code duplication for other functions (e.g. ancestral sampling) that
2995 // make use of conditional probability vectors, but would not benefit from subtree caching
2996 {
2997 if (!tcc) {
2998 return;
2999 }
3000
3001 long alphabetDimension = theFilter->GetDimension(),
3002 siteCount = theFilter->GetPatternCount();
3003
3004 for (long nodeID = 0; nodeID < flatTree.lLength; nodeID++) {
3005 hyFloat * conditionals = iNodeCache +(nodeID * siteCount) * alphabetDimension;
3006 long currentTCCIndex = siteCount * nodeID,
3007 currentTCCBit = currentTCCIndex % _HY_BITMASK_WIDTH_;
3008
3009 currentTCCIndex /= _HY_BITMASK_WIDTH_;
3010 for (long siteID = 0; siteID < siteCount; siteID++, conditionals += alphabetDimension) {
3011 if (siteID && (tcc->list_data[currentTCCIndex] & bitMaskArray.masks[currentTCCBit]) > 0) {
3012 for (long k = 0; k < alphabetDimension; k++) {
3013 conditionals[k] = conditionals[k-alphabetDimension];
3014 }
3015 }
3016 if (++currentTCCBit == _HY_BITMASK_WIDTH_) {
3017 currentTCCBit = 0;
3018 currentTCCIndex ++;
3019 }
3020
3021 }
3022 }
3023 }
3024
3025
3026
3027
3028 /*----------------------------------------------------------------------------------------------------------*/
3029
GetNodeFromFlatIndex(long index) const3030 const _CalcNode* _TheTree::GetNodeFromFlatIndex(long index) const {
3031 return index < flatLeaves.lLength ? (((_CalcNode**) flatCLeaves.list_data)[index]):
3032 (((_CalcNode**) flatTree.list_data) [index - flatLeaves.lLength]);
3033 }
3034
3035 /*----------------------------------------------------------------------------------------------------------*/
3036
ComputeLLWithBranchCache(_SimpleList & siteOrdering,long brID,hyFloat * cache,_DataSetFilter const * theFilter,long siteFrom,long siteTo,long catID,hyFloat * storageVec)3037 hyFloat _TheTree::ComputeLLWithBranchCache (
3038 _SimpleList& siteOrdering,
3039 long brID,
3040 hyFloat* cache,
3041 _DataSetFilter const* theFilter,
3042 long siteFrom,
3043 long siteTo,
3044 long catID,
3045 hyFloat* storageVec
3046 )
3047 {
3048 auto bookkeeping = [&siteOrdering, &storageVec, &theFilter] (const long siteID, const hyFloat accumulator, hyFloat& correction, hyFloat& result) -> void {
3049
3050 long direct_index = siteOrdering.list_data[siteID];
3051
3052 if (storageVec) {
3053 storageVec [direct_index] = accumulator;
3054 } else {
3055 if (accumulator <= 0.0) {
3056 //fprintf (stderr, "ZERO TERM AT SITE %ld (direct %ld) EVAL %ld\n",siteID,direct_index, likeFuncEvalCallCount);
3057 /*for (long s = 0; s < theFilter->NumberSpecies(); s++) {
3058 fprintf (stderr, "%s", theFilter->RetrieveState(direct_index, s).get_str());
3059 }
3060 fprintf (stderr, "\n");*/
3061 throw (1L+direct_index);
3062 }
3063
3064 hyFloat term;
3065 long site_frequency = theFilter->theFrequencies.get(direct_index);
3066 if ( site_frequency > 1L) {
3067 term = log(accumulator) * site_frequency - correction;
3068 } else {
3069 term = log(accumulator) - correction;
3070 }
3071
3072 /*if (likeFuncEvalCallCount == 15098) {
3073 fprintf (stderr, "CACHE, %ld, %ld, %20.15lg, %20.15lg, %20.15lg\n", likeFuncEvalCallCount, siteID, accumulator, correction, term);
3074 }*/
3075
3076 hyFloat temp_sum = result + term;
3077 correction = (temp_sum - result) - term;
3078 result = temp_sum;
3079 //result += log(accumulator) * theFilter->theFrequencies [siteOrdering.list_data[siteID]];
3080 }
3081 };
3082
3083 const unsigned long alphabetDimension = theFilter->GetDimension(),
3084 siteCount = theFilter->GetPatternCount();
3085
3086 if (siteTo > siteCount) {
3087 siteTo = siteCount;
3088 }
3089
3090 hyFloat const * branchConditionals = cache + siteFrom * alphabetDimension;
3091 hyFloat const * rootConditionals = branchConditionals + siteCount * alphabetDimension;
3092 hyFloat result = 0.0,
3093 correction = 0.0;
3094
3095
3096 //printf ("ComputeLLWithBranchCache @ %d catID = %d branchID = %d\n", likeFuncEvalCallCount, catID, brID);
3097
3098 _CalcNode const *givenTreeNode = GetNodeFromFlatIndex (brID);
3099
3100 hyFloat const * transitionMatrix = givenTreeNode->GetCompExp(catID)->theData;
3101
3102
3103 // cases by alphabet dimension
3104 /*if (likeFuncEvalCallCount == 15098) {
3105 fprintf (stderr, "\nBRANCH ID %ld (%s)\n", brID, givenTreeNode->GetName()->get_str());
3106 for (long e = 0; e < 16; e++) {
3107 fprintf (stderr, "%ld => %lg\n", transitionMatrix[e]);
3108 }
3109 }*/
3110
3111 try {
3112 switch (alphabetDimension) {
3113 /****
3114
3115 NUCLEOTIDES
3116
3117 ****/
3118 case 4UL: {
3119 #ifdef _SLKP_USE_AVX_INTRINSICS
3120 __m256d tmatrix_transpose [4] = {
3121 (__m256d) {transitionMatrix[0],transitionMatrix[4],transitionMatrix[8],transitionMatrix[12]},
3122 (__m256d) {transitionMatrix[1],transitionMatrix[5],transitionMatrix[9],transitionMatrix[13]},
3123 (__m256d) {transitionMatrix[2],transitionMatrix[6],transitionMatrix[10],transitionMatrix[14]},
3124 (__m256d) {transitionMatrix[3],transitionMatrix[7],transitionMatrix[11],transitionMatrix[15]}
3125 };
3126 #endif
3127
3128 for (unsigned long siteID = siteFrom; siteID < siteTo; siteID++) {
3129 hyFloat accumulator = 0.;
3130 #ifdef _SLKP_USE_AVX_INTRINSICS
3131 __m256d root_c = _mm256_loadu_pd (rootConditionals),
3132 probs = _mm256_loadu_pd (theProbs),
3133 b_cond0 = _mm256_set1_pd(branchConditionals[0]),
3134 b_cond1 = _mm256_set1_pd(branchConditionals[1]),
3135 b_cond2 = _mm256_set1_pd(branchConditionals[2]),
3136 b_cond3 = _mm256_set1_pd(branchConditionals[3]),
3137 s01 = _mm256_add_pd ( _mm256_mul_pd (b_cond0, tmatrix_transpose[0]), _mm256_mul_pd (b_cond1, tmatrix_transpose[1])),
3138 s23 = _mm256_add_pd ( _mm256_mul_pd (b_cond2, tmatrix_transpose[2]), _mm256_mul_pd (b_cond3, tmatrix_transpose[3]));
3139 accumulator = _avx_sum_4(_mm256_mul_pd (_mm256_mul_pd (root_c, probs), _mm256_add_pd (s01,s23)));
3140
3141 #else
3142 accumulator = rootConditionals[0] * theProbs[0] *
3143 (branchConditionals[0] * transitionMatrix[0] + branchConditionals[1] * transitionMatrix[1] + branchConditionals[2] * transitionMatrix[2] + branchConditionals[3] * transitionMatrix[3]) +
3144 rootConditionals[1] * theProbs[1] *
3145 (branchConditionals[0] * transitionMatrix[4] + branchConditionals[1] * transitionMatrix[5] + branchConditionals[2] * transitionMatrix[6] + branchConditionals[3] * transitionMatrix[7]) +
3146 rootConditionals[2] * theProbs[2] *
3147 (branchConditionals[0] * transitionMatrix[8] + branchConditionals[1] * transitionMatrix[9] + branchConditionals[2] * transitionMatrix[10] + branchConditionals[3] * transitionMatrix[11]) +
3148 rootConditionals[3] * theProbs[3] *
3149 (branchConditionals[0] * transitionMatrix[12] + branchConditionals[1] * transitionMatrix[13] + branchConditionals[2] * transitionMatrix[14] + branchConditionals[3] * transitionMatrix[15]);
3150 #endif
3151 rootConditionals += 4UL;
3152 branchConditionals += 4UL;
3153 bookkeeping (siteID, accumulator, correction, result);
3154 } // siteID
3155 }
3156 break;
3157 /****
3158
3159 AMINOACIDS
3160
3161 ****/
3162 case 20UL: {
3163 for (unsigned long siteID = siteFrom; siteID < siteTo; siteID++) {
3164 hyFloat accumulator = 0.;
3165 #ifdef _SLKP_USE_AVX_INTRINSICS
3166 __m256d bc_vector[5] = {_mm256_loadu_pd(branchConditionals),
3167 _mm256_loadu_pd(branchConditionals+4UL),
3168 _mm256_loadu_pd(branchConditionals+8UL),
3169 _mm256_loadu_pd(branchConditionals+12UL),
3170 _mm256_loadu_pd(branchConditionals+16UL)};
3171
3172
3173 hyFloat const * tm = transitionMatrix;
3174
3175 for (unsigned long p = 0UL; p < 20UL; p++, rootConditionals++) {
3176
3177 __m256d t_matrix[5] = { _mm256_loadu_pd(tm),
3178 _mm256_loadu_pd(tm+4UL),
3179 _mm256_loadu_pd(tm+8UL),
3180 _mm256_loadu_pd(tm+12UL),
3181 _mm256_loadu_pd(tm+16UL)};
3182
3183
3184 t_matrix[0] = _mm256_mul_pd(t_matrix[0], bc_vector[0]);
3185 t_matrix[1] = _mm256_mul_pd(t_matrix[1], bc_vector[1]);
3186 t_matrix[2] = _mm256_mul_pd(t_matrix[2], bc_vector[2]);
3187 t_matrix[3] = _mm256_mul_pd(t_matrix[3], bc_vector[3]);
3188 t_matrix[4] = _mm256_mul_pd(t_matrix[4], bc_vector[4]);
3189
3190 t_matrix[0] = _mm256_add_pd (t_matrix[0],t_matrix[1]);
3191 t_matrix[1] = _mm256_add_pd (t_matrix[2],t_matrix[3]);
3192 t_matrix[3] = _mm256_add_pd (t_matrix[0],t_matrix[1]);
3193
3194 tm += 20UL;
3195
3196 accumulator += *rootConditionals * theProbs[p] * _avx_sum_4(_mm256_add_pd (t_matrix[3],t_matrix[4]));
3197 }
3198 #else // _SLKP_USE_AVX_INTRINSICS
3199 unsigned long rmx = 0UL;
3200 for (unsigned long p = 0UL; p < 20UL; p++,rootConditionals++) {
3201 hyFloat r2 = 0.;
3202
3203 for (unsigned long c = 0UL; c < 20UL; c+=4UL, rmx +=4UL) {
3204 r2 += (branchConditionals[c] * transitionMatrix[rmx] +
3205 branchConditionals[c+1] * transitionMatrix[rmx+1]) +
3206 (branchConditionals[c+2] * transitionMatrix[rmx+2] +
3207 branchConditionals[c+3] * transitionMatrix[rmx+3]);
3208 }
3209
3210 accumulator += *rootConditionals * theProbs[p] * r2;
3211 }
3212 #endif // _SLKP_USE_AVX_INTRINSICS
3213 branchConditionals += 20UL;
3214 bookkeeping (siteID, accumulator, correction, result);
3215
3216 } // siteID
3217
3218 } // case 20
3219 break;
3220 /****
3221
3222 CODONS
3223
3224 ****/
3225 case 60UL:
3226 case 61UL:
3227 case 62UL:
3228 case 63UL: {
3229 for (unsigned long siteID = siteFrom; siteID < siteTo; siteID++) {
3230 hyFloat accumulator = 0.;
3231
3232 unsigned long rmx = 0UL;
3233 for (unsigned long p = 0UL; p < alphabetDimension; p++,rootConditionals++) {
3234 hyFloat r2 = 0.;
3235 unsigned long c = 0UL;
3236
3237 #ifdef _SLKP_USE_AVX_INTRINSICS
3238
3239 __m256d sum256 = _mm256_setzero_pd ();
3240
3241 for (; c < 60UL; c+=12UL, rmx +=12UL) {
3242
3243 __m256d branches0 = _mm256_loadu_pd (branchConditionals+c),
3244 branches1 = _mm256_loadu_pd (branchConditionals+c+4),
3245 branches2 = _mm256_loadu_pd (branchConditionals+c+8),
3246 matrix0 = _mm256_loadu_pd (transitionMatrix+rmx),
3247 matrix1 = _mm256_loadu_pd (transitionMatrix+rmx+4),
3248 matrix2 = _mm256_loadu_pd (transitionMatrix+rmx+8);
3249
3250 branches0 = _mm256_mul_pd(branches0, matrix0);
3251 branches1 = _mm256_mul_pd(branches1, matrix1);
3252 branches2 = _mm256_mul_pd(branches2, matrix2);
3253
3254 branches0 = _mm256_add_pd (branches0,branches2);
3255 sum256 = _mm256_add_pd (branches0,_mm256_add_pd (sum256, branches1));
3256 }
3257
3258 r2 = _avx_sum_4(sum256);
3259
3260 #else // _SLKP_USE_AVX_INTRINSICS
3261 for (; c < 60UL; c+=4UL, rmx +=4UL) {
3262 r2 += (branchConditionals[c] * transitionMatrix[rmx] +
3263 branchConditionals[c+1] * transitionMatrix[rmx+1]) +
3264 (branchConditionals[c+2] * transitionMatrix[rmx+2] +
3265 branchConditionals[c+3] * transitionMatrix[rmx+3]);
3266 }
3267 #endif
3268
3269 for (; c < alphabetDimension; c++, rmx ++) {
3270 r2 += branchConditionals[c] * transitionMatrix[rmx];
3271 }
3272
3273 accumulator += *rootConditionals * theProbs[p] * r2;
3274 }
3275
3276 branchConditionals += alphabetDimension;
3277 bookkeeping (siteID, accumulator, correction, result);
3278
3279 }
3280 } // cases 60-63
3281 break;
3282 default: { // valid alphabetDimension >= 2
3283
3284 if (alphabetDimension % 2) { // odd
3285 unsigned long alphabetDimension_minus1 = alphabetDimension-1;
3286 for (unsigned long siteID = siteFrom; siteID < siteTo; siteID++) {
3287 hyFloat accumulator = 0.;
3288
3289 unsigned long rmx = 0UL;
3290 for (unsigned long p = 0UL; p < alphabetDimension; p++,rootConditionals++) {
3291 hyFloat r2 = 0.;
3292
3293 for (unsigned long c = 0UL; c < alphabetDimension_minus1; c+=2UL, rmx +=2UL) {
3294 r2 += branchConditionals[c] * transitionMatrix[rmx] +
3295 branchConditionals[c+1] * transitionMatrix[rmx+1];
3296 }
3297
3298 r2 += branchConditionals[alphabetDimension_minus1] * transitionMatrix[rmx++];
3299
3300 accumulator += *rootConditionals * theProbs[p] * r2;
3301 }
3302
3303 branchConditionals += alphabetDimension;
3304 bookkeeping (siteID, accumulator, correction, result);
3305
3306 }
3307 } else {
3308 for (unsigned long siteID = siteFrom; siteID < siteTo; siteID++) {
3309 hyFloat accumulator = 0.;
3310
3311 unsigned long rmx = 0UL;
3312 for (unsigned long p = 0UL; p < alphabetDimension; p++,rootConditionals++) {
3313 hyFloat r2 = 0.;
3314
3315 for (unsigned long c = 0UL; c < alphabetDimension; c+=2UL, rmx +=2UL) {
3316 r2 += branchConditionals[c] * transitionMatrix[rmx] +
3317 branchConditionals[c+1] * transitionMatrix[rmx+1];
3318 }
3319
3320 accumulator += *rootConditionals * theProbs[p] * r2;
3321 }
3322
3323 branchConditionals += alphabetDimension;
3324 bookkeeping (siteID, accumulator, correction, result);
3325
3326 }
3327 }
3328
3329 } // default
3330 } // switch (alphabetDimension)
3331 } catch (long site) {
3332 #pragma omp critical
3333 {
3334 hy_global::ReportWarning (_String("Site ") & _String(site) & " evaluated to a 0 probability in ComputeLLWithBranchCache");
3335 }
3336 return -INFINITY;
3337 }
3338 return result;
3339 }
3340
3341 /*----------------------------------------------------------------------------------------------------------*/
3342
ComputeTwoSequenceLikelihood(_SimpleList & siteOrdering,_DataSetFilter const * theFilter,long * lNodeFlags,_Vector * lNodeResolutions,long siteFrom,long siteTo,long catID,hyFloat * storageVec)3343 hyFloat _TheTree::ComputeTwoSequenceLikelihood
3344 (
3345 _SimpleList & siteOrdering,
3346 _DataSetFilter const* theFilter,
3347 long * lNodeFlags,
3348 _Vector* lNodeResolutions,
3349 long siteFrom,
3350 long siteTo,
3351 long catID,
3352 hyFloat* storageVec
3353 )
3354 // the updateNodes flags the nodes (leaves followed by inodes in the same order as flatLeaves and flatNodes)
3355 // that must be recomputed
3356 {
3357 // process the leaves first
3358
3359 long alphabetDimension = theFilter->GetDimension(),
3360 siteCount = theFilter->GetPatternCount(),
3361 alphabetDimensionmod4 = alphabetDimension-alphabetDimension%4;
3362
3363 _CalcNode *theNode = ((_CalcNode*) flatCLeaves (0));
3364 hyFloat * _hprestrict_ transitionMatrix
3365 = theNode->GetCompExp(catID)->theData,
3366 result = 0.;
3367
3368 if (siteTo > siteCount) {
3369 siteTo = siteCount;
3370 }
3371
3372 for (long siteID = siteFrom; siteID < siteTo; siteID++) {
3373 hyFloat *tMatrix = transitionMatrix,
3374 sum = 0.;
3375
3376 long siteState1 = lNodeFlags[siteOrdering.list_data[siteID]],
3377 siteState2 = lNodeFlags[siteCount + siteOrdering.list_data[siteID]];
3378
3379 if (siteState1 >= 0)
3380 // a single character state; sweep down the appropriate column
3381 {
3382 if (siteState2 >= 0) { // both completely resolved;
3383 sum = tMatrix[siteState1*alphabetDimension + siteState2];
3384 } else { // first resolved, second is not
3385 hyFloat* childVector = lNodeResolutions->theData + (-siteState2-1) * alphabetDimension;
3386 tMatrix += siteState1*alphabetDimension;
3387 if (alphabetDimension == 4) { // special case for nuc data
3388 sum = tMatrix[0] * childVector[0] +
3389 tMatrix[1] * childVector[1] +
3390 tMatrix[2] * childVector[2] +
3391 tMatrix[3] * childVector[3];
3392 } else {
3393 for (long c = 0; c < alphabetDimensionmod4; c+=4) // 4 - unroll the loop
3394 sum += tMatrix[c] * childVector[c] +
3395 tMatrix[c+1] * childVector[c+1] +
3396 tMatrix[c+2] * childVector[c+2] +
3397 tMatrix[c+3] * childVector[c+3];
3398
3399 for (long c = alphabetDimensionmod4; c < alphabetDimension; c++) {
3400 sum += tMatrix[c] * childVector[c];
3401 }
3402 }
3403 }
3404 sum *= theProbs[siteState1];
3405 } else {
3406 if (siteState2 >=0 ) { // second resolved, but not the first
3407 hyFloat* childVector = lNodeResolutions->theData + (-siteState1-1) * alphabetDimension;
3408 tMatrix += siteState2;
3409 if (alphabetDimension == 4) { // special case for nuc data
3410 sum = tMatrix[0] * childVector[0] * theProbs[0]+
3411 tMatrix[4] * childVector[1] * theProbs[1]+
3412 tMatrix[8] * childVector[2] * theProbs[2]+
3413 tMatrix[12] * childVector[3] * theProbs[3];
3414
3415 } else {
3416 for (long c = 0; c < alphabetDimensionmod4; c+=4, tMatrix += 4*alphabetDimension) // 4 - unroll the loop
3417 sum += tMatrix[0] * childVector[c] * theProbs[c]+
3418 tMatrix[alphabetDimension] * childVector[c+1] * theProbs[c+1]+
3419 tMatrix[alphabetDimension+alphabetDimension] * childVector[c+2] * theProbs[c+2]+
3420 tMatrix[alphabetDimension+alphabetDimension+alphabetDimension] * childVector[c+3] * theProbs[c+3];
3421
3422 for (long c = alphabetDimensionmod4; c < alphabetDimension; c++, tMatrix += alphabetDimension) {
3423 sum += tMatrix[0] * childVector[c] * theProbs[c];
3424 }
3425 }
3426 } else
3427 // both unresolved
3428 {
3429 hyFloat *childVector1 = lNodeResolutions->theData + (-siteState1-1) * alphabetDimension,
3430 *childVector2 = lNodeResolutions->theData + (-siteState2-1) * alphabetDimension;
3431
3432 if (alphabetDimension == 4) { // special case for nuc data
3433 sum = (tMatrix[0] * childVector2[0] + tMatrix[1] * childVector2[1] + tMatrix[2] * childVector2[2] + tMatrix[3] * childVector2[3]) * childVector1[0] * theProbs[0]+
3434 (tMatrix[4] * childVector2[0] + tMatrix[5] * childVector2[1] + tMatrix[6] * childVector2[2] + tMatrix[7] * childVector2[3]) * childVector1[1] * theProbs[1]+
3435 (tMatrix[8] * childVector2[0] + tMatrix[9] * childVector2[1] + tMatrix[10] * childVector2[2] + tMatrix[11] * childVector2[3]) * childVector1[2] * theProbs[2] +
3436 (tMatrix[12] * childVector2[0] + tMatrix[13] * childVector2[1] + tMatrix[14] * childVector2[2] + tMatrix[15] * childVector2[3]) * childVector1[3] * theProbs[3];
3437
3438 } else {
3439 for (long r = 0; r < alphabetDimension; r++) { // 4 - unroll the loop
3440 if (childVector1[r] > 0.0) {
3441 hyFloat sum2 = 0.0;
3442 for (long c = 0; c < alphabetDimensionmod4; c+=4, tMatrix += 4) // 4 - unroll the loop
3443 sum2 += tMatrix[0] * childVector2[c] +
3444 tMatrix[1] * childVector2[c+1]+
3445 tMatrix[2] * childVector2[c+2]+
3446 tMatrix[3] * childVector2[c+3];
3447
3448 for (long c = alphabetDimensionmod4; c < alphabetDimension; c++, tMatrix ++) {
3449 sum2 += tMatrix[0] * childVector2[c];
3450 }
3451
3452 sum += sum2 * childVector1[r] * theProbs[r];
3453 } else {
3454 tMatrix += alphabetDimension;
3455 }
3456 }
3457 }
3458 }
3459 }
3460 if (storageVec) {
3461 storageVec [siteOrdering.list_data[siteID]] = sum;
3462 } else {
3463 if (sum <= 0.0) {
3464 return -INFINITY;
3465 } else {
3466 //printf ("%d: %g\n", siteID, sum);
3467 result += log(sum) * theFilter->theFrequencies.get (siteOrdering.list_data[siteID]);
3468 }
3469 }
3470 }
3471
3472 return result;
3473 }
3474
3475
3476
3477 //_______________________________________________________________________________________________
3478
SampleAncestorsBySequence(_DataSetFilter const * dsf,_SimpleList const & siteOrdering,node<long> * currentNode,_AVLListX const * nodeToIndex,hyFloat const * iNodeCache,_List & result,_SimpleList * parentStates,_List & expandedSiteMap,hyFloat const * catAssignments,long catCount)3479 void _TheTree::SampleAncestorsBySequence (_DataSetFilter const* dsf, _SimpleList const& siteOrdering, node<long>* currentNode, _AVLListX const* nodeToIndex, hyFloat const* iNodeCache,
3480 _List& result, _SimpleList* parentStates, _List& expandedSiteMap, hyFloat const* catAssignments, long catCount)
3481
3482 // must be called initially with the root node
3483
3484
3485 // dsf: the filter to sample from
3486 // siteOrdering: the map from cache ordering to actual pattern ordering
3487 // currentNode: the node index to sample for
3488 // nodeToIndex: an AVL that maps the address of an internal node pointed to by node<long> to its order in the tree postorder traversal
3489 // iNodeCache: internal node likelihood caches
3490 // results: the list that will store sampled strings
3491 // parentStates: sampled states for the parent of the current node
3492 // expandedSiteMap: a list of simple lists giving site indices for each unique column pattern in the alignment
3493 // catAssignments: a vector assigning a (partition specific) rate category to each site (nil if no rate variation)
3494 // catCount: the number of rate classes
3495
3496 // this needs to be updated to deal with traversal caches!
3497 {
3498 long childrenCount = currentNode->get_num_nodes();
3499
3500 if (childrenCount) {
3501 long siteCount = dsf->GetPatternCount (),
3502 alphabetDimension = dsf->GetDimension (),
3503 nodeIndex = nodeToIndex->GetXtra (nodeToIndex->Find ((BaseRef)currentNode)),
3504 unitLength = dsf->GetUnitLength(),
3505 catBlockShifter = catAssignments?(dsf->GetPatternCount()*GetINodeCount()):0;
3506
3507
3508 _CalcNode * currentTreeNode = ((_CalcNode*) flatTree (nodeIndex));
3509 _SimpleList sampledStates (dsf->GetSiteCountInUnits (), 0, 0);
3510
3511 hyFloat const * transitionMatrix = (catAssignments|| !parentStates)?nil:currentTreeNode->GetCompExp()->theData;
3512 hyFloat const * conditionals = catAssignments?nil:(iNodeCache + nodeIndex * siteCount * alphabetDimension);
3513 hyFloat * cache = new hyFloat [alphabetDimension];
3514
3515 for (long pattern = 0; pattern < siteCount; pattern++) {
3516 _SimpleList* patternMap = (_SimpleList*) expandedSiteMap (siteOrdering.list_data[pattern]);
3517 if (catAssignments) {
3518 long localCatID = catAssignments[siteOrdering.list_data[pattern]];
3519 if (parentStates) {
3520 transitionMatrix = currentTreeNode->GetCompExp(localCatID)->theData;
3521 }
3522
3523 conditionals = iNodeCache + localCatID*alphabetDimension*catBlockShifter + (pattern + nodeIndex * siteCount) * alphabetDimension;
3524 }
3525
3526 for (long site = 0; site < patternMap->lLength; site++) {
3527 long siteID = patternMap->list_data[site];
3528
3529 hyFloat randVal = genrand_real2(),
3530 totalSum = 0.;
3531
3532 hyFloat const * matrixRow;
3533
3534 if (parentStates == nil) {
3535 matrixRow = theProbs;
3536 } else {
3537 matrixRow = transitionMatrix + parentStates->list_data[siteID] * alphabetDimension;
3538 }
3539
3540 for (long i = 0; i<alphabetDimension; i++) {
3541 totalSum += (cache[i] = matrixRow[i]*conditionals[i]);
3542 }
3543
3544 randVal *= totalSum;
3545 totalSum = 0.0;
3546 long sampledChar = -1;
3547 while (totalSum < randVal) {
3548 sampledChar ++;
3549 totalSum += cache[sampledChar];
3550 }
3551
3552 sampledStates.list_data[siteID] = sampledChar;
3553 }
3554
3555 if (catAssignments == nil) {
3556 conditionals += alphabetDimension;
3557 }
3558 }
3559
3560 delete [] cache;
3561
3562 _SimpleList conversion;
3563 _AVLListXL conversionAVL (&conversion);
3564
3565 _StringBuffer * sampledSequence = new _StringBuffer (siteCount*unitLength);
3566 _String letterValue ((unsigned long) unitLength);
3567 for (long charIndexer = 0; charIndexer < sampledStates.countitems(); charIndexer++) {
3568 dsf->ConvertCodeToLettersBuffered (dsf->CorrectCode(sampledStates.list_data[charIndexer]), unitLength, letterValue, &conversionAVL);
3569 (*sampledSequence) << letterValue;
3570 }
3571 sampledSequence->TrimSpace();
3572 result.AppendNewInstance(sampledSequence);
3573 //printf ("%d: %s\n", nodeIndex, sampledSequence->sData);
3574
3575 for (long child = 1L; child <= childrenCount; child ++) {
3576 SampleAncestorsBySequence (dsf, siteOrdering, currentNode->go_down(child), nodeToIndex, iNodeCache, result, &sampledStates, expandedSiteMap, catAssignments, catCount);
3577 }
3578 }
3579 }
3580
3581 //_______________________________________________________________________________________________
3582
RecoverAncestralSequences(_DataSetFilter const * dsf,_SimpleList const & siteOrdering,_List const & expandedSiteMap,hyFloat * iNodeCache,hyFloat const * catAssignments,long catCount,long * lNodeFlags,_Vector * lNodeResolutions,bool alsoDoLeaves)3583 _List* _TheTree::RecoverAncestralSequences (_DataSetFilter const* dsf,
3584 _SimpleList const& siteOrdering,
3585 _List const& expandedSiteMap,
3586 hyFloat * iNodeCache,
3587 hyFloat const* catAssignments,
3588 long catCount,
3589 long* lNodeFlags,
3590 _Vector * lNodeResolutions,
3591 bool alsoDoLeaves
3592 )
3593
3594
3595 // dsf: the filter to sample from
3596 // siteOrdering: the map from cache ordering to actual pattern ordering
3597 // expandedSiteMap: a list of simple lists giving site indices for each unique column pattern in the alignment
3598 // iNodeCache: internal node likelihood caches
3599 // catAssignments: a vector assigning a (partition specific) rate category to each site
3600 // catCount: the number of rate classes
3601 // alsoDoLeaves: if true, also return ML reconstruction of observed (or partially observed) sequences
3602 {
3603 long patternCount = dsf->GetPatternCount (),
3604 alphabetDimension = dsf->GetDimension (),
3605 unitLength = dsf->GetUnitLength (),
3606 iNodeCount = GetINodeCount (),
3607 leafCount = GetLeafCount (),
3608 siteCount = dsf->GetSiteCountInUnits (),
3609 allNodeCount = 0,
3610 stateCacheDim = (alsoDoLeaves? (iNodeCount + leafCount): (iNodeCount));
3611
3612 long *stateCache = new long [patternCount*(iNodeCount-1)*alphabetDimension],
3613 *leafBuffer = new long [(alsoDoLeaves?leafCount*patternCount:1)*alphabetDimension],
3614 *initiaStateCache = stateCache;
3615
3616 // a Patterns x Int-Nodes x CharStates integer table
3617 // with the best character assignment for node i given that its parent state is j for a given site
3618
3619 hyFloat *buffer = new hyFloat [alphabetDimension];
3620 // iNodeCache will be OVERWRITTEN with conditional pair (i,j) conditional likelihoods
3621
3622
3623 _SimpleList taggedInternals (iNodeCount, 0, 0),
3624 postToIn;
3625
3626 MapPostOrderToInOrderTraversal (postToIn);
3627 // all nodes except the root
3628
3629 allNodeCount = iNodeCount + leafCount - 1;
3630
3631 for (long nodeID = 0; nodeID < allNodeCount; nodeID++) {
3632 long parent_index = flatParents.get (nodeID),
3633 node_index = nodeID;
3634
3635 bool is_leaf = nodeID < flatLeaves.countitems();
3636
3637
3638 if (!is_leaf) {
3639 node_index -= flatLeaves.countitems();
3640 AddBranchToForcedRecomputeList (node_index);
3641 }
3642
3643 hyFloat * parentConditionals = iNodeCache + parent_index * alphabetDimension * patternCount;
3644
3645 if (taggedInternals.get(parent_index) == 0L) {
3646 // mark the parent for update and clear its conditionals if needed
3647 taggedInternals[parent_index] = 1L;
3648 InitializeArray(parentConditionals, patternCount*alphabetDimension, 1.);
3649 }
3650
3651 _CalcNode * tree_node_object = is_leaf? ((_CalcNode*) flatCLeaves (node_index)):((_CalcNode*) flatTree (node_index));
3652 hyFloat const* transition_matrix = nil;
3653
3654 if (!catAssignments) {
3655 _Matrix* comp_exp = tree_node_object->GetCompExp();
3656 if (!comp_exp) {
3657 hy_global::HandleApplicationError(_String ("Internal error in ") & __PRETTY_FUNCTION__ & ". Transition matrix not computed for " & *tree_node_object->GetName());
3658 delete [] stateCache; delete [] leafBuffer; delete [] buffer;
3659 return nil;
3660 }
3661 transition_matrix = comp_exp->theData;
3662 }
3663
3664 // this will need to be toggled on a per site basis
3665 hyFloat * childVector;
3666
3667 if (!is_leaf) {
3668 childVector = iNodeCache + (node_index * patternCount) * alphabetDimension;
3669 }
3670
3671 for (long siteID = 0; siteID < patternCount; siteID++, parentConditionals += alphabetDimension) {
3672 if (catAssignments) {
3673 transition_matrix = tree_node_object->GetCompExp(catAssignments[siteOrdering.list_data[siteID]])->theData;
3674 }
3675
3676 hyFloat const *tMatrix = transition_matrix;
3677 if (is_leaf) {
3678 long siteState = lNodeFlags[node_index*patternCount + siteOrdering.list_data[siteID]] ;
3679 if (siteState >= 0L) { // a fully resolved leaf
3680 tMatrix += siteState;
3681 for (long k = 0; k < alphabetDimension; k++, tMatrix += alphabetDimension) {
3682 parentConditionals[k] *= *tMatrix;
3683 }
3684 if (alsoDoLeaves) {
3685 InitializeArray (leafBuffer, alphabetDimension, (const long)siteState);
3686 leafBuffer += alphabetDimension;
3687 }
3688
3689 continue;
3690 } else {// an ambiguous leaf
3691 childVector = lNodeResolutions->theData + (-siteState-1L) * alphabetDimension;
3692 }
3693
3694 }
3695
3696 // now repopulate this vector as necessary -- if we are here this means
3697 // that the subtree below has been completely processed,
3698 // the i-th cell of childVector contains the likelihood of the _optimal_
3699 // assignment in the subtree below given that the character at the current
3700 // node is i.
3701
3702 // hence, given parent state 'p', we optimize
3703 // max_i pr (p->i) childVector [i] and store it in the p cell of vector childVector
3704
3705 hyFloat overallMax = 0.0;
3706
3707 long *stateBuffer = is_leaf?leafBuffer:stateCache;
3708
3709 // check for degeneracy
3710
3711 bool completely_unresolved = ArrayAll (childVector, alphabetDimension, [] (hyFloat x, unsigned long) {return x == 1.;});
3712
3713 if (completely_unresolved) {
3714 InitializeArray(stateBuffer, alphabetDimension, -1L);
3715 } else {
3716 for (long p = 0L; p < alphabetDimension; p++) {
3717 hyFloat max_lik = 0.;
3718 long max_idx = 0L;
3719
3720 for (long c = 0L; c < alphabetDimension; c++) {
3721 hyFloat thisV = tMatrix[c] * childVector[c];
3722 if (thisV > max_lik) {
3723 max_lik = thisV;
3724 max_idx = c;
3725 }
3726 }
3727
3728 stateBuffer [p] = max_idx;
3729 buffer [p] = max_lik;
3730
3731 if (max_lik > overallMax) {
3732 overallMax = max_lik;
3733 }
3734
3735 tMatrix += alphabetDimension;
3736 }
3737
3738 if (overallMax > 0.0 && overallMax < _lfScalingFactorThreshold) {
3739 for (long k = 0L; k < alphabetDimension; k++) {
3740 buffer[k] *= _lfScalerUpwards;
3741 }
3742 }
3743
3744 // buffer[p] now contains the maximum likelihood of the tree
3745 // from this point forward given that parent state is p
3746 // and stateBuffer[p] stores the maximizing assignment
3747 // for this node
3748
3749 for (long k = 0; k < alphabetDimension; k++) {
3750 if (stateBuffer[k] >= 0L) {
3751 parentConditionals[k] *= buffer[k];
3752 }
3753 }
3754 }
3755
3756 if (is_leaf) {
3757 if (alsoDoLeaves) {
3758 leafBuffer += alphabetDimension;
3759 }
3760 } else {
3761 stateCache += alphabetDimension;
3762 }
3763
3764 childVector += alphabetDimension;
3765 }
3766 }
3767
3768 _List *result = new _List;
3769 for (long k = 0; k < stateCacheDim; k++) {
3770 result->AppendNewInstance (new _String((unsigned long)siteCount*unitLength));
3771 }
3772
3773 hyFloat * _hprestrict_ rootConditionals = iNodeCache + alphabetDimension * ((iNodeCount-1) * patternCount);
3774 _SimpleList parentStates (stateCacheDim,0,0),
3775 conversion;
3776
3777 stateCache -= patternCount*(iNodeCount-1)*alphabetDimension;
3778 if (alsoDoLeaves) {
3779 leafBuffer -= patternCount*leafCount*alphabetDimension;
3780 }
3781
3782 _AVLListXL conversionAVL (&conversion);
3783 _String codeBuffer ((unsigned long)unitLength);
3784
3785 for (long siteID = 0; siteID < patternCount; siteID++, rootConditionals += alphabetDimension) {
3786 hyFloat max_lik = 0.;
3787 long max_idx = 0;
3788
3789 long howManyOnes = 0;
3790 for (long k = 0; k < alphabetDimension; k++) {
3791 howManyOnes += rootConditionals[k]==1.;
3792 }
3793
3794 _SimpleList const* patternMap = (_SimpleList const*) expandedSiteMap.GetItem(siteOrdering.list_data[siteID]);
3795
3796 if (howManyOnes != alphabetDimension) {
3797 for (long c = 0; c < alphabetDimension; c++) {
3798 hyFloat thisV = theProbs[c] * rootConditionals[c];
3799 if (thisV > max_lik) {
3800 max_lik = thisV;
3801 max_idx = c;
3802 }
3803 }
3804
3805 parentStates.list_data[iNodeCount-1] = max_idx;
3806 for (long nodeID = iNodeCount-2; nodeID >=0 ; nodeID--) {
3807 long parentState = parentStates.list_data[flatParents.list_data [nodeID+flatLeaves.lLength]];
3808 if (parentState == -1) {
3809 parentStates.list_data[nodeID] = -1;
3810 } else {
3811 parentStates.list_data[nodeID] = stateCache[(patternCount*nodeID+siteID)*alphabetDimension + parentState];
3812 }
3813 }
3814 if (alsoDoLeaves)
3815 for (long nodeID = 0; nodeID <leafCount ; nodeID++) {
3816 long parentState = parentStates.list_data[flatParents.list_data [nodeID]];
3817 if (parentState == -1) {
3818 parentStates.list_data[nodeID+iNodeCount] = -1;
3819 } else {
3820 parentStates.list_data[nodeID+iNodeCount] = leafBuffer[(patternCount*nodeID+siteID)*alphabetDimension + parentState];
3821 }
3822 }
3823 } else {
3824 parentStates.Populate(stateCacheDim,-1,0);
3825 }
3826
3827 for (long nodeID = 0; nodeID < stateCacheDim ; nodeID++) {
3828 dsf->ConvertCodeToLettersBuffered (dsf->CorrectCode(parentStates.list_data[nodeID]), unitLength, codeBuffer, &conversionAVL);
3829 _String *sequence = (_String*) (*result)(nodeID<iNodeCount?postToIn.list_data[nodeID]:nodeID);
3830
3831 for (long site = 0; site < patternMap->lLength; site++) {
3832 unsigned long offset = patternMap->list_data[site]*unitLength;
3833 for (long charS = 0; charS < unitLength; charS ++) {
3834 sequence->set_char (offset + charS, codeBuffer.char_at(charS));
3835 }
3836 }
3837 }
3838 }
3839
3840 delete [] initiaStateCache;
3841 delete [] leafBuffer;
3842 delete [] buffer;
3843
3844 return result;
3845 }
3846
3847 //_______________________________________________________________________________________________
3848
SetupCategoryMapsForNodes(_List & containerVariables,_SimpleList & classCounter,_SimpleList & multipliers)3849 void _TheTree::SetupCategoryMapsForNodes (_List& containerVariables, _SimpleList& classCounter, _SimpleList& multipliers)
3850 {
3851 _TreeIterator ti (this, _HY_TREE_TRAVERSAL_POSTORDER);
3852 while (_CalcNode* iterator = ti.Next()) {
3853 iterator->SetupCategoryMap (containerVariables,classCounter,multipliers);
3854 }
3855 }
3856
3857
3858
3859 //_______________________________________________________________________________________________
3860
Process3TaxonNumericFilter(_DataSetFilterNumeric * dsf,long catID)3861 hyFloat _TheTree::Process3TaxonNumericFilter (_DataSetFilterNumeric* dsf, long catID)
3862 {
3863
3864 hyFloat *l0 = dsf->probabilityVectors.theData +
3865 dsf->categoryShifter * catID + dsf->theNodeMap.list_data[0]*dsf->shifter,
3866 *l1 = dsf->probabilityVectors.theData +
3867 dsf->categoryShifter * catID + dsf->theNodeMap.list_data[1]*dsf->shifter,
3868 *l2 = dsf->probabilityVectors.theData +
3869 dsf->categoryShifter * catID + dsf->theNodeMap.list_data[2]*dsf->shifter,
3870 * matrix0 = ((_CalcNode*)(LocateVar(theRoot->get_node(1)->in_object)))->GetCompExp(catID)->theData,
3871 * matrix1 = ((_CalcNode*)(LocateVar(theRoot->get_node(2)->in_object)))->GetCompExp(catID)->theData,
3872 * matrix2 = ((_CalcNode*)(LocateVar(theRoot->get_node(3)->in_object)))->GetCompExp(catID)->theData,
3873 overallResult = 0.;
3874
3875 long patternCount = dsf->GetPatternCount();
3876
3877 hyFloat currentAccumulator = 1.;
3878
3879 for (long patternIndex = 0; patternIndex < patternCount; patternIndex ++, l0+=4, l1+=4, l2+=4) {
3880 hyFloat rp0 = l0[0] * matrix0[0]+ l0[1] * matrix0[1] + l0[2] * matrix0[2] + l0[3] * matrix0[3];
3881 hyFloat rp1 = l0[0] * matrix0[4]+ l0[1] * matrix0[5] + l0[2] * matrix0[6] + l0[3] * matrix0[7];
3882 hyFloat rp2 = l0[0] * matrix0[8]+ l0[1] * matrix0[9] + l0[2] * matrix0[10] + l0[3] * matrix0[11];
3883 hyFloat rp3 = l0[0] * matrix0[12]+ l0[1] * matrix0[13] + l0[2] * matrix0[14] + l0[3] * matrix0[15];
3884
3885 rp0 *= l1[0] * matrix1[0] + l1[1] * matrix1[1] + l1[2] * matrix1[2] + l1[3] * matrix1[3];
3886 rp1 *= l1[0] * matrix1[4] + l1[1] * matrix1[5] + l1[2] * matrix1[6] + l1[3] * matrix1[7];
3887 rp2 *= l1[0] * matrix1[8] + l1[1] * matrix1[9] + l1[2] * matrix1[10] + l1[3] * matrix1[11];
3888 rp3 *= l1[0] * matrix1[12]+ l1[1] * matrix1[13] + l1[2] * matrix1[14] + l1[3] * matrix1[15];
3889
3890 rp0 *= l2[0] * matrix2[0] + l2[1] * matrix2[1] + l2[2] * matrix2[2] + l2[3] * matrix2[3];
3891 rp1 *= l2[0] * matrix2[4] + l2[1] * matrix2[5] + l2[2] * matrix2[6] + l2[3] * matrix2[7];
3892 rp2 *= l2[0] * matrix2[8] + l2[1] * matrix2[9] + l2[2] * matrix2[10] + l2[3] * matrix2[11];
3893 rp3 *= l2[0] * matrix2[12]+ l2[1] * matrix2[13] + l2[2] * matrix2[14] + l2[3] * matrix2[15];
3894
3895 hyFloat result = theProbs[0]*rp0+
3896 theProbs[1]*rp1+
3897 theProbs[2]*rp2+
3898 theProbs[3]*rp3;
3899
3900
3901 if (result<=0.0) {
3902 return -INFINITY;
3903 }
3904
3905 long patternFreq = dsf->theFrequencies[patternIndex];
3906 for (long freqIterator = 0; freqIterator < patternFreq; freqIterator++) {
3907 hyFloat tryMultiplication = currentAccumulator*result;
3908 if (tryMultiplication > 1.e-300) {
3909 currentAccumulator = tryMultiplication;
3910 } else {
3911 overallResult += myLog (currentAccumulator);
3912 currentAccumulator = result;
3913 }
3914 }
3915 }
3916 return overallResult + myLog (currentAccumulator);
3917 }
3918
3919 //_______________________________________________________________________________________________
3920
_RemoveNodeList(_SimpleList const & clean_indices)3921 void _TheTree::_RemoveNodeList (_SimpleList const& clean_indices) {
3922 if (compExp) {
3923 DeleteAndZeroObject(compExp);
3924 }
3925 clean_indices.Each([] (long var_idx, unsigned long) -> void {
3926 DeleteVariable(var_idx, true);
3927 });
3928 }
3929