StaticModel.cc 21.7 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>

23
#include <deque>
sebastien's avatar
sebastien committed
24
#include <algorithm>
25
#include <iterator>
sebastien's avatar
sebastien committed
26
27
#include <functional>

28
#ifdef DEBUG
29
30
31
32
33
// For select2nd()
# ifdef __GNUC__
#  include <ext/functional>
using namespace __gnu_cxx;
# endif
34
#endif
sebastien's avatar
sebastien committed
35

36
37
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
38
39
#include <boost/graph/strong_components.hpp>
#include <boost/graph/topological_sort.hpp>
40

sebastien's avatar
sebastien committed
41
#include "StaticModel.hh"
42
#include "MinimumFeedbackSet.hh"
sebastien's avatar
sebastien committed
43

44
45
using namespace boost;

sebastien's avatar
sebastien committed
46
47
StaticModel::StaticModel(SymbolTable &symbol_table_arg,
                         NumericalConstants &num_constants_arg) :
48
  ModelTree(symbol_table_arg, num_constants_arg)
sebastien's avatar
sebastien committed
49
50
51
52
{
}

void
53
StaticModel::writeStaticMFile(ostream &output, const string &func_name) const
sebastien's avatar
sebastien committed
54
55
{
  // Writing comments and function definition command
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
  output << "function [residual, g1, g2] = " << func_name << "(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
         << "residual = zeros( " << equations.size() << ", 1);" << endl
         << endl
         << "%" << endl
         << "% Model equations" << endl
         << "%" << endl
         << endl;

  writeModelLocalVariables(output, oMatlabStaticModel);

  writeTemporaryTerms(temporary_terms, output, oMatlabStaticModel);

  writeModelEquations(output, oMatlabStaticModel);

  output << "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;
sebastien's avatar
sebastien committed
86
87
88
89
90
91

  // 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;
92
      int symb_id = it->first.second;
sebastien's avatar
sebastien committed
93
94
      NodeID d1 = it->second;

95
      output << "  g1(" << eq+1 << "," << symbol_table.getTypeSpecificID(symb_id)+1 << ")=";
96
97
      d1->writeOutput(output, oMatlabStaticModel, temporary_terms);
      output << ";" << endl;
sebastien's avatar
sebastien committed
98
99
    }

100
101
102
103
104
105
106
107
108
109
110
111
  output << "  if ~isreal(g1)" << endl
         << "    g1 = real(g1)+2*imag(g1);" << endl
         << "  end" << endl
         << "end" << endl
         << "if nargout >= 3," << endl
         << "%" << endl
         << "% Hessian matrix" << endl
         << "%" << endl
         << endl;

  int g2ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
  if (second_derivatives.size())
112
    {
113
      output << "  v2 = zeros(" << NNZDerivatives[1] << ",3);" << endl;
114

115
116
117
118
      // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
      int k = 0; // Keep the line of a 2nd derivative in v2
      for (second_derivatives_type::const_iterator it = second_derivatives.begin();
           it != second_derivatives.end(); it++)
119
        {
120
121
122
123
124
125
126
          int eq = it->first.first;
          int symb_id1 = it->first.second.first;
          int symb_id2 = it->first.second.second;
          NodeID d2 = it->second;

          int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
          int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
sebastien's avatar
sebastien committed
127

128
129
          int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
          int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
sebastien's avatar
sebastien committed
130

131
132
133
134
135
          output << "v2(" << k+1 << ",1)=" << eq + 1 << ";" << endl
                 << "v2(" << k+1 << ",2)=" << col_nb + 1 << ";" << endl
                 << "v2(" << k+1 << ",3)=";
          d2->writeOutput(output, oMatlabStaticModel, temporary_terms);
          output << ";" << endl;
sebastien's avatar
sebastien committed
136
137

          k++;
138
139
140
141
142
143
144
145
146

          // Treating symetric elements
          if (symb_id1 != symb_id2)
            {
              output << "v2(" << k+1 << ",1)=" << eq + 1 << ";" << endl
                     << "v2(" << k+1 << ",2)=" << col_nb_sym + 1 << ";" << endl
                     << "v2(" << k+1 << ",3)=v2(" << k << ",3);" << endl;
              k++;
            }
147
        }
sebastien's avatar
sebastien committed
148

149
      output << "  g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << equations.size() << "," << g2ncols << ");" << endl;
