DataTree.cc 19.3 KB
Newer Older
1
/*
2
 * Copyright (C) 2003-2018 Dynare Team
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * This file is part of Dynare.
 *
 * Dynare is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Dynare is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 */

20
#include <cstdlib>
21
#include <cassert>
22
#include <iostream>
23
24
25
#include <regex>

#include <boost/filesystem.hpp>
26
27
28

#include "DataTree.hh"

29
30
void
DataTree::initConstants()
31
{
32
33
34
  Zero = AddNonNegativeConstant("0");
  One = AddNonNegativeConstant("1");
  Two = AddNonNegativeConstant("2");
35
36

  MinusOne = AddUMinus(One);
37

38
39
  NaN = AddNonNegativeConstant("NaN");
  Infinity = AddNonNegativeConstant("Inf");
40
  MinusInfinity = AddUMinus(Infinity);
sebastien's avatar
sebastien committed
41

42
  Pi = AddNonNegativeConstant("3.141592653589793");
43
44
}

45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
DataTree::DataTree(SymbolTable &symbol_table_arg,
                   NumericalConstants &num_constants_arg,
                   ExternalFunctionsTable &external_functions_table_arg,
                   bool is_dynamic_arg) :
  symbol_table{symbol_table_arg},
  num_constants{num_constants_arg},
  external_functions_table{external_functions_table_arg},
  is_dynamic {is_dynamic_arg}
{
  initConstants();
}

DataTree::DataTree(const DataTree &d) :
  symbol_table {d.symbol_table},
  num_constants {d.num_constants},
  external_functions_table {d.external_functions_table},
  is_dynamic {d.is_dynamic},
  local_variables_vector {d.local_variables_vector}
{
  // Constants must be initialized first because they are used in some Add* methods
  initConstants();

  for (const auto & it : d.node_list)
68
    it->clone(*this);
69
70
71
72

  assert(node_list.size() == d.node_list.size());

  for (const auto & it : d.local_variables_table)
73
    local_variables_table[it.first] = it.second->clone(*this);
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
}

DataTree &
DataTree::operator=(const DataTree &d)
{
  assert (&symbol_table == &d.symbol_table);
  assert (&num_constants == &d.num_constants);
  assert (&external_functions_table == &d.external_functions_table);
  assert (is_dynamic == d.is_dynamic);

  num_const_node_map.clear();
  variable_node_map.clear();
  unary_op_node_map.clear();
  binary_op_node_map.clear();
  trinary_op_node_map.clear();
  external_function_node_map.clear();
  var_expectation_node_map.clear();
  pac_expectation_node_map.clear();
  first_deriv_external_function_node_map.clear();
  second_deriv_external_function_node_map.clear();

  node_list.clear();

  // Constants must be initialized first because they are used in some Add* methods
  initConstants();

  for (const auto & it : d.node_list)
101
    it->clone(*this);
102
103
104
105
106
107

  assert(node_list.size() == d.node_list.size());

  local_variables_vector = d.local_variables_vector;

  for (const auto & it : d.local_variables_table)
108
    local_variables_table[it.first] = it.second->clone(*this);
109
110
111
112

  return *this;
}

113
expr_t
114
DataTree::AddNonNegativeConstant(const string &value)
115
{
116
  int id = num_constants.AddNonNegativeConstant(value);
117

118
  auto it = num_const_node_map.find(id);
119
120
  if (it != num_const_node_map.end())
    return it->second;
121
122
123
124
125
126

  auto sp = make_unique<NumConstNode>(*this, node_list.size(), id);
  auto p = sp.get();
  node_list.push_back(move(sp));
  num_const_node_map[id] = p;
  return p;
127
128
}

sebastien's avatar
sebastien committed
129
VariableNode *
130
DataTree::AddVariable(int symb_id, int lag)
131
{
132
133
134
135
136
137
  if (lag != 0 && !is_dynamic)
    {
      cerr << "Leads/lags not authorized in this DataTree" << endl;
      exit(EXIT_FAILURE);
    }

138
  auto it = variable_node_map.find({ symb_id, lag });
139
140
  if (it != variable_node_map.end())
    return it->second;
141
142
143
144
145
146

  auto sp = make_unique<VariableNode>(*this, node_list.size(), symb_id, lag);
  auto p = sp.get();
  node_list.push_back(move(sp));
  variable_node_map[{ symb_id, lag }] = p;
  return p;
147
}
sebastien's avatar
sebastien committed
148

