StaticModel.cc 15.2 KB
Newer Older
sebastien's avatar
sebastien committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*
 * Copyright (C) 2003-2009 Dynare Team
 *
 * 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
#include <cassert>

sebastien's avatar
sebastien committed
23
24
25
26
27
28
29
30
#include <algorithm>
#include <functional>

/*
#include <ext/functional>
using namespace __gnu_cxx;
*/

31
32
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
33

sebastien's avatar
sebastien committed
34
35
#include "StaticModel.hh"

36
37
using namespace boost;

sebastien's avatar
sebastien committed
38
39
StaticModel::StaticModel(SymbolTable &symbol_table_arg,
                         NumericalConstants &num_constants_arg) :
40
  ModelTree(symbol_table_arg, num_constants_arg)
sebastien's avatar
sebastien committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
{
}

void
StaticModel::writeStaticMFile(const string &static_basename) const
{
  string filename = static_basename + ".m";

  ofstream mStaticModelFile;
  mStaticModelFile.open(filename.c_str(), ios::out | ios::binary);
  if (!mStaticModelFile.is_open())
    {
      cerr << "Error: Can't open file " << filename << " for writing" << endl;
      exit(EXIT_FAILURE);
    }
  // Writing comments and function definition command
  mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "(y, x, params)" << endl
                   << "%" << endl
                   << "% Status : Computes static model for Dynare" << endl
                   << "%" << endl
                   << "% Warning : this file is generated automatically by Dynare" << endl
                   << "%           from model file (.mod)" << endl << endl;

  writeStaticModel(mStaticModelFile);

  mStaticModelFile.close();
}

void
StaticModel::writeStaticCFile(const string &static_basename) const
{
  string filename = static_basename + ".c";

  ofstream mStaticModelFile;
  mStaticModelFile.open(filename.c_str(), ios::out | ios::binary);
  if (!mStaticModelFile.is_open())
    {
      cerr << "Error: Can't open file " << filename << " for writing" << endl;
      exit(EXIT_FAILURE);
    }
  mStaticModelFile << "/*" << endl
                   << " * " << filename << " : Computes static model for Dynare" << endl
                   << " * Warning : this file is generated automatically by Dynare" << endl
                   << " *           from model file (.mod)" << endl
                   << endl
                   << " */" << endl
                   << "#include <math.h>" << endl
                   << "#include \"mex.h\"" << endl;

  // Writing the function Static
  writeStaticModel(mStaticModelFile);

  // Writing the gateway routine
  mStaticModelFile << "/* The gateway routine */" << endl
                   << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
                   << "{" << endl
                   << "  double *y, *x, *params;" << endl
                   << "  double *residual, *g1;" << endl
                   << endl
                   << "  /* Create a pointer to the input matrix y. */" << endl
                   << "  y = mxGetPr(prhs[0]);" << endl
                   << endl
                   << "  /* Create a pointer to the input matrix x. */" << endl
                   << "  x = mxGetPr(prhs[1]);" << endl
                   << endl
                   << "  /* Create a pointer to the input matrix params. */" << endl
                   << "  params = mxGetPr(prhs[2]);" << endl
                   << endl
                   << "  residual = NULL;" << endl
                   << "  if (nlhs >= 1)" << endl
                   << "  {" << endl
                   << "      /* Set the output pointer to the output matrix residual. */" << endl
                   << "      plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
                   << "     /* Create a C pointer to a copy of the output matrix residual. */" << endl
                   << "     residual = mxGetPr(plhs[0]);" << endl
                   << "  }" << endl
                   << endl
                   << "  g1 = NULL;" << endl
                   << "  if (nlhs >= 2)" << endl
                   << "  {" << endl
                   << "      /* Set the output pointer to the output matrix g1. */" << endl
                   << "      plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl
                   << "      /* Create a C pointer to a copy of the output matrix g1. */" << endl
                   << "      g1 = mxGetPr(plhs[1]);" << endl
                   << "  }" << endl
                   << endl
                   << "  /* Call the C Static. */" << endl
                   << "  Static(y, x, params, residual, g1);" << endl
                   << "}" << endl;

  mStaticModelFile.close();
}

