DataTree.cc 19 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
DataTree::DataTree(SymbolTable &symbol_table_arg,
                   NumericalConstants &num_constants_arg,
31
                   ExternalFunctionsTable &external_functions_table_arg) :
32 33
  symbol_table(symbol_table_arg),
  num_constants(num_constants_arg),
34
  external_functions_table(external_functions_table_arg)
35
{
36 37 38
  Zero = AddNonNegativeConstant("0");
  One = AddNonNegativeConstant("1");
  Two = AddNonNegativeConstant("2");
39 40

  MinusOne = AddUMinus(One);
41

42 43
  NaN = AddNonNegativeConstant("NaN");
  Infinity = AddNonNegativeConstant("Inf");
44
  MinusInfinity = AddUMinus(Infinity);
sebastien's avatar
sebastien committed
45

46
  Pi = AddNonNegativeConstant("3.141592653589793");
47 48
}

49
DataTree::~DataTree() = default;
50

51
expr_t
52
DataTree::AddNonNegativeConstant(const string &value)
53
{
54
  int id = num_constants.AddNonNegativeConstant(value);
55

56
  auto it = num_const_node_map.find(id);
57 58
  if (it != num_const_node_map.end())
    return it->second;
59 60 61 62 63 64

  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;
65 66
}

sebastien's avatar
sebastien committed
67 68
VariableNode *
DataTree::AddVariableInternal(int symb_id, int lag)
69
{
70
  auto it = variable_node_map.find({ symb_id, lag });
71 72
  if (it != variable_node_map.end())
    return it->second;
73 74 75 76 77 78

  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;
79
}
sebastien's avatar
sebastien committed
80

81 82 83
bool
DataTree::ParamUsedWithLeadLagInternal() const
{
84
  for (const auto & it : variable_node_map)
85
    if (symbol_table.getType(it.first.first) == SymbolType::parameter && it.first.second != 0)
86 87 88 89
      return true;
  return false;
}

sebastien's avatar
sebastien committed
90 91
VariableNode *
DataTree::AddVariable(int symb_id, int lag)
sebastien's avatar
sebastien committed
92
{
93
  assert(lag == 0);
sebastien's avatar
sebastien committed
94
  return AddVariableInternal(symb_id, lag);
sebastien's avatar
sebastien committed
95 96
}

97 98
expr_t
DataTree::AddPlus(expr_t iArg1, expr_t iArg2)
99 100 101 102
{
  if (iArg1 != Zero && iArg2 != Zero)
    {
      // Simplify x+(-y) in x-y
103
      auto *uarg2 = dynamic_cast<UnaryOpNode *>(iArg2);
104
      if (uarg2 != nullptr && uarg2->get_op_code() == UnaryOpcode::uminus)
sebastien's avatar
sebastien committed
105
        return AddMinus(iArg1, uarg2->get_arg());
106 107 108 109 110

      // To treat commutativity of "+"
      // Nodes iArg1 and iArg2 are sorted by index
      if (iArg1->idx > iArg2->idx)
        {
111
          expr_t tmp = iArg1;
112 113 114
          iArg1 = iArg2;
          iArg2 = tmp;
        }
115
      return AddBinaryOp(iArg1, BinaryOpcode::plus, iArg2);
116 117 118 119 120 121 122 123 124
    }
  else if (iArg1 != Zero)
    return iArg1;
  else if (iArg2 != Zero)
    return iArg2;
  else
    return Zero;
}

125 126
expr_t
DataTree::AddMinus(expr_t iArg1, expr_t iArg2)
127
{
128
  if (iArg2 == Zero)
129
    return iArg1;
130 131

  if (iArg1 == Zero)
132
    return AddUMinus(iArg2);
133 134

  if (iArg1 == iArg2)
135
    return Zero;
136

137
  return AddBinaryOp(iArg1, BinaryOpcode::minus, iArg2);
138 139
}

140 141
expr_t
DataTree::AddUMinus(expr_t iArg1)
142 143 144 145
{
  if (iArg1 != Zero)
    {
      // Simplify -(-x) in x
146
      auto *uarg = dynamic_cast<UnaryOpNode *>(iArg1);
147
      if (uarg != nullptr && uarg->get_op_code() == UnaryOpcode::uminus)
sebastien's avatar
sebastien committed
148
        return uarg->get_arg();
149

150
      return AddUnaryOp(UnaryOpcode::uminus, iArg1);
151 152 153 154 155
    }
  else
    return Zero;
}