sebastien's avatar
sebastien committed
150
    }
151
152
153
154
  else // Either hessian is all zero, or we didn't compute it
    output << "  g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");" << endl;

  output << "end;" << endl; // Close the if nargout >= 3 statement
sebastien's avatar
sebastien committed
155
156
157
}

void
158
StaticModel::writeStaticFile(const string &basename, bool block) const
sebastien's avatar
sebastien committed
159
{
160
161
162
163
164
  string filename = basename + "_static.m";

  ofstream output;
  output.open(filename.c_str(), ios::out | ios::binary);
  if (!output.is_open())
sebastien's avatar
sebastien committed
165
    {
166
167
      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
      exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
168
    }
169

170
  if (block)
171
172
173
    writeStaticBlockMFSFile(output, basename + "_static");
  else
    writeStaticMFile(output, basename + "_static");
174
175

  output.close();
sebastien's avatar
sebastien committed
176
177
178
}

void
179
StaticModel::computingPass(bool block, bool hessian, bool no_tmp_terms)
sebastien's avatar
sebastien committed
180
{
181

182

183
  // Compute derivatives w.r. to all endogenous
184
  set<int> vars;
185
186
  for(int i = 0; i < symbol_table.endo_nbr(); i++)
    vars.insert(symbol_table.getID(eEndogenous, i));
sebastien's avatar
sebastien committed
187
188

  // Launch computations
189
190
191
192
193
194
195
196
197
  cout << "Computing static model derivatives:" << endl
       << " - order 1" << endl;
  computeJacobian(vars);

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

199
  if (block)
sebastien's avatar
sebastien committed
200
    {
201
202
203
      computeNormalization();
      computeSortedBlockDecomposition();
      computeMFS();
204
      computeSortedRecursive();
205
      computeBlockMFSJacobian();
sebastien's avatar
sebastien committed
206
    }
207
208
  else if (!no_tmp_terms)
    computeTemporaryTerms(true);
sebastien's avatar
sebastien committed
209
}
210
211
212
213

int
StaticModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException)
{
214
215
  if (symbol_table.getType(symb_id) == eEndogenous)
    return symb_id;
sebastien's avatar
sebastien committed
216
  else
sebastien's avatar
sebastien committed
217
    return -1;
218
}
219
220

void
221
StaticModel::computeNormalization()
222
{
223
  const int n = equation_number();
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242

  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++)
243
        add_edge(i + n, it->first, g);
244
245
246
    }

  // Compute maximum cardinality matching
247
  vector<int> mate_map(2*n);
248

249
#if 1
250
  bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
251
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
#else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
  fill(mate_map.begin(), mate_map.end(), graph_traits<BipartiteGraph>::null_vertex());

  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;

      int j = natural_endo2eqs.find(i)->second;

      put(&mate_map[0], i, n+j);
      put(&mate_map[0], n+j, i);
    }

  edmonds_augmenting_path_finder<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type> augmentor(g, &mate_map[0], get(vertex_index, g));
  bool not_maximum_yet = true;
  while(not_maximum_yet)
    {
      not_maximum_yet = augmentor.augment_matching();
    }
  augmentor.get_current_matching(&mate_map[0]);

  bool check = maximum_cardinality_matching_verifier<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(vertex_index, g));
#endif
278
279
280

  assert(check);

sebastien's avatar
sebastien committed
281
  // Check if all variables are normalized
282
  vector<int>::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
sebastien's avatar
sebastien committed
283
284
285
286
287
288
289
290
  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);
    }

291
#ifdef DEBUG
292
  for(int i = 0; i < n; i++)
sebastien's avatar
sebastien committed
293
    cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
294
295
         << " matched with equation " << (mate_map[i]-n+1) << endl;
#endif
sebastien's avatar
sebastien committed
296
297

  // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
298
  endo2eq.resize(equation_number());
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
  transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));

#ifdef DEBUG
  multimap<int, int> natural_endo2eqs;
  computeNormalizedEquations(natural_endo2eqs);

  int n1 = 0, n2 = 0;

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

      n1++;

      pair<multimap<int, int>::const_iterator, multimap<int, int>::const_iterator> x = natural_endo2eqs.equal_range(i);