void
StaticModel::writeStaticModel(ostream &StaticOutput) const
{
  ostringstream model_output;    // Used for storing model equations
  ostringstream jacobian_output; // Used for storing jacobian equations
  ostringstream hessian_output;

  ExprNodeOutputType output_type = (mode == eDLLMode ? oCStaticModel : oMatlabStaticModel);

  writeModelLocalVariables(model_output, output_type);

sebastien's avatar
sebastien committed
145
  writeTemporaryTerms(temporary_terms, model_output, output_type);
sebastien's avatar
sebastien committed
146
147
148
149
150
151
152
153

  writeModelEquations(model_output, output_type);

  // Write Jacobian w.r. to endogenous only
  for (first_derivatives_type::const_iterator it = first_derivatives.begin();
       it != first_derivatives.end(); it++)
    {
      int eq = it->first.first;
154
      int symb_id = it->first.second;
sebastien's avatar
sebastien committed
155
156
      NodeID d1 = it->second;

157
158
      ostringstream g1;
      g1 << "  g1";
159
      jacobianHelper(g1, eq, symbol_table.getTypeSpecificID(symb_id), output_type);
sebastien's avatar
sebastien committed
160

161
162
163
      jacobian_output << g1.str() << "=" << g1.str() << "+";
      d1->writeOutput(jacobian_output, output_type, temporary_terms);
      jacobian_output << ";" << endl;
sebastien's avatar
sebastien committed
164
165
    }

166
  // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
sebastien's avatar
sebastien committed
167
  int k = 0; // Keep the line of a 2nd derivative in v2
168
169
170
171
  for (second_derivatives_type::const_iterator it = second_derivatives.begin();
       it != second_derivatives.end(); it++)
    {
      int eq = it->first.first;
172
173
      int symb_id1 = it->first.second.first;
      int symb_id2 = it->first.second.second;
174
175
176
177
178
179
180
181
      NodeID d2 = it->second;

      int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
      int tsid2 = symbol_table.getTypeSpecificID(symb_id2);

      int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
      int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;

sebastien's avatar
sebastien committed
182
      hessian_output << "v2";
183
      hessianHelper(hessian_output, k, 0, output_type);
sebastien's avatar
sebastien committed
184
185
186
      hessian_output << "=" << eq + 1 << ";" << endl;

      hessian_output << "v2";
187
      hessianHelper(hessian_output, k, 1, output_type);
sebastien's avatar
sebastien committed
188
189
190
      hessian_output << "=" << col_nb + 1 << ";" << endl;

      hessian_output << "v2";
191
      hessianHelper(hessian_output, k, 2, output_type);
sebastien's avatar
sebastien committed
192
      hessian_output << "=";
193
194
195
      d2->writeOutput(hessian_output, output_type, temporary_terms);
      hessian_output << ";" << endl;

sebastien's avatar
sebastien committed
196
197
      k++;

198
199
200
      // Treating symetric elements
      if (symb_id1 != symb_id2)
        {
sebastien's avatar
sebastien committed
201
          hessian_output << "v2";
202
          hessianHelper(hessian_output, k, 0, output_type);
sebastien's avatar
sebastien committed
203
204
205
          hessian_output << "=" << eq + 1 << ";" << endl;

          hessian_output << "v2";
206
          hessianHelper(hessian_output, k, 1, output_type);
sebastien's avatar
sebastien committed
207
208
209
          hessian_output << "=" << col_nb_sym + 1 << ";" << endl;

          hessian_output << "v2";
210
          hessianHelper(hessian_output, k, 2, output_type);
sebastien's avatar
sebastien committed
211
          hessian_output << "=v2"; 
212
          hessianHelper(hessian_output, k-1, 2, output_type);
sebastien's avatar
sebastien committed
213
214
215
          hessian_output << ";" << endl;

          k++;
216
217
        }
    }
sebastien's avatar
sebastien committed
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

  // Writing ouputs
  if (mode != eDLLMode)
    {
      StaticOutput << "residual = zeros( " << equations.size() << ", 1);" << endl << endl
                   << "%" << endl
                   << "% Model equations" << endl
                   << "%" << endl
                   << endl
                   << model_output.str()
                   << "if ~isreal(residual)" << endl
                   << "  residual = real(residual)+imag(residual).^2;" << endl
                   << "end" << endl
                   << "if nargout >= 2," << endl
                   << "  g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");" << endl
                   << endl
                   << "%" << endl
                   << "% Jacobian matrix" << endl
                   << "%" << endl
                   << endl
                   << jacobian_output.str()
                   << "  if ~isreal(g1)" << endl
                   << "    g1 = real(g1)+2*imag(g1);" << endl
                   << "  end" << endl
                   << "end" << endl;
243

sebastien's avatar
sebastien committed
244
245
246
247
248
249
250
      // Initialize g2 matrix
      StaticOutput << "if nargout >= 3," << endl
                   << "%" << endl
                   << "% Hessian matrix" << endl
                   << "%" << endl
                   << endl;
      int ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
251
      if (second_derivatives.size())
sebastien's avatar
sebastien committed
252
253
254
255
256
257
        StaticOutput << "  v2 = zeros(" << NNZDerivatives[1] << ",3);" << endl
                     << hessian_output.str()
                     << "  g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << equations.size() << "," << ncols << ");" << endl;
      else // Either hessian is all zero, or we didn't compute it
        StaticOutput << "  g2 = sparse([],[],[]," << equations.size() << "," << ncols << ");" << endl;
      StaticOutput << "end;" << endl;