156 157
expr_t
DataTree::AddTimes(expr_t iArg1, expr_t iArg2)
158 159 160 161 162 163 164 165 166 167 168
{
  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)
        {
169
          expr_t tmp = iArg1;
170 171 172
          iArg1 = iArg2;
          iArg2 = tmp;
        }
173
      return AddBinaryOp(iArg1, BinaryOpcode::times, iArg2);
174 175 176 177 178 179 180 181 182 183 184
    }
  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;
}

185
expr_t
186
DataTree::AddDivide(expr_t iArg1, expr_t iArg2) noexcept(false)
187
{
188
  if (iArg2 == One)
189
    return iArg1;
190 191 192

  // This test should be before the next two, otherwise 0/0 won't be rejected
  if (iArg2 == Zero)
193
    {
194
      cerr << "ERROR: Division by zero!" << endl;
195
      throw DivisionByZeroException();
196
    }
197 198 199 200 201 202 203

  if (iArg1 == Zero)
    return Zero;

  if (iArg1 == iArg2)
    return One;

204
  return AddBinaryOp(iArg1, BinaryOpcode::divide, iArg2);
205 206
}

207 208
expr_t
DataTree::AddLess(expr_t iArg1, expr_t iArg2)
209
{
210
  return AddBinaryOp(iArg1, BinaryOpcode::less, iArg2);
211 212
}

213 214
expr_t
DataTree::AddGreater(expr_t iArg1, expr_t iArg2)
215
{
216
  return AddBinaryOp(iArg1, BinaryOpcode::greater, iArg2);
217 218
}

219 220
expr_t
DataTree::AddLessEqual(expr_t iArg1, expr_t iArg2)
221
{
222
  return AddBinaryOp(iArg1, BinaryOpcode::lessEqual, iArg2);
223 224
}

225 226
expr_t
DataTree::AddGreaterEqual(expr_t iArg1, expr_t iArg2)
227
{
228
  return AddBinaryOp(iArg1, BinaryOpcode::greaterEqual, iArg2);
229 230
}

231 232
expr_t
DataTree::AddEqualEqual(expr_t iArg1, expr_t iArg2)
233
{
234
  return AddBinaryOp(iArg1, BinaryOpcode::equalEqual, iArg2);
235 236
}

237 238
expr_t
DataTree::AddDifferent(expr_t iArg1, expr_t iArg2)
239
{
240
  return AddBinaryOp(iArg1, BinaryOpcode::different, iArg2);
241 242
}

243 244
expr_t
DataTree::AddPower(expr_t iArg1, expr_t iArg2)
245
{
246
  if (iArg1 != Zero && iArg2 != Zero && iArg1 != One && iArg2 != One)
247
    return AddBinaryOp(iArg1, BinaryOpcode::power, iArg2);
248 249
  else if (iArg1 == One)
    return One;
250 251 252 253 254 255 256 257
  else if (iArg2 == One)
    return iArg1;
  else if (iArg2 == Zero)
    return One;
  else
    return Zero;
}

258 259 260 261
expr_t
DataTree::AddPowerDeriv(expr_t iArg1, expr_t iArg2, int powerDerivOrder)
{
  assert(powerDerivOrder > 0);
262
  return AddBinaryOp(iArg1, BinaryOpcode::powerDeriv, iArg2, powerDerivOrder);
263 264
}

265 266 267
expr_t
DataTree::AddDiff(expr_t iArg1)
{
268
  return AddUnaryOp(UnaryOpcode::diff, iArg1);
269 270
}

271 272 273
expr_t
DataTree::AddAdl(expr_t iArg1, const string &name, const vector<int> &lags)
{
274
  return AddUnaryOp(UnaryOpcode::adl, iArg1, 0, 0, 0, string(name), lags);
275 276
}

277 278
expr_t
DataTree::AddExp(expr_t iArg1)
279 280
{
  if (iArg1 != Zero)
281
    return AddUnaryOp(UnaryOpcode::exp, iArg1);
282 283 284 285
  else
    return One;
}

