DataTree.cc 20.7 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
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();
}

57
DataTree::~DataTree() = default;
58

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
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)
    it->cloneDynamic(*this);

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

  for (const auto & it : d.local_variables_table)
    local_variables_table[it.first] = it.second->cloneDynamic(*this);
}

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)
    it->cloneDynamic(*this);

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

  local_variables_vector = d.local_variables_vector;

  for (const auto & it : d.local_variables_table)
    local_variables_table[it.first] = it.second->cloneDynamic(*this);

  return *this;
}

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

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

  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;
129
130
}

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

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

  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;
149
}
sebastien's avatar
sebastien committed
150

151
152
153
bool
DataTree::ParamUsedWithLeadLagInternal() const
{
154
  for (const auto & it : variable_node_map)
155
    if (symbol_table.getType(it.first.first) == SymbolType::parameter && it.first.second != 0)
156
157
158
159
      return true;
  return false;
}

160
161
expr_t
DataTree::AddPlus(expr_t iArg1, expr_t iArg2)
162
163
164
165
{
  if (iArg1 != Zero && iArg2 != Zero)
    {
      // Simplify x+(-y) in x-y
166
      auto *uarg2 = dynamic_cast<UnaryOpNode *>(iArg2);
167
      if (uarg2 != nullptr && uarg2->get_op_code() == UnaryOpcode::uminus)
sebastien's avatar
sebastien committed
168
        return AddMinus(iArg1, uarg2->get_arg());
169
170
171
172
173

      // To treat commutativity of "+"
      // Nodes iArg1 and iArg2 are sorted by index
      if (iArg1->idx > iArg2->idx)
        {
174
          expr_t tmp = iArg1;
175
176
177
          iArg1 = iArg2;
          iArg2 = tmp;
        }
178
      return AddBinaryOp(iArg1, BinaryOpcode::plus, iArg2);
179
180
181
182
183
184
185
186
187
    }
  else if (iArg1 != Zero)
    return iArg1;
  else if (iArg2 != Zero)
    return iArg2;
  else
    return Zero;
}

188
189
expr_t
DataTree::AddMinus(expr_t iArg1, expr_t iArg2)
190
{
191
  if (iArg2 == Zero)
192
    return iArg1;
193
194

  if (iArg1 == Zero)
195
    return AddUMinus(iArg2);
196
197

  if (iArg1 == iArg2)
198
    return Zero;
199

200
  return AddBinaryOp(iArg1, BinaryOpcode::minus, iArg2);
201
202
}

203
204
expr_t
DataTree::AddUMinus(expr_t iArg1)
205
206
207
208
{
  if (iArg1 != Zero)
    {
      // Simplify -(-x) in x
209
      auto *uarg = dynamic_cast<UnaryOpNode *>(iArg1);
210
      if (uarg != nullptr && uarg->get_op_code() == UnaryOpcode::uminus)
sebastien's avatar
sebastien committed
211
        return uarg->get_arg();
212

213
      return AddUnaryOp(UnaryOpcode::uminus, iArg1);
214
215
216
217
218
    }
  else
    return Zero;
}

219
220
expr_t
DataTree::AddTimes(expr_t iArg1, expr_t iArg2)
221
222
223
224
225
226
227
228
229
230
231
{
  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)
        {
232
          expr_t tmp = iArg1;
233
234
235
          iArg1 = iArg2;
          iArg2 = tmp;
        }
236
      return AddBinaryOp(iArg1, BinaryOpcode::times, iArg2);
237
238
239
240
241
242
243
244
245
246
247
    }
  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;
}

248
expr_t
249
DataTree::AddDivide(expr_t iArg1, expr_t iArg2) noexcept(false)
250
{
251
  if (iArg2 == One)
252
    return iArg1;
253
254
255

  // This test should be before the next two, otherwise 0/0 won't be rejected
  if (iArg2 == Zero)
256
    {
257
      cerr << "ERROR: Division by zero!" << endl;
258
      throw DivisionByZeroException();
259
    }
260
261
262
263
264
265
266

  if (iArg1 == Zero)
    return Zero;

  if (iArg1 == iArg2)
    return One;

267
  return AddBinaryOp(iArg1, BinaryOpcode::divide, iArg2);
268
269
}

270
271
expr_t
DataTree::AddLess(expr_t iArg1, expr_t iArg2)
272
{
273
  return AddBinaryOp(iArg1, BinaryOpcode::less, iArg2);
274
275
}