149
150
151
152
153
154
155
156
157
158
159
160
VariableNode *
DataTree::getVariable(int symb_id, int lag) const
{
  auto it = variable_node_map.find({ symb_id, lag });
  if (it == variable_node_map.end())
    {
      cerr << "DataTree::getVariable: unknown variable node for symb_id=" << symb_id << " and lag=" << lag << endl;
      exit(EXIT_FAILURE);
    }
  return it->second;
}

161
162
163
bool
DataTree::ParamUsedWithLeadLagInternal() const
{
164
  for (const auto & it : variable_node_map)
165
    if (symbol_table.getType(it.first.first) == SymbolType::parameter && it.first.second != 0)
166
167
168
169
      return true;
  return false;
}

170
171
expr_t
DataTree::AddPlus(expr_t iArg1, expr_t iArg2)
172
173
174
{
  if (iArg1 != Zero && iArg2 != Zero)
    {
175
      // Simplify x+(-y) in x-y
176
      auto uarg2 = dynamic_cast<UnaryOpNode *>(iArg2);
177
178
      if (uarg2 != nullptr && uarg2->op_code == UnaryOpcode::uminus)
        return AddMinus(iArg1, uarg2->arg);
179

180
181
      // Simplify (-x)+y in y-x
      auto uarg1 = dynamic_cast<UnaryOpNode *>(iArg1);
182
183
      if (uarg1 != nullptr && uarg1->op_code == UnaryOpcode::uminus)
        return AddMinus(iArg2, uarg1->arg);
184
185
186
187
188
189

      // To treat commutativity of "+"
      // Nodes iArg1 and iArg2 are sorted by index
      if (iArg1->idx > iArg2->idx)
        swap(iArg1, iArg2);

190
      return AddBinaryOp(iArg1, BinaryOpcode::plus, iArg2);
191
192
193
194
195
196
197
198
199
    }
  else if (iArg1 != Zero)
    return iArg1;
  else if (iArg2 != Zero)
    return iArg2;
  else
    return Zero;
}

200
201
expr_t
DataTree::AddMinus(expr_t iArg1, expr_t iArg2)
202
{
203
  if (iArg2 == Zero)
204
    return iArg1;
205
206

  if (iArg1 == Zero)
207
    return AddUMinus(iArg2);
208
209

  if (iArg1 == iArg2)
210
    return Zero;
211

212
  return AddBinaryOp(iArg1, BinaryOpcode::minus, iArg2);
213
214
}

215
216
expr_t
DataTree::AddUMinus(expr_t iArg1)
217
218
219
220
{
  if (iArg1 != Zero)
    {
      // Simplify -(-x) in x
221
      auto *uarg = dynamic_cast<UnaryOpNode *>(iArg1);
222
223
      if (uarg != nullptr && uarg->op_code == UnaryOpcode::uminus)
        return uarg->arg;
224

225
      return AddUnaryOp(UnaryOpcode::uminus, iArg1);
226
227
228
229
230
    }
  else
    return Zero;
}

231
232
expr_t
DataTree::AddTimes(expr_t iArg1, expr_t iArg2)
233
234
235
236
237
238
239
240
241
242
{
  if (iArg1 == MinusOne)
    return AddUMinus(iArg2);
  else if (iArg2 == MinusOne)
    return AddUMinus(iArg1);
  else if (iArg1 != Zero && iArg1 != One && iArg2 != Zero && iArg2 != One)
    {
      // To treat commutativity of "*"
      // Nodes iArg1 and iArg2 are sorted by index
      if (iArg1->idx > iArg2->idx)
243
        swap(iArg1, iArg2);
244
      return AddBinaryOp(iArg1, BinaryOpcode::times, iArg2);
245
246
247
248
249
250
251
252
253
254
255
    }
  else if (iArg1 != Zero && iArg1 != One && iArg2 == One)
    return iArg1;
  else if (iArg2 != Zero && iArg2 != One && iArg1 == One)
    return iArg2;
  else if (iArg2 == One && iArg1 == One)
    return One;
  else
    return Zero;
}

