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 <float.h>
42 #include "defines.h"
43 #include "formula.h"
44
45
46 #include "parser.h"
47 #include "batchlan.h"
48 #include "function_templates.h"
49 #include "global_things.h"
50 #include "global_object_lists.h"
51 #include "polynoml.h"
52
53 using namespace hy_global;
54 using namespace hyphy_global_objects;
55
56 //Constants
57 hyFloat const sqrtPi = 1.77245385090551603,
58 twoOverSqrtPi = 2./sqrtPi;
59
60 _Formula * current_formula_being_computed = nil;
61
62 //__________________________________________________________________________________
63
64
_Formula(void)65 _Formula::_Formula (void) {
66 Initialize();
67 }
68
69 //__________________________________________________________________________________
70
_Formula(HBLObjectRef p,bool is_a_var)71 _Formula::_Formula (HBLObjectRef p, bool is_a_var) {
72 theTree = nil;
73 resultCache = nil;
74 recursion_calls = nil;
75 call_count = 0UL;
76 simpleExpressionStatus = nil;
77
78 if (!is_a_var) {
79 theFormula.AppendNewInstance (new _Operation (p));
80 } else {
81 _Variable* v = (_Variable*)p;
82 theFormula.AppendNewInstance (new _Operation (true,*v->GetName(),v->IsGlobal(), nil));
83 }
84 }
85
86 //__________________________________________________________________________________
_Formula(_String const & s,_VariableContainer const * theParent,_String * reportErrors)87 _Formula::_Formula (_String const &s, _VariableContainer const* theParent, _String* reportErrors) {
88 simpleExpressionStatus = nil;
89 ParseFormula (s, theParent, reportErrors);
90 }
91
92 //__________________________________________________________________________________
93
_Formula(const _Polynomial * source)94 _Formula::_Formula(const _Polynomial* source) {
95 Initialize();
96 _PolynomialData * poly_terms = source->GetTheTerms();
97
98 _List poly_vars;
99 for (long i = 0L; i < source->GetNoVariables(); i++) {
100 poly_vars << source->GetIthVariable(i);
101 }
102
103 for (long i = 0L; i < poly_terms->NumberOfTerms(); i++) {
104 hyFloat c = poly_terms->GetCoeff(i);
105 /*theFormula.AppendNewInstance (new _Operation (new _Constant (poly_terms->GetCoeff(i))));*/
106 bool started_term = false;
107 if (!CheckEqual(c, 1.0)) {
108 theFormula.AppendNewInstance (new _Operation (new _Constant (poly_terms->GetCoeff(i))));
109 started_term = true;
110 }
111 long *exponents = poly_terms->GetTerm(i);
112 poly_vars.ForEach ([this, exponents, &started_term] (BaseRefConst v, unsigned long idx) -> void {
113 if (exponents[idx] != 0L) {
114 this->theFormula.AppendNewInstance (new _Operation (*((_Variable const*)v)));
115 if (exponents[idx] != 1L) {
116 this->theFormula.AppendNewInstance (new _Operation (new _Constant (exponents[idx])));
117 this->theFormula.AppendNewInstance (new _Operation (HY_OP_CODE_POWER, 2));
118 }
119 if (started_term) {
120 this->theFormula.AppendNewInstance (new _Operation (HY_OP_CODE_MUL, 2));
121 } else {
122 started_term = true;
123 }
124 }
125 });
126 if (!started_term) {
127 this->theFormula.AppendNewInstance (new _Operation (new _Constant (1.0)));
128 }
129 if (i) {
130 this->theFormula.AppendNewInstance (new _Operation (HY_OP_CODE_ADD, 2));
131 }
132 }
133 }
134
135 //__________________________________________________________________________________
Initialize(void)136 void _Formula::Initialize (void) {
137 theTree = nil;
138 resultCache = nil;
139 call_count = 0UL;
140 recursion_calls = nil;
141 simpleExpressionStatus = nil;
142 }
143
144 //__________________________________________________________________________________
_Formula(_Formula const & rhs)145 _Formula::_Formula (_Formula const& rhs) {
146 *this = rhs;
147 }
148
149 //__________________________________________________________________________________
operator =(_Formula const & rhs)150 const _Formula& _Formula::operator = (_Formula const& rhs) {
151 if (this != &rhs) {
152 Initialize();
153 Duplicate (&rhs);
154 }
155 return *this;
156 }
157
158 //__________________________________________________________________________________
Duplicate(_Formula const * f_cast)159 void _Formula::Duplicate (_Formula const * f_cast) {
160
161 theFormula.Duplicate (& f_cast->theFormula);
162 theStack.theStack.Duplicate(& f_cast->theStack.theStack);
163 call_count = f_cast->call_count;
164 if (f_cast->recursion_calls) {
165 recursion_calls = (HBLObjectRef)f_cast->recursion_calls->makeDynamic();
166 } else {
167 recursion_calls = nil;
168 }
169 if (f_cast->theTree) {
170 theTree = f_cast->theTree->duplicate_tree();
171 } else {
172 theTree = nil;
173 }
174 simpleExpressionStatus = nil;
175 if (f_cast->resultCache) {
176 resultCache = (_List*)f_cast->resultCache->makeDynamic();
177 } else {
178 resultCache = nil;
179 }
180 }
181
182 //__________________________________________________________________________________
DuplicateReference(const _Formula * f)183 void _Formula::DuplicateReference (const _Formula* f) {
184 for (unsigned long i=0; i<f->theFormula.countitems(); i++) {
185 _Operation *ith_term = f->GetIthTerm(i);
186 // -2 is used to code for value substitutions (x__)
187 if (ith_term->IsValueSubstitution()) {
188 theFormula.AppendNewInstance(new _Operation ((HBLObjectRef)LocateVar (ith_term->GetAVariable())->Compute()->makeDynamic()));
189 } else {
190 theFormula && ith_term;
191 }
192 }
193 }
194
195 //__________________________________________________________________________________
196
ParseAndCompute(_String const & expression,bool use_exceptions,long type,_hyExecutionContext * context)197 HBLObjectRef _Formula::ParseAndCompute (_String const& expression, bool use_exceptions, long type, _hyExecutionContext * context) {
198
199 _String error_message;
200 _Formula f;
201
202 long parse_result = f.ParseFormula (expression, context ? context -> GetContext() : nil , &error_message);
203
204 try {
205 if (error_message.nonempty()) {
206 throw (_String ("Failed to parse ") & expression.Enquote () & " with the following error: " & error_message);
207 }
208 if (parse_result != HY_FORMULA_EXPRESSION) {
209 throw (expression.Enquote () & " did not parse to a simple expression");
210 }
211 if (f.IsEmpty ()) {
212 throw (expression.Enquote () & " parsed to an empty expression");
213 }
214 if (!(type == HY_ANY_OBJECT || f.ObjectClass() == type)) {
215 // TODO SLKP 20170704: ObjectClass will compute the expression with current values which may fail
216 throw (expression.Enquote () & " did not evaluate to a " & FetchObjectNameFromType (type));
217 }
218 } catch (_String const & e) {
219 if (use_exceptions) {
220 throw (e);
221 } else {
222 HandleApplicationError (e);
223 return nil;
224 }
225 }
226
227 HBLObjectRef result = f.Compute();
228 result->AddAReference();
229 return result;
230 }
231
232
233 //__________________________________________________________________________________
makeDynamic(void) const234 BaseRef _Formula::makeDynamic (void) const {
235 _Formula * res = new _Formula;
236 res->Duplicate(this);
237 return (BaseRef)res;
238 }
239 //__________________________________________________________________________________
240
~_Formula(void)241 _Formula::~_Formula (void) {
242 Clear();
243 }
244
245 //__________________________________________________________________________________
Clear(void)246 void _Formula::Clear (void) {
247 if (theTree) {
248 theTree->delete_tree();
249 delete theTree;
250 }
251 theTree = nil;
252 if (resultCache) {
253 DeleteAndZeroObject(resultCache);
254 }
255
256 if (simpleExpressionStatus){
257 delete [] simpleExpressionStatus;
258 simpleExpressionStatus = nil;
259 }
260
261 theFormula.Clear();
262 if (recursion_calls) {
263 delete (recursion_calls);
264 }
265
266 // theStack.Clear();
267 }
268 //__________________________________________________________________________________
toRPN(_hyFormulaStringConversionMode mode,_List * matched_names)269 _StringBuffer const _Formula::toRPN (_hyFormulaStringConversionMode mode, _List* matched_names) {
270 _StringBuffer r;
271 if (theFormula.countitems()) {
272 SubtreeToString (r,nil,0,nil,GetIthTerm(0UL), mode);
273 for (unsigned long k=1UL; k<theFormula.countitems(); k++) {
274 r <<'|';
275 SubtreeToString (r,nil,0,nil,GetIthTerm(k), mode);
276 }
277 }
278 return r;
279 }
280 //__________________________________________________________________________________
toStr(_hyFormulaStringConversionMode mode,_List * matched_names,bool drop_tree)281 BaseRef _Formula::toStr (_hyFormulaStringConversionMode mode, _List* matched_names, bool drop_tree) {
282 ConvertToTree(false);
283
284 _StringBuffer * result = new _StringBuffer (64UL);
285
286 long savepd = print_digit_specification;
287 print_digit_specification = 0L;
288
289 if (theTree) { // there is something to do
290 SubtreeToString (*result, theTree, -1, matched_names, nil, mode);
291 } else {
292 if (theFormula.countitems()) {
293 (*result) << "RPN:" << toRPN (mode, matched_names);
294 }
295 }
296
297 print_digit_specification = savepd;
298 result->TrimSpace();
299 if (theTree && drop_tree) {
300 theTree->delete_tree();
301 delete theTree;
302 theTree = nil;
303 }
304
305 return result;
306 }
307 //__________________________________________________________________________________
DuplicateFormula(node<long> * src,_Formula & tgt) const308 node<long>* _Formula::DuplicateFormula (node<long>* src, _Formula& tgt) const {
309
310 node<long>* copied = new node<long>;
311
312 for (long k=1L; k<=src->get_num_nodes(); k++) {
313 copied->add_node (*DuplicateFormula (src->go_down (k), tgt));
314 }
315
316 copied->in_object = tgt.theFormula.lLength;
317 tgt.theFormula && theFormula.GetItem(src->in_object);
318
319 return copied;
320 }
321
322 //__________________________________________________________________________________
Differentiate(_String const & var_name,bool bail,bool convert_from_tree)323 _Formula* _Formula::Differentiate (_String const & var_name, bool bail, bool convert_from_tree) {
324
325 _Variable * dx = FetchVar(LocateVarByName(var_name));
326
327 if (!dx) { // create the dX variable if it's not here ye
328 return new _Formula (new _Constant (0.0));
329 }
330
331 long dx_id = dx->get_index();
332
333 ConvertToTree ();
334
335 //printf ("\n **** Diff %s on %s\n\n", _String ((_String*)toStr(kFormulaStringConversionNormal)).get_str(), var_name.get_str());
336
337 _SimpleList var_refs = PopulateAndSort ([&] (_AVLList & parameter_list) -> void {
338 this->ScanFForVariables (parameter_list, true, true, true);
339 });
340
341 if (var_refs.empty()) { // handle the case of constant expressions here
342 return new _Formula (new _Constant (0.0));
343 }
344
345 _Formula* res = new _Formula ();
346 _Formula ** dydx = new _Formula* [var_refs.countitems()] {0};// stores precomputed derivatives for all the
347
348 auto dydx_cleanup = [&] () -> void {
349 for (unsigned long k=0UL; k < var_refs.countitems(); k++) {
350 if (dydx[k]) {
351 delete dydx[k];
352 }
353 }
354 delete [] dydx;
355 };
356
357 node<long>* dTree = nil;
358
359 try {
360 for (unsigned long k=0UL; k < var_refs.countitems(); k++) {
361 _Variable* thisVar = LocateVar (var_refs.GetElement(k));
362 _Formula * dYdX;
363
364 if (thisVar->IsIndependent()) {
365 dYdX = new _Formula ((*thisVar->GetName() == var_name)?new _Constant (1.0):new _Constant (0.0));
366 } else {
367 dYdX = thisVar->varFormula->Differentiate (var_name, bail, false);
368 if (dYdX->IsEmpty()) {
369 delete dYdX;
370 dydx_cleanup ();
371 return res;
372 }
373 }
374 dYdX->ConvertToTree();
375 dydx [k] = dYdX;
376 }
377
378 // SortLists (&varRefs, &dydx);
379 // this is already sorted coming from PopulateAndSort
380
381
382 if (!(dTree = InternalDifferentiate (theTree, dx_id, var_refs, dydx, *res))) {
383 throw (_String ("Differentiation of ") & _String((_String*)toStr(kFormulaStringConversionNormal)) & " failed.");
384 }
385 } catch (_String const &e) {
386 dydx_cleanup ();
387
388 if (bail) {
389 HandleApplicationError (e);
390 res->Clear();
391 return res;
392 } else {
393 delete res;
394 return nil;
395 }
396 }
397
398 dydx_cleanup ();
399
400 //res->theFormula.AppendNewInstance (new _Operation(new _Constant (0.0))) ; // why is this here?
401 res->theTree = dTree;
402 res->InternalSimplify (dTree);
403 //printf ("%s\n", _String ((_String*)res->theFormula.toStr(kFormulaStringConversionNormal)).get_str());
404
405 if (convert_from_tree) {
406 res->ConvertFromTree ();
407 }
408 //printf ("%s\n", _String ((_String*)res->theFormula.toStr(kFormulaStringConversionNormal)).get_str());
409
410 // try polynomial simplification
411 _Polynomial* is_poly = (_Polynomial*)res->ConstructPolynomial();
412 if (is_poly) {
413 _Formula * simplified_polynomial = new _Formula (is_poly);
414 //printf ("%s\n", _String ((_String*)simplified_polynomial->toStr(kFormulaStringConversionNormal)).get_str());
415 delete res;
416
417 //printf ("\n RESULT : %s \n", _String ((_String*)simplified_polynomial->toStr(kFormulaStringConversionNormal)).get_str());
418
419 return simplified_polynomial;
420 } else {
421 //printf ("%s\n", _String ((_String*)res->toStr(kFormulaStringConversionNormal)).get_str());
422 }
423
424 return res;
425
426 }
427
428 //__________________________________________________________________________________
InternalDifferentiate(node<long> * currentSubExpression,long varID,_SimpleList const & varRefs,_Formula * const * dydx,_Formula & tgt)429 node<long>* _Formula::InternalDifferentiate (node<long>* currentSubExpression, long varID, _SimpleList const & varRefs, _Formula * const * dydx, _Formula& tgt) {
430 _Operation * op = GetIthTerm(currentSubExpression->in_object);
431
432 if (op->theData!=-1) {
433 long k = varRefs.BinaryFind (op->GetAVariable());
434 if (k<0L) {
435 return nil;
436 }
437
438 _Formula const* dYdX = dydx[k];
439 return dYdX->DuplicateFormula (dYdX->theTree, tgt);
440 }
441
442 if (op->theNumber) {
443 _Formula src (new _Constant (0.0));
444 src.ConvertToTree ();
445 return src.DuplicateFormula (src.theTree, tgt);
446 }
447
448 node<long>* newNode = new node<long>;
449 node <long> * created_nodes [7]{nil};
450 created_nodes [0] = newNode;
451
452 try {
453 switch (op->opCode) {
454 case HY_OP_CODE_MUL: {
455 /**
456 d (X*Y) = X*dY + Y*dX
457 */
458
459 created_nodes [1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
460
461 if (!created_nodes [1]) {
462 throw (2);
463 }
464
465 node<long>* YtimesDX = new node<long>;
466 created_nodes[3] = YtimesDX;
467 YtimesDX->add_node (*created_nodes [1],*DuplicateFormula (currentSubExpression->go_down(2),tgt));
468 YtimesDX->in_object = tgt.theFormula.AppendNewInstance (new _Operation (HY_OP_CODE_MUL,2)) - 1;
469
470
471 created_nodes [2] = InternalDifferentiate (currentSubExpression->go_down(2), varID, varRefs, dydx, tgt);
472 if (!created_nodes [2]) {
473 throw (3);
474 }
475 node<long>* XtimesDY = new node<long>;
476 XtimesDY->add_node (*created_nodes [2],*DuplicateFormula (currentSubExpression->go_down(1),tgt));
477 XtimesDY->in_object = tgt.theFormula.AppendNewInstance (new _Operation (HY_OP_CODE_MUL,2)) - 1;
478 newNode->add_node (*YtimesDX,*XtimesDY);
479 newNode->in_object = tgt.theFormula. AppendNewInstance (new _Operation (HY_OP_CODE_ADD,2)) - 1;
480
481 return newNode;
482 }
483 break;
484
485 case HY_OP_CODE_ADD: // +
486 case HY_OP_CODE_SUB: { // -
487
488 // TODO SLKP 20170426: clean up the rest of the function
489
490 created_nodes[1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
491
492 if (!created_nodes[1]) {
493 throw (1);
494 }
495
496 long isUnary = (currentSubExpression->get_num_nodes()==1);
497
498 if (!isUnary) {
499 created_nodes[2] = InternalDifferentiate (currentSubExpression->go_down(2), varID, varRefs, dydx, tgt);
500 if (!created_nodes[2]) {
501 throw (2);
502 }
503 }
504
505
506 newNode->add_node (*created_nodes[1]);
507 if (!isUnary) {
508 newNode->add_node (*created_nodes[2]);
509 }
510 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (*op)) - 1;
511 return newNode;
512 }
513 break;
514
515 case HY_OP_CODE_DIV: { // /
516
517 /*
518 d (X/Y) = (y * dX - x * dY ) * Y^{-2}
519 */
520 created_nodes [1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
521
522 if (!created_nodes [1]) {
523 throw (2);
524 }
525
526
527 node<long>* yDx = new node<long>;
528 created_nodes[3] = yDx;
529 yDx->add_node (*created_nodes [1],*DuplicateFormula (currentSubExpression->go_down(2),tgt)); // dX * Y
530 yDx->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL,2))-1;
531
532 created_nodes [2] = InternalDifferentiate (currentSubExpression->go_down(2), varID, varRefs, dydx, tgt);
533 if (!created_nodes [2]) {
534 throw (4);
535 }
536
537 node<long>* y_raise_m2 = new node<long>;
538 node<long>* yDx_minus_xDy = new node<long>;
539 node<long>* xDy = new node<long>;
540 node<long>* m2 = new node<long>;
541
542 xDy->add_node (*created_nodes [2],*DuplicateFormula (currentSubExpression->go_down(1),tgt)); // dY * X
543 xDy->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL,2))-1;
544
545 yDx_minus_xDy->add_node (*yDx,*xDy);
546 yDx_minus_xDy->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_SUB,2))-1; // numerator
547
548 y_raise_m2->add_node (*DuplicateFormula (currentSubExpression->go_down(2),tgt), *m2);
549 m2->in_object = tgt.theFormula.AppendNewInstance(new _Operation (new _Constant (-2.)))-1;
550 y_raise_m2->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_POWER,2))-1; // y^{-2}
551
552 newNode->add_node(*yDx_minus_xDy, *y_raise_m2);
553 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL,2)) - 1;
554
555 return newNode;
556 }
557 break;
558
559 case HY_OP_CODE_ARCTAN: { // Arctan
560
561 // Arctax (y)' = y' / (1+y^2)
562
563 created_nodes [1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
564
565 if (!created_nodes[1]) {
566 throw (1);
567 }
568
569 node<long>* one = new node<long>,
570 * two = new node<long>,
571 * denom = new node<long>,
572 * y_squared = new node<long>;
573
574 y_squared->add_node(*DuplicateFormula(currentSubExpression->go_down(1), tgt), *two);
575 two->in_object = tgt.theFormula.AppendNewInstance(new _Operation (new _Constant (2.))) - 1;
576 y_squared->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_POWER,2)) - 1;
577 denom->add_node(*y_squared, *one);
578 one->in_object = tgt.theFormula.AppendNewInstance(new _Operation (new _Constant (1.))) - 1;
579 denom->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_ADD,2)) - 1;
580 newNode->add_node(*created_nodes[1], *denom);
581 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_DIV,2)) - 1;
582
583 return newNode;
584 }
585 break;
586
587 case HY_OP_CODE_COS: {
588 // d Cos (f[x]) = - f'[x] * Sin (x)
589
590 created_nodes[1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
591
592 if (! created_nodes[1] ) {
593 throw (1);
594 }
595
596 node <long> * func_fx = new node <long>;
597 func_fx->add_node(*DuplicateFormula(currentSubExpression->go_down(1), tgt));
598 func_fx->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_SIN,1)) - 1;
599 node <long> * product_term = new node <long>;
600 product_term->add_node (*created_nodes[1], *func_fx);
601 product_term->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL,2)) - 1;
602 newNode->add_node(*product_term);
603 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_SUB,1)) - 1;
604
605 return newNode;
606 }
607 break;
608
609 case HY_OP_CODE_ERF: { // Erf
610
611 // d ERF (f [X]) = Exp(-f[X]^2) / 2Pi
612
613 created_nodes[1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
614
615 if (!created_nodes[1]) {
616 throw (1);
617 }
618
619 node <long> * f_squared = new node<long>,
620 * two = new node<long>;
621
622 f_squared->add_node(*DuplicateFormula(currentSubExpression->go_down(1), tgt), *two);
623 two->in_object = tgt.theFormula.AppendNewInstance(new _Operation (new _Constant (2.))) - 1L;
624 f_squared->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_POWER, 2)) - 1L;
625
626 node <long> * minus_f_squared = new node<long>;
627 minus_f_squared->add_node(*f_squared);
628 minus_f_squared->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_SUB, 1)) - 1L;
629
630 node <long> * exp_minus_f_squared = new node<long>;
631 exp_minus_f_squared->add_node(*minus_f_squared);
632 exp_minus_f_squared->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_EXP, 1)) - 1L;
633
634
635 node <long> * two_pi = new node<long>,
636 * div_by_two_pi = new node<long>;
637
638 div_by_two_pi->add_node(*exp_minus_f_squared, *two_pi);
639 two_pi->in_object = tgt.theFormula.AppendNewInstance(new _Operation (new _Constant (twoOverSqrtPi))) - 1L;
640 div_by_two_pi ->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_DIV, 2)) - 1L;
641
642 newNode->add_node(*created_nodes[1], *div_by_two_pi);
643 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL, 2)) - 1L;
644
645 return newNode;
646 }
647 break;
648
649 case HY_OP_CODE_EXP: // HY_OP_CODE_EXP
650 case HY_OP_CODE_SIN: { // HY_OP_CODE_SIN
651
652 // d Exp (f[X]) = f'[X] Exp (f[X])
653 // d Sin (f[X]) = f'[X] Cos (f[X])
654
655 created_nodes[1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
656
657 if (! created_nodes[1] ) {
658 throw (1);
659 }
660
661 node <long> * func_fx = new node <long>;
662 func_fx->add_node(*DuplicateFormula(currentSubExpression->go_down(1), tgt));
663 func_fx->in_object = tgt.theFormula.AppendNewInstance(new _Operation (op->opCode==HY_OP_CODE_SIN ? HY_OP_CODE_COS : HY_OP_CODE_EXP,1)) - 1;
664 newNode->add_node (*created_nodes[1], *func_fx);
665 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL,2)) - 1;
666
667 return newNode;
668 }
669 break;
670
671 case HY_OP_CODE_LOG: {
672 // Log (f(x))' = f'(x) * f(x)^{-1}
673 created_nodes[1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
674
675 if (!created_nodes[1]) {
676 throw (2);
677 }
678
679 node <long> * n1 = new node<long>, * n2 = new node<long>;
680 n2->add_node (*DuplicateFormula (currentSubExpression->go_down(1),tgt), *n1); // f(x) , {-1}
681 n1->in_object = tgt.theFormula.AppendNewInstance (new _Operation (new _Constant (-1))) - 1; // push -1
682 n2->in_object = tgt.theFormula.AppendNewInstance (new _Operation (HY_OP_CODE_POWER,2)) - 1; // push '^'
683 newNode->add_node (*created_nodes[1], *n2); // f'(x) , f(x) (^) {-1}
684 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL,2))-1; // push *
685
686 return newNode;
687 }
688 break;
689
690 case HY_OP_CODE_SQRT: {
691
692 // d (Sqrt (f[x])) = 0.5 * f'[x] / Sqrt (f[x]))
693 created_nodes[1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
694
695 if (!created_nodes[1]) {
696 throw (1);
697 }
698
699 node <long>* half = new node<long>;
700 node <long>* half_times_dfx = new node<long>;
701
702 half_times_dfx->add_node (*created_nodes[1],*half);
703 half->in_object = tgt.theFormula.AppendNewInstance(new _Operation (new _Constant (0.5))) - 1;
704 half_times_dfx->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL, 2)) - 1;
705
706 node <long> * sqrt_fx = new node<long>;
707 sqrt_fx->add_node (*DuplicateFormula (currentSubExpression->go_down(1),tgt));
708 sqrt_fx->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_SQRT, 1)) - 1;
709
710 newNode->add_node (*half_times_dfx,*sqrt_fx);
711 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_DIV, 2)) - 1;
712 return newNode;
713 }
714 break;
715
716 case HY_OP_CODE_TAN: {
717 // Tan (f[X]) = f'[X] / Cos ^ 2 (f[X])
718
719 created_nodes[1] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
720
721 if (!created_nodes[1]) {
722 throw (1);
723 }
724
725 node <long>* two = new node<long>;
726 node <long>* cos_f = new node<long>;
727 node <long>* cos_f_squared = new node<long>;
728
729 cos_f->add_node(*DuplicateFormula(currentSubExpression->go_down(1), tgt));
730 cos_f->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_COS, 1)) - 1;
731 cos_f_squared->add_node(*cos_f, *two);
732 two->in_object = tgt.theFormula.AppendNewInstance(new _Operation (new _Constant (2.))) - 1;
733 cos_f_squared->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_POWER, 2)) - 1;
734
735 newNode->add_node(*created_nodes[1], *cos_f_squared);
736 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_DIV, 2)) - 1;
737
738 return newNode;
739 }
740 break;
741
742 case HY_OP_CODE_POWER: // ^
743 // f[x]^g[x] (g'[x] Log[f[x]] + f'[x]g[x] * (f[x]^{-1}))
744 {
745 node <long> * f_raise_g = new node <long>;
746 f_raise_g->add_node (*DuplicateFormula (currentSubExpression->go_down(1),tgt),*DuplicateFormula (currentSubExpression->go_down(2),tgt));
747 f_raise_g->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_POWER, 2)) - 1;
748 created_nodes[1] = f_raise_g;
749 created_nodes[2] = InternalDifferentiate (currentSubExpression->go_down(2), varID, varRefs, dydx, tgt);
750 if (!created_nodes[2]) {
751 throw (3);
752 }
753 node <long> * log_f = new node <long>;
754 created_nodes[3] = log_f;
755 log_f->add_node (*DuplicateFormula (currentSubExpression->go_down(1),tgt));
756 log_f->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_LOG, 1)) - 1;
757 node <long> * dg_times_log_fx = new node <long>;
758 created_nodes[4] = dg_times_log_fx;
759 dg_times_log_fx->add_node (*created_nodes[2], *log_f);
760 dg_times_log_fx->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL, 2)) - 1;
761 created_nodes[5] = InternalDifferentiate (currentSubExpression->go_down(1), varID, varRefs, dydx, tgt);
762 if (!created_nodes[5]) {
763 throw (6);
764 }
765 node <long> * df_times_g = new node <long>;
766 df_times_g->add_node (*created_nodes[5],*DuplicateFormula (currentSubExpression->go_down(2),tgt));
767 df_times_g->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL, 2)) - 1;
768 node <long> * inverse_f = new node <long>;
769 node <long> * minus1 = new node<long>;
770 inverse_f->add_node (*DuplicateFormula (currentSubExpression->go_down(1),tgt),*minus1);
771 minus1->in_object = tgt.theFormula.AppendNewInstance(new _Operation (new _Constant (-1.))) - 1;
772 inverse_f->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_POWER, 2)) - 1;
773 node <long> * product_term = new node<long>; // f'[x]g[x] [*] (f[x]^{-1}
774 product_term->add_node (*df_times_g, *inverse_f);
775 product_term->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL, 2)) - 1;
776 node <long> * sum_term = new node<long>; // g'[x] Log[f[x]] <+> f'[x]g[x] * (f[x]^{-1})
777 sum_term->add_node (*dg_times_log_fx, *product_term);
778 sum_term->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_ADD, 2)) - 1;
779 newNode->add_node (*f_raise_g, *sum_term);
780 newNode->in_object = tgt.theFormula.AppendNewInstance(new _Operation (HY_OP_CODE_MUL, 2)) - 1;
781
782
783 return newNode;
784 }
785 }
786 } catch (int cleanup_length) {
787 for (long i = 0L; i < cleanup_length; i++) {
788 if (created_nodes [i]) {
789 created_nodes[i]->delete_tree(true);
790 }
791 }
792 return nil;
793 }
794 delete (newNode);
795 return nil;
796 }
797
798
799
800 //__________________________________________________________________________________
InternalSimplify(node<long> * top_node)801 bool _Formula::InternalSimplify (node<long>* top_node) {
802 // returns true if the subexpression at
803 // and below startnode is constant
804 _Operation* op = GetIthTerm(top_node->get_data());
805 long n_children = top_node->get_num_nodes();
806
807 if (n_children == 0L) {
808 return !op->IsAVariable();
809 }
810
811
812 bool all_constant = true,
813 left_constant = true,
814 right_constant = (n_children>1L);
815
816 long prune_this_child = -1;
817
818 hyFloat evaluated_to = 0.0;
819 HBLObjectRef replace_with = nil;
820
821
822 //printf ("InternalSimplify %x\n", startNode);
823
824 for (unsigned long k=1UL; k<=n_children; k++) {
825 if (k==1UL) {
826 left_constant = InternalSimplify (top_node->go_down(k));
827 } else if (k==2UL) {
828 right_constant = InternalSimplify (top_node->go_down(k));
829 } else {
830 if (!InternalSimplify (top_node->go_down(k))) {
831 all_constant = false;
832 }
833 }
834 }
835
836 all_constant = all_constant && left_constant && (n_children==1 || right_constant);
837
838 if (op->opCode > HY_OP_CODE_NONE) {
839 if (all_constant) { // this executes the subxpression starting at the current node
840 _Stack scrap;
841 for (unsigned long k=1UL; k<=n_children; k++) {
842 ((_Operation*)theFormula (top_node->go_down(k)->get_data()))->Execute (scrap);
843 }
844 op->Execute (scrap);
845 replace_with = (HBLObjectRef)scrap.Pop();//->makeDynamic();
846 } else {
847 if (left_constant||right_constant) {
848
849 HBLObjectRef constant_value = ((_Operation*)theFormula (top_node->go_down(left_constant?1:2)->get_data()))->GetANumber();
850
851
852 if (constant_value->ObjectClass() == NUMBER) {
853 evaluated_to = constant_value->Value();
854
855 switch (op->opCode) {
856 case HY_OP_CODE_MUL: { // *
857 if (CheckEqual (evaluated_to,0.0)) { // *0 => 0
858 //printf ("*0\n");
859 replace_with = new _Constant (0.0);
860 break;
861 }
862 if (CheckEqual (evaluated_to,1.0)) { // x*1 => x
863 //printf ("*1\n");
864 prune_this_child = left_constant?1:2;
865 break;
866 }
867 }
868 break;
869
870 case HY_OP_CODE_ADD: { // +
871 if (CheckEqual (evaluated_to,0.0)) { // x+0 => x
872 //printf ("+0\n");
873 prune_this_child = left_constant?1:2;
874 }
875 break;
876 }
877
878 case HY_OP_CODE_SUB: { // x-0 => x
879 // 0-x => -x
880
881 if (CheckEqual (evaluated_to,0.0)) {
882 //printf ("-0\n");
883 prune_this_child = left_constant? -2 : 2;
884 }
885 break;
886 }
887
888 case HY_OP_CODE_DIV: { // /
889 if (left_constant&&CheckEqual (evaluated_to,0.0)) { // 0/x => 0
890 replace_with = new _Constant (0.0);
891 //printf ("0/\n");
892
893 break;
894 }
895 if (right_constant&&CheckEqual (evaluated_to,1.0)) { // x/1 => x
896 //printf ("/1\n");
897 prune_this_child = 2;
898 break;
899 }
900 }
901 break;
902
903 case HY_OP_CODE_POWER: { // ^
904 if (left_constant&&CheckEqual (evaluated_to,1.0)) { // 1^x => 1
905 //printf ("1^\n");
906 replace_with = new _Constant (1.0);
907 break;
908 }
909 if (right_constant&&CheckEqual (evaluated_to,1.0)) { // x^1 => ?
910 //printf ("^1\n");
911 prune_this_child = 1;
912 break;
913 }
914 }
915 break;
916 }
917 }
918 }
919 }
920 }
921
922 if (replace_with) {
923 all_constant = true;
924 for (int k=1; k <= n_children; k++) {
925 top_node->go_down(k)->delete_tree(true);
926 }
927 top_node->kill_all_nodes();
928 top_node->in_object = theFormula.lLength;
929 theFormula < (new _Operation(replace_with));
930 } else {
931
932 if (prune_this_child !=- 1L) {
933 if (prune_this_child > 0L) {
934 top_node->go_down(prune_this_child)->delete_tree(true);
935 top_node->kill_node (prune_this_child);
936 node <long>* pruned_tree = top_node->go_down(1);
937
938 top_node->kill_all_nodes();
939
940 for (unsigned long k=1; k<=pruned_tree->get_num_nodes(); k++) {
941 top_node->add_node (*pruned_tree->go_down(k));
942 }
943 top_node->in_object = pruned_tree->in_object;
944 delete (pruned_tree);
945
946 } else { // 0-? => -?
947 top_node->go_down(1)->delete_tree(true);
948 top_node->kill_node(1);
949 op->SetTerms(1);
950 //startNode->kill_node(1);
951 }
952 }
953 }
954
955 // check to see if this op may be of the kind F (F^-1 (X)), like Log (Exp (..)))
956
957 if (!all_constant && top_node->get_num_nodes() == 1) {
958 long check_op_code = GetIthTerm(top_node->get_data())->opCode;
959 node <long>* child_node = top_node->go_down(1);
960 if (child_node->get_num_nodes() == 1) {
961 long child_op_code = GetIthTerm(child_node->get_data())->opCode;
962 bool cancel_pair = _Operation::AreOpsInverse (child_op_code, check_op_code);
963 if (cancel_pair) {
964 node<long> * remaining_trunk = child_node->go_down(1);
965 delete (top_node->go_down(1));
966 top_node->in_object = remaining_trunk->in_object;
967 top_node->kill_all_nodes();
968 for (long k = 1; k <= remaining_trunk->get_num_nodes(); k++) {
969 top_node->add_node (*remaining_trunk->go_down (k));
970 }
971 delete remaining_trunk;
972
973 }
974 }
975 }
976
977 return all_constant;
978 }
979
980
981 //__________________________________________________________________________________
SubtreeToString(_StringBuffer & result,node<long> * top_node,int op_level,_List * match_names,_Operation * this_node_op,_hyFormulaStringConversionMode mode)982 void _Formula::SubtreeToString (_StringBuffer & result, node<long>* top_node, int op_level, _List* match_names, _Operation* this_node_op, _hyFormulaStringConversionMode mode) {
983
984 if (!this_node_op) {
985 if (!top_node) {
986 return;
987 }
988 this_node_op = GetIthTerm (top_node->get_data());
989 }
990
991 // decide what to do about this operation
992
993 if (this_node_op->IsAVariable(false) || this_node_op->IsValueSubstitution()) {
994 // this operation is just a variable - add ident to string and return
995 long node_variable_index = this_node_op->GetAVariable();
996 if (node_variable_index < 0L) {
997 return;
998 }
999
1000 if (mode != kFormulaStringConversionNormal) {
1001
1002 _Variable *node_variable = LocateVar(node_variable_index);
1003
1004 if (mode == kFormulaStringConversionSubstiteValues) {
1005 if (hy_x_variable && (node_variable->get_index()==hy_x_variable->get_index())) {
1006 result << hy_x_variable->GetName();
1007 return;
1008 }
1009 }
1010
1011 HBLObjectRef replace_with = node_variable->Compute();
1012
1013 if (replace_with->ObjectClass () == NUMBER) {
1014 if (mode == kFormulaStringConversionReportRanges) {
1015 result << node_variable->GetName()
1016 << '[';
1017 result.AppendNewInstance(new _String (replace_with->Value()));
1018 result << ':';
1019 result.AppendNewInstance(new _String (node_variable->GetLowerBound()));
1020 result << '-';
1021 result.AppendNewInstance(new _String (node_variable->GetUpperBound()));
1022 result << ']';
1023
1024 } else {
1025 result.AppendNewInstance(new _String (replace_with->Value()));
1026 }
1027 } else if (replace_with->ObjectClass () == STRING) {
1028 result.AppendNewInstance((_String*)replace_with->toStr());
1029 if (this_node_op->IsValueSubstitution()) {
1030 result << "__";
1031 }
1032 } else {
1033 result << node_variable->GetName();
1034 if (this_node_op->IsValueSubstitution()) {
1035 result << "__";
1036 }
1037 }
1038 } else {
1039 _String * variable_name = LocateVar(node_variable_index)->GetName();
1040 if (match_names) {
1041
1042 long f = ((_List*)(*match_names)(0))->FindObject (variable_name);
1043
1044 if (f == kNotFound) {
1045 result<<variable_name;
1046 } else {
1047 result << (_String*)((_List*)(*match_names)(1))->GetItem(f);
1048 }
1049 } else {
1050 result<<variable_name;
1051 }
1052 if (this_node_op->IsValueSubstitution()) {
1053 result << "__";
1054 }
1055 }
1056 return;
1057 }
1058
1059 long node_op_count = this_node_op->GetNoTerms();
1060 if (node_op_count > 0) {
1061 // a built-in operation or a function call
1062 // check if it's a built-in binary operation
1063
1064 long f = _Operation::BinOpCode(this_node_op->GetCode());
1065
1066 if (f != kNotFound) {
1067 // indeed - a binary operation is what we have. check if need to wrap the return in parentheses
1068
1069 if (!top_node || top_node->get_num_nodes() == 2 ) {
1070 unsigned char op_precedence = opPrecedence(f),
1071 op_precedence_right = op_precedence;
1072
1073 if (associativeOps.Find(f) == kNotFound) {
1074 op_precedence_right ++;
1075 }
1076 if (op_level >= 0) { // need to worry about op's precedence
1077 bool need_parens = op_precedence < op_level;
1078
1079 if (need_parens && top_node ) { // put parentheses around the return of this expression
1080 result<<'(';
1081 SubtreeToString (result, top_node->go_down(1),op_precedence,match_names, nil, mode);
1082 result <<&this_node_op->GetCode();
1083 SubtreeToString (result, top_node->go_down(2),op_precedence_right,match_names, nil, mode);
1084 result<<')';
1085 return;
1086 }
1087 }
1088
1089 if (top_node) {
1090 SubtreeToString (result, top_node->go_down(1),op_precedence,match_names,nil, mode);
1091 }
1092 result << &this_node_op->GetCode();
1093 if (top_node) {
1094 SubtreeToString (result, top_node->go_down(2),op_precedence_right,match_names,nil, mode);
1095 }
1096 return;
1097 } else { // mixed binary-unary operation
1098 result<<&this_node_op->GetCode();
1099 if (top_node) {
1100 result<<'(';
1101 SubtreeToString (result, top_node->go_down(1),opPrecedence(f),match_names,nil, mode);
1102 result<<')';
1103 }
1104 return;
1105 }
1106 } else {
1107 if (this_node_op->TheCode() != HY_OP_CODE_MACCESS) {
1108 result<<&this_node_op->GetCode();
1109 if (top_node) {
1110 result<<'(';
1111 SubtreeToString (result, top_node->go_down(1),-1,match_names,nil, mode);
1112 for (long k=2UL; k <= node_op_count; k++) {
1113 result << ',';
1114 SubtreeToString (result, top_node->go_down(k),-1,match_names,nil, mode);
1115 }
1116 result<<')';
1117 }
1118 } else { // matrix element access - treat specially
1119 if (top_node) {
1120 SubtreeToString (result, top_node->go_down(1),-1,match_names,nil,mode);
1121 for (long k=2; k<=node_op_count; k++) {
1122 result<<'[';
1123 SubtreeToString (result, top_node->go_down(k),-1,match_names,nil, mode);
1124 result<<']';
1125 }
1126 }
1127 }
1128 }
1129 return;
1130 }
1131 if (node_op_count < 0L) {
1132 // a user-defined function
1133 long func_id = this_node_op->UserFunctionID();
1134 result<< & GetBFFunctionNameByIndex(func_id);
1135 if (top_node) {
1136 result<<'(';
1137 long argument_count = GetBFFunctionArgumentCount(func_id);
1138 SubtreeToString (result, top_node->go_down(1),-1,match_names,nil,mode);
1139 for (long k=2L; k<=argument_count; k++) {
1140 result<<',';
1141 SubtreeToString (result, top_node->go_down(k),-1,match_names,nil,mode);
1142 }
1143 result<<')';
1144 }
1145 return;
1146 }
1147
1148 HBLObjectRef op_data = this_node_op->GetANumber();
1149 if (op_data) {
1150 _String* conv = (_String*)op_data->toStr();
1151 if (op_data->ObjectClass()==STRING) {
1152 (result <<'"').AppendNewInstance (conv) << '"';
1153 } else {
1154 if (op_data->ObjectClass() == NUMBER && op_data->Value() < 0.0) {
1155 (result <<'(').AppendNewInstance (conv) << ')';
1156 } else {
1157 result.AppendNewInstance (conv);
1158 }
1159 }
1160 } else {
1161 result << "<null>";
1162 }
1163 }
1164
1165 //__________________________________________________________________________________
IsEmpty(void) const1166 bool _Formula::IsEmpty(void) const {
1167 return theFormula.empty();
1168 }
1169
1170 //__________________________________________________________________________________
Newton(_Formula & derivative,_Variable * unknown,hyFloat target_value,hyFloat left,hyFloat right)1171 hyFloat _Formula::Newton(_Formula& derivative, _Variable* unknown, hyFloat target_value, hyFloat left, hyFloat right) {
1172 // find a root of the formulaic expression, using Newton's method, given the derivative and a bracketed root.
1173 // will alternate between bisections and Newton iterations based on what is fatser
1174 // check that there is indeed a sign change on the interval
1175
1176 auto set_and_compute = [&] (hyFloat x) -> hyFloat {
1177 unknown->SetValue(x);
1178 return Compute()->Value()-target_value;
1179 };
1180
1181 hyFloat func_left, func_right, // keeps track of function values in the current interval, [left,right]
1182 root_guess, func_root_guess,
1183 lastCorrection = 100.,
1184 newCorrection;
1185
1186 func_left = set_and_compute (left);
1187 if (func_left==0.0) {
1188 return left;
1189 }
1190 unknown->SetValue(right);
1191 func_right = set_and_compute (right);
1192 if (func_right==0.0) {
1193 return right;
1194 }
1195
1196 if (func_left*func_right>0.0) { // bracket fail
1197 ReportWarning (_String((_String*)toStr(kFormulaStringConversionReportRanges)) & " = "&_String(target_value)&" has no (or multiple) roots in ["&_String(left)&",Inf); " & func_left & "-" & func_right);
1198 return left;
1199 }
1200 // else all is good we can start the machine
1201 bool useNewton = false;
1202
1203 root_guess = (right+left) * .5;
1204
1205 for (unsigned long iterCount = 0L; fabs(right-left)/MAX(left,right) > kMachineEpsilon*10. && iterCount < 200UL; iterCount++) {
1206 func_root_guess = set_and_compute (root_guess);
1207 if (func_root_guess == 0.) {
1208 return root_guess;
1209 }
1210 // get the correction term from the derivative
1211 hyFloat df_dx = derivative.Compute()->Value(),
1212 adjusted_root_guess;
1213
1214 useNewton = true;
1215 if (df_dx==0.0) {
1216 useNewton = false;
1217 } else {
1218 newCorrection = -func_root_guess/df_dx;
1219
1220 if (fabs(newCorrection/func_root_guess)<kMachineEpsilon*2. || fabs(newCorrection)<kMachineEpsilon*2.) { // correction too small - the root has been found
1221 return root_guess;
1222 }
1223
1224 if (fabs(newCorrection/lastCorrection)>4.) { // Newton correction too large - revert to bisection
1225 useNewton = false;
1226 }
1227
1228 adjusted_root_guess = root_guess +newCorrection;
1229 if (adjusted_root_guess<=left || adjusted_root_guess >=right) {
1230 useNewton = false;
1231 } else {
1232 lastCorrection = newCorrection;
1233 }
1234 }
1235
1236 if (useNewton) {
1237 root_guess = adjusted_root_guess;
1238 } else {
1239 if (func_root_guess==0.0) {
1240 return root_guess;
1241 }
1242 if (func_root_guess*func_left > 0.0) { // move to the right half
1243 func_left = func_root_guess;
1244 left = root_guess;
1245 } else {
1246 right = root_guess;
1247 }
1248 root_guess = (right+left) * .5;
1249 }
1250 }
1251 return root_guess;
1252 }
1253
1254
1255 //__________________________________________________________________________________
Brent(_Variable * unknown,hyFloat a,hyFloat b,hyFloat tol,_List * store_evals,hyFloat rhs)1256 hyFloat _Formula::Brent(_Variable* unknown, hyFloat a, hyFloat b, hyFloat tol, _List* store_evals, hyFloat rhs) {
1257 // find a root of the formulaic expression, using Brent's method and a bracketed root.
1258 // check that there is indeed a sign change on the interval
1259
1260 auto set_and_compute = [&] (hyFloat x) -> hyFloat {
1261 unknown->SetValue(x);
1262 hyFloat fx = Compute()->Value()-rhs;
1263
1264 if (store_evals) {
1265 store_evals->AppendNewInstance(new _Constant (x));
1266 store_evals->AppendNewInstance(new _Constant (fx));
1267 }
1268
1269 return fx;
1270 };
1271
1272 hyFloat fa = 0.0,fb = 0.0,fc,d = b-a,e = b-a ,min1,min2,xm,p,q,r,s,tol1,
1273 c = b;
1274
1275 min1 = unknown->GetLowerBound();
1276 min2 = unknown->GetUpperBound();
1277
1278 long it = 0;
1279
1280 if (a>b) { // autobracket to the left
1281 fb = set_and_compute (b);
1282
1283 if (b<0.00001 && b>-0.00001) {
1284 a = b-0.0001;
1285 } else {
1286 a = b-fabs(b)*0.1;
1287 }
1288
1289 if (a<min1) {
1290 a = min1+0.5*(b-min1);
1291 }
1292
1293 fa = set_and_compute (a);
1294
1295 for (long k=0L; k<50L && fb*fa >= 0.0; k++) {
1296
1297 d = (b-a)*GOLDEN_RATIO;
1298 b = a;
1299 a -= d;
1300
1301 if (a<min1) {
1302 if (b>min1) {
1303 a = min1;
1304 } else {
1305 break;
1306 }
1307 }
1308
1309 fb = fa;
1310 fa = set_and_compute (a);
1311 }
1312 } else if (CheckEqual (a,b)) { // autobracket to the right
1313 fa = set_and_compute (a);
1314
1315 a = b;
1316
1317 if ((b<0.00001)&&(b>-0.00001)) {
1318 b = b+0.0001;
1319 } else {
1320 b = b+fabs(b)*0.1;
1321 }
1322
1323 if (b>min2) {
1324 b = a+0.5*(min2-a);
1325 }
1326
1327 fb = set_and_compute (b);
1328
1329 for (long k=0L; k<50L && fb*fa >= 0.0; k++) {
1330
1331 d = (b-a)*GOLDEN_RATIO;
1332 a = b;
1333 b += d;
1334
1335 if (b>min2) {
1336 if (a<min2) {
1337 b = min2;
1338 } else {
1339 break;
1340 }
1341 }
1342
1343 fa = fb;
1344 fb = set_and_compute (b);
1345 }
1346 }
1347
1348
1349 if (fa == 0.0) {
1350 fa = set_and_compute (a);
1351 if (fa == 0.0) {
1352 return a;
1353 }
1354 }
1355
1356 if (fb == 0.0) {
1357 fb = set_and_compute (b);
1358 if (fb == 0.0) {
1359 return b;
1360 }
1361 }
1362
1363 if (fa*fb<0.0) {
1364 fc = fb;
1365 c = b;
1366
1367 for (it = 0; it < MAX_BRENT_ITERATES; it++) {
1368 if (fb*fc>0.0) {
1369 fc = fa;
1370 c = a;
1371 e = d = b-a;
1372 }
1373
1374 if (fabs (fc) < fabs (fb)) {
1375 a = b;
1376 b = c;
1377 c = a;
1378 fa = fb;
1379 fb = fc;
1380 fc = fa;
1381 }
1382
1383 tol1 = 2.*fabs(b)*kMachineEpsilon+.5*tol;
1384
1385 xm = .5*(c-b);
1386
1387 if (fabs(xm)<=tol1 || fb == 0.0) {
1388 return b;
1389 }
1390
1391 if (fabs(e)>=tol1 && fabs (fa) > fabs (fb)) {
1392 s = fb/fa;
1393 if (a==c) {
1394 p = 2.*xm*s;
1395 q = 1.-s;
1396 } else {
1397 q = fa/fc;
1398 r = fb/fc;
1399 p = s*(2.*xm*q*(q-r)-(b-a)*(r-1.0));
1400 q = (q-1.)*(r-1.)*(s-1.);
1401 }
1402 if (p>0.0) {
1403 q = -q;
1404 }
1405
1406 if (p<0.0) {
1407 p = -p;
1408 }
1409
1410 min1 = 3.*xm*q-fabs(tol1*q);
1411 min2 = fabs (e*q);
1412 if (2.*p < (min1<min2?min1:min2)) {
1413 e = d;
1414 d = p/q;
1415 } else {
1416 d = xm;
1417 e = d;
1418 }
1419 } else {
1420 d = xm;
1421 e = d;
1422 }
1423 a = b;
1424 fa = fb;
1425 if (fabs(d)>tol1) {
1426 b+=d;
1427 } else {
1428 if (xm > 0.) {
1429 b += fabs (tol1);
1430 } else {
1431 b -= fabs (tol1);
1432 }
1433 }
1434 fb = set_and_compute (b);
1435 }
1436 }
1437
1438
1439 /*for (long i = 0; i < theFormula.lLength; i++) {
1440 _Operation *op_i = GetIthTerm(i);
1441 printf ("%ld: %s\n", i+1, _String((_String*)op_i->toStr()).sData);
1442 }*/
1443
1444 _String msg ((_String*)toStr(kFormulaStringConversionSubstiteValues));
1445 msg = msg & "=" & rhs;
1446 if (it < MAX_BRENT_ITERATES) {
1447 msg = msg & " has no (or multiple) roots in ["&_String(a)&","&_String(b)&"]";
1448 } else {
1449 msg = msg & " failed to converge to sufficient precision in " & MAX_BRENT_ITERATES &" iterates.";
1450 }
1451 ReportWarning (msg);
1452 return b;
1453 }
1454 //__________________________________________________________________________________
Newton(_Formula & derivative,hyFloat target_value,hyFloat left,hyFloat max_right,_Variable * unknown)1455 hyFloat _Formula::Newton(_Formula& derivative, hyFloat target_value, hyFloat left, hyFloat max_right, _Variable* unknown) {
1456 // given a monotone function and a left bracket bound, found the right bracket bound and solve
1457
1458 // check that there is indeed a sign change on the interval
1459 unknown->SetValue(left);
1460 hyFloat t1 = Compute()->Value(), right = left, t2, step = 1.0;
1461
1462 if (max_right-left < step * 100.) {
1463 step = (max_right-left) * 0.01;
1464 }
1465 if (step==0.0) {
1466 return left;
1467 }
1468 do {
1469 right += step;
1470 if (right>max_right) { // function doesn't seem to have a root
1471 ReportWarning (_String((_String*)toStr(kFormulaStringConversionSubstiteValues))&"="&_String(target_value)&" has no (or multiple) roots in ["&_String(left)&","&_String(right)&")");
1472 return left;
1473 }
1474 unknown->SetValue(right);
1475 t2 = Compute()->Value();
1476 step*=2;
1477 if (right+step>max_right)
1478 if (right<max_right) {
1479 step = max_right-right;
1480 }
1481 } while ((target_value-t1)*(target_value-t2)>0);
1482 return Newton (derivative, unknown, target_value, left, right);
1483
1484 }
1485
1486 //__________________________________________________________________________________
Newton(_Variable * unknown,hyFloat target_value,hyFloat x_min,hyFloat left,hyFloat right)1487 hyFloat _Formula::Newton(_Variable* unknown, hyFloat target_value, hyFloat x_min, hyFloat left, hyFloat right) {
1488 // check that there is indeed a sign change on the interval
1489 hyFloat t1,t2,t3,t4,t5,lastCorrection = 100, newCorrection;
1490 _String msg;
1491 t1 =Integral(unknown, x_min, left)-target_value;
1492 if (t1==0.0) {
1493 return left;
1494 }
1495 t2 = t1+Integral(unknown, left, right);
1496 if (t2==0.0) {
1497 return right;
1498 }
1499 if (t1*t2>0.0) {
1500 ReportWarning (_String((_String*)toStr(kFormulaStringConversionSubstiteValues))&"="&_String(target_value)&" has no (or multiple) roots in ["&_String(left)&","&_String(right)&")");
1501 return left;
1502 }
1503 // else all is good we can start the machine
1504 bool useNewton = false;
1505
1506 t3 = (right + left) * 0.5;
1507
1508 while (right-left>1e-6) { // stuff to do
1509 if (!useNewton) {
1510 t3 = (right+left) * 0.5;
1511 }
1512 unknown->SetValue(t3);
1513 t4 = Integral(unknown, x_min, t3)-target_value;
1514 // get the correction term from the derivative
1515 unknown->SetValue(t3);
1516 t5 = Compute()->Value();
1517 useNewton = true;
1518 if (t5==0.0) {
1519 useNewton = false;
1520 } else {
1521 newCorrection = -t4/t5;
1522 if (fabs(newCorrection)<1e-5) { // correction too small - the root has been found
1523 return t3;
1524 }
1525 if (fabs(newCorrection/lastCorrection)>4) { // Newton correction too large - revert to bisection
1526 useNewton = false;
1527 }
1528 t5 = t3+newCorrection;
1529 if ((t5<=left)||(t5>=right)) {
1530 useNewton = false;
1531 } else {
1532 lastCorrection = newCorrection;
1533 }
1534 }
1535 if (useNewton) {
1536 t3 = t5;
1537 } else {
1538 t4 = Integral(unknown, x_min, t3)-target_value;
1539 if (t4==0.0) {
1540 return t3;
1541 }
1542 if (t4*t1 >0) {
1543 t1 = t4;
1544 left = t3;
1545 } else {
1546 right = t3;
1547 }
1548 }
1549
1550 }
1551 return t3;
1552 }
1553
1554 //__________________________________________________________________________________
Newton(_Variable * unknown,hyFloat target_value,hyFloat x_min,hyFloat left)1555 hyFloat _Formula::Newton( _Variable* unknown, hyFloat target_value,hyFloat x_min, hyFloat left) {
1556 // given a monotone function and a left bracket bound, found the right bracket bound and solve
1557 // check that there is indeed a sign change on the interval
1558 hyFloat t1 = Integral(unknown, x_min, left), right = left, t2, step = 1.0;
1559 do {
1560 right += step;
1561 t2 = Integral(unknown, right-step, right);
1562 step*=2;
1563 if (right>=1e10) { // function doesn't seem to have a root
1564 ReportWarning (_String((_String*)toStr(kFormulaStringConversionSubstiteValues))&"="&_String(target_value)&" has no (or multiple) roots in ["&_String(left)&",Inf)");
1565 return 0.0;
1566 }
1567 } while ((target_value-t1)*(target_value-t2-t1)>=0);
1568 return Newton ( unknown, target_value,x_min, left, right);
1569
1570 }
1571
1572 //__________________________________________________________________________________
1573
Integral(_Variable * dx,hyFloat left,hyFloat right,bool infinite)1574 hyFloat _Formula::Integral(_Variable* dx, hyFloat left, hyFloat right, bool infinite) {
1575 // uses Romberg's intergation method
1576 if (infinite) { // tweak "right" here
1577 hyFloat value = 1.0, step = 1.0, right1= -1.;
1578 right = left;
1579 while (value>1e-8) {
1580 right+=step;
1581 dx->SetValue(right);
1582 value = fabs(Compute()->Value());
1583 if ( value < 1e-4 && right1 < 0.) { // SLKP 20181205 why is this here?
1584 right1 = right;
1585 }
1586 step *= 2;
1587 if (step > 100000) { // integrand decreasing too slowly
1588 HandleApplicationError (_String(*(_String*)toStr(kFormulaStringConversionNormal)).Enquote() & " decreases too slowly to be integrated over an infinite interval");
1589 return 0.0;
1590 }
1591 }
1592
1593
1594 if (right1<right-kMachineEpsilon) {
1595 return Integral(dx,left,right1,false)+Integral(dx,right1,right,false);
1596 } else {
1597 return Integral(dx,left,right1,false);
1598 }
1599 }
1600
1601 hyFloat precision_factor = hy_env::EnvVariableGetNumber(hy_env::integration_precision_factor);
1602 long max_iterations = hy_env::EnvVariableGetNumber(hy_env::integration_maximum_iterations);
1603
1604 hyFloat ss = 0.,
1605 dss,
1606 *s = new hyFloat [max_iterations],
1607 *h = new hyFloat [max_iterations+1L];
1608
1609
1610 h[0]=1.0;
1611
1612 long interpolate_steps = 5L,
1613 stack_depth = 0L;
1614
1615 _SimpleList fvidx_aux,
1616 changing_vars,
1617 idx_map;
1618
1619
1620 _AVLList fvidx (&fvidx_aux);
1621
1622 hyFloat * ic = new hyFloat[interpolate_steps],
1623 * id = new hyFloat[interpolate_steps];
1624
1625 _SimpleFormulaDatum
1626 * stack = nil,
1627 * vvals = nil;
1628
1629
1630 if (AmISimple (stack_depth,fvidx)) {
1631 stack = new _SimpleFormulaDatum [stack_depth];
1632
1633 ConvertToSimple (fvidx);
1634 if (fvidx.countitems()) { // could be a constant expression with no variable dependancies
1635 vvals = new _SimpleFormulaDatum [fvidx.countitems()];
1636 //fvidx.ReorderList();
1637 long dx_index = dx->get_index();
1638 for (long vi = 0; vi < fvidx_aux.countitems(); vi++) {
1639 _Variable* variable_in_expression = LocateVar (fvidx_aux.Element(vi));
1640 if (variable_in_expression->CheckFForDependence (dx_index,true)) {
1641 changing_vars << fvidx_aux.Element(vi);
1642 idx_map << vi;
1643 }
1644 vvals[vi].value = variable_in_expression->Compute()->Value();
1645 }
1646
1647 long depends_on_dx = fvidx.FindLong(dx_index);
1648
1649 if (depends_on_dx >= 0) {
1650 changing_vars.InsertElement ((BaseRef)dx_index,0,false,false);
1651 idx_map.InsertElement ((BaseRef)fvidx.FindLong(dx_index),0,false,false);
1652 }
1653 } else {
1654 vvals = new _SimpleFormulaDatum [1L];
1655 }
1656 } else {
1657 stack_depth = -1L;
1658 }
1659
1660 for (long step = 0L; step < max_iterations; step ++) {
1661 //printf ("%d\n", step);
1662 if (stack_depth >=0) { // compiled
1663 s[step] = TrapezoidLevelKSimple(*this, dx, left, right, step+1L, stack, vvals,changing_vars,idx_map);
1664 } else {
1665 s[step] = TrapezoidLevelK(*this, dx, left, right, step+1L);
1666 }
1667 if (step >= 4L) {
1668 ss = InterpolateValue(&h[step-4L],&s[step-4L],interpolate_steps,ic,id,0.0, dss);
1669 if (fabs(dss)<= precision_factor * fabs(ss)) {
1670
1671 BatchDeleteArray (s,h,ic,id);
1672 if (stack_depth >=0L ) {
1673 ConvertFromSimple(fvidx);
1674 BatchDeleteArray (stack, vvals);
1675 }
1676 return ss;
1677 }
1678 }
1679 h[step+1] = h[step]/9.0;
1680 }
1681
1682 if (stack_depth >=0) {
1683 ConvertFromSimple(fvidx);
1684 BatchDeleteArray (stack, vvals);
1685 }
1686
1687 ReportWarning (_String("Integral of ")&*_String((_String*)toStr(kFormulaStringConversionNormal)) & " over ["&_String(left)&","&_String(right)&"] converges slowly, loss of precision may occur. Change either INTEGRATION_PRECISION_FACTOR or INTEGRATION_MAX_ITERATES");
1688 BatchDeleteArray (s,h,ic,id);
1689 return ss;
1690 }
1691 //__________________________________________________________________________________
InterpolateValue(hyFloat * theX,hyFloat * theY,long n,hyFloat * c,hyFloat * d,hyFloat x,hyFloat & err)1692 hyFloat InterpolateValue (hyFloat* theX, hyFloat* theY, long n, hyFloat *c , hyFloat *d, hyFloat x, hyFloat& err) {
1693 // Neville's algoruthm for polynomial interpolation (Numerical Recipes' rawinterp)
1694 hyFloat y,
1695 den,
1696 dif = 1e10,
1697 dift,
1698 ho,
1699 hp,
1700 w;
1701
1702 long ns = 0L;
1703
1704 for (long i=0L; i<n; i++) {
1705 dift = fabs(x-theX[i]);
1706 if (dift<dif) {
1707 ns = i;
1708 dif = dift;
1709 }
1710 c[i] = d[i] = theY[i];
1711 }
1712
1713 y = theY[ns--];
1714
1715 for (long m=1L; m<n; m++) {
1716 for (long i=0L; i<=n-m-1; i++) {
1717 ho = theX[i]-x;
1718 hp = theX[i+m]-x;
1719 w = c[i+1]-d[i];
1720 den = w/(ho-hp);
1721 d[i] = hp*den;
1722 c[i] = ho*den;
1723 }
1724 err = 2*ns< (n-m)? c[ns+1]: d[ns--];
1725 y += err;
1726 }
1727
1728 return y;
1729 }
1730
1731
1732
1733 //__________________________________________________________________________________
1734
TrapezoidLevelKSimple(_Formula & f,_Variable * xvar,hyFloat left,hyFloat right,long k,_SimpleFormulaDatum * stack,_SimpleFormulaDatum * values,_SimpleList & changingVars,_SimpleList & varToStack)1735 hyFloat TrapezoidLevelKSimple (_Formula&f, _Variable* xvar, hyFloat left, hyFloat right, long k, _SimpleFormulaDatum * stack, _SimpleFormulaDatum* values, _SimpleList& changingVars, _SimpleList& varToStack) {
1736
1737 hyFloat x,
1738 tnm,
1739 sum,
1740 del,
1741 ddel;
1742
1743 static hyFloat s;
1744
1745 //_Constant dummy;
1746
1747 auto set_and_compute = [&] (hyFloat x) -> hyFloat {
1748 if (changingVars.countitems() == 1L) {
1749 values[varToStack.get(0)].value = x;
1750 } else {
1751 xvar->SetNumericValue (x);
1752 changingVars.Each ([&] (long v, unsigned long vi) -> void {
1753 values[varToStack.get(vi)].value = LocateVar(v)->Compute()->Value();
1754 });
1755 }
1756 return f.ComputeSimple(stack, values);
1757 };
1758
1759
1760 if (k==1) {
1761 s = set_and_compute ((left+right)*0.5);
1762 return s;
1763 }
1764
1765 long it =1;
1766
1767 for (long j=1L; j<k-1; j++) {
1768 it *= 3L;
1769 }
1770
1771 tnm = it;
1772 del = (right-left)/(3.0*tnm);
1773 ddel = del+del;
1774 x = left+del*.5;
1775 sum = 0.;
1776 for (long j=1L; j<=it; j++, x+=del) {
1777 sum += set_and_compute (x);
1778 x+=ddel;
1779 sum += set_and_compute (x);
1780 }
1781 s = (s+(right-left)*sum/tnm)/3.0;
1782 return s;
1783 }
1784
1785 //__________________________________________________________________________________
TrapezoidLevelK(_Formula & f,_Variable * xvar,hyFloat left,hyFloat right,long k)1786 hyFloat TrapezoidLevelK (_Formula&f, _Variable* xvar, hyFloat left, hyFloat right, long k) {
1787 hyFloat x,
1788 tnm,
1789 sum,
1790 del,
1791 ddel;
1792
1793 static hyFloat s;
1794
1795 long it = 1L;
1796
1797 if (k==1L) {
1798 xvar -> SetNumericValue((left + right) * 0.5);
1799 return s = f.Compute()->Value();
1800 }
1801
1802 for (long j = 1L; j < k-1; j++) {
1803 it *= 3L;
1804 }
1805
1806
1807 tnm = it;
1808 del = (right-left)/(3.0*tnm);
1809 ddel = del+del;
1810 x = left+del*.5;
1811
1812 sum = 0.0;
1813
1814 for (long j=1L; j<=it; j++, x+=del) {
1815 xvar->SetNumericValue(x);
1816 sum += f.Compute()->Value();
1817 x+=ddel;
1818 xvar->SetNumericValue(x);
1819 sum += f.Compute()->Value();
1820 }
1821 s = (s+(right-left)*sum/tnm)/3.0;
1822 return s;
1823 }
1824
1825 //__________________________________________________________________________________
MeanIntegral(_Variable * dx,hyFloat left,hyFloat right,bool infinite)1826 hyFloat _Formula::MeanIntegral(_Variable* dx, hyFloat left, hyFloat right, bool infinite) {
1827 _Formula newF;
1828 newF.Duplicate(this);
1829 newF.theFormula < new _Operation (true, *(dx->GetName())) < new _Operation (HY_OP_CODE_MUL, 2);
1830 return newF.Integral (dx,left,right, infinite);
1831 }
1832
1833 //__________________________________________________________________________________
NumberOperations(void) const1834 long _Formula::NumberOperations(void) const {
1835 // number of operations in the formula
1836 return theFormula.countitems();
1837 }
1838
1839 //__________________________________________________________________________________
1840
ExtractMatrixExpArguments(_List * storage)1841 long _Formula::ExtractMatrixExpArguments (_List* storage) {
1842 long count = 0L;
1843
1844 if (resultCache && resultCache->empty() == false) {
1845 long cacheID = 0L;
1846 bool cacheUpdated = false;
1847 // whether or not a cached result was used
1848
1849
1850 for (unsigned long i=0UL; i + 1UL <theFormula.countitems(); i++) {
1851 _Operation* this_op = GetIthTerm(i),
1852 * next_op = GetIthTerm(i + 1UL);
1853
1854 if (! cacheUpdated && next_op->CanResultsBeCached(this_op)) {
1855 /*
1856 StringToConsole("\n----\n");
1857 ObjectToConsole(FetchVar(LocateVarByName("k"))->Compute());
1858 NLToConsole();*/
1859
1860 _Stack temp;
1861 this_op->Execute (temp);
1862
1863
1864 _Matrix *currentArg = (_Matrix*)temp.Pop(true),
1865 *cachedArg = (_Matrix*)((HBLObjectRef)(*resultCache)(cacheID)),
1866 *diff = nil;
1867
1868 if (cachedArg->ObjectClass() == MATRIX) {
1869 diff = (_Matrix*)cachedArg->SubObj(currentArg, nil);
1870 }
1871
1872 if (diff && diff->MaxElement() <= 1e-12) {
1873 cacheID += 2;
1874 i ++;
1875 } else {
1876 cacheUpdated = true;
1877 cacheID++;
1878 if (next_op->CanResultsBeCached(this_op, true)) {
1879 storage->AppendNewInstance(currentArg);
1880 count ++;
1881 }
1882 }
1883 DeleteObject (diff);
1884 continue;
1885 }
1886 if (cacheUpdated) {
1887 cacheID++;
1888 cacheUpdated = false;
1889 }
1890 }
1891 }
1892
1893 return count;
1894 }
1895
1896 //__________________________________________________________________________________
Dereference(bool ignore_context,_hyExecutionContext * theContext)1897 _Variable * _Formula::Dereference (bool ignore_context, _hyExecutionContext* theContext) {
1898 _Variable * result = nil;
1899 HBLObjectRef computedValue = Compute (0, theContext->GetContext(), nil, theContext->GetErrorBuffer());
1900 if (computedValue && computedValue->ObjectClass() == STRING) {
1901 result = (_Variable*)((_FString*)computedValue)->Dereference(ignore_context, theContext, true);
1902 }
1903
1904 if (!result) {
1905 theContext->ReportError( (_String ("Failed to dereference '") & _String ((_String*)toStr(kFormulaStringConversionNormal)) & "' in the " & (ignore_context ? "global" : "local") & " context"));
1906 }
1907
1908 return result;
1909 }
1910
1911 //unsigned long ticker = 0UL;
1912
1913 //__________________________________________________________________________________
Compute(long startAt,_VariableContainer const * nameSpace,_List * additionalCacheArguments,_String * errMsg,long valid_type,bool can_cache)1914 HBLObjectRef _Formula::Compute (long startAt, _VariableContainer const * nameSpace, _List* additionalCacheArguments, _String* errMsg, long valid_type, bool can_cache) {
1915 // compute the value of the formula
1916 // TODO SLKP 20170925 Needs code review
1917 _Stack * scrap_here;
1918 current_formula_being_computed = this;
1919 if (theFormula.empty()) {
1920 theStack.theStack.Clear();
1921 theStack.Push (new _MathObject, false);
1922 scrap_here = &theStack;
1923 } else {
1924 bool wellDone = true;
1925
1926
1927 if (call_count++) {
1928 scrap_here = new _Stack;
1929 } else {
1930 scrap_here = &theStack;
1931 if (startAt == 0) {
1932 theStack.Reset();
1933 }
1934 }
1935
1936 /*ticker++;
1937 if (ticker >= 1462440) {
1938 printf ("\n_Formula::Compute (%x, %d) %ld terms, stack depth %ld\n", this, ticker, theFormula.lLength, theStack.theStack.lLength);
1939 }*/
1940
1941 const unsigned long term_count = NumberOperations();
1942
1943 if (startAt == 0L && resultCache && !resultCache->empty()) {
1944 long cacheID = 0L;
1945 // where in the cache are we currently looking
1946 bool cacheUpdated = false;
1947 // whether or not a cached result was used
1948
1949 for (unsigned long i=0; i< term_count ; i++) {
1950 _Operation* thisOp = ItemAt (i);
1951 if ( i + 1UL < term_count ) {
1952 _Operation* nextOp = ItemAt (i+1UL);
1953
1954 if (! cacheUpdated && nextOp->CanResultsBeCached(thisOp)) {
1955 if (!thisOp->Execute(*scrap_here,nameSpace, errMsg)) {
1956 wellDone = false;
1957 break;
1958 }
1959
1960 _Matrix *currentArg = (_Matrix*)scrap_here->Pop(false),
1961 *cachedArg = (_Matrix*)((HBLObjectRef)(*resultCache)(cacheID)),
1962 *diff = nil;
1963
1964 if (cachedArg->ObjectClass() == MATRIX) {
1965 diff = (_Matrix*)cachedArg->SubObj(currentArg, nil);
1966 }
1967
1968 bool no_difference = diff && diff->MaxElement() <= 1e-12;
1969
1970 if (no_difference || (additionalCacheArguments && !additionalCacheArguments->empty() && nextOp->CanResultsBeCached(thisOp,true))) {
1971 DeleteObject (scrap_here->Pop ());
1972 if (no_difference) {
1973 scrap_here->Push ((HBLObjectRef)(*resultCache)(cacheID+1));
1974 } else {
1975
1976 scrap_here->Push ((HBLObjectRef)additionalCacheArguments->GetItem (0));
1977 resultCache->Replace(cacheID,scrap_here->Pop(false),true);
1978 resultCache->Replace(cacheID+1,(HBLObjectRef)additionalCacheArguments->GetItem (0),false);
1979 additionalCacheArguments->Delete (0, false);
1980 //printf ("_Formula::Compute additional arguments %ld\n", additionalCacheArguments->lLength);
1981 }
1982 cacheID += 2;
1983 i ++;
1984 //printf ("Used formula cache %s\n", _String((_String*)nextOp->toStr()).sData);
1985 } else {
1986 cacheUpdated = true;
1987 resultCache->Replace(cacheID++,scrap_here->Pop(false),true);
1988 //printf ("Updated formula cache %s\n", _String((_String*)nextOp->toStr()).sData);
1989 }
1990 DeleteObject (diff);
1991 continue;
1992 }
1993 }
1994 if (!thisOp->Execute(*scrap_here,nameSpace, errMsg, can_cache && call_count == 1)) { // does this always get executed?
1995 wellDone = false;
1996 break;
1997 }
1998 if (cacheUpdated) {
1999 resultCache->Replace(cacheID++,scrap_here->Pop(false),true);
2000 cacheUpdated = false;
2001 }
2002 }
2003 } else {
2004
2005 for (unsigned long i=startAt; i< term_count; i++) {
2006 if (!ItemAt (i)->Execute(*scrap_here, nameSpace, errMsg, can_cache && call_count == 1)) {
2007 wellDone = false;
2008 break;
2009 }
2010
2011 }
2012 }
2013 if (scrap_here->StackDepth() != 1L || !wellDone) {
2014 _String errorText = _String ("'") & _String((_String*)toStr(kFormulaStringConversionNormal)) & _String("' evaluated with errors ");
2015 if (errMsg && errMsg->nonempty()) {
2016 errorText = errorText & " " & errMsg->Enquote('(',')');
2017 }
2018
2019 if (scrap_here->StackDepth() > 1 && wellDone) {
2020 errorText = errorText & " Unconsumed values on the stack";
2021 for (long stack_id = scrap_here->StackDepth()-1; stack_id >= 0; stack_id --) {
2022 errorText = errorText & "\n[" & (stack_id+1) & "]------------------\n" & (_String*) scrap_here->Pop(false)->toStr();
2023 }
2024 errorText & "\n------------------\n";
2025 }
2026
2027 current_formula_being_computed = nil;
2028
2029 if (errMsg) {
2030 *errMsg = *errMsg & errorText;
2031 }
2032 else {
2033 HandleApplicationError (errorText);
2034 }
2035 scrap_here->Reset();
2036 scrap_here->Push (new _Constant (0.0), false);
2037 }
2038 }
2039
2040
2041 HBLObjectRef return_value = scrap_here->Pop(false);
2042
2043 if (theFormula.lLength) {
2044 DeleteObject (recursion_calls);
2045 if (--call_count) {
2046 recursion_calls = return_value;
2047 return_value->AddAReference();
2048 if (scrap_here != &theStack) {
2049 delete scrap_here;
2050 }
2051
2052 } else {
2053 recursion_calls = nil;
2054 }
2055 }
2056 current_formula_being_computed = nil;
2057 return valid_type == HY_ANY_OBJECT ? return_value : ((return_value->ObjectClass() & valid_type) ? return_value : nil);
2058 }
2059
2060 //__________________________________________________________________________________
CheckSimpleTerm(HBLObjectRef thisObj)2061 bool _Formula::CheckSimpleTerm (HBLObjectRef thisObj) {
2062 if (thisObj) {
2063 long oc = thisObj->ObjectClass();
2064 if (oc != NUMBER) {
2065 if (oc ==MATRIX) {
2066 _Matrix * mv = (_Matrix*)thisObj->Compute();
2067 if (mv->is_dense() && mv->IsIndependent ()) {
2068 return true;
2069 }
2070 }
2071 } else {
2072 return true;
2073 }
2074 }
2075 return false;
2076 }
2077
2078 //__________________________________________________________________________________
operator +(const _Formula & operand2)2079 _Formula const _Formula::operator+ (const _Formula& operand2) {
2080 return *PatchFormulasTogether (*this,operand2,HY_OP_CODE_ADD);
2081 }
2082
2083 //__________________________________________________________________________________
operator -(const _Formula & operand2)2084 _Formula const _Formula::operator- (const _Formula& operand2) {
2085 return *PatchFormulasTogether (*this,operand2,HY_OP_CODE_SUB);
2086 }
2087
2088 //__________________________________________________________________________________
operator *(const _Formula & operand2)2089 _Formula const _Formula::operator* (const _Formula& operand2) {
2090 return *PatchFormulasTogether (*this,operand2,HY_OP_CODE_MUL);
2091 }
2092
2093 //__________________________________________________________________________________
operator /(const _Formula & operand2)2094 _Formula const _Formula::operator/ (const _Formula& operand2) {
2095 return *PatchFormulasTogether (*this,operand2,HY_OP_CODE_DIV);
2096 }
2097
2098 //__________________________________________________________________________________
operator ^(const _Formula & operand2)2099 _Formula const _Formula::operator^ (const _Formula& operand2) {
2100 return *PatchFormulasTogether (*this,operand2,HY_OP_CODE_POWER);
2101 }
2102
2103 //__________________________________________________________________________________
PatchFormulasTogether(const _Formula & op1,const _Formula & op2,const char op_code)2104 _Formula* _Formula::PatchFormulasTogether (const _Formula& op1, const _Formula& op2, const char op_code) {
2105 _Formula * result = new _Formula;
2106 result->DuplicateReference(&op1);
2107 result->DuplicateReference(&op2);
2108 result->theFormula.AppendNewInstance(new _Operation (op_code, 2));
2109 return result;
2110 }
2111
2112 //__________________________________________________________________________________
PatchFormulasTogether(const _Formula & op1,HBLObjectRef op2,const char op_code)2113 _Formula* _Formula::PatchFormulasTogether (const _Formula& op1, HBLObjectRef op2, const char op_code) {
2114 _Formula * result = new _Formula;
2115 result->DuplicateReference(&op1);
2116 result->theFormula.AppendNewInstance(new _Operation (op2));
2117 result->theFormula.AppendNewInstance(new _Operation (op_code, 2));
2118 return result;
2119 }
2120
2121
2122 //__________________________________________________________________________________
ConvertMatrixArgumentsToSimpleOrComplexForm(bool make_complex)2123 void _Formula::ConvertMatrixArgumentsToSimpleOrComplexForm (bool make_complex) {
2124 unsigned long term_count = NumberOperations();
2125
2126 if (make_complex) {
2127 if (resultCache) {
2128 DeleteObject (resultCache);
2129 resultCache = nil;
2130 }
2131 } else {
2132 if (!resultCache) {
2133 resultCache = new _List();
2134 for (unsigned long i = 1UL ; i< term_count ; i++) {
2135 if ( ItemAt(i)->CanResultsBeCached(ItemAt (i-1))) {
2136 resultCache->AppendNewInstance(new _MathObject());
2137 resultCache->AppendNewInstance(new _MathObject());
2138 }
2139 }
2140 }
2141 }
2142
2143 for (unsigned long i = 0UL ; i< term_count ; i++) {
2144 _Operation* this_op = ItemAt (i);
2145 _Matrix * op_matrix = nil;
2146
2147 if (this_op->theNumber) {
2148 if (this_op->theNumber->ObjectClass() == MATRIX) {
2149 op_matrix = (_Matrix*)this_op->theNumber;
2150 }
2151 } else {
2152 if (this_op->theData >= 0L) {
2153 _Variable* thisVar = LocateVar (this_op->theData);
2154 if (thisVar->ObjectClass() == MATRIX) {
2155 op_matrix = (_Matrix*)thisVar->GetValue();
2156 }
2157 }
2158 }
2159
2160 if (op_matrix) {
2161 if (make_complex) {
2162 op_matrix->MakeMeGeneral();
2163 } else {
2164 op_matrix->MakeMeSimple();
2165 }
2166 }
2167 }
2168 }
2169
2170 //__________________________________________________________________________________
StackDepth(long from,long to) const2171 long _Formula::StackDepth (long from, long to) const {
2172 _SimpleList::NormalizeCoordinates(from, to, NumberOperations());
2173 long result = 0L;
2174
2175 for ( long i = from; i <= to; i++) {
2176 result += ItemAt(i)->StackDepth ();
2177 }
2178
2179 return result;
2180
2181 }
2182
2183
2184 //__________________________________________________________________________________
AmISimple(long & stack_depth,_AVLList & variable_index)2185 bool _Formula::AmISimple (long& stack_depth, _AVLList& variable_index) {
2186 if (theFormula.empty()) {
2187 return true;
2188 }
2189
2190 long loc_depth = 0L;
2191
2192 for (unsigned long i=0UL; i<theFormula.countitems(); i++) {
2193 _Operation* this_op =ItemAt (i);
2194 loc_depth++;
2195 if ( this_op->theData<-2 || this_op->numberOfTerms<0) {
2196 if (this_op->theData < -2 && i == 0UL) {
2197 variable_index.InsertNumber (this_op->GetAVariable());
2198 continue;
2199 }
2200 return false;
2201 }
2202
2203 if (this_op->theNumber) {
2204 if (this_op->theNumber->ObjectClass() != NUMBER) {
2205 return false;
2206 }
2207 } else {
2208 if (this_op->theData >= 0) {
2209 _Variable* this_var = LocateVar (this_op->theData);
2210 if (this_var->ObjectClass()!=NUMBER) {
2211 HBLObjectRef cv = this_var->GetValue();
2212 if (!CheckSimpleTerm (cv)) {
2213 return false;
2214 }
2215 }
2216 variable_index.InsertNumber (this_op->GetAVariable());
2217 } else {
2218 long op_code = this_op->TheCode();
2219
2220 if (simpleOperationCodes.Find(op_code)==kNotFound) {
2221 return false;
2222 } else if ((op_code == HY_OP_CODE_MACCESS || op_code == HY_OP_CODE_MCOORD || op_code == HY_OP_CODE_MUL) && this_op->GetNoTerms() != 2) {
2223 return false;
2224 }
2225
2226 loc_depth -= this_op->GetNoTerms();
2227 }
2228 }
2229 if (loc_depth>stack_depth) {
2230 stack_depth = loc_depth;
2231 } else if (loc_depth==0L) {
2232 HandleApplicationError (_String("Invalid formula (no return value) passed to ") & __PRETTY_FUNCTION__ & " :" & _String ((_String*)toStr(kFormulaStringConversionNormal)).Enquote());
2233 return false;
2234 }
2235 }
2236 return true;
2237 }
2238
2239 //__________________________________________________________________________________
IsArrayAccess(void)2240 bool _Formula::IsArrayAccess (void){
2241 if (!theFormula.empty()) {
2242 return (ItemAt (theFormula.countitems()-1)->TheCode() == HY_OP_CODE_MACCESS);
2243 }
2244 return false;
2245 }
2246
2247 //__________________________________________________________________________________
ConvertToSimple(_AVLList & variable_index)2248 bool _Formula::ConvertToSimple (_AVLList& variable_index) {
2249 bool has_volatiles = false;
2250 if (simpleExpressionStatus) {
2251 delete [] simpleExpressionStatus;
2252 simpleExpressionStatus = nil;
2253 }
2254 if (!theFormula.empty()) {
2255
2256 simpleExpressionStatus = new long [theFormula.countitems()];
2257
2258 for (unsigned long i=0UL; i<theFormula.countitems(); i++) {
2259 _Operation* this_op = ItemAt (i);
2260 simpleExpressionStatus[i] = -4L;
2261 if (this_op->theNumber) {
2262 simpleExpressionStatus[i] = -1L;
2263 } else if (this_op->theData >= 0) {
2264 this_op->theData = variable_index.FindLong (this_op->theData);
2265 simpleExpressionStatus[i] = this_op->theData;
2266 } else if (this_op->opCode == HY_OP_CODE_SUB && this_op->numberOfTerms == 1) {
2267 this_op->opCode = (long)MinusNumber;
2268 } else {
2269 if (this_op->opCode == HY_OP_CODE_MACCESS) {
2270 this_op->numberOfTerms = -2;
2271 } else {
2272 if (this_op->opCode == HY_OP_CODE_MCOORD) {
2273 this_op->numberOfTerms = -3;
2274 }
2275 }
2276 if (this_op->opCode == HY_OP_CODE_RANDOM || this_op->opCode == HY_OP_CODE_TIME)
2277 has_volatiles = true;
2278
2279 if (this_op->numberOfTerms == 2) {
2280 switch (this_op->opCode) {
2281 case HY_OP_CODE_ADD:
2282 case HY_OP_CODE_SUB:
2283 case HY_OP_CODE_MUL:
2284 case HY_OP_CODE_DIV:
2285 simpleExpressionStatus[i] = -10000L - this_op->opCode;
2286 }
2287 }
2288 this_op->opCode = simpleOperationFunctions(simpleOperationCodes.Find(this_op->opCode));
2289 }
2290
2291 }
2292 }
2293 return has_volatiles;
2294 }
2295
2296 //__________________________________________________________________________________
ConvertFromSimple(_AVLList const & variableIndex)2297 void _Formula::ConvertFromSimple (_AVLList const& variableIndex) {
2298 delete [] simpleExpressionStatus;
2299 simpleExpressionStatus = nil;
2300 ConvertFromSimpleList (*variableIndex.dataList);
2301 }
2302
2303 //__________________________________________________________________________________
ConvertFromSimpleList(_SimpleList const & variableIndex)2304 void _Formula::ConvertFromSimpleList (_SimpleList const& variableIndex) {
2305 if (theFormula.empty()) {
2306 return;
2307 }
2308
2309 for (unsigned long i=0UL; i<theFormula.countitems(); i++) {
2310 _Operation* this_op = ItemAt (i);
2311 if (this_op->theNumber) {
2312 continue;
2313 } else {
2314 if (this_op->theData>-1) {
2315 this_op->theData = variableIndex.get (this_op->theData);
2316 } else if (this_op->opCode == (long)MinusNumber) {
2317 this_op->opCode = HY_OP_CODE_SUB;
2318 } else {
2319 if (this_op->opCode == (long)FastMxAccess) {
2320 this_op->numberOfTerms = 2;
2321 } else {
2322 if (this_op->opCode == (long)FastMxWrite) {
2323 this_op->numberOfTerms = 2;
2324 }
2325 }
2326 this_op->opCode = simpleOperationCodes(simpleOperationFunctions.Find(this_op->opCode));
2327 }
2328 }
2329 }
2330 }
2331
2332 //__________________________________________________________________________________
ComputeSimple(_SimpleFormulaDatum * stack,_SimpleFormulaDatum * varValues)2333 hyFloat _Formula::ComputeSimple (_SimpleFormulaDatum* stack, _SimpleFormulaDatum* varValues) {
2334 if (theFormula.nonempty()) {
2335 long stackTop = 0;
2336 unsigned long upper_bound = NumberOperations();
2337
2338 for (unsigned long i=0UL; i<upper_bound; i++) {
2339
2340 if (simpleExpressionStatus[i] == -1L) {
2341 stack[stackTop++].value = ItemAt (i)->theNumber->Value();
2342 /*if (thisOp->theNumber) {
2343 stack[stackTop++].value = thisOp->theNumber->Value();
2344 continue;*/
2345 } else {
2346 if (simpleExpressionStatus[i] >= 0) {
2347 stack[stackTop++] = varValues[simpleExpressionStatus[i]];
2348 } else {
2349 if (simpleExpressionStatus[i] <= -10000L) {
2350 stackTop--;
2351 switch (-simpleExpressionStatus[i] - 10000L) {
2352 case HY_OP_CODE_ADD:
2353 stack[stackTop-1].value = stack[stackTop-1].value + stack[stackTop].value;
2354 break;
2355 case HY_OP_CODE_SUB:
2356 stack[stackTop-1].value = stack[stackTop-1].value - stack[stackTop].value;
2357 break;
2358 case HY_OP_CODE_MUL:
2359 stack[stackTop-1].value = stack[stackTop-1].value * stack[stackTop].value;
2360 break;
2361 case HY_OP_CODE_DIV:
2362 stack[stackTop-1].value = stack[stackTop-1].value / stack[stackTop].value;
2363 break;
2364 default:
2365 HandleApplicationError ("Internal error in _Formula::ComputeSimple - unsupported shortcut operation.)", true);
2366 return 0.0;
2367 }
2368 } else {
2369 _Operation const* thisOp = ItemAt (i);
2370 stackTop--;
2371 if (thisOp->numberOfTerms==2) {
2372 hyFloat (*theFunc) (hyFloat, hyFloat);
2373 theFunc = (hyFloat(*)(hyFloat,hyFloat))thisOp->opCode;
2374 if (stackTop<1L) {
2375 HandleApplicationError ("Internal error in _Formula::ComputeSimple - stack underflow.)", true);
2376 return 0.0;
2377 }
2378 stack[stackTop-1].value = (*theFunc)(stack[stackTop-1].value,stack[stackTop].value);
2379 } else {
2380 switch (thisOp->numberOfTerms) {
2381 case -2 : {
2382 hyFloat (*theFunc) (hyPointer,hyFloat);
2383 theFunc = (hyFloat(*)(hyPointer,hyFloat))thisOp->opCode;
2384 if (stackTop<1L) {
2385 HandleApplicationError ("Internal error in _Formula::ComputeSimple - stack underflow.)", true);
2386 return 0.0;
2387 }
2388 stack[stackTop-1].value = (*theFunc)(stack[stackTop-1].reference,stack[stackTop].value);
2389 break;
2390 }
2391 case -3 : {
2392 void (*theFunc) (hyPointer,hyFloat,hyFloat);
2393 theFunc = (void(*)(hyPointer,hyFloat,hyFloat))thisOp->opCode;
2394 if (stackTop != 2 || i != theFormula.lLength - 1) {
2395 HandleApplicationError ("Internal error in _Formula::ComputeSimple - stack underflow or MCoord command is not the last one.)", true);
2396
2397 return 0.0;
2398 }
2399 //stackTop = 0;
2400 // value, reference, index
2401 (*theFunc)(stack[1].reference,stack[2].value, stack[0].value);
2402 break;
2403 }
2404 default: {
2405 hyFloat (*theFunc) (hyFloat);
2406 theFunc = (hyFloat(*)(hyFloat))thisOp->opCode;
2407 stack[stackTop].value = (*theFunc)(stack[stackTop].value);
2408 ++stackTop;
2409 }
2410 }
2411
2412 }
2413 }
2414 }
2415 }
2416 }
2417 return stack[0].value;
2418 }
2419 return 0.0;
2420 }
2421
2422 //__________________________________________________________________________________
EqualFormula(_Formula * f)2423 bool _Formula::EqualFormula (_Formula* f) {
2424 if (theFormula.countitems() == f->theFormula.countitems()) {
2425 unsigned long const upper_bound = NumberOperations();
2426
2427 for (unsigned long i=0UL; i<upper_bound; i++) {
2428 if (ItemAt (i) -> EqualOp (f->ItemAt(i)) == false) {
2429 return false;
2430 }
2431 }
2432 return true;
2433 }
2434 return false;
2435 }
2436
2437 //__________________________________________________________________________________
ConstructPolynomial(void)2438 HBLObjectRef _Formula::ConstructPolynomial (void) {
2439 theStack.Reset();
2440 bool wellDone = true;
2441 _String errMsg;
2442
2443 for (unsigned long i=0UL; wellDone && i<theFormula.countitems(); i++) {
2444 wellDone = ItemAt (i)->ExecutePolynomial(theStack, nil, &errMsg);
2445 /*if (wellDone) {
2446 printf ("%ld (%ld)\n", i, theStack.StackDepth());
2447 for (long k = 0; k < theStack.StackDepth(); k++) {
2448 printf ("\t%s\n", _String((_String*)theStack.Peek (k)->toStr()).get_str());
2449 }
2450 }*/
2451 }
2452
2453 if (theStack.StackDepth() == 1 && wellDone) {
2454 return theStack.Pop(false);
2455 }
2456
2457 return nil;
2458 }
2459
2460 //__________________________________________________________________________________
HasChanged(bool ingoreCats)2461 bool _Formula::HasChanged (bool ingoreCats) {
2462 unsigned long const upper_bound = NumberOperations();
2463
2464 for (unsigned long i=0UL; i<upper_bound; i++) {
2465 long data_id;
2466 _Operation * this_op = ItemAt (i);
2467
2468 if (this_op->IsAVariable()) {
2469 data_id = this_op->GetAVariable();
2470 if (data_id>=0) {
2471 if (((_Variable*)(((BaseRef*)(variablePtrs.list_data))[data_id]))->HasChanged(ingoreCats)) {
2472 return true;
2473 }
2474 } else if (this_op->theNumber->HasChanged()) {
2475 return true;
2476 }
2477 } else if (this_op->opCode == HY_OP_CODE_BRANCHLENGTH||this_op->opCode == HY_OP_CODE_RANDOM||this_op->opCode == HY_OP_CODE_TIME)
2478 // time, random or branch length
2479 {
2480 return true;
2481 } else if (this_op->numberOfTerms<0L) {
2482 data_id = -this_op->numberOfTerms-2;
2483 if (IsBFFunctionIndexValid (data_id)) {
2484 if (GetBFFunctionType (data_id) == kBLFunctionSkipUpdate) {
2485 continue;
2486 }
2487 }
2488 return true;
2489 }
2490 }
2491 return false;
2492 }
2493
2494
2495 //__________________________________________________________________________________
2496
ScanFormulaForHBLFunctions(_AVLListX & collection,bool recursive,bool simplify)2497 void _Formula::ScanFormulaForHBLFunctions (_AVLListX& collection , bool recursive, bool simplify) {
2498
2499
2500 auto handle_function_id = [&collection, recursive] (const long hbl_id) -> void {
2501 if (IsBFFunctionIndexValid(hbl_id)) {
2502 _String function_name = GetBFFunctionNameByIndex(hbl_id);
2503
2504 if (collection.Find(&function_name) < 0) {
2505 collection.Insert (new _String (function_name), HY_BL_HBL_FUNCTION, false, false);
2506 if (recursive) {
2507 GetBFFunctionBody(hbl_id).BuildListOfDependancies(collection, true);
2508 }
2509 }
2510 }
2511 };
2512
2513 ConvertToTree(false);
2514
2515 if (theTree) {
2516
2517 if (simplify) {
2518 InternalSimplify(theTree);
2519 }
2520 node_iterator<long> ni (theTree, _HY_TREE_TRAVERSAL_PREORDER);
2521
2522 while (node<long>* iterator = ni.Next()) {
2523 _Operation *this_op = GetIthTerm(iterator->get_data());
2524
2525 long hbl_id = -1L;
2526
2527 if (this_op -> IsHBLFunctionCall()) {
2528 hbl_id = this_op -> GetHBLFunctionID();
2529 } else {
2530 if (this_op->opCode == HY_OP_CODE_CALL) {
2531 node <long>* function_to_call = iterator->go_down (1);
2532 _Operation * function_to_call_value = GetIthTerm (function_to_call->get_data());
2533
2534 if (function_to_call->get_num_nodes() == 0 && function_to_call_value -> IsConstantOfType(STRING)) {
2535 hbl_id = FindBFFunctionName (((_FString*)function_to_call_value->theNumber->Compute())->get_str());
2536 } else {
2537 ReportWarning ("Cannot export Call function arguments which are run-time dependent");
2538 }
2539 } else {
2540 if (this_op->opCode == HY_OP_CODE_MACCESS) { // handle AVL iterators
2541 if (this_op->GetNoTerms() == 3) { // [][]
2542 if (iterator->go_down(2)->get_num_nodes() == 0 && iterator->go_down(3)->get_num_nodes() == 0) {
2543 _Operation* bracket_1 = GetIthTerm (iterator->go_down(2)->get_data());
2544 _Operation* bracket_2 = GetIthTerm (iterator->go_down(3)->get_data());
2545 if (bracket_1->IsConstantOfType(STRING) && bracket_2->IsConstantOfType(STRING)) {
2546 handle_function_id (FindBFFunctionName (((_FString*)bracket_1->theNumber->Compute())->get_str()));
2547 handle_function_id (FindBFFunctionName (((_FString*)bracket_2->theNumber->Compute())->get_str()));
2548 continue;
2549 }
2550 }
2551 ReportWarning ("Potentially missed dependence on a function in [][]; arguments are run-time dependent");
2552 }
2553 }
2554 }
2555 }
2556
2557 handle_function_id (hbl_id);
2558
2559
2560
2561 }
2562 }
2563 }
2564
2565 //__________________________________________________________________________________
HasChangedSimple(_SimpleList & variableIndex)2566 bool _Formula::HasChangedSimple (_SimpleList& variableIndex) {
2567 unsigned long const upper_bound = NumberOperations();
2568
2569 for (unsigned long i=0UL; i<upper_bound; i++) {
2570 _Operation * this_op = ItemAt (i);
2571 if (this_op->theNumber) {
2572 continue;
2573 } else if (this_op->theData >= 0) {
2574 if (((_Variable*)(((BaseRef*)(variablePtrs.list_data))[variableIndex.list_data[this_op->theData]]))->HasChanged(false)) {
2575 return true;
2576 }
2577 } else {
2578 if (this_op->opCode == (long) RandomNumber) {
2579 return true;
2580 }
2581 }
2582 }
2583 return false;
2584 }
2585
2586 //__________________________________________________________________________________
ScanFForVariables(_AVLList & l,bool includeGlobals,bool includeAll,bool includeCategs,bool skipMatrixAssignments,_AVLListX * tagger,long weight) const2587 void _Formula::ScanFForVariables (_AVLList&l, bool includeGlobals, bool includeAll, bool includeCategs, bool skipMatrixAssignments, _AVLListX* tagger, long weight) const {
2588 unsigned long const upper_bound = NumberOperations();
2589
2590 for (unsigned long i=0UL; i<upper_bound; i++) {
2591 _Operation * this_op = ItemAt (i);
2592 if (this_op->IsAVariable()) {
2593 if (!includeGlobals)
2594 // This change was part of a commit that introduced an optimizer bug (suspected location:
2595 // src/core/batchlan2.cpp:2220). This change is suspicious as well (removed and undocumented condition).
2596 //if ((((_Variable*)LocateVar(theObj->GetAVariable()))->IsGlobal())||
2597 // (((_Variable*)LocateVar(theObj->GetIndex()))->ObjectClass()!=NUMBER)) {
2598 if (((_Variable*)LocateVar(this_op->GetAVariable()))->IsGlobal()) {
2599 continue;
2600 }
2601
2602 long f = this_op->GetAVariable();
2603
2604 if (f>=0) {
2605 _Variable * v = LocateVar(f);
2606
2607 if (v->IsCategory()&&includeCategs) {
2608 v->ScanForVariables (l,includeGlobals,tagger, weight);
2609 }
2610
2611 if(includeAll || v->ObjectClass()==NUMBER) {
2612 l.Insert ((BaseRef)f);
2613 if (tagger) {
2614 tagger -> UpdateValue((BaseRef)f, weight, 0);
2615 }
2616 }
2617
2618 if (skipMatrixAssignments) {
2619 if (v->ObjectClass()!=MATRIX || !this_op->AssignmentVariable()) {
2620 v->ScanForVariables(l,includeGlobals,tagger, weight);
2621 }
2622 } else if (!v->IsIndependent()) {
2623 v->ScanForVariables(l,includeGlobals,tagger);
2624 }
2625 } else if (this_op->theNumber)
2626 if (this_op->theNumber->ObjectClass()==MATRIX) {
2627 ((_Matrix*)this_op->theNumber)->ScanForVariables(l,includeGlobals,tagger, weight);
2628 }
2629 }
2630 }
2631 }
2632
2633 //__________________________________________________________________________________
ScanFForType(_SimpleList & l,int type)2634 void _Formula::ScanFForType (_SimpleList &l, int type) {
2635 unsigned long const upper_bound = NumberOperations();
2636
2637 for (unsigned long i=0UL; i<upper_bound; i++) {
2638 _Operation * this_op = ItemAt (i);
2639 if (this_op->IsAVariable()) {
2640 long f = this_op->GetAVariable();
2641
2642 if (f>=0) {
2643 _Variable * v = LocateVar(f);
2644
2645 if(v->ObjectClass()==type) {
2646 l << f;
2647 }
2648
2649 }
2650 }
2651 }
2652 }
2653
2654 //__________________________________________________________________________________
CheckFForDependence(long varID,bool checkAll)2655 bool _Formula::CheckFForDependence (long varID, bool checkAll) {
2656 unsigned long const upper_bound = NumberOperations();
2657
2658 for (unsigned long i=0UL; i<upper_bound; i++) {
2659 _Operation * this_op = ItemAt (i);
2660 if (this_op->IsAVariable()) {
2661 long f = this_op->GetAVariable();
2662 if (f>=0) {
2663 if (f == varID) {
2664 return true;
2665 }
2666 if (checkAll) {
2667 if (LocateVar(f)->CheckFForDependence(varID)) {
2668 return true;
2669 }
2670 }
2671 }
2672 }
2673 }
2674 return false;
2675 }
2676
2677 //__________________________________________________________________________________
CheckFForDependence(_AVLList const & indices,bool checkAll)2678 bool _Formula::CheckFForDependence (_AVLList const& indices, bool checkAll) {
2679 unsigned long const upper_bound = NumberOperations();
2680
2681 for (unsigned long i=0UL; i<upper_bound; i++) {
2682 _Operation * this_op = ItemAt (i);
2683 if (this_op->IsAVariable()) {
2684 long f = this_op->GetAVariable();
2685 if (f>=0) {
2686 if (indices.FindLong (f) >= 0) {
2687 return true;
2688 }
2689 if (checkAll) {
2690 if (LocateVar(f)->CheckFForDependence(indices)) {
2691 return true;
2692 }
2693 }
2694 }
2695 }
2696 }
2697 return false;
2698 }
2699
2700
2701 //__________________________________________________________________________________
LocalizeFormula(_Formula & ref,_String & parentName,_SimpleList & iv,_SimpleList & iiv,_SimpleList & dv,_SimpleList & idv)2702 void _Formula::LocalizeFormula (_Formula& ref, _String& parentName, _SimpleList& iv, _SimpleList& iiv, _SimpleList& dv, _SimpleList& idv) {
2703 unsigned long const upper_bound = ref.NumberOperations();
2704
2705 for (unsigned long i=0UL; i<upper_bound; i++) {
2706 _Operation * this_op = ref.ItemAt (i);
2707 if (this_op->IsAVariable()) {
2708 long vIndex = this_op->GetAVariable();
2709 _Variable* theV = LocateVar (vIndex);
2710 if (theV->IsGlobal()) {
2711 theFormula&& ref.theFormula(i);
2712 continue;
2713 }
2714 if (theV->IsContainer()) {
2715 continue;
2716 }
2717 _String localized_name = parentName&"."&*(theV->GetName());
2718
2719 _Variable * localized_var = CheckReceptacle (&localized_name, kEmptyString, false, false);
2720
2721 if (theV->IsIndependent()) {
2722 iv<<localized_var->get_index();
2723 iiv<<vIndex;
2724 } else {
2725 dv<<localized_var->get_index();
2726 idv<<vIndex;
2727
2728 }
2729 theFormula.AppendNewInstance( new _Operation (true, localized_name));
2730 } else {
2731 theFormula&& ref.theFormula(i);
2732 }
2733 }
2734 }
2735
2736 //__________________________________________________________________________________
DependsOnVariable(long idx)2737 bool _Formula::DependsOnVariable (long idx) {
2738 unsigned long const upper_bound = NumberOperations();
2739
2740 for (unsigned long i=0UL; i<upper_bound; i++) {
2741 _Operation * this_op = ItemAt (i);
2742 if (this_op->IsAVariable() && this_op->GetAVariable() == idx) {
2743 return true;
2744 }
2745 }
2746
2747 return false;
2748 }
2749
2750 //__________________________________________________________________________________
GetIthTerm(long idx) const2751 _Operation* _Formula::GetIthTerm (long idx) const {
2752 if (idx >= 0 && idx < theFormula.lLength) {
2753 return (_Operation*)theFormula.GetItem(idx);
2754 }
2755
2756 return nil;
2757 }
2758
2759 //__________________________________________________________________________________
ItemAt(long idx) const2760 _Operation* _Formula::ItemAt (long idx) const {
2761 return (_Operation*)theFormula.GetItem(idx);
2762 }
2763
2764 //__________________________________________________________________________________
IsAConstant(void)2765 bool _Formula::IsAConstant (void) {
2766
2767 unsigned long const upper_bound = NumberOperations();
2768
2769 for (unsigned long i=0UL; i<upper_bound; i++) {
2770 if ( ItemAt (i)->IsAVariable()) {
2771 return false;
2772 }
2773 }
2774 return true;
2775
2776 }
2777
2778 //__________________________________________________________________________________
IsConstant(bool strict)2779 bool _Formula::IsConstant (bool strict) {
2780 unsigned long const upper_bound = NumberOperations();
2781
2782 for (unsigned long i=0UL; i<upper_bound; i++) {
2783 if (ItemAt (i)->IsConstant(strict) == false) {
2784 return false;
2785 }
2786 }
2787
2788 return true;
2789 }
2790
2791 //__________________________________________________________________________________
PushTerm(BaseRef object)2792 void _Formula::PushTerm (BaseRef object) {
2793 _List * test_list;
2794 if ((test_list = dynamic_cast<_List*>(object))) {
2795 theFormula << *test_list;
2796 } else {
2797 theFormula << object;
2798 }
2799 }
2800
2801 //__________________________________________________________________________________
SimplifyConstants(void)2802 void _Formula::SimplifyConstants (void){
2803 ConvertToTree ();
2804 InternalSimplify (theTree);
2805 ConvertFromTree ();
2806 }
2807
2808 //__________________________________________________________________________________
GetTheMatrix(void)2809 HBLObjectRef _Formula::GetTheMatrix (void) {
2810 if (theFormula.countitems()==1) {
2811 _Operation* firstOp = ItemAt (0);
2812 HBLObjectRef ret = firstOp->GetANumber();
2813 if (ret&&(ret->ObjectClass()==MATRIX)) {
2814 return ret;
2815 } else {
2816 if (firstOp->theData!=-1) {
2817 _Variable* firstVar = LocateVar (firstOp->GetAVariable());
2818 ret = firstVar->GetValue();
2819 if (ret&&(ret->ObjectClass()==MATRIX)) {
2820 return ret;
2821 }
2822 }
2823 }
2824 }
2825 return nil;
2826 }
2827
2828 //__________________________________________________________________________________
ObjectClass(void)2829 long _Formula::ObjectClass (void) {
2830 if (theStack.StackDepth()) {
2831 return ((HBLObjectRef)theStack.theStack.list_data[0])->ObjectClass();
2832 }
2833
2834 HBLObjectRef res = Compute();
2835
2836 if (res) {
2837 return res->ObjectClass();
2838 }
2839
2840 return HY_UNDEFINED;
2841 }
2842
2843
2844
2845 //__________________________________________________________________________________
ParseFormula(_String const & s,_VariableContainer const * theParent,_String * reportErrors)2846 long _Formula::ParseFormula (_String const &s, _VariableContainer const* theParent, _String* reportErrors) {
2847 theTree = nil;
2848 resultCache = nil;
2849 recursion_calls = nil;
2850 call_count = 0UL;
2851
2852 _FormulaParsingContext fpc (reportErrors, theParent);
2853
2854 _String formula_copy (s);
2855
2856 long return_value = Parse (this, formula_copy, fpc, nil);
2857
2858 if (return_value != HY_FORMULA_EXPRESSION) {
2859 Clear();
2860 }
2861
2862 return return_value;
2863
2864 }
2865
2866 //__________________________________________________________________________________
ConvertToTree(bool err_msg)2867 void _Formula::ConvertToTree (bool err_msg) {
2868 if (!theTree && theFormula.countitems()) { // work to do
2869 _SimpleList nodeStack;
2870
2871 unsigned long const upper_bound = NumberOperations();
2872 theTree = nil;
2873 try {
2874
2875 for (unsigned long i=0UL; i<upper_bound; i++) {
2876
2877 _Operation* currentOp = ItemAt(i);
2878
2879 if (currentOp->theNumber || currentOp->theData >= 0L || currentOp->theData <= -2L) { // a data bit
2880 node<long>* leafNode = new node<long>;
2881 leafNode->init(i);
2882 nodeStack<<(long)leafNode;
2883 } else { // an operation
2884 long nTerms = currentOp->GetNoTerms();
2885 if (nTerms<0L) {
2886 nTerms = GetBFFunctionArgumentCount(currentOp->opCode);
2887 }
2888
2889 if (nTerms>nodeStack.lLength) {
2890 throw (_String ("Insufficient number of arguments for a call to ") & _String ((_String*)currentOp->toStr()) & " while converting " & toRPN(kFormulaStringConversionNormal).Enquote() & " to a parse tree");
2891 }
2892
2893 node<long>* operationNode = new node<long>;
2894 operationNode->init(i);
2895 for (long j=0; j<nTerms; j++) {
2896 operationNode->prepend_node(*((node<long>*)nodeStack.Pop()));
2897 }
2898 nodeStack<<(long)operationNode;
2899 }
2900 }
2901 if (nodeStack.lLength!=1) {
2902 throw ((_String)"The expression '" & toRPN (kFormulaStringConversionNormal) & "' has " & (long)nodeStack.lLength & " terms left on the stack after evaluation");
2903 } else {
2904 theTree = (node<long>*)nodeStack(0);
2905 }
2906 } catch (_String const& e) {
2907 if (err_msg) {
2908 HandleApplicationError(e);
2909 }
2910 nodeStack.Each ([this] (long v, unsigned long) -> void {
2911 node<long>* operationNode = (node<long>*)v;
2912 operationNode->delete_tree();
2913 delete (operationNode);
2914 });
2915 theTree = nil;
2916 }
2917 }
2918 }
2919 //__________________________________________________________________________________
ConvertFromTree(void)2920 void _Formula::ConvertFromTree (void) {
2921 if (theTree) { // work to do
2922 _SimpleList termOrder;
2923 node_iterator<long> ni (theTree, _HY_TREE_TRAVERSAL_POSTORDER);
2924
2925 while (node<long>* iterator = ni.Next()) {
2926 termOrder<<iterator->get_data();
2927 }
2928
2929 if (termOrder.lLength!=theFormula.lLength) { // something has changed
2930 _List newFormula;
2931 for (long i=0; i<termOrder.lLength; i++) {
2932 newFormula<<theFormula(termOrder(i));
2933 }
2934 theFormula.Clear();
2935 theFormula.Duplicate(&newFormula);
2936 //ConvertToTree();
2937 }
2938 theTree->delete_tree();
2939 delete (theTree);
2940 theTree = nil;
2941 }
2942 }
2943
2944