276
277
expr_t
DataTree::AddGreater(expr_t iArg1, expr_t iArg2)
278
{
279
  return AddBinaryOp(iArg1, BinaryOpcode::greater, iArg2);
280
281
}

282
283
expr_t
DataTree::AddLessEqual(expr_t iArg1, expr_t iArg2)
284
{
285
  return AddBinaryOp(iArg1, BinaryOpcode::lessEqual, iArg2);
286
287
}

288
289
expr_t
DataTree::AddGreaterEqual(expr_t iArg1, expr_t iArg2)
290
{
291
  return AddBinaryOp(iArg1, BinaryOpcode::greaterEqual, iArg2);
292
293
}

294
295
expr_t
DataTree::AddEqualEqual(expr_t iArg1, expr_t iArg2)
296
{
297
  return AddBinaryOp(iArg1, BinaryOpcode::equalEqual, iArg2);
298
299
}

300
301
expr_t
DataTree::AddDifferent(expr_t iArg1, expr_t iArg2)
302
{
303
  return AddBinaryOp(iArg1, BinaryOpcode::different, iArg2);
304
305
}

306
307
expr_t
DataTree::AddPower(expr_t iArg1, expr_t iArg2)
308
{
309
  if (iArg1 != Zero && iArg2 != Zero && iArg1 != One && iArg2 != One)
310
    return AddBinaryOp(iArg1, BinaryOpcode::power, iArg2);
311
312
  else if (iArg1 == One)
    return One;
313
314
315
316
317
318
319
320
  else if (iArg2 == One)
    return iArg1;
  else if (iArg2 == Zero)
    return One;
  else
    return Zero;
}

321
322
323
324
expr_t
DataTree::AddPowerDeriv(expr_t iArg1, expr_t iArg2, int powerDerivOrder)
{
  assert(powerDerivOrder > 0);
325
  return AddBinaryOp(iArg1, BinaryOpcode::powerDeriv, iArg2, powerDerivOrder);
326
327
}

328
329
330
expr_t
DataTree::AddDiff(expr_t iArg1)
{
331
  return AddUnaryOp(UnaryOpcode::diff, iArg1);
332
333
}

334
335
336
expr_t
DataTree::AddAdl(expr_t iArg1, const string &name, const vector<int> &lags)
{
337
  return AddUnaryOp(UnaryOpcode::adl, iArg1, 0, 0, 0, string(name), lags);
338
339
}

340
341
expr_t
DataTree::AddExp(expr_t iArg1)
342
343
{
  if (iArg1 != Zero)
344
    return AddUnaryOp(UnaryOpcode::exp, iArg1);
345
346
347
348
  else
    return One;
}

349
350
expr_t
DataTree::AddLog(expr_t iArg1)
351
352
{
  if (iArg1 != Zero && iArg1 != One)
353
    return AddUnaryOp(UnaryOpcode::log, iArg1);
354
355
356
357
  else if (iArg1 == One)
    return Zero;
  else
    {
358
      cerr << "ERROR: log(0) not defined!" << endl;
359
      exit(EXIT_FAILURE);
360
361
362
    }
}

363
364
expr_t
DataTree::AddLog10(expr_t iArg1)
365
366
{
  if (iArg1 != Zero && iArg1 != One)
367
    return AddUnaryOp(UnaryOpcode::log10, iArg1);
368
369
370
371
  else if (iArg1 == One)
    return Zero;
  else
    {
372
      cerr << "ERROR: log10(0) not defined!" << endl;
373
      exit(EXIT_FAILURE);
374
375
376
    }
}

377
378
expr_t
DataTree::AddCos(expr_t iArg1)
379
380
{
  if (iArg1 != Zero)
381
    return AddUnaryOp(UnaryOpcode::cos, iArg1);
382
383
384
385
  else
    return One;
}

386
387
expr_t
DataTree::AddSin(expr_t iArg1)
388
389
{
  if (iArg1 != Zero)
390
    return AddUnaryOp(UnaryOpcode::sin, iArg1);
391
392
393
394
  else
    return Zero;
}

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

404
405
expr_t
DataTree::AddAcos(expr_t iArg1)
406
407
{
  if (iArg1 != One)
408
    return AddUnaryOp(UnaryOpcode::acos, iArg1);
409
410
411
412
  else
    return Zero;
}