256
expr_t
257
DataTree::AddDivide(expr_t iArg1, expr_t iArg2) noexcept(false)
258
{
259
  if (iArg2 == One)
260
    return iArg1;
261
262
263

  // This test should be before the next two, otherwise 0/0 won't be rejected
  if (iArg2 == Zero)
264
    {
265
      cerr << "ERROR: Division by zero!" << endl;
266
      throw DivisionByZeroException();
267
    }
268
269
270
271
272
273
274

  if (iArg1 == Zero)
    return Zero;

  if (iArg1 == iArg2)
    return One;

275
  return AddBinaryOp(iArg1, BinaryOpcode::divide, iArg2);
276
277
}

278
279
expr_t
DataTree::AddLess(expr_t iArg1, expr_t iArg2)
280
{
281
  return AddBinaryOp(iArg1, BinaryOpcode::less, iArg2);
282
283
}

284
285
expr_t
DataTree::AddGreater(expr_t iArg1, expr_t iArg2)
286
{
287
  return AddBinaryOp(iArg1, BinaryOpcode::greater, iArg2);
288
289
}

290
291
expr_t
DataTree::AddLessEqual(expr_t iArg1, expr_t iArg2)
292
{
293
  return AddBinaryOp(iArg1, BinaryOpcode::lessEqual, iArg2);
294
295
}

296
297
expr_t
DataTree::AddGreaterEqual(expr_t iArg1, expr_t iArg2)
298
{
299
  return AddBinaryOp(iArg1, BinaryOpcode::greaterEqual, iArg2);
300
301
}

302
303
expr_t
DataTree::AddEqualEqual(expr_t iArg1, expr_t iArg2)
304
{
305
  return AddBinaryOp(iArg1, BinaryOpcode::equalEqual, iArg2);
306
307
}

308
309
expr_t
DataTree::AddDifferent(expr_t iArg1, expr_t iArg2)
310
{
311
  return AddBinaryOp(iArg1, BinaryOpcode::different, iArg2);
312
313
}

314
315
expr_t
DataTree::AddPower(expr_t iArg1, expr_t iArg2)
316
{
317
  if (iArg1 != Zero && iArg2 != Zero && iArg1 != One && iArg2 != One)
318
    return AddBinaryOp(iArg1, BinaryOpcode::power, iArg2);
319
320
  else if (iArg1 == One)
    return One;
321
322
323
324
325
326
327
328
  else if (iArg2 == One)
    return iArg1;
  else if (iArg2 == Zero)
    return One;
  else
    return Zero;
}

329
330
331
332
expr_t
DataTree::AddPowerDeriv(expr_t iArg1, expr_t iArg2, int powerDerivOrder)
{
  assert(powerDerivOrder > 0);
333
  return AddBinaryOp(iArg1, BinaryOpcode::powerDeriv, iArg2, powerDerivOrder);
334
335
}

336
337
338
expr_t
DataTree::AddDiff(expr_t iArg1)
{
339
  return AddUnaryOp(UnaryOpcode::diff, iArg1);
340
341
}

342
343
344
expr_t
DataTree::AddAdl(expr_t iArg1, const string &name, const vector<int> &lags)
{
345
  return AddUnaryOp(UnaryOpcode::adl, iArg1, 0, 0, 0, string(name), lags);
346
347
}

348
349
expr_t
DataTree::AddExp(expr_t iArg1)
350
351
{
  if (iArg1 != Zero)
352
    return AddUnaryOp(UnaryOpcode::exp, iArg1);
353
354
355
356
  else
    return One;
}