sebastien's avatar
sebastien committed
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
299
300
    }
  else
    {
      StaticOutput << "void Static(double *y, double *x, double *params, double *residual, double *g1)" << endl
                   << "{" << endl
                   << "  double lhs, rhs;" << endl
        // Writing residual equations
                   << "  /* Residual equations */" << endl
                   << "  if (residual == NULL)" << endl
                   << "    return;" << endl
                   << "  else" << endl
                   << "    {" << endl
                   << model_output.str()
        // Writing Jacobian
                   << "     /* Jacobian for endogenous variables without lag */" << endl
                   << "     if (g1 == NULL)" << endl
                   << "       return;" << endl
                   << "     else" << endl
                   << "       {" << endl
                   << jacobian_output.str()
                   << "       }" << endl
                   << "    }" << endl
                   << "}" << endl << endl;
    }
}

void
StaticModel::writeStaticFile(const string &basename) const
{
  switch (mode)
    {
    case eStandardMode:
    case eSparseDLLMode:
    case eSparseMode:
      writeStaticMFile(basename + "_static");
      break;
    case eDLLMode:
      writeStaticCFile(basename + "_static");
      break;
    }
}

void
301
StaticModel::computingPass(bool hessian, bool no_tmp_terms)
sebastien's avatar
sebastien committed
302
{
303
  // Compute derivatives w.r. to all endogenous
304
  set<int> vars;
305
306
  for(int i = 0; i < symbol_table.endo_nbr(); i++)
    vars.insert(symbol_table.getID(eEndogenous, i));
sebastien's avatar
sebastien committed
307
308

  // Launch computations
309
310
311
312
313
314
315
316
317
  cout << "Computing static model derivatives:" << endl
       << " - order 1" << endl;
  computeJacobian(vars);

  if (hessian)
    {
      cout << " - order 2" << endl;
      computeHessian(vars);
    }
sebastien's avatar
sebastien committed
318
319

  if (!no_tmp_terms)
320
    computeTemporaryTerms();
sebastien's avatar
sebastien committed
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339

  /*
  vector<int> endo2eq(equation_number());
  computeNormalization(endo2eq);

  multimap<int, int> natural_endo2eqs;
  computeNormalizedEquations(natural_endo2eqs);

  for(int i = 0; i < symbol_table.endo_nbr(); i++)
    {
      if (natural_endo2eqs.count(i) == 0)
        continue;

      pair<multimap<int, int>::const_iterator, multimap<int, int>::const_iterator> x = natural_endo2eqs.equal_range(i);
      if (find_if(x.first, x.second, compose1(bind2nd(equal_to<int>(), endo2eq[i]), select2nd<multimap<int, int>::value_type>())) == x.second)
        cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
             << " not used." << endl;
    }
  */
sebastien's avatar
sebastien committed
340
}
341
342
343
344