413
414
expr_t
DataTree::AddAsin(expr_t iArg1)
415
416
{
  if (iArg1 != Zero)
417
    return AddUnaryOp(UnaryOpcode::asin, iArg1);
418
419
420
421
  else
    return Zero;
}

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

431
432
expr_t
DataTree::AddCosh(expr_t iArg1)
433
434
{
  if (iArg1 != Zero)
435
    return AddUnaryOp(UnaryOpcode::cosh, iArg1);
436
437
438
439
  else
    return One;
}

440
441
expr_t
DataTree::AddSinh(expr_t iArg1)
442
443
{
  if (iArg1 != Zero)
444
    return AddUnaryOp(UnaryOpcode::sinh, iArg1);
445
446
447
448
  else
    return Zero;
}

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

458
459
expr_t
DataTree::AddAcosh(expr_t iArg1)
460
461
{
  if (iArg1 != One)
462
    return AddUnaryOp(UnaryOpcode::acosh, iArg1);
463
464
465
466
  else
    return Zero;
}

467
468
expr_t
DataTree::AddAsinh(expr_t iArg1)
469
470
{
  if (iArg1 != Zero)
471
    return AddUnaryOp(UnaryOpcode::asinh, iArg1);
472
473
474
475
  else
    return Zero;
}

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

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

494
495
496
497
498
499
500
501
expr_t
DataTree::AddAbs(expr_t iArg1)
{
  if (iArg1 == Zero)
    return Zero;
  if (iArg1 == One)
    return One;
  else
502
    return AddUnaryOp(UnaryOpcode::abs, iArg1);
503
504
505
506
507
508
509
510
511
512
}

expr_t
DataTree::AddSign(expr_t iArg1)
{
  if (iArg1 == Zero)
    return Zero;
  if (iArg1 == One)
    return One;
  else
513
    return AddUnaryOp(UnaryOpcode::sign, iArg1);
514
515
}

516
517
expr_t
DataTree::AddErf(expr_t iArg1)
518
519
{
  if (iArg1 != Zero)
520
    return AddUnaryOp(UnaryOpcode::erf, iArg1);
521
522
523
524
  else
    return Zero;
}

525
526
expr_t
DataTree::AddMax(expr_t iArg1, expr_t iArg2)
527
{
528
  return AddBinaryOp(iArg1, BinaryOpcode::max, iArg2);
529
530
}

531
532
expr_t
DataTree::AddMin(expr_t iArg1, expr_t iArg2)
533
{
534
  return AddBinaryOp(iArg1, BinaryOpcode::min, iArg2);
535
536
}

537
538
expr_t
DataTree::AddNormcdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
539
{
540
  return AddTrinaryOp(iArg1, TrinaryOpcode::normcdf, iArg2, iArg3);
541
542
}

543
544
expr_t
DataTree::AddNormpdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
545
{
546
  return AddTrinaryOp(iArg1, TrinaryOpcode::normpdf, iArg2, iArg3);
547
548
}

549
550
expr_t
DataTree::AddSteadyState(expr_t iArg1)
551
{
552
  return AddUnaryOp(UnaryOpcode::steadyState, iArg1);
553
554
}

555
556
557
expr_t
DataTree::AddSteadyStateParamDeriv(expr_t iArg1, int param_symb_id)
{
558
  return AddUnaryOp(UnaryOpcode::steadyStateParamDeriv, iArg1, 0, param_symb_id);
559
560
561
562
563
}

expr_t
DataTree::AddSteadyStateParam2ndDeriv(expr_t iArg1, int param1_symb_id, int param2_symb_id)
{
564
  return AddUnaryOp(UnaryOpcode::steadyStateParam2ndDeriv, iArg1, 0, param1_symb_id, param2_symb_id);
565
566
}

567
568
expr_t
DataTree::AddExpectation(int iArg1, expr_t iArg2)
569
{
570
  return AddUnaryOp(UnaryOpcode::expectation, iArg2, iArg1);
571
572
}

Houtan Bastani's avatar
Houtan Bastani committed
573
expr_t
574
DataTree::AddVarExpectation(const string &model_name)
Houtan Bastani's avatar
Houtan Bastani committed
575
{
576
  auto it = var_expectation_node_map.find(model_name);
577
578
579
  if (it != var_expectation_node_map.end())
    return it->second;

580
581
582
583
584
  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
585
586
}