357
358
expr_t
DataTree::AddLog(expr_t iArg1)
359
360
{
  if (iArg1 != Zero && iArg1 != One)
361
    return AddUnaryOp(UnaryOpcode::log, iArg1);
362
363
364
365
  else if (iArg1 == One)
    return Zero;
  else
    {
366
      cerr << "ERROR: log(0) not defined!" << endl;
367
      exit(EXIT_FAILURE);
368
369
370
    }
}

371
372
expr_t
DataTree::AddLog10(expr_t iArg1)
373
374
{
  if (iArg1 != Zero && iArg1 != One)
375
    return AddUnaryOp(UnaryOpcode::log10, iArg1);
376
377
378
379
  else if (iArg1 == One)
    return Zero;
  else
    {
380
      cerr << "ERROR: log10(0) not defined!" << endl;
381
      exit(EXIT_FAILURE);
382
383
384
    }
}

385
386
expr_t
DataTree::AddCos(expr_t iArg1)
387
388
{
  if (iArg1 != Zero)
389
    return AddUnaryOp(UnaryOpcode::cos, iArg1);
390
391
392
393
  else
    return One;
}

394
395
expr_t
DataTree::AddSin(expr_t iArg1)
396
397
{
  if (iArg1 != Zero)
398
    return AddUnaryOp(UnaryOpcode::sin, iArg1);
399
400
401
402
  else
    return Zero;
}

403
404
expr_t
DataTree::AddTan(expr_t iArg1)
405
406
{
  if (iArg1 != Zero)
407
    return AddUnaryOp(UnaryOpcode::tan, iArg1);
408
409
410
411
  else
    return Zero;
}

412
413
expr_t
DataTree::AddAcos(expr_t iArg1)
414
415
{
  if (iArg1 != One)
416
    return AddUnaryOp(UnaryOpcode::acos, iArg1);
417
418
419
420
  else
    return Zero;
}

421
422
expr_t
DataTree::AddAsin(expr_t iArg1)
423
424
{
  if (iArg1 != Zero)
425
    return AddUnaryOp(UnaryOpcode::asin, iArg1);
426
427
428
429
  else
    return Zero;
}

430
431
expr_t
DataTree::AddAtan(expr_t iArg1)
432
433
{
  if (iArg1 != Zero)
434
    return AddUnaryOp(UnaryOpcode::atan, iArg1);
435
436
437
438
  else
    return Zero;
}

439
440
expr_t
DataTree::AddCosh(expr_t iArg1)
441
442
{
  if (iArg1 != Zero)
443
    return AddUnaryOp(UnaryOpcode::cosh, iArg1);
444
445
446
447
  else
    return One;
}

448
449
expr_t
DataTree::AddSinh(expr_t iArg1)
450
451
{
  if (iArg1 != Zero)
452
    return AddUnaryOp(UnaryOpcode::sinh, iArg1);
453
454
455
456
  else
    return Zero;
}

457
458
expr_t
DataTree::AddTanh(expr_t iArg1)
459
460
{
  if (iArg1 != Zero)
461
    return AddUnaryOp(UnaryOpcode::tanh, iArg1);
462
463
464
465
  else
    return Zero;
}

466
467
expr_t
DataTree::AddAcosh(expr_t iArg1)
468
469
{
  if (iArg1 != One)
470
    return AddUnaryOp(UnaryOpcode::acosh, iArg1);
471
472
473
474
  else
    return Zero;
}

475
476
expr_t
DataTree::AddAsinh(expr_t iArg1)
477
478
{
  if (iArg1 != Zero)
479
    return AddUnaryOp(UnaryOpcode::asinh, iArg1);
480
481
482
483
  else
    return Zero;
}

484
485
expr_t
DataTree::AddAtanh(expr_t iArg1)
486
487
{
  if (iArg1 != Zero)
488
    return AddUnaryOp(UnaryOpcode::atanh, iArg1);
489
490
491
492
  else
    return Zero;
}

493
494
expr_t
DataTree::AddSqrt(expr_t iArg1)
495
496
{
  if (iArg1 != Zero)
497
    return AddUnaryOp(UnaryOpcode::sqrt, iArg1);
498
499
500
501
  else
    return Zero;
}