286 287
expr_t
DataTree::AddLog(expr_t iArg1)
288 289
{
  if (iArg1 != Zero && iArg1 != One)
290
    return AddUnaryOp(UnaryOpcode::log, iArg1);
291 292 293 294
  else if (iArg1 == One)
    return Zero;
  else
    {
295
      cerr << "ERROR: log(0) not defined!" << endl;
296
      exit(EXIT_FAILURE);
297 298 299
    }
}

300 301
expr_t
DataTree::AddLog10(expr_t iArg1)
302 303
{
  if (iArg1 != Zero && iArg1 != One)
304
    return AddUnaryOp(UnaryOpcode::log10, iArg1);
305 306 307 308
  else if (iArg1 == One)
    return Zero;
  else
    {
309
      cerr << "ERROR: log10(0) not defined!" << endl;
310
      exit(EXIT_FAILURE);
311 312 313
    }
}

314 315
expr_t
DataTree::AddCos(expr_t iArg1)
316 317
{
  if (iArg1 != Zero)
318
    return AddUnaryOp(UnaryOpcode::cos, iArg1);
319 320 321 322
  else
    return One;
}

323 324
expr_t
DataTree::AddSin(expr_t iArg1)
325 326
{
  if (iArg1 != Zero)
327
    return AddUnaryOp(UnaryOpcode::sin, iArg1);
328 329 330 331
  else
    return Zero;
}

332 333
expr_t
DataTree::AddTan(expr_t iArg1)
334 335
{
  if (iArg1 != Zero)
336
    return AddUnaryOp(UnaryOpcode::tan, iArg1);
337 338 339 340
  else
    return Zero;
}

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

350 351
expr_t
DataTree::AddAsin(expr_t iArg1)
352 353
{
  if (iArg1 != Zero)
354
    return AddUnaryOp(UnaryOpcode::asin, iArg1);
355 356 357 358
  else
    return Zero;
}

359 360
expr_t
DataTree::AddAtan(expr_t iArg1)
361 362
{
  if (iArg1 != Zero)
363
    return AddUnaryOp(UnaryOpcode::atan, iArg1);
364 365 366 367
  else
    return Zero;
}

368 369
expr_t
DataTree::AddCosh(expr_t iArg1)
370 371
{
  if (iArg1 != Zero)
372
    return AddUnaryOp(UnaryOpcode::cosh, iArg1);
373 374 375 376
  else
    return One;
}

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

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

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

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

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

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

431 432 433 434 435 436 437 438
expr_t
DataTree::AddAbs(expr_t iArg1)
{
  if (iArg1 == Zero)
    return Zero;
  if (iArg1 == One)
    return One;
  else
439
    return AddUnaryOp(UnaryOpcode::abs, iArg1);
440 441 442 443 444 445 446 447 448 449
}

expr_t
DataTree::AddSign(expr_t iArg1)
{
  if (iArg1 == Zero)
    return Zero;
  if (iArg1 == One)
    return One;
  else
450
    return AddUnaryOp(UnaryOpcode::sign, iArg1);
451 452
}

453 454
expr_t
DataTree::AddErf(expr_t iArg1)
455 456
{
  if (iArg1 != Zero)
457
    return AddUnaryOp(UnaryOpcode::erf, iArg1);
458 459 460 461
  else
    return Zero;
}

462 463
expr_t
DataTree::AddMax(expr_t iArg1, expr_t iArg2)
464
{
465
  return AddBinaryOp(iArg1, BinaryOpcode::max, iArg2);
466 467
}

468 469
expr_t
DataTree::AddMin(expr_t iArg1, expr_t iArg2)
470
{
471
  return AddBinaryOp(iArg1, BinaryOpcode::min, iArg2);
472 473
}

474 475
expr_t
DataTree::AddNormcdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
476
{
477
  return AddTrinaryOp(iArg1, TrinaryOpcode::normcdf, iArg2, iArg3);
478 479
}

480 481
expr_t
DataTree::AddNormpdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
482
{
483
  return AddTrinaryOp(iArg1, TrinaryOpcode::normpdf, iArg2, iArg3);
484 485
}

486 487
expr_t
DataTree::AddSteadyState(expr_t iArg1)
488
{
489
  return AddUnaryOp(UnaryOpcode::steadyState, iArg1);
490 491
}