315
      if (find_if(x.first, x.second, compose1(bind2nd(equal_to<int>(), endo2eq[i]), select2nd<multimap<int, int>::value_type>())) == x.second)
316
317
318
319
320
321
322
323
        cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
             << " not used." << endl;
      else
        n2++;
    }

  cout << "Used " << n2 << " natural normalizations out of " << n1 << ", for a total of " << n << " equations." << endl;
#endif
sebastien's avatar
sebastien committed
324
325
326
}

void
327
StaticModel::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
sebastien's avatar
sebastien committed
328
329
330
331
332
333
334
335
336
337
338
339
340
{
  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);
341
      if (endo.find(make_pair(symbol_table.getTypeSpecificID(symb_id), 0)) != endo.end())
sebastien's avatar
sebastien committed
342
343
        continue;

344
      endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i));
sebastien's avatar
sebastien committed
345
346
      cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl;
    }
347
}
348
349
350
351
352
353

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

void
356
StaticModel::computeSortedBlockDecomposition()
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
{
  const int n = equation_number();

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

  // Compute graph representation of static model
  typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
  DirectedGraph g(n);

  set<pair<int, int> > endo;
  for(int i = 0; i < n; i++)
    {
      endo.clear();
      equations[endo2eq[i]]->collectEndogenous(endo);
      for(set<pair<int, int> >::const_iterator it = endo.begin();
          it != endo.end(); it++)
373
        add_edge(it->first, i, g);
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
    }

  // Compute strongly connected components
  vector<int> endo2block(n);
  int m = strong_components(g, &endo2block[0]);

  // Create directed acyclic graph associated to the strongly connected components
  DirectedGraph dag(m);
  graph_traits<DirectedGraph>::edge_iterator ei, ei_end;
  for(tie(ei, ei_end) = edges(g); ei != ei_end; ++ei)
    {
      int s = endo2block[source(*ei, g)];
      int t = endo2block[target(*ei, g)];
      if (s != t)
        add_edge(s, t, dag);
    }

  // Compute topological sort of DAG (ordered list of unordered SCC)
  deque<int> ordered2unordered;
  topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order
  // Construct mapping from unordered SCC to ordered SCC
  vector<int> unordered2ordered(m);
  for(int i = 0; i < m; i++)
    unordered2ordered[ordered2unordered[i]] = i;

  // Fill in data structure representing blocks
  blocks.clear();
  blocks.resize(m);
  for(int i = 0; i < n; i++)
    blocks[unordered2ordered[endo2block[i]]].insert(i);

405
#ifdef DEBUG
406
407
408
  cout << "Found " << m << " blocks" << endl;
  for(int i = 0; i < m; i++)
    cout << " Block " << i << " of size " << blocks[i].size() << endl;
409
#endif
410
411
412
}

void
413
StaticModel::computeMFS()
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
{
  const int n = equation_number();
  assert((int) endo2eq.size() == n);

  const int nblocks = blocks.size();
  blocksMFS.clear();
  blocksMFS.resize(nblocks);

  // Iterate over blocks
  for(int b = 0; b < nblocks; b++)
    {
      // Construct subgraph for MFS computation, where vertex number is position in the block
      int p = blocks[b].size();
      MFS::AdjacencyList_type g(p);

      // Construct v_index and v_index1 properties, and a mapping between type specific IDs and vertex descriptors
      property_map<MFS::AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, g);
      property_map<MFS::AdjacencyList_type, vertex_index1_t>::type v_index1 = get(vertex_index1, g);
      map<int, graph_traits<MFS::AdjacencyList_type>::vertex_descriptor> tsid2vertex;
      int j = 0;
      for(set<int>::const_iterator it = blocks[b].begin(); it != blocks[b].end(); ++it)
        {
          tsid2vertex[*it] = vertex(j, g);
          put(v_index, vertex(j, g), *it);
          put(v_index1, vertex(j, g), *it);
          j++;
        }

      // Add edges, loop over endogenous in the block
      set<pair<int, int> > endo;
      for(set<int>::const_iterator it = blocks[b].begin(); it != blocks[b].end(); ++it)
        {
          endo.clear();

          // Test if associated equation is in normalized form, and compute set of endogenous appearing in it
          ExprNode *lhs = equations[endo2eq[*it]]->get_arg1();
          VariableNode *lhs_var = dynamic_cast<VariableNode *>(lhs);
          if (lhs_var == NULL || lhs_var->get_symb_id() != symbol_table.getID(eEndogenous, *it))
            lhs->collectEndogenous(endo); // Only collect endogenous of LHS if not normalized form
          ExprNode *rhs = equations[endo2eq[*it]]->get_arg2();
          rhs->collectEndogenous(endo);

          for(set<pair<int, int> >::const_iterator it2 = endo.begin();
              it2 != endo.end(); ++it2)
458
459
            if (blocks[b].find(it2->first) != blocks[b].end()) // Add edge only if vertex member of this block
              add_edge(tsid2vertex[it2->first], tsid2vertex[*it], g);
460
461
462
463
464
        }

      // Compute minimum feedback set
      MFS::Minimal_set_of_feedback_vertex(blocksMFS[b], g);

465
#ifdef DEBUG
466
      cout << "Block " << b << ": " << blocksMFS[b].size() << "/" << blocks[b].size() << " in MFS" << endl;
467
#endif
468
469
    }
}
470
471