502
503
504
505
506
507
508
509
expr_t
DataTree::AddAbs(expr_t iArg1)
{
  if (iArg1 == Zero)
    return Zero;
  if (iArg1 == One)
    return One;
  else
510
    return AddUnaryOp(UnaryOpcode::abs, iArg1);
511
512
513
514
515
516
517
518
519
520
}

expr_t
DataTree::AddSign(expr_t iArg1)
{
  if (iArg1 == Zero)
    return Zero;
  if (iArg1 == One)
    return One;
  else
521
    return AddUnaryOp(UnaryOpcode::sign, iArg1);
522
523
}

524
525
expr_t
DataTree::AddErf(expr_t iArg1)
526
527
{
  if (iArg1 != Zero)
528
    return AddUnaryOp(UnaryOpcode::erf, iArg1);
529
530
531
532
  else
    return Zero;
}

533
534
expr_t
DataTree::AddMax(expr_t iArg1, expr_t iArg2)
535
{
536
  return AddBinaryOp(iArg1, BinaryOpcode::max, iArg2);
537
538
}

539
540
expr_t
DataTree::AddMin(expr_t iArg1, expr_t iArg2)
541
{
542
  return AddBinaryOp(iArg1, BinaryOpcode::min, iArg2);
543
544
}

545
546
expr_t
DataTree::AddNormcdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
547
{
548
  return AddTrinaryOp(iArg1, TrinaryOpcode::normcdf, iArg2, iArg3);
549
550
}

551
552
expr_t
DataTree::AddNormpdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
553
{
554
  return AddTrinaryOp(iArg1, TrinaryOpcode::normpdf, iArg2, iArg3);
555
556
}

557
558
expr_t
DataTree::AddSteadyState(expr_t iArg1)
559
{
560
  return AddUnaryOp(UnaryOpcode::steadyState, iArg1);
561
562
}

563
564
565
expr_t
DataTree::AddSteadyStateParamDeriv(expr_t iArg1, int param_symb_id)
{
566
  return AddUnaryOp(UnaryOpcode::steadyStateParamDeriv, iArg1, 0, param_symb_id);
567
568
569
570
571
}

expr_t
DataTree::AddSteadyStateParam2ndDeriv(expr_t iArg1, int param1_symb_id, int param2_symb_id)
{
572
  return AddUnaryOp(UnaryOpcode::steadyStateParam2ndDeriv, iArg1, 0, param1_symb_id, param2_symb_id);
573
574
}

575
576
expr_t
DataTree::AddExpectation(int iArg1, expr_t iArg2)
577
{
578
  return AddUnaryOp(UnaryOpcode::expectation, iArg2, iArg1);
579
580
}

Houtan Bastani's avatar
Houtan Bastani committed
581
expr_t
582
DataTree::AddVarExpectation(const string &model_name)
Houtan Bastani's avatar
Houtan Bastani committed
583
{
584
  auto it = var_expectation_node_map.find(model_name);
585
586
587
  if (it != var_expectation_node_map.end())
    return it->second;

588
589
590
591
592
  auto sp = make_unique<VarExpectationNode>(*this, node_list.size(), model_name);
  auto p = sp.get();
  node_list.push_back(move(sp));
  var_expectation_node_map[model_name] = p;
  return p;
Houtan Bastani's avatar
Houtan Bastani committed
593
594
}

595
expr_t
Houtan Bastani's avatar
Houtan Bastani committed
596
DataTree::AddPacExpectation(const string &model_name)
597
{
598
  auto it = pac_expectation_node_map.find(model_name);
599
600
  if (it != pac_expectation_node_map.end())
    return it->second;
601

602
603
604
605
606
  auto sp = make_unique<PacExpectationNode>(*this, node_list.size(), model_name);
  auto p = sp.get();
  node_list.push_back(move(sp));
  pac_expectation_node_map[model_name] = p;
  return p;
607
608
}

609
610
expr_t
DataTree::AddEqual(expr_t iArg1, expr_t iArg2)
611
{
612
  return AddBinaryOp(iArg1, BinaryOpcode::equal, iArg2);
613
614
615
}