int
StaticModel::computeDerivID(int symb_id, int lag)
{
345
346
347
  if (symbol_table.getType(symb_id) == eEndogenous)
    return symb_id;
  else
348
349
350
351
352
353
    return -1;
}

int
StaticModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException)
{
354
355
  if (symbol_table.getType(symb_id) == eEndogenous)
    return symb_id;
sebastien's avatar
sebastien committed
356
  else
357
    throw UnknownDerivIDException();
358
}
359
360

void
sebastien's avatar
sebastien committed
361
StaticModel::computeNormalization(vector<int> &endo_to_eq) const
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
{
  int n = equation_number();

  assert(n == symbol_table.endo_nbr());

  typedef adjacency_list<vecS, vecS, undirectedS> BipartiteGraph;

  /*
    Vertices 0 to n-1 are for endogenous (using type specific ID)
    Vertices n to 2*n-1 are for equations (using equation no.)
  */
  BipartiteGraph g(2 * n);

  // Fill in the graph
  set<pair<int, int> > endo;
  for(int i = 0; i < n; i++)
    {
      endo.clear();
      equations[i]->collectEndogenous(endo);
      for(set<pair<int, int> >::const_iterator it = endo.begin();
          it != endo.end(); it++)
        add_edge(i + n, symbol_table.getTypeSpecificID(it->first), g);
    }

  // Compute maximum cardinality matching
sebastien's avatar
sebastien committed
387
388
  typedef vector<graph_traits<BipartiteGraph>::vertex_descriptor> mate_map_t;
  mate_map_t mate_map(2*n);
389
390
391
392
393

  bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);

  assert(check);

sebastien's avatar
sebastien committed
394
395
396
397
398
399
400
401
402
403
  // Check if all variables are normalized
  mate_map_t::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
  if (it != mate_map.begin() + n)
    {
      cerr << "ERROR: Could not normalize static model. Variable "
           << symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin()))
           << " is not in the maximum cardinality matching." << endl;
      exit(EXIT_FAILURE);
    }

404
  for(int i = 0; i < n; i++)
sebastien's avatar
sebastien committed
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
    cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
           << " matched with equation " << (mate_map[i]-n+1) << endl;

  assert((int) endo_to_eq.size() == n);

  // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
  transform(mate_map.begin(), mate_map.begin() + n, endo_to_eq.begin(), bind2nd(minus<int>(), n));
}

void
StaticModel::computeNormalizedEquations(multimap<int, int> &endo_to_eqs) const
{
  for(int i = 0; i < equation_number(); i++)
    {
      VariableNode *lhs = dynamic_cast<VariableNode *>(equations[i]->get_arg1());
      if (lhs == NULL)
        continue;

      int symb_id = lhs->get_symb_id();
      if (symbol_table.getType(symb_id) != eEndogenous)
        continue;

      set<pair<int, int> > endo;
      equations[i]->get_arg2()->collectEndogenous(endo);
      if (endo.find(make_pair(symb_id, 0)) != endo.end())
        continue;

      endo_to_eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i));
      cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl;
    }
435
}
436
437
438
439
440
441

void
StaticModel::writeLatexFile(const string &basename) const
{
  writeLatexModelFile(basename + "_static.tex", oLatexStaticModel);
}