void
472
StaticModel::computeSortedRecursive()
473
{
474
475
476
477
478
  const int nblocks = blocks.size();
  blocksRecursive.clear();
  blocksRecursive.resize(nblocks);

  for(int b = 0; b < nblocks; b++)
479
    {
480
481
482
      // Construct the set of recursive vars
      // The index in this vector will be the vertex descriptor in the graph
      vector<int> recurs_vars;
483
484
      set_difference(blocks[b].begin(), blocks[b].end(),
                     blocksMFS[b].begin(), blocksMFS[b].end(),
485
                     back_inserter(recurs_vars));
486

487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
      // Construct graph representation of recursive vars
      typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
      DirectedGraph dag(recurs_vars.size());
      set<pair<int, int> > endo;
      for(int i = 0; i < (int) recurs_vars.size(); i++)
        {
          endo.clear();
          equations[endo2eq[recurs_vars[i]]]->get_arg2()->collectEndogenous(endo);
          for(set<pair<int, int> >::const_iterator it = endo.begin();
              it != endo.end(); it++)
            {
              vector<int>::const_iterator it2 = find(recurs_vars.begin(), recurs_vars.end(), it->first);
              if (it2 != recurs_vars.end())
                {
                  int source_vertex = it2 - recurs_vars.begin();
                  add_edge(source_vertex, i, dag);
                }
            }
        }
      // Compute topological sort
      deque<int> ordered_recurs_vertices;
      topological_sort(dag, front_inserter(ordered_recurs_vertices)); // We use a front inserter because topological_sort returns the inverse order

      // Construct the final order
      for(deque<int>::const_iterator it = ordered_recurs_vertices.begin();
          it != ordered_recurs_vertices.end(); it++)
        blocksRecursive[b].push_back(recurs_vars[*it]);
    }
}

void
StaticModel::computeBlockMFSJacobian()
{
  blocksMFSJacobian.clear();
  for(int b = 0; b < (int) blocks.size(); b++)
    {
      // Create the map of recursive vars to their normalized equation
      map<int, NodeID> recursive2eq;
      for(vector<int>::const_iterator it = blocksRecursive[b].begin();
          it != blocksRecursive[b].end(); it++)
        recursive2eq[symbol_table.getID(eEndogenous, *it)] = equations[endo2eq[*it]];
528
529
530
531
532

      for(set<int>::const_iterator it = blocksMFS[b].begin();
          it != blocksMFS[b].end(); it++)
        {
          int eq_no = endo2eq[*it];
533
534
535
536
          for(set<int>::const_iterator it2 = blocksMFS[b].begin();
              it2 != blocksMFS[b].end(); it2++)
            {
              int deriv_id = symbol_table.getID(eEndogenous, *it2);
537
              NodeID d = equations[eq_no]->getChainRuleDerivative(deriv_id, recursive2eq);
538
539
540
              if (d != Zero)
                blocksMFSJacobian[make_pair(eq_no, *it2)] = d;
            }
541
542
543
544
545
        }
    }
}

