ModelTree.cc 9.27 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
#include <iostream>
23
#include <fstream>
24
25
26
27
28

#include "ModelTree.hh"

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

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

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

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];
sebastien's avatar
sebastien committed
215
216
      NodeID lhs = eq_node->get_arg1();
      NodeID rhs = eq_node->get_arg2();
217

218
219
220
221
222
223
224
225
226
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
      // Test if the right hand side of the equation is empty.
      double vrhs = 1.0;
      try
        {
          vrhs = rhs->eval(eval_context_type());
        }
      catch(ExprNode::EvalException &e)
        {
        }

      if (vrhs!=0)// The right hand side of the equation is not empty ==> residual=lhs-rhs;
        {
          output << "lhs =";
          lhs->writeOutput(output, output_type, temporary_terms);
          output << ";" << endl;

          output << "rhs =";
          rhs->writeOutput(output, output_type, temporary_terms);
          output << ";" << endl;

          output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type) 
                 << eq + ARRAY_SUBSCRIPT_OFFSET(output_type) 
                 << RIGHT_ARRAY_SUBSCRIPT(output_type) 
                 << "= lhs-rhs;" << endl;
        }
      else// The right hand side of the equation is empty ==> residual=lhs;
        {
          output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type) 
                 << eq + ARRAY_SUBSCRIPT_OFFSET(output_type) 
                 << RIGHT_ARRAY_SUBSCRIPT(output_type) 
                 << " = "; 
          lhs->writeOutput(output, output_type, temporary_terms);
          output << ";" << endl;      
        }
252
253
    }
}
254

255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
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();
}

299
void
sebastien's avatar
sebastien committed
300
ModelTree::addEquation(NodeID eq)
301
{
sebastien's avatar
sebastien committed
302
  BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
303
  assert(beq != NULL && beq->get_op_code() == oEqual);
304
305
306

  equations.push_back(beq);
}
307
308
309
310
311
312

void
ModelTree::addEquationTags(int i, const string &key, const string &value)
{
  equation_tags.push_back(make_pair(i, make_pair(key, value)));
}
sebastien's avatar
sebastien committed
313
314
315
316
317
318
319
320
321

void
ModelTree::addAuxEquation(NodeID eq)
{
  BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
  assert(beq != NULL && beq->get_op_code() == oEqual);

  aux_equations.push_back(beq);
}