ModelTree.cc 8.82 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
#include <cassert>
22
23
24
25
26
27
#include <iostream>

#include "ModelTree.hh"

ModelTree::ModelTree(SymbolTable &symbol_table_arg,
                     NumericalConstants &num_constants_arg) :
28
  DataTree(symbol_table_arg, num_constants_arg),
sebastien's avatar
sebastien committed
29
  mode(eStandardMode)
30
{
31
32
  for(int i=0; i < 3; i++)
    NNZDerivatives[i] = 0;
33
34
35
36
}

int
ModelTree::equation_number() const
37
38
39
{
  return(equations.size());
}
40
41
42
43

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

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

68
69
70
71
72
void
ModelTree::computeHessian(const set<int> &vars)
{
  for (first_derivatives_type::const_iterator it = first_derivatives.begin();
       it != first_derivatives.end(); it++)
73
    {
74
75
76
77
78
79
80
      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++)
81
        {
82
83
84
85
86
87
88
89
          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;
90
91
92
93
	  if (var2 == var1)
	    ++NNZDerivatives[1];
	  else
	    NNZDerivatives[1] += 2;
94
95
        }
    }
96
}
97

98
99
100
101
102
void
ModelTree::computeThirdDerivatives(const set<int> &vars)
{
  for (second_derivatives_type::const_iterator it = second_derivatives.begin();
       it != second_derivatives.end(); it++)
103
    {
104
105
106
107
108
109
110
111
112
113
114
      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++)
115
        {
116
117
118
119
120
121
122
123
          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;
124
125
126
127
128
129
	  if (var3 == var2 && var2 == var1)
	    ++NNZDerivatives[2];
	  else if (var3 == var2 || var2 == var1)
	    NNZDerivatives[2] += 3;
	  else
	    NNZDerivatives[2] += 6;
130
131
132
133
134
        }
    }
}

void
135
ModelTree::computeTemporaryTerms()
136
137
138
139
140
141
{
  map<NodeID, int> reference_count;
  temporary_terms.clear();

  bool is_matlab = (mode != eDLLMode);

142
143
  for (vector<BinaryOpNode *>::iterator it = equations.begin();
       it != equations.end(); it++)
144
145
    (*it)->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);

146
147
  for (first_derivatives_type::iterator it = first_derivatives.begin();
       it != first_derivatives.end(); it++)
148
149
    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);

150
151
152
  for (second_derivatives_type::iterator it = second_derivatives.begin();
       it != second_derivatives.end(); it++)
    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
153

154
155
156
  for (third_derivatives_type::iterator it = third_derivatives.begin();
       it != third_derivatives.end(); it++)
    it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
157
158
159
}

void
sebastien's avatar
sebastien committed
160
161
ModelTree::writeTemporaryTerms(const temporary_terms_type &tt, ostream &output,
                               ExprNodeOutputType output_type) const
162
{
sebastien's avatar
sebastien committed
163
  // Local var used to keep track of temp nodes already written
164
  temporary_terms_type tt2;
165

166
  if (tt.size() > 0 && (IS_C(output_type)))
sebastien's avatar
sebastien committed
167
    output << "double" << endl;
168

sebastien's avatar
sebastien committed
169
170
  for (temporary_terms_type::const_iterator it = tt.begin();
       it != tt.end(); it++)
171
    {
172
      if (IS_C(output_type) && it != tt.begin())
173
        output << "," << endl;
174

sebastien's avatar
sebastien committed
175
      (*it)->writeOutput(output, output_type, tt);
176
      output << " = ";
177

178
      (*it)->writeOutput(output, output_type, tt2);
179

180
181
      // Insert current node into tt2
      tt2.insert(*it);
182

183
      if (IS_MATLAB(output_type))
184
185
        output << ";" << endl;
    }
186
  if (IS_C(output_type))
187
188
    output << ";" << endl;
}
189
190
191

void
ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type) const
192
193
194
195
196
197
{
  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;
198

199
      if (IS_C(output_type))
200
        output << "double ";
201

sebastien's avatar
sebastien committed
202
      output << symbol_table.getName(id) << " = ";
203
204
205
206
207
      // Use an empty set for the temporary terms
      value->writeOutput(output, output_type, temporary_terms_type());
      output << ";" << endl;
    }
}
208
209
210

void
ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const
211
212
213
214
{
  for (int eq = 0; eq < (int) equations.size(); eq++)
    {
      BinaryOpNode *eq_node = equations[eq];
215

sebastien's avatar
sebastien committed
216
      NodeID lhs = eq_node->get_arg1();
217
218
219
      output << "lhs =";
      lhs->writeOutput(output, output_type, temporary_terms);
      output << ";" << endl;
220

sebastien's avatar
sebastien committed
221
      NodeID rhs = eq_node->get_arg2();
222
223
224
      output << "rhs =";
      rhs->writeOutput(output, output_type, temporary_terms);
      output << ";" << endl;
225

226
      output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type) << "= lhs-rhs;" << endl;
227
228
    }
}
229

230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
void
ModelTree::writeLatexModelFile(const string &filename, ExprNodeOutputType output_type) const
{
  ofstream output;
  output.open(filename.c_str(), ios::out | ios::binary);
  if (!output.is_open())
    {
      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
      exit(EXIT_FAILURE);
    }

  output << "\\documentclass[10pt,a4paper]{article}" << endl
         << "\\usepackage[landscape]{geometry}" << endl
         << "\\usepackage{fullpage}" << endl
         << "\\begin{document}" << endl
         << "\\footnotesize" << endl;

  // Write model local variables
  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;

      output << "\\begin{equation*}" << endl
             << symbol_table.getName(id) << " = ";
      // Use an empty set for the temporary terms
      value->writeOutput(output, output_type, temporary_terms_type());
      output << endl << "\\end{equation*}" << endl;
    }

  for (int eq = 0; eq < (int) equations.size(); eq++)
    {
      output << "\\begin{equation}" << endl
             << "% Equation " << eq+1 << endl;
      equations[eq]->writeOutput(output, output_type, temporary_terms_type());
      output << endl << "\\end{equation}" << endl;
    }

  output << "\\end{document}" << endl;

  output.close();
}

274
void
sebastien's avatar
sebastien committed
275
ModelTree::addEquation(NodeID eq)
276
{
sebastien's avatar
sebastien committed
277
  BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
278
  assert(beq != NULL && beq->get_op_code() == oEqual);
279
280
281
282
283

  equations.push_back(beq);
}

void
284
ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
285
{
286
287
  output << LEFT_ARRAY_SUBSCRIPT(output_type);
  if (IS_MATLAB(output_type))
288
289
290
    output << eq_nb + 1 << ", " << col_nb + 1;
  else
    output << eq_nb + col_nb * equations.size();
291
  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
292
}
293
294
295
296
297
298
299
300
301
302
303

void
ModelTree::hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
{
  output << LEFT_ARRAY_SUBSCRIPT(output_type);
  if (IS_MATLAB(output_type))
    output << row_nb + 1 << ", " << col_nb + 1;
  else
    output << row_nb + col_nb * NNZDerivatives[1];
  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
}