void
546
StaticModel::writeOutput(ostream &output, bool block) const
547
{
548
  if (!block)
549
    return;
550

551
  output << "M_.blocksMFS = cell(" << blocksMFS.size() << ", 1);" << endl;
552
553
  for(int b = 0; b < (int) blocks.size(); b++)
    {
554
555
      output << "M_.blocksMFS{" << b+1 << "} = [ ";
      transform(blocksMFS[b].begin(), blocksMFS[b].end(), ostream_iterator<int>(output, "; "), bind2nd(plus<int>(), 1));
556
557
558
      output << "];" << endl;
    }
}
559

560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
void
StaticModel::writeLocalVars(ostream &output, NodeID expr, set<int> &local_var_written) const
{
  set<int> expr_local_var;
  expr->collectModelLocalVariables(expr_local_var);
  vector<int> new_local_var;
  set_difference(expr_local_var.begin(), expr_local_var.end(),
                 local_var_written.begin(), local_var_written.end(),
                 back_inserter(new_local_var));

  for(vector<int>::const_iterator it = new_local_var.begin();
      it != new_local_var.end(); it++)
    {
      output << symbol_table.getName(*it) << " = ";
      map<int, NodeID>::const_iterator it2 = local_variables_table.find(*it);
      it2->second->writeOutput(output, oMatlabStaticModel, temporary_terms_type());
      output << ";" << endl;
      local_var_written.insert(*it);
    }
}

581
582
583
void
StaticModel::writeStaticBlockMFSFile(ostream &output, const string &func_name) const
{
584
  output << "function [residual, g1, y] = " << func_name << "(nblock, y, x, params)" << endl
585
586
587
588
         << "  switch nblock" << endl;

  for(int b = 0; b < (int) blocks.size(); b++)
    {
589
590
      set<int> local_var;

591
592
      output << "    case " << b+1 << endl
             << "      % Variables not in minimum feedback set" << endl;
593
594
      for(vector<int>::const_iterator it = blocksRecursive[b].begin();
          it != blocksRecursive[b].end(); it++)
595
        {
596
597
598
599
600
          NodeID eq = equations[endo2eq[*it]];

          writeLocalVars(output, eq, local_var);

          eq->writeOutput(output, oMatlabStaticModel, temporary_terms_type());
601
602
603
604
605
606
607
608
609
610
          output << ";" << endl;
        }

      output << "      % Model residuals" << endl
             << "residual = zeros(" << blocksMFS[b].size() << ",1);" << endl;

      int i = 1;
      for (set<int>::const_iterator it = blocksMFS[b].begin();
           it != blocksMFS[b].end(); it++)
        {
611
612
613
          BinaryOpNode *eq = equations[endo2eq[*it]];

          writeLocalVars(output, eq, local_var);
614

615
          output << "residual(" << i << ")=(";
616

617
          NodeID lhs = eq->get_arg1();
618
619
620
          lhs->writeOutput(output, oMatlabStaticModel, temporary_terms_type());
          output << ")-(";

621
          NodeID rhs = eq->get_arg2();
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
          rhs->writeOutput(output, oMatlabStaticModel, temporary_terms_type());
          output << ");" << endl;

          i++;
        }

      output << "      % Jacobian matrix" << endl
             << "g1 = zeros(" << blocksMFS[b].size() << ", " << blocksMFS[b].size() << ");" << endl;
      i = 1;
      for (set<int>::const_iterator it = blocksMFS[b].begin();
           it != blocksMFS[b].end(); it++)
        {
          int eq_no = endo2eq[*it];
          int j = 1;
          for (set<int>::const_iterator it2 = blocksMFS[b].begin();
               it2 != blocksMFS[b].end(); it2++)
            {
              map<pair<int, int>, NodeID>::const_iterator itd = blocksMFSJacobian.find(make_pair(eq_no, *it2));
              if (itd != blocksMFSJacobian.end())
                {
                  output << "g1(" << i << "," << j << ")=";
                  itd->second->writeOutput(output, oMatlabStaticModel, temporary_terms_type());
                  output << ";" << endl;
                }
              j++;
            }
          i++;
        }
    }
651
652
653

  output << "  end" << endl
         << "end" << endl;
654
}
sebastien's avatar
sebastien committed
655
656
657
658
659
660
661
662
663
664

void
StaticModel::writeAuxVarInitval(ostream &output) const
{
  for(int i = 0; i < (int) aux_equations.size(); i++)
    {
      dynamic_cast<ExprNode *>(aux_equations[i])->writeOutput(output);
      output << ";" << endl;
    }
}