void
616
DataTree::AddLocalVariable(int symb_id, expr_t value) noexcept(false)
617
{
618
  assert(symbol_table.getType(symb_id) == SymbolType::modelLocalVariable);
sebastien's avatar
sebastien committed
619

620
  // Throw an exception if symbol already declared
621
  auto it = local_variables_table.find(symb_id);
622
  if (it != local_variables_table.end())
623
    throw LocalVariableException(symbol_table.getName(symb_id));
624

625
  local_variables_table[symb_id] = value;
626
  local_variables_vector.push_back(symb_id);
627
628
}

629
630
expr_t
DataTree::AddExternalFunction(int symb_id, const vector<expr_t> &arguments)
631
{
632
  assert(symbol_table.getType(symb_id) == SymbolType::externalFunction);
633

634
  auto it = external_function_node_map.find({ arguments, symb_id });
635
636
637
  if (it != external_function_node_map.end())
    return it->second;

638
639
640
641
642
  auto sp = make_unique<ExternalFunctionNode>(*this, node_list.size(), symb_id, arguments);
  auto p = sp.get();
  node_list.push_back(move(sp));
  external_function_node_map[{ arguments, symb_id }] = p;
  return p;
643
}
sebastien's avatar
sebastien committed
644

645
expr_t
646
DataTree::AddFirstDerivExternalFunction(int top_level_symb_id, const vector<expr_t> &arguments, int input_index)
647
{
648
  assert(symbol_table.getType(top_level_symb_id) == SymbolType::externalFunction);
649

650
  auto it
651
    = first_deriv_external_function_node_map.find({ arguments, input_index, top_level_symb_id });
652
653
654
  if (it != first_deriv_external_function_node_map.end())
    return it->second;

655
656
657
658
659
  auto sp = make_unique<FirstDerivExternalFunctionNode>(*this, node_list.size(), top_level_symb_id, arguments, input_index);
  auto p = sp.get();
  node_list.push_back(move(sp));
  first_deriv_external_function_node_map[{ arguments, input_index, top_level_symb_id }] = p;
  return p;
660
}
661

662
expr_t
663
DataTree::AddSecondDerivExternalFunction(int top_level_symb_id, const vector<expr_t> &arguments, int input_index1, int input_index2)
664
{
665
  assert(symbol_table.getType(top_level_symb_id) == SymbolType::externalFunction);
666

667
  auto it
668
    = second_deriv_external_function_node_map.find({ arguments, input_index1, input_index2,
669
          top_level_symb_id });
670
671
672
  if (it != second_deriv_external_function_node_map.end())
    return it->second;

673
674
675
676
677
  auto sp = make_unique<SecondDerivExternalFunctionNode>(*this, node_list.size(), top_level_symb_id, arguments, input_index1, input_index2);
  auto p = sp.get();
  node_list.push_back(move(sp));
  second_deriv_external_function_node_map[{ arguments, input_index1, input_index2, top_level_symb_id }] = p;
  return p;
678
}
sebastien's avatar
sebastien committed
679

680
681
682
bool
DataTree::isSymbolUsed(int symb_id) const
{
683
684
  for (const auto & it : variable_node_map)
    if (it.first.first == symb_id)
685
686
687
688
689
690
691
692
      return true;

  if (local_variables_table.find(symb_id) != local_variables_table.end())
    return true;

  return false;
}

693
int
694
DataTree::getDerivID(int symb_id, int lag) const noexcept(false)
695
696
697
698
{
  throw UnknownDerivIDException();
}

699
SymbolType
700
DataTree::getTypeByDerivID(int deriv_id) const noexcept(false)
701
702
703
704
705
{
  throw UnknownDerivIDException();
}

int
706
DataTree::getLagByDerivID(int deriv_id) const noexcept(false)
707
708
709
710
711
{
  throw UnknownDerivIDException();
}

int
712
DataTree::getSymbIDByDerivID(int deriv_id) const noexcept(false)
713
714
715
716
717
718
719
720
721
{
  throw UnknownDerivIDException();
}