587
expr_t
Houtan Bastani's avatar
Houtan Bastani committed
588
DataTree::AddPacExpectation(const string &model_name)
589
{
590
  auto it = pac_expectation_node_map.find(model_name);
591
592
  if (it != pac_expectation_node_map.end())
    return it->second;
593

594
595
596
597
598
  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;
599
600
}

601
602
expr_t
DataTree::AddEqual(expr_t iArg1, expr_t iArg2)
603
{
604
  return AddBinaryOp(iArg1, BinaryOpcode::equal, iArg2);
605
606
607
}

void
608
DataTree::AddLocalVariable(int symb_id, expr_t value) noexcept(false)
609
{
610
  assert(symbol_table.getType(symb_id) == SymbolType::modelLocalVariable);
sebastien's avatar
sebastien committed
611

612
  // Throw an exception if symbol already declared
613
  auto it = local_variables_table.find(symb_id);
614
  if (it != local_variables_table.end())
615
    throw LocalVariableException(symbol_table.getName(symb_id));
616

617
  local_variables_table[symb_id] = value;
618
  local_variables_vector.push_back(symb_id);
619
620
}

621
622
expr_t
DataTree::AddExternalFunction(int symb_id, const vector<expr_t> &arguments)
623
{
624
  assert(symbol_table.getType(symb_id) == SymbolType::externalFunction);
625

626
  auto it = external_function_node_map.find({ arguments, symb_id });
627
628
629
  if (it != external_function_node_map.end())
    return it->second;

630
631
632
633
634
  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;
635
}
sebastien's avatar
sebastien committed
636

637
expr_t
638
DataTree::AddFirstDerivExternalFunction(int top_level_symb_id, const vector<expr_t> &arguments, int input_index)
639
{
640
  assert(symbol_table.getType(top_level_symb_id) == SymbolType::externalFunction);
641

642
  auto it
643
    = first_deriv_external_function_node_map.find({ arguments, input_index, top_level_symb_id });
644
645
646
  if (it != first_deriv_external_function_node_map.end())
    return it->second;

647
648
649
650
651
  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;
652
}
653

654
expr_t
655
DataTree::AddSecondDerivExternalFunction(int top_level_symb_id, const vector<expr_t> &arguments, int input_index1, int input_index2)
656
{
657
  assert(symbol_table.getType(top_level_symb_id) == SymbolType::externalFunction);
658

659
  auto it
660
    = second_deriv_external_function_node_map.find({ arguments, input_index1, input_index2,
661
          top_level_symb_id });
662
663
664
  if (it != second_deriv_external_function_node_map.end())
    return it->second;

665
666
667
668
669
  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;
670
}
sebastien's avatar
sebastien committed
671

672
673
674
bool
DataTree::isSymbolUsed(int symb_id) const
{
675
676
  for (const auto & it : variable_node_map)
    if (it.first.first == symb_id)
677
678
679
680
681
682
683
684
      return true;

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

  return false;
}

685
int
686
DataTree::getDerivID(int symb_id, int lag) const noexcept(false)
687
688
689
690
{
  throw UnknownDerivIDException();
}

691
SymbolType
692
DataTree::getTypeByDerivID(int deriv_id) const noexcept(false)
693
694
695
696
697
{
  throw UnknownDerivIDException();
}

int
698
DataTree::getLagByDerivID(int deriv_id) const noexcept(false)
699
700
701
702
703
{
  throw UnknownDerivIDException();
}

int
704
DataTree::getSymbIDByDerivID(int deriv_id) const noexcept(false)
705
706
707
708
709
710
711
712
713
{
  throw UnknownDerivIDException();
}

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

714
int
715
DataTree::getDynJacobianCol(int deriv_id) const noexcept(false)
716
717
718
{
  throw UnknownDerivIDException();
}
sebastien's avatar
sebastien committed
719
720

bool
Sébastien Villemot's avatar
Sébastien Villemot committed
721
DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
sebastien's avatar
sebastien committed
722
{
723
  for (const auto & it : unary_op_node_map)
724
    if (get<1>(it.first) == opcode)
Sébastien Villemot's avatar
Sébastien Villemot committed
725
726
727
728
729
730
731
732
      return true;

  return false;
}

bool
DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
{
733
  for (const auto & it : binary_op_node_map)
734
    if (get<2>(it.first) == opcode)
Sébastien Villemot's avatar
Sébastien Villemot committed
735
736
737
738
739
740
741
742
      return true;

  return false;
}