492 493 494
expr_t
DataTree::AddSteadyStateParamDeriv(expr_t iArg1, int param_symb_id)
{
495
  return AddUnaryOp(UnaryOpcode::steadyStateParamDeriv, iArg1, 0, param_symb_id);
496 497 498 499 500
}

expr_t
DataTree::AddSteadyStateParam2ndDeriv(expr_t iArg1, int param1_symb_id, int param2_symb_id)
{
501
  return AddUnaryOp(UnaryOpcode::steadyStateParam2ndDeriv, iArg1, 0, param1_symb_id, param2_symb_id);
502 503
}

504 505
expr_t
DataTree::AddExpectation(int iArg1, expr_t iArg2)
506
{
507
  return AddUnaryOp(UnaryOpcode::expectation, iArg2, iArg1);
508 509
}

Houtan Bastani's avatar
Houtan Bastani committed
510
expr_t
511
DataTree::AddVarExpectation(const string &model_name)
Houtan Bastani's avatar
Houtan Bastani committed
512
{
513
  auto it = var_expectation_node_map.find(model_name);
514 515 516
  if (it != var_expectation_node_map.end())
    return it->second;

517 518 519 520 521
  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
522 523
}

524
expr_t
Houtan Bastani's avatar
Houtan Bastani committed
525
DataTree::AddPacExpectation(const string &model_name)
526
{
527
  auto it = pac_expectation_node_map.find(model_name);
528 529
  if (it != pac_expectation_node_map.end())
    return it->second;
530

531 532 533 534 535
  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;
536 537
}

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

void
545
DataTree::AddLocalVariable(int symb_id, expr_t value) noexcept(false)
546
{
547
  assert(symbol_table.getType(symb_id) == SymbolType::modelLocalVariable);
sebastien's avatar
sebastien committed
548

549
  // Throw an exception if symbol already declared
550
  auto it = local_variables_table.find(symb_id);
551
  if (it != local_variables_table.end())
552
    throw LocalVariableException(symbol_table.getName(symb_id));
553

554
  local_variables_table[symb_id] = value;
555
  local_variables_vector.push_back(symb_id);
556 557
}

558 559
expr_t
DataTree::AddExternalFunction(int symb_id, const vector<expr_t> &arguments)
560
{
561
  assert(symbol_table.getType(symb_id) == SymbolType::externalFunction);
562

563
  auto it = external_function_node_map.find({ arguments, symb_id });
564 565 566
  if (it != external_function_node_map.end())
    return it->second;

567 568 569 570 571
  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;
572
}
sebastien's avatar
sebastien committed
573

574
expr_t
575
DataTree::AddFirstDerivExternalFunction(int top_level_symb_id, const vector<expr_t> &arguments, int input_index)
576
{
577
  assert(symbol_table.getType(top_level_symb_id) == SymbolType::externalFunction);
578

579
  auto it
580
    = first_deriv_external_function_node_map.find({ arguments, input_index, top_level_symb_id });
581 582 583
  if (it != first_deriv_external_function_node_map.end())
    return it->second;

584 585 586 587 588
  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;
589
}
590

591
expr_t
592
DataTree::AddSecondDerivExternalFunction(int top_level_symb_id, const vector<expr_t> &arguments, int input_index1, int input_index2)
593
{
594
  assert(symbol_table.getType(top_level_symb_id) == SymbolType::externalFunction);
595

596
  auto it
597
    = second_deriv_external_function_node_map.find({ arguments, input_index1, input_index2,
598
          top_level_symb_id });
599 600 601
  if (it != second_deriv_external_function_node_map.end())
    return it->second;

602 603 604 605 606
  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;
607
}
sebastien's avatar
sebastien committed
608

609 610 611
bool
DataTree::isSymbolUsed(int symb_id) const
{
612 613
  for (const auto & it : variable_node_map)
    if (it.first.first == symb_id)
614 615 616 617 618 619 620 621
      return true;

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

  return false;
}

622
int
623
DataTree::getDerivID(int symb_id, int lag) const noexcept(false)
624 625 626 627
{
  throw UnknownDerivIDException();
}

