ModelTree.cc 6.8 KB
Newer Older
1
/*
sebastien's avatar
trunk:    
sebastien committed
2
 * Copyright (C) 2003-2009 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
22
23
24
25
26
#include <iostream>

#include "ModelTree.hh"

ModelTree::ModelTree(SymbolTable &symbol_table_arg,
                     NumericalConstants &num_constants_arg) :
27
  DataTree(symbol_table_arg, num_constants_arg),
sebastien's avatar
sebastien committed
28
  mode(eStandardMode)
29
30
31
32
33
{
}

int
ModelTree::equation_number() const
34
35
36
{
  return(equations.size());
}
37
38
39
40

void
ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
                           ExprNodeOutputType output_type,
sebastien's avatar
sebastien committed
41
                           const temporary_terms_type &temporary_terms) const
42
{
43
  first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag)));
44
45
46
47
48
  if (it != first_derivatives.end())
    (it->second)->writeOutput(output, output_type, temporary_terms);
  else
    output << 0;
}
49
50

void
51
ModelTree::computeJacobian(const set<int> &vars)
52
{
53
54
  for(set<int>::const_iterator it = vars.begin();
      it != vars.end(); it++)
55
    for (int eq = 0; eq < (int) equations.size(); eq++)
56
      {
57
        NodeID d1 = equations[eq]->getDerivative(*it);
58
59
        if (d1 == Zero)
          continue;
60
        first_derivatives[make_pair(eq, *it)] = d1;
61
      }
62
}
63

64
65
66
67
68
void
ModelTree::computeHessian(const set<int> &vars)
{
  for (first_derivatives_type::const_iterator it = first_derivatives.begin();
       it != first_derivatives.end(); it++)
69
    {
70
71
72
73
74
75
76
      int eq = it->first.first;
      int var1 = it->first.second;
      NodeID d1 = it->second;

      // Store only second derivatives with var2 <= var1
      for(set<int>::const_iterator it2 = vars.begin();
          it2 != vars.end(); it2++)
77
        {
78
79
80
81
82
83
84
85
          int var2 = *it2;
          if (var2 > var1)
            continue;

          NodeID d2 = d1->getDerivative(var2);
          if (d2 == Zero)
            continue;
          second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
86
87
        }
    }
88
}
89

90
91
92
93
94
void
ModelTree::computeThirdDerivatives(const set<int> &vars)
{
  for (second_derivatives_type::const_iterator it = second_derivatives.begin();
       it != second_derivatives.end(); it++)
95
    {
96
97
98
99
100
101
102
103
104
105
106
      int eq = it->first.first;

      int var1 = it->first.second.first;
      int var2 = it->first.second.second;
      // By construction, var2 <= var1

      NodeID d2 = it->second;

      // Store only third derivatives such that var3 <= var2 <= var1
      for(set<int>::const_iterator it2 = vars.begin();
          it2 != vars.end(); it2++)
107
        {
108
109
110
111
112
113
114
115
          int var3 = *it2;
          if (var3 > var2)
            continue;

          NodeID d3 = d2->getDerivative(var3);
          if (d3 == Zero)
            continue;
          third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
116
117
118
119
120
        }
    }
}

void
121
ModelTree::computeTemporaryTerms()
122
123
124
125
126
127
{
  map<NodeID, int> reference_count;
  temporary_terms.clear();

  bool is_matlab = (mode != eDLLMode);

128
129
  for (vector<BinaryOpNode *>::iterator it = equations.begin();
       it != equations.end(); it++)
130
131
    (*it)->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);

132
133
  for (first_derivatives_type::iterator it = first_derivatives.begin();
       it != first_derivatives.end(); it++)
134
135
    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);

136
137
138
  for (second_derivatives_type::iterator it = second_derivatives.begin();
       it != second_derivatives.end(); it++)
    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
139

140
141
142
  for (third_derivatives_type::iterator it = third_derivatives.begin();
       it != third_derivatives.end(); it++)
    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
143
144
145
}

void
sebastien's avatar
sebastien committed
146
147
ModelTree::writeTemporaryTerms(const temporary_terms_type &tt, ostream &output,
                               ExprNodeOutputType output_type) const
148
{
sebastien's avatar
sebastien committed
149
  // Local var used to keep track of temp nodes already written
150
  temporary_terms_type tt2;
151

sebastien's avatar
sebastien committed
152
153
  if (tt.size() > 0 && (!OFFSET(output_type)))
    output << "double" << endl;
154

sebastien's avatar
sebastien committed
155
156
  for (temporary_terms_type::const_iterator it = tt.begin();
       it != tt.end(); it++)
157
    {
sebastien's avatar
sebastien committed
158
      if (!OFFSET(output_type) && it != tt.begin())
159
        output << "," << endl;
160

sebastien's avatar
sebastien committed
161
      (*it)->writeOutput(output, output_type, tt);
162
      output << " = ";
163

164
      (*it)->writeOutput(output, output_type, tt2);
165

166
167
      // Insert current node into tt2
      tt2.insert(*it);
168

169
170
171
172
173
174
      if (OFFSET(output_type))
        output << ";" << endl;
    }
  if (!OFFSET(output_type))
    output << ";" << endl;
}
175
176
177

void
ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type) const
178
179
180
181
182
183
{
  for (map<int, NodeID>::const_iterator it = local_variables_table.begin();
       it != local_variables_table.end(); it++)
    {
      int id = it->first;
      NodeID value = it->second;
184

185
186
      if (!OFFSET(output_type))
        output << "double ";
187

sebastien's avatar
sebastien committed
188
      output << symbol_table.getName(id) << " = ";
189
190
191
192
193
      // Use an empty set for the temporary terms
      value->writeOutput(output, output_type, temporary_terms_type());
      output << ";" << endl;
    }
}
194
195
196

void
ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const
197
198
199
200
{
  for (int eq = 0; eq < (int) equations.size(); eq++)
    {
      BinaryOpNode *eq_node = equations[eq];
201

sebastien's avatar
sebastien committed
202
      NodeID lhs = eq_node->get_arg1();
203
204
205
      output << "lhs =";
      lhs->writeOutput(output, output_type, temporary_terms);
      output << ";" << endl;
206

sebastien's avatar
sebastien committed
207
      NodeID rhs = eq_node->get_arg2();
208
209
210
      output << "rhs =";
      rhs->writeOutput(output, output_type, temporary_terms);
      output << ";" << endl;
211

212
213
214
      output << "residual" << LPAR(output_type) << eq + OFFSET(output_type) << RPAR(output_type) << "= lhs-rhs;" << endl;
    }
}
215
216

void
sebastien's avatar
sebastien committed
217
ModelTree::addEquation(NodeID eq)
218
{
sebastien's avatar
sebastien committed
219
  BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
220

sebastien's avatar
sebastien committed
221
  if (beq == NULL || beq->get_op_code() != oEqual)
222
    {
sebastien's avatar
sebastien committed
223
      cerr << "ModelTree::addEquation: you didn't provide an equal node!" << endl;
224
      exit(EXIT_FAILURE);
225
226
227
228
229
230
231
    }

  equations.push_back(beq);
}

void
ModelTree::matrixHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
232
233
234
235
236
237
238
239
{
  output << LPAR(output_type);
  if (OFFSET(output_type))
    output << eq_nb + 1 << ", " << col_nb + 1;
  else
    output << eq_nb + col_nb * equations.size();
  output << RPAR(output_type);
}