bool
DataTree::isTrinaryOpUsed(TrinaryOpcode opcode) const
{
743
  for (const auto & it : trinary_op_node_map)
744
    if (get<3>(it.first) == opcode)
sebastien's avatar
sebastien committed
745
746
747
748
      return true;

  return false;
}
749
750
751
752

bool
DataTree::isExternalFunctionUsed(int symb_id) const
{
753
754
  for (const auto & it : external_function_node_map)
    if (it.first.second == symb_id)
755
756
757
758
759
760
761
762
      return true;

  return false;
}

bool
DataTree::isFirstDerivExternalFunctionUsed(int symb_id) const
{
763
  for (const auto & it : first_deriv_external_function_node_map)
764
    if (get<2>(it.first) == symb_id)
765
766
767
768
769
770
771
772
      return true;

  return false;
}

bool
DataTree::isSecondDerivExternalFunctionUsed(int symb_id) const
{
773
  for (const auto & it : second_deriv_external_function_node_map)
774
    if (get<3>(it.first) == symb_id)
775
776
777
778
      return true;

  return false;
}
Sébastien Villemot's avatar
Sébastien Villemot committed
779
780
781
782
783

int
DataTree::minLagForSymbol(int symb_id) const
{
  int r = 0;
784
785
786
  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
787
788
  return r;
}
789
790
791
792

void
DataTree::writePowerDerivCHeader(ostream &output) const
{
793
  if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
794
795
796
797
    output << "double getPowerDeriv(double, double, int);" << endl;
}

void
798
DataTree::writePowerDeriv(ostream &output) const
799
{
800
  if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
801
802
803
804
805
    output << "/*" << endl
           << " * The k-th derivative of x^p" << endl
           << " */" << endl
           << "double getPowerDeriv(double x, double p, int k)" << endl
           << "{" << endl
806
807
808
           << "#ifdef _MSC_VER" << endl
           << "# define nearbyint(x) (fabs((x)-floor(x)) < fabs((x)-ceil(x)) ? floor(x) : ceil(x))" << endl
           << "#endif" << endl
809
           << "  if ( fabs(x) < " << near_zero << " && p > 0 && k > p && fabs(p-nearbyint(p)) < " << near_zero << " )" << endl
810
811
812
813
814
815
816
817
818
819
820
           << "    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;
}
821
822
823
824

void
DataTree::writeNormcdfCHeader(ostream &output) const
{
825
#if defined(_WIN32) || defined(__CYGWIN32__)
826
  if (isTrinaryOpUsed(TrinaryOpcode::normcdf))
827
828
829
830
831
832
833
834
835
    output << "#ifdef _MSC_VER" << endl
           << "double normcdf(double);" << endl
           << "#endif" << endl;
#endif
}

void
DataTree::writeNormcdf(ostream &output) const
{
836
#if defined(_WIN32) || defined(__CYGWIN32__)
837
  if (isTrinaryOpUsed(TrinaryOpcode::normcdf))
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
    output << endl
           << "#ifdef _MSC_VER" << endl
           << "/*" << endl
           << " * Define normcdf for MSVC compiler" << endl
           << " */" << endl
           << "double normcdf(double x)" << endl
           << "{" << endl
           << "#if _MSC_VER >= 1700" << endl
           << "  return 0.5 * erfc(-x * M_SQRT1_2);" << endl
           << "#else" << endl
           << "  // From http://www.johndcook.com/blog/cpp_phi" << endl
           << "  double a1 =  0.254829592;" << endl
           << "  double a2 = -0.284496736;" << endl
           << "  double a3 =  1.421413741;" << endl
           << "  double a4 = -1.453152027;" << endl
           << "  double a5 =  1.061405429;" << endl
           << "  double p  =  0.3275911;" << endl
           << "  int sign = (x < 0) ? -1 : 1;" << endl
           << "  x = fabs(x)/sqrt(2.0);" << endl
           << "  // From the Handbook of Mathematical Functions by Abramowitz and Stegun, formula 7.1.26" << endl
           << "  double t = 1.0/(1.0 + p*x);" << endl
           << "  double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);" << endl
           << "  return 0.5*(1.0 + sign*y);" << endl
           << "#endif" << endl
           << "}" << endl
           << "#endif" << endl;
#endif
}
866
867
868
869
870
871
872
873
874

string
DataTree::packageDir(const string &package)
{
  regex pat{"\\."};
  string dirname = "+" + regex_replace(package, pat, "/+");
  boost::filesystem::create_directories(dirname);
  return dirname;
}