628
SymbolType
629
DataTree::getTypeByDerivID(int deriv_id) const noexcept(false)
630 631 632 633 634
{
  throw UnknownDerivIDException();
}

int
635
DataTree::getLagByDerivID(int deriv_id) const noexcept(false)
636 637 638 639 640
{
  throw UnknownDerivIDException();
}

int
641
DataTree::getSymbIDByDerivID(int deriv_id) const noexcept(false)
642 643 644 645 646 647 648 649 650
{
  throw UnknownDerivIDException();
}

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

651
int
652
DataTree::getDynJacobianCol(int deriv_id) const noexcept(false)
653 654 655
{
  throw UnknownDerivIDException();
}
sebastien's avatar
sebastien committed
656 657

bool
Sébastien Villemot's avatar
Sébastien Villemot committed
658
DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
sebastien's avatar
sebastien committed
659
{
660
  for (const auto & it : unary_op_node_map)
661
    if (get<1>(it.first) == opcode)
Sébastien Villemot's avatar
Sébastien Villemot committed
662 663 664 665 666 667 668 669
      return true;

  return false;
}

bool
DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
{
670
  for (const auto & it : binary_op_node_map)
671
    if (get<2>(it.first) == opcode)
Sébastien Villemot's avatar
Sébastien Villemot committed
672 673 674 675 676 677 678 679
      return true;

  return false;
}

bool
DataTree::isTrinaryOpUsed(TrinaryOpcode opcode) const
{
680
  for (const auto & it : trinary_op_node_map)
681
    if (get<3>(it.first) == opcode)
sebastien's avatar
sebastien committed
682 683 684 685
      return true;

  return false;
}
686 687 688 689

bool
DataTree::isExternalFunctionUsed(int symb_id) const
{
690 691
  for (const auto & it : external_function_node_map)
    if (it.first.second == symb_id)
692 693 694 695 696 697 698 699
      return true;

  return false;
}

bool
DataTree::isFirstDerivExternalFunctionUsed(int symb_id) const
{
700
  for (const auto & it : first_deriv_external_function_node_map)
701
    if (get<2>(it.first) == symb_id)
702 703 704 705 706 707 708 709
      return true;

  return false;
}

bool
DataTree::isSecondDerivExternalFunctionUsed(int symb_id) const
{
710
  for (const auto & it : second_deriv_external_function_node_map)
711
    if (get<3>(it.first) == symb_id)
712 713 714 715
      return true;

  return false;
}
Sébastien Villemot's avatar
Sébastien Villemot committed
716 717 718 719 720

int
DataTree::minLagForSymbol(int symb_id) const
{
  int r = 0;
721 722 723
  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
724 725
  return r;
}
726 727 728 729

void
DataTree::writePowerDerivCHeader(ostream &output) const
{
730
  if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
731 732 733 734
    output << "double getPowerDeriv(double, double, int);" << endl;
}

void
735
DataTree::writePowerDeriv(ostream &output) const
736
{
737
  if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
738 739 740 741 742
    output << "/*" << endl
           << " * The k-th derivative of x^p" << endl
           << " */" << endl
           << "double getPowerDeriv(double x, double p, int k)" << endl
           << "{" << endl
743 744 745
           << "#ifdef _MSC_VER" << endl
           << "# define nearbyint(x) (fabs((x)-floor(x)) < fabs((x)-ceil(x)) ? floor(x) : ceil(x))" << endl
           << "#endif" << endl
746
           << "  if ( fabs(x) < " << near_zero << " && p > 0 && k > p && fabs(p-nearbyint(p)) < " << near_zero << " )" << endl
747 748 749 750 751 752 753 754 755 756 757
           << "    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;
}
758 759 760 761 762

void
DataTree::writeNormcdfCHeader(ostream &output) const
{
#if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
763
  if (isTrinaryOpUsed(TrinaryOpcode::normcdf))
764 765 766 767 768 769 770 771 772 773
    output << "#ifdef _MSC_VER" << endl
           << "double normcdf(double);" << endl
           << "#endif" << endl;
#endif
}

void
DataTree::writeNormcdf(ostream &output) const
{
#if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
774
  if (isTrinaryOpUsed(TrinaryOpcode::normcdf))
775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802
    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
}
803 804 805 806 807 808 809 810 811

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