void
DataTree::addAllParamDerivId(set<int> &deriv_id_set)
{
}

722
int
723
DataTree::getDynJacobianCol(int deriv_id) const noexcept(false)
724
725
726
{
  throw UnknownDerivIDException();
}
sebastien's avatar
sebastien committed
727
728

bool
Sébastien Villemot's avatar
Sébastien Villemot committed
729
DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
sebastien's avatar
sebastien committed
730
{
731
  for (const auto & it : unary_op_node_map)
732
    if (get<1>(it.first) == opcode)
Sébastien Villemot's avatar
Sébastien Villemot committed
733
734
735
736
737
738
739
740
      return true;

  return false;
}

bool
DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
{
741
  for (const auto & it : binary_op_node_map)
742
    if (get<2>(it.first) == opcode)
Sébastien Villemot's avatar
Sébastien Villemot committed
743
744
745
746
747
748
749
750
      return true;

  return false;
}

bool
DataTree::isTrinaryOpUsed(TrinaryOpcode opcode) const
{
751
  for (const auto & it : trinary_op_node_map)
752
    if (get<3>(it.first) == opcode)
sebastien's avatar
sebastien committed
753
754
755
756
      return true;

  return false;
}
757
758
759
760

bool
DataTree::isExternalFunctionUsed(int symb_id) const
{
761
762
  for (const auto & it : external_function_node_map)
    if (it.first.second == symb_id)
763
764
765
766
767
768
769
770
      return true;

  return false;
}

bool
DataTree::isFirstDerivExternalFunctionUsed(int symb_id) const
{
771
  for (const auto & it : first_deriv_external_function_node_map)
772
    if (get<2>(it.first) == symb_id)
773
774
775
776
777
778
779
780
      return true;

  return false;
}

bool
DataTree::isSecondDerivExternalFunctionUsed(int symb_id) const
{
781
  for (const auto & it : second_deriv_external_function_node_map)
782
    if (get<3>(it.first) == symb_id)
783
784
785
786
      return true;

  return false;
}
Sébastien Villemot's avatar
Sébastien Villemot committed
787
788
789
790
791

int
DataTree::minLagForSymbol(int symb_id) const
{
  int r = 0;
792
793
794
  for (const auto & it : variable_node_map)
    if (it.first.first == symb_id && it.first.second < r)
      r = it.first.second;
Sébastien Villemot's avatar
Sébastien Villemot committed
795
796
  return r;
}
797
798
799
800

void
DataTree::writePowerDerivCHeader(ostream &output) const
{
801
  if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
802
803
804
805
    output << "double getPowerDeriv(double, double, int);" << endl;
}

void
806
DataTree::writePowerDeriv(ostream &output) const
807
{
808
  if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
809
810
811
812
813
    output << "/*" << endl
           << " * The k-th derivative of x^p" << endl
           << " */" << endl
           << "double getPowerDeriv(double x, double p, int k)" << endl
           << "{" << endl
814
815
816
           << "#ifdef _MSC_VER" << endl
           << "# define nearbyint(x) (fabs((x)-floor(x)) < fabs((x)-ceil(x)) ? floor(x) : ceil(x))" << endl
           << "#endif" << endl
817
           << "  if ( fabs(x) < " << near_zero << " && p > 0 && k > p && fabs(p-nearbyint(p)) < " << near_zero << " )" << endl
818
819
820
821
822
823
824
825
826
827
828
           << "    return 0.0;" << endl
           << "  else" << endl
           << "    {" << endl
           << "      int i = 0;" << endl
           << "      double dxp = pow(x, p-k);" << endl
           << "      for (; i<k; i++)" << endl
           << "        dxp *= p--;" << endl
           << "      return dxp;" << endl
           << "    }" << endl
           << "}" << endl;
}
829

830
831
832
833
834
835
836
837
string
DataTree::packageDir(const string &package)
{
  regex pat{"\\."};
  string dirname = "+" + regex_replace(package, pat, "/+");
  boost::filesystem::create_directories(dirname);
  return dirname;
}