1 /*
2 * Copyright © 2003-2020 Dynare Team
3 *
4 * This file is part of Dynare.
5 *
6 * Dynare is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * Dynare is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 #include "ModelTree.hh"
21 #include "MinimumFeedbackSet.hh"
22 #pragma GCC diagnostic push
23 #pragma GCC diagnostic ignored "-Wold-style-cast"
24 #pragma GCC diagnostic ignored "-Wsign-compare"
25 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
26 #include <boost/graph/adjacency_list.hpp>
27 #include <boost/graph/max_cardinality_matching.hpp>
28 #include <boost/graph/strong_components.hpp>
29 #include <boost/graph/topological_sort.hpp>
30 #pragma GCC diagnostic pop
31
32 #ifdef __APPLE__
33 # include <mach-o/dyld.h>
34 #endif
35
36 #include <regex>
37
38 using namespace MFS;
39
40 void
copyHelper(const ModelTree & m)41 ModelTree::copyHelper(const ModelTree &m)
42 {
43 auto f = [this](expr_t e) { return e->clone(*this); };
44
45 // Equations
46 for (const auto &it : m.equations)
47 equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
48 for (const auto &it : m.aux_equations)
49 aux_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
50
51 auto convert_deriv_map = [f](map<vector<int>, expr_t> dm)
52 {
53 map<vector<int>, expr_t> dm2;
54 for (const auto &it : dm)
55 dm2.emplace(it.first, f(it.second));
56 return dm2;
57 };
58
59 // Derivatives
60 for (const auto &it : m.derivatives)
61 derivatives.push_back(convert_deriv_map(it));
62 for (const auto &it : m.params_derivatives)
63 params_derivatives[it.first] = convert_deriv_map(it.second);
64
65 auto convert_temporary_terms_t = [f](temporary_terms_t tt)
66 {
67 temporary_terms_t tt2;
68 for (const auto &it : tt)
69 tt2.insert(f(it));
70 return tt2;
71 };
72
73 // Temporary terms
74 for (const auto &it : m.temporary_terms)
75 temporary_terms.insert(f(it));
76 for (const auto &it : m.temporary_terms_mlv)
77 temporary_terms_mlv[dynamic_cast<VariableNode *>(f(it.first))] = f(it.second);
78 for (const auto &it : m.temporary_terms_derivatives)
79 temporary_terms_derivatives.push_back(convert_temporary_terms_t(it));
80 for (const auto &it : m.temporary_terms_idxs)
81 temporary_terms_idxs[f(it.first)] = it.second;
82 for (const auto &it : m.params_derivs_temporary_terms)
83 params_derivs_temporary_terms[it.first] = convert_temporary_terms_t(it.second);
84 for (const auto &it : m.params_derivs_temporary_terms_idxs)
85 params_derivs_temporary_terms_idxs[f(it.first)] = it.second;
86
87 // Other stuff
88 for (const auto &it : m.trend_symbols_map)
89 trend_symbols_map[it.first] = f(it.second);
90 for (const auto &it : m.nonstationary_symbols_map)
91 nonstationary_symbols_map[it.first] = {it.second.first, f(it.second.second)};
92 }
93
ModelTree(SymbolTable & symbol_table_arg,NumericalConstants & num_constants_arg,ExternalFunctionsTable & external_functions_table_arg,bool is_dynamic_arg)94 ModelTree::ModelTree(SymbolTable &symbol_table_arg,
95 NumericalConstants &num_constants_arg,
96 ExternalFunctionsTable &external_functions_table_arg,
97 bool is_dynamic_arg) :
98 DataTree{symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg},
99 derivatives(4),
100 NNZDerivatives(4, 0),
101 temporary_terms_derivatives(4)
102 {
103 }
104
ModelTree(const ModelTree & m)105 ModelTree::ModelTree(const ModelTree &m) :
106 DataTree{m},
107 user_set_add_flags{m.user_set_add_flags},
108 user_set_subst_flags{m.user_set_subst_flags},
109 user_set_add_libs{m.user_set_add_libs},
110 user_set_subst_libs{m.user_set_subst_libs},
111 user_set_compiler{m.user_set_compiler},
112 equations_lineno{m.equations_lineno},
113 equation_tags{m.equation_tags},
114 equation_tags_xref{m.equation_tags_xref},
115 computed_derivs_order{m.computed_derivs_order},
116 NNZDerivatives{m.NNZDerivatives},
117 equation_reordered{m.equation_reordered},
118 variable_reordered{m.variable_reordered},
119 inv_equation_reordered{m.inv_equation_reordered},
120 inv_variable_reordered{m.inv_variable_reordered},
121 is_equation_linear{m.is_equation_linear},
122 endo2eq{m.endo2eq},
123 epilogue{m.epilogue},
124 prologue{m.prologue},
125 block_lag_lead{m.block_lag_lead},
126 cutoff{m.cutoff},
127 mfs{m.mfs}
128 {
129 copyHelper(m);
130 }
131
132 ModelTree &
operator =(const ModelTree & m)133 ModelTree::operator=(const ModelTree &m)
134 {
135 DataTree::operator=(m);
136
137 equations.clear();
138 equations_lineno = m.equations_lineno;
139 aux_equations.clear();
140 equation_tags = m.equation_tags;
141 equation_tags_xref = m.equation_tags_xref;
142 computed_derivs_order = m.computed_derivs_order;
143 NNZDerivatives = m.NNZDerivatives;
144
145 derivatives.clear();
146 params_derivatives.clear();
147
148 temporary_terms.clear();
149 temporary_terms_mlv.clear();
150 temporary_terms_derivatives.clear();
151 params_derivs_temporary_terms.clear();
152 params_derivs_temporary_terms_idxs.clear();
153
154 trend_symbols_map.clear();
155 nonstationary_symbols_map.clear();
156
157 equation_reordered = m.equation_reordered;
158 variable_reordered = m.variable_reordered;
159 inv_equation_reordered = m.inv_equation_reordered;
160 inv_variable_reordered = m.inv_variable_reordered;
161 is_equation_linear = m.is_equation_linear;
162 endo2eq = m.endo2eq;
163 epilogue = m.epilogue;
164 prologue = m.prologue;
165 block_lag_lead = m.block_lag_lead;
166 cutoff = m.cutoff;
167 mfs = m.mfs;
168
169 user_set_add_flags = m.user_set_add_flags;
170 user_set_subst_flags = m.user_set_subst_flags;
171 user_set_add_libs = m.user_set_add_libs;
172 user_set_subst_libs = m.user_set_subst_libs;
173 user_set_compiler = m.user_set_compiler;
174
175 copyHelper(m);
176
177 return *this;
178 }
179
180 bool
computeNormalization(const jacob_map_t & contemporaneous_jacobian,bool verbose)181 ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
182 {
183 const int n = equations.size();
184
185 assert(n == symbol_table.endo_nbr());
186
187 using BipartiteGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>;
188
189 /*
190 Vertices 0 to n-1 are for endogenous (using type specific ID)
191 Vertices n to 2*n-1 are for equations (using equation no.)
192 */
193 BipartiteGraph g(2 * n);
194
195 // Fill in the graph
196 set<pair<int, int>> endo;
197
198 for (const auto &it : contemporaneous_jacobian)
199 add_edge(it.first.first + n, it.first.second, g);
200
201 // Compute maximum cardinality matching
202 vector<int> mate_map(2*n);
203
204 #if 1
205 bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
206 #else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
207 fill(mate_map.begin(), mate_map.end(), boost::graph_traits<BipartiteGraph>::null_vertex());
208
209 auto natural_endo2eqs = computeNormalizedEquations();
210
211 for (int i = 0; i < symbol_table.endo_nbr(); i++)
212 {
213 if (natural_endo2eqs.count(i) == 0)
214 continue;
215
216 int j = natural_endo2eqs.find(i)->second;
217
218 put(&mate_map[0], i, n+j);
219 put(&mate_map[0], n+j, i);
220 }
221
222 boost::edmonds_augmenting_path_finder<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type> augmentor(g, &mate_map[0], get(boost::vertex_index, g));
223 while (augmentor.augment_matching())
224 {
225 };
226
227 augmentor.get_current_matching(&mate_map[0]);
228
229 bool check = boost::maximum_cardinality_matching_verifier<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(boost::vertex_index, g));
230 #endif
231
232 assert(check);
233
234 #ifdef DEBUG
235 for (int i = 0; i < n; i++)
236 cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
237 << " matched with equation " << (mate_map[i]-n+1) << endl;
238 #endif
239
240 // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
241 endo2eq.resize(equations.size());
242 transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), [=](int i) { return i-n; });
243
244 #ifdef DEBUG
245 auto natural_endo2eqs = computeNormalizedEquations(natural_endo2eqs);
246
247 int n1 = 0, n2 = 0;
248
249 for (int i = 0; i < symbol_table.endo_nbr(); i++)
250 {
251 if (natural_endo2eqs.count(i) == 0)
252 continue;
253
254 n1++;
255
256 auto x = natural_endo2eqs.equal_range(i);
257 if (find_if(x.first, x.second, [=](auto y) { return y.second == endo2eq[i]; }) == x.second)
258 cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, i))
259 << " not used." << endl;
260 else
261 n2++;
262 }
263
264 cout << "Used " << n2 << " natural normalizations out of " << n1 << ", for a total of " << n << " equations." << endl;
265 #endif
266
267 // Check if all variables are normalized
268 if (auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
269 it != mate_map.begin() + n)
270 {
271 if (verbose)
272 cerr << "ERROR: Could not normalize the model. Variable "
273 << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
274 << " is not in the maximum cardinality matching." << endl;
275 check = false;
276 }
277 return check;
278 }
279
280 void
computeNonSingularNormalization(jacob_map_t & contemporaneous_jacobian,double cutoff,jacob_map_t & static_jacobian,dynamic_jacob_map_t & dynamic_jacobian)281 ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian)
282 {
283 bool check = false;
284
285 cout << "Normalizing the model..." << endl;
286
287 int n = equations.size();
288
289 // compute the maximum value of each row of the contemporaneous Jacobian matrix
290 //jacob_map normalized_contemporaneous_jacobian;
291 jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
292 vector<double> max_val(n, 0.0);
293 for (const auto &it : contemporaneous_jacobian)
294 if (fabs(it.second) > max_val[it.first.first])
295 max_val[it.first.first] = fabs(it.second);
296
297 for (auto &iter : normalized_contemporaneous_jacobian)
298 iter.second /= max_val[iter.first.first];
299
300 //We start with the highest value of the cutoff and try to normalize the model
301 double current_cutoff = 0.99999999;
302
303 int suppressed = 0;
304 while (!check && current_cutoff > 1e-19)
305 {
306 jacob_map_t tmp_normalized_contemporaneous_jacobian;
307 int suppress = 0;
308 for (auto &iter : normalized_contemporaneous_jacobian)
309 if (fabs(iter.second) > max(current_cutoff, cutoff))
310 tmp_normalized_contemporaneous_jacobian[{ iter.first.first, iter.first.second }] = iter.second;
311 else
312 suppress++;
313
314 if (suppress != suppressed)
315 check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false);
316 suppressed = suppress;
317 if (!check)
318 {
319 current_cutoff /= 2;
320 // In this last case try to normalize with the complete jacobian
321 if (current_cutoff <= 1e-19)
322 check = computeNormalization(normalized_contemporaneous_jacobian, false);
323 }
324 }
325
326 if (!check)
327 {
328 cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
329 //if no non-singular normalization can be found, try to find a normalization even with a potential singularity
330 jacob_map_t tmp_normalized_contemporaneous_jacobian;
331 set<pair<int, int>> endo;
332 for (int i = 0; i < n; i++)
333 {
334 endo.clear();
335 equations[i]->collectEndogenous(endo);
336 for (const auto &it : endo)
337 tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
338 }
339 check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
340 if (check)
341 {
342 // Update the jacobian matrix
343 for (const auto &[key, ignore] : tmp_normalized_contemporaneous_jacobian)
344 {
345 if (static_jacobian.find({ key.first, key.second }) == static_jacobian.end())
346 static_jacobian[{ key.first, key.second }] = 0;
347 if (dynamic_jacobian.find({ 0, key.first, key.second }) == dynamic_jacobian.end())
348 dynamic_jacobian[{ 0, key.first, key.second }] = nullptr;
349 if (contemporaneous_jacobian.find({ key.first, key.second }) == contemporaneous_jacobian.end())
350 contemporaneous_jacobian[{ key.first, key.second }] = 0;
351 try
352 {
353 if (derivatives[1].find({ key.first, getDerivID(symbol_table.getID(SymbolType::endogenous, key.second), 0) }) == derivatives[1].end())
354 derivatives[1][{ key.first, getDerivID(symbol_table.getID(SymbolType::endogenous, key.second), 0) }] = Zero;
355 }
356 catch (DataTree::UnknownDerivIDException &e)
357 {
358 cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, key.second))
359 << " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block." << endl;
360 exit(EXIT_FAILURE);
361 }
362 }
363 }
364 }
365
366 if (!check)
367 {
368 cerr << "No normalization could be computed. Aborting." << endl;
369 exit(EXIT_FAILURE);
370 }
371 }
372
373 multimap<int, int>
computeNormalizedEquations() const374 ModelTree::computeNormalizedEquations() const
375 {
376 multimap<int, int> endo2eqs;
377 for (size_t i = 0; i < equations.size(); i++)
378 {
379 auto lhs = dynamic_cast<VariableNode *>(equations[i]->arg1);
380 if (!lhs)
381 continue;
382
383 int symb_id = lhs->symb_id;
384 if (symbol_table.getType(symb_id) != SymbolType::endogenous)
385 continue;
386
387 set<pair<int, int>> endo;
388 equations[i]->arg2->collectEndogenous(endo);
389 if (endo.find({ symbol_table.getTypeSpecificID(symb_id), 0 }) != endo.end())
390 continue;
391
392 endo2eqs.emplace(symbol_table.getTypeSpecificID(symb_id), static_cast<int>(i));
393 cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << i+1 << endl;
394 }
395 return endo2eqs;
396 }
397
398 void
evaluateAndReduceJacobian(const eval_context_t & eval_context,jacob_map_t & contemporaneous_jacobian,jacob_map_t & static_jacobian,dynamic_jacob_map_t & dynamic_jacobian,double cutoff,bool verbose)399 ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose)
400 {
401 int nb_elements_contemparenous_Jacobian = 0;
402 set<vector<int>> jacobian_elements_to_delete;
403 for (const auto &[indices, d1] : derivatives[1])
404 {
405 int deriv_id = indices[1];
406 if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
407 {
408 int eq = indices[0];
409 int symb = getSymbIDByDerivID(deriv_id);
410 int var = symbol_table.getTypeSpecificID(symb);
411 int lag = getLagByDerivID(deriv_id);
412 double val = 0;
413 try
414 {
415 val = d1->eval(eval_context);
416 }
417 catch (ExprNode::EvalExternalFunctionException &e)
418 {
419 val = 1;
420 }
421 catch (ExprNode::EvalException &e)
422 {
423 cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " (line " << equations_lineno[eq] << ") and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
424 d1->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModelSparse, temporary_terms, {});
425 cerr << endl;
426 exit(EXIT_FAILURE);
427 }
428 if (fabs(val) < cutoff)
429 {
430 if (verbose)
431 cout << "the coefficient related to variable " << var << " with lag " << lag << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")" << endl;
432 jacobian_elements_to_delete.insert({ eq, deriv_id });
433 }
434 else
435 {
436 if (lag == 0)
437 {
438 nb_elements_contemparenous_Jacobian++;
439 contemporaneous_jacobian[{ eq, var }] = val;
440 }
441 if (static_jacobian.find({ eq, var }) != static_jacobian.end())
442 static_jacobian[{ eq, var }] += val;
443 else
444 static_jacobian[{ eq, var }] = val;
445 dynamic_jacobian[{ lag, eq, var }] = d1;
446 }
447 }
448 }
449
450 // Get rid of the elements of the Jacobian matrix below the cutoff
451 for (const auto &it : jacobian_elements_to_delete)
452 derivatives[1].erase(it);
453
454 if (jacobian_elements_to_delete.size() > 0)
455 {
456 cout << jacobian_elements_to_delete.size() << " elements among " << derivatives[1].size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded" << endl
457 << "The contemporaneous incidence matrix has " << nb_elements_contemparenous_Jacobian << " elements" << endl;
458 }
459 }
460
461 vector<pair<int, int>>
select_non_linear_equations_and_variables(vector<bool> is_equation_linear,const dynamic_jacob_map_t & dynamic_jacobian,vector<int> & equation_reordered,vector<int> & variable_reordered,vector<int> & inv_equation_reordered,vector<int> & inv_variable_reordered,lag_lead_vector_t & equation_lag_lead,lag_lead_vector_t & variable_lag_lead,vector<unsigned int> & n_static,vector<unsigned int> & n_forward,vector<unsigned int> & n_backward,vector<unsigned int> & n_mixed)462 ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_linear, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered,
463 vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered,
464 lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead,
465 vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed)
466 {
467 vector<int> eq2endo(equations.size(), 0);
468 /*equation_reordered.resize(equations.size());
469 variable_reordered.resize(equations.size());*/
470 unsigned int num = 0;
471 for (auto it : endo2eq)
472 if (!is_equation_linear[it])
473 num++;
474 vector<int> endo2block = vector<int>(endo2eq.size(), 1);
475 vector<pair<set<int>, pair<set<int>, vector<int>>>> components_set(num);
476 int i = 0, j = 0;
477 for (auto it : endo2eq)
478 if (!is_equation_linear[it])
479 {
480 equation_reordered[i] = it;
481 variable_reordered[i] = j;
482 endo2block[j] = 0;
483 components_set[endo2block[j]].first.insert(i);
484 i++;
485 j++;
486 }
487 getVariableLeadLagByBlock(dynamic_jacobian, endo2block, endo2block.size(), equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
488 n_static = vector<unsigned int>(endo2eq.size(), 0);
489 n_forward = vector<unsigned int>(endo2eq.size(), 0);
490 n_backward = vector<unsigned int>(endo2eq.size(), 0);
491 n_mixed = vector<unsigned int>(endo2eq.size(), 0);
492 for (unsigned int i = 0; i < endo2eq.size(); i++)
493 {
494 if (variable_lag_lead[variable_reordered[i]].first != 0 && variable_lag_lead[variable_reordered[i]].second != 0)
495 n_mixed[i]++;
496 else if (variable_lag_lead[variable_reordered[i]].first == 0 && variable_lag_lead[variable_reordered[i]].second != 0)
497 n_forward[i]++;
498 else if (variable_lag_lead[variable_reordered[i]].first != 0 && variable_lag_lead[variable_reordered[i]].second == 0)
499 n_backward[i]++;
500 else if (variable_lag_lead[variable_reordered[i]].first == 0 && variable_lag_lead[variable_reordered[i]].second == 0)
501 n_static[i]++;
502 }
503 cout.flush();
504 int nb_endo = is_equation_linear.size();
505 vector<pair<int, int>> blocks(1, {i, i});
506 inv_equation_reordered.resize(nb_endo);
507 inv_variable_reordered.resize(nb_endo);
508 for (int i = 0; i < nb_endo; i++)
509 {
510 inv_variable_reordered[variable_reordered[i]] = i;
511 inv_equation_reordered[equation_reordered[i]] = i;
512 }
513 return blocks;
514 }
515
516 bool
computeNaturalNormalization()517 ModelTree::computeNaturalNormalization()
518 {
519 bool bool_result = true;
520 set<pair<int, int>> result;
521 endo2eq.resize(equations.size());
522 for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
523 if (!is_equation_linear[eq])
524 {
525 BinaryOpNode *eq_node = equations[eq];
526 expr_t lhs = eq_node->arg1;
527 result.clear();
528 lhs->collectDynamicVariables(SymbolType::endogenous, result);
529 if (result.size() == 1 && result.begin()->second == 0)
530 {
531 //check if the endogenous variable has not been already used in an other match !
532 if (find(endo2eq.begin(), endo2eq.end(), result.begin()->first) == endo2eq.end())
533 endo2eq[result.begin()->first] = eq;
534 else
535 {
536 bool_result = false;
537 break;
538 }
539 }
540 }
541 return bool_result;
542 }
543
544 void
computePrologueAndEpilogue(const jacob_map_t & static_jacobian_arg,vector<int> & equation_reordered,vector<int> & variable_reordered)545 ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered)
546 {
547 vector<int> eq2endo(equations.size(), 0);
548 equation_reordered.resize(equations.size());
549 variable_reordered.resize(equations.size());
550 int n = equations.size();
551 vector<bool> IM(n*n);
552 int i = 0;
553 for (auto it : endo2eq)
554 {
555 eq2endo[it] = i;
556 equation_reordered[i] = i;
557 variable_reordered[it] = i;
558 i++;
559 }
560 if (cutoff == 0)
561 {
562 set<pair<int, int>> endo;
563 for (int i = 0; i < n; i++)
564 {
565 endo.clear();
566 equations[i]->collectEndogenous(endo);
567 for (const auto &it : endo)
568 IM[i * n + endo2eq[it.first]] = true;
569 }
570 }
571 else
572 for (const auto &it : static_jacobian_arg)
573 IM[it.first.first * n + endo2eq[it.first.second]] = true;
574 bool something_has_been_done = true;
575 prologue = 0;
576 int k = 0;
577 // Find the prologue equations and place first the AR(1) shock equations first
578 while (something_has_been_done)
579 {
580 int tmp_prologue = prologue;
581 something_has_been_done = false;
582 for (int i = prologue; i < n; i++)
583 {
584 int nze = 0;
585 for (int j = tmp_prologue; j < n; j++)
586 if (IM[i * n + j])
587 {
588 nze++;
589 k = j;
590 }
591 if (nze == 1)
592 {
593 for (int j = 0; j < n; j++)
594 {
595 bool tmp_bool = IM[tmp_prologue * n + j];
596 IM[tmp_prologue * n + j] = IM[i * n + j];
597 IM[i * n + j] = tmp_bool;
598 }
599 int tmp = equation_reordered[tmp_prologue];
600 equation_reordered[tmp_prologue] = equation_reordered[i];
601 equation_reordered[i] = tmp;
602 for (int j = 0; j < n; j++)
603 {
604 bool tmp_bool = IM[j * n + tmp_prologue];
605 IM[j * n + tmp_prologue] = IM[j * n + k];
606 IM[j * n + k] = tmp_bool;
607 }
608 tmp = variable_reordered[tmp_prologue];
609 variable_reordered[tmp_prologue] = variable_reordered[k];
610 variable_reordered[k] = tmp;
611 tmp_prologue++;
612 something_has_been_done = true;
613 }
614 }
615 prologue = tmp_prologue;
616 }
617
618 something_has_been_done = true;
619 epilogue = 0;
620 // Find the epilogue equations
621 while (something_has_been_done)
622 {
623 int tmp_epilogue = epilogue;
624 something_has_been_done = false;
625 for (int i = prologue; i < n - static_cast<int>(epilogue); i++)
626 {
627 int nze = 0;
628 for (int j = prologue; j < n - tmp_epilogue; j++)
629 if (IM[j * n + i])
630 {
631 nze++;
632 k = j;
633 }
634 if (nze == 1)
635 {
636 for (int j = 0; j < n; j++)
637 {
638 bool tmp_bool = IM[(n - 1 - tmp_epilogue) * n + j];
639 IM[(n - 1 - tmp_epilogue) * n + j] = IM[k * n + j];
640 IM[k * n + j] = tmp_bool;
641 }
642 int tmp = equation_reordered[n - 1 - tmp_epilogue];
643 equation_reordered[n - 1 - tmp_epilogue] = equation_reordered[k];
644 equation_reordered[k] = tmp;
645 for (int j = 0; j < n; j++)
646 {
647 bool tmp_bool = IM[j * n + n - 1 - tmp_epilogue];
648 IM[j * n + n - 1 - tmp_epilogue] = IM[j * n + i];
649 IM[j * n + i] = tmp_bool;
650 }
651 tmp = variable_reordered[n - 1 - tmp_epilogue];
652 variable_reordered[n - 1 - tmp_epilogue] = variable_reordered[i];
653 variable_reordered[i] = tmp;
654 tmp_epilogue++;
655 something_has_been_done = true;
656 }
657 }
658 epilogue = tmp_epilogue;
659 }
660 }
661
662 equation_type_and_normalized_equation_t
equationTypeDetermination(const map<tuple<int,int,int>,expr_t> & first_order_endo_derivatives,const vector<int> & Index_Var_IM,const vector<int> & Index_Equ_IM,int mfs) const663 ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
664 {
665 expr_t lhs;
666 BinaryOpNode *eq_node;
667 EquationType Equation_Simulation_Type;
668 equation_type_and_normalized_equation_t V_Equation_Simulation_Type(equations.size());
669 for (unsigned int i = 0; i < equations.size(); i++)
670 {
671 int eq = Index_Equ_IM[i];
672 int var = Index_Var_IM[i];
673 eq_node = equations[eq];
674 lhs = eq_node->arg1;
675 Equation_Simulation_Type = E_SOLVE;
676 auto derivative = first_order_endo_derivatives.find({ eq, var, 0 });
677 pair<bool, expr_t> res;
678 if (derivative != first_order_endo_derivatives.end())
679 {
680 set<pair<int, int>> result;
681 derivative->second->collectEndogenous(result);
682 auto d_endo_variable = result.find({ var, 0 });
683 //Determine whether the equation could be evaluated rather than to be solved
684 if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
685 Equation_Simulation_Type = E_EVALUATE;
686 else
687 {
688 vector<tuple<int, expr_t, expr_t>> List_of_Op_RHS;
689 res = equations[eq]->normalizeEquation(var, List_of_Op_RHS);
690 if (mfs == 2)
691 {
692 if (d_endo_variable == result.end() && res.second)
693 Equation_Simulation_Type = E_EVALUATE_S;
694 }
695 else if (mfs == 3)
696 {
697 if (res.second) // The equation could be solved analytically
698 Equation_Simulation_Type = E_EVALUATE_S;
699 }
700 }
701 }
702 V_Equation_Simulation_Type[eq] = { Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second) };
703 }
704 return V_Equation_Simulation_Type;
705 }
706
707 void
getVariableLeadLagByBlock(const dynamic_jacob_map_t & dynamic_jacobian,const vector<int> & components_set,int nb_blck_sim,lag_lead_vector_t & equation_lead_lag,lag_lead_vector_t & variable_lead_lag,const vector<int> & equation_reordered,const vector<int> & variable_reordered) const708 ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const
709 {
710 int nb_endo = symbol_table.endo_nbr();
711 variable_lead_lag = lag_lead_vector_t(nb_endo, { 0, 0 });
712 equation_lead_lag = lag_lead_vector_t(nb_endo, { 0, 0 });
713 vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
714 for (int i = 0; i < nb_endo; i++)
715 {
716 if (i < static_cast<int>(prologue))
717 {
718 variable_blck[variable_reordered[i]] = i;
719 equation_blck[equation_reordered[i]] = i;
720 }
721 else if (i < static_cast<int>(components_set.size() + prologue))
722 {
723 variable_blck[variable_reordered[i]] = components_set[i-prologue] + prologue;
724 equation_blck[equation_reordered[i]] = components_set[i-prologue] + prologue;
725 }
726 else
727 {
728 variable_blck[variable_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
729 equation_blck[equation_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
730 }
731 }
732 for (const auto &it : dynamic_jacobian)
733 {
734 auto [lag, j_1, i_1] = it.first;
735 if (variable_blck[i_1] == equation_blck[j_1])
736 {
737 if (lag > variable_lead_lag[i_1].second)
738 variable_lead_lag[i_1] = { variable_lead_lag[i_1].first, lag };
739 if (lag < -variable_lead_lag[i_1].first)
740 variable_lead_lag[i_1] = { -lag, variable_lead_lag[i_1].second };
741 if (lag > equation_lead_lag[j_1].second)
742 equation_lead_lag[j_1] = { equation_lead_lag[j_1].first, lag };
743 if (lag < -equation_lead_lag[j_1].first)
744 equation_lead_lag[j_1] = { -lag, equation_lead_lag[j_1].second };
745 }
746 }
747 }
748
749 void
computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t & static_jacobian,const dynamic_jacob_map_t & dynamic_jacobian,vector<int> & equation_reordered,vector<int> & variable_reordered,vector<pair<int,int>> & blocks,const equation_type_and_normalized_equation_t & Equation_Type,bool verbose_,bool select_feedback_variable,int mfs,vector<int> & inv_equation_reordered,vector<int> & inv_variable_reordered,lag_lead_vector_t & equation_lag_lead,lag_lead_vector_t & variable_lag_lead,vector<unsigned int> & n_static,vector<unsigned int> & n_forward,vector<unsigned int> & n_backward,vector<unsigned int> & n_mixed) const750 ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const
751 {
752 int nb_var = variable_reordered.size();
753 int n = nb_var - prologue - epilogue;
754
755 AdjacencyList_t G2(n);
756
757 // It is necessary to manually initialize vertex_index property since this graph uses listS and not vecS as underlying vertex container
758 auto v_index = get(boost::vertex_index, G2);
759 for (int i = 0; i < n; i++)
760 put(v_index, vertex(i, G2), i);
761
762 vector<int> reverse_equation_reordered(nb_var), reverse_variable_reordered(nb_var);
763
764 for (int i = 0; i < nb_var; i++)
765 {
766 reverse_equation_reordered[equation_reordered[i]] = i;
767 reverse_variable_reordered[variable_reordered[i]] = i;
768 }
769 jacob_map_t tmp_normalized_contemporaneous_jacobian;
770 if (cutoff == 0)
771 {
772 set<pair<int, int>> endo;
773 for (int i = 0; i < nb_var; i++)
774 {
775 endo.clear();
776 equations[i]->collectEndogenous(endo);
777 for (const auto &it : endo)
778 tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
779 }
780 }
781 else
782 tmp_normalized_contemporaneous_jacobian = static_jacobian;
783 for (const auto &[key, value] : tmp_normalized_contemporaneous_jacobian)
784 if (reverse_equation_reordered[key.first] >= static_cast<int>(prologue) && reverse_equation_reordered[key.first] < static_cast<int>(nb_var - epilogue)
785 && reverse_variable_reordered[key.second] >= static_cast<int>(prologue) && reverse_variable_reordered[key.second] < static_cast<int>(nb_var - epilogue)
786 && key.first != endo2eq[key.second])
787 add_edge(vertex(reverse_equation_reordered[endo2eq[key.second]]-prologue, G2),
788 vertex(reverse_equation_reordered[key.first]-prologue, G2),
789 G2);
790
791 vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
792 boost::iterator_property_map<int *, boost::property_map<AdjacencyList_t, boost::vertex_index_t>::type, int, int &> endo2block_map(&endo2block[0], get(boost::vertex_index, G2));
793
794 // Compute strongly connected components
795 int num = strong_components(G2, endo2block_map);
796
797 blocks = vector<pair<int, int>>(num, { 0, 0 });
798
799 // Create directed acyclic graph associated to the strongly connected components
800 using DirectedGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS>;
801 DirectedGraph dag(num);
802
803 for (unsigned int i = 0; i < num_vertices(G2); i++)
804 {
805 AdjacencyList_t::out_edge_iterator it_out, out_end;
806 AdjacencyList_t::vertex_descriptor vi = vertex(i, G2);
807 for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
808 {
809 int t_b = endo2block_map[target(*it_out, G2)];
810 int s_b = endo2block_map[source(*it_out, G2)];
811 if (s_b != t_b)
812 add_edge(s_b, t_b, dag);
813 }
814 }
815
816 // Compute topological sort of DAG (ordered list of unordered SCC)
817 deque<int> ordered2unordered;
818 topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order
819
820 // Construct mapping from unordered SCC to ordered SCC
821 vector<int> unordered2ordered(num);
822 for (int i = 0; i < num; i++)
823 unordered2ordered[ordered2unordered[i]] = i;
824
825 //This vector contains for each block:
826 // - first set = equations belonging to the block,
827 // - second set = the feeback variables,
828 // - third vector = the reordered non-feedback variables.
829 vector<tuple<set<int>, set<int>, vector<int>>> components_set(num);
830 for (unsigned int i = 0; i < endo2block.size(); i++)
831 {
832 endo2block[i] = unordered2ordered[endo2block[i]];
833 blocks[endo2block[i]].first++;
834 get<0>(components_set[endo2block[i]]).insert(i);
835 }
836
837 getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
838
839 vector<int> tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
840 int order = prologue;
841 //Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
842 if (select_feedback_variable)
843 {
844 for (int i = 0; i < n; i++)
845 if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
846 || variable_lag_lead[variable_reordered[i+prologue]].second > 0
847 || variable_lag_lead[variable_reordered[i+prologue]].first > 0
848 || equation_lag_lead[equation_reordered[i+prologue]].second > 0
849 || equation_lag_lead[equation_reordered[i+prologue]].first > 0
850 || mfs == 0)
851 add_edge(vertex(i, G2), vertex(i, G2), G2);
852 }
853 else
854 for (int i = 0; i < n; i++)
855 if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
856 add_edge(vertex(i, G2), vertex(i, G2), G2);
857
858 //Determines the dynamic structure of each equation
859 n_static = vector<unsigned int>(prologue+num+epilogue, 0);
860 n_forward = vector<unsigned int>(prologue+num+epilogue, 0);
861 n_backward = vector<unsigned int>(prologue+num+epilogue, 0);
862 n_mixed = vector<unsigned int>(prologue+num+epilogue, 0);
863
864 for (int i = 0; i < static_cast<int>(prologue); i++)
865 if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
866 n_mixed[i]++;
867 else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
868 n_forward[i]++;
869 else if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
870 n_backward[i]++;
871 else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
872 n_static[i]++;
873
874 //For each block, the minimum set of feedback variable is computed
875 // and the non-feedback variables are reordered to get
876 // a sub-recursive block without feedback variables
877
878 for (int i = 0; i < num; i++)
879 {
880 AdjacencyList_t G = extract_subgraph(G2, get<0>(components_set[i]));
881 set<int> feed_back_vertices;
882 AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
883 auto v_index = get(boost::vertex_index, G);
884 get<1>(components_set[i]) = feed_back_vertices;
885 blocks[i].second = feed_back_vertices.size();
886 vector<int> Reordered_Vertice;
887 Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
888
889 //First we have the recursive equations conditional on feedback variables
890 for (int j = 0; j < 4; j++)
891 for (int its : Reordered_Vertice)
892 {
893 bool something_done = false;
894 if (j == 2 && variable_lag_lead[tmp_variable_reordered[its+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its+prologue]].second != 0)
895 {
896 n_mixed[prologue+i]++;
897 something_done = true;
898 }
899 else if (j == 3 && variable_lag_lead[tmp_variable_reordered[its+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its+prologue]].second != 0)
900 {
901 n_forward[prologue+i]++;
902 something_done = true;
903 }
904 else if (j == 1 && variable_lag_lead[tmp_variable_reordered[its+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its+prologue]].second == 0)
905 {
906 n_backward[prologue+i]++;
907 something_done = true;
908 }
909 else if (j == 0 && variable_lag_lead[tmp_variable_reordered[its+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[its+prologue]].second == 0)
910 {
911 n_static[prologue+i]++;
912 something_done = true;
913 }
914 if (something_done)
915 {
916 equation_reordered[order] = tmp_equation_reordered[its+prologue];
917 variable_reordered[order] = tmp_variable_reordered[its+prologue];
918 order++;
919 }
920 }
921
922 get<2>(components_set[i]) = Reordered_Vertice;
923 //Second we have the equations related to the feedback variables
924 for (int j = 0; j < 4; j++)
925 for (int feed_back_vertice : feed_back_vertices)
926 {
927 bool something_done = false;
928 if (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0)
929 {
930 n_mixed[prologue+i]++;
931 something_done = true;
932 }
933 else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0)
934 {
935 n_forward[prologue+i]++;
936 something_done = true;
937 }
938 else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0)
939 {
940 n_backward[prologue+i]++;
941 something_done = true;
942 }
943 else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0)
944 {
945 n_static[prologue+i]++;
946 something_done = true;
947 }
948 if (something_done)
949 {
950 equation_reordered[order] = tmp_equation_reordered[v_index[vertex(feed_back_vertice, G)]+prologue];
951 variable_reordered[order] = tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue];
952 order++;
953 }
954 }
955 }
956
957 for (int i = 0; i < static_cast<int>(epilogue); i++)
958 if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
959 n_mixed[prologue+num+i]++;
960 else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
961 n_forward[prologue+num+i]++;
962 else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
963 n_backward[prologue+num+i]++;
964 else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
965 n_static[prologue+num+i]++;
966
967 inv_equation_reordered.resize(nb_var);
968 inv_variable_reordered.resize(nb_var);
969 for (int i = 0; i < nb_var; i++)
970 {
971 inv_variable_reordered[variable_reordered[i]] = i;
972 inv_equation_reordered[equation_reordered[i]] = i;
973 }
974 }
975
976 void
printBlockDecomposition(const vector<pair<int,int>> & blocks) const977 ModelTree::printBlockDecomposition(const vector<pair<int, int>> &blocks) const
978 {
979 int largest_block = 0,
980 Nb_SimulBlocks = 0,
981 Nb_feedback_variable = 0;
982 unsigned int Nb_TotalBlocks = getNbBlocks();
983 for (unsigned int block = 0; block < Nb_TotalBlocks; block++)
984 {
985 BlockSimulationType simulation_type = getBlockSimulationType(block);
986 if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
987 {
988 Nb_SimulBlocks++;
989 int size = getBlockSize(block);
990 if (size > largest_block)
991 {
992 largest_block = size;
993 Nb_feedback_variable = getBlockMfs(block);
994 }
995 }
996 }
997
998 int Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks;
999 cout << Nb_TotalBlocks << " block(s) found:" << endl
1000 << " " << Nb_RecursBlocks << " recursive block(s) and " << Nb_SimulBlocks << " simultaneous block(s)." << endl
1001 << " the largest simultaneous block has " << largest_block << " equation(s)" << endl
1002 << " and " << Nb_feedback_variable << " feedback variable(s)." << endl;
1003 }
1004
1005 block_type_firstequation_size_mfs_t
reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t & dynamic_jacobian,vector<pair<int,int>> & blocks,const equation_type_and_normalized_equation_t & Equation_Type,const vector<int> & variable_reordered,const vector<int> & equation_reordered,vector<unsigned int> & n_static,vector<unsigned int> & n_forward,vector<unsigned int> & n_backward,vector<unsigned int> & n_mixed,vector<tuple<int,int,int,int>> & block_col_type,bool linear_decomposition)1006 ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<tuple<int, int, int, int>> &block_col_type, bool linear_decomposition)
1007 {
1008 int i = 0;
1009 int count_equ = 0, blck_count_simult = 0;
1010 int Blck_Size, MFS_Size;
1011 int Lead, Lag;
1012 block_type_firstequation_size_mfs_t block_type_size_mfs;
1013 BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
1014 int eq = 0;
1015 unsigned int l_n_static = 0, l_n_forward = 0, l_n_backward = 0, l_n_mixed = 0;
1016 for (i = 0; i < static_cast<int>(prologue+blocks.size()+epilogue); i++)
1017 {
1018 int first_count_equ = count_equ;
1019 if (i < static_cast<int>(prologue))
1020 {
1021 Blck_Size = 1;
1022 MFS_Size = 1;
1023 }
1024 else if (i < static_cast<int>(prologue+blocks.size()))
1025 {
1026 Blck_Size = blocks[blck_count_simult].first;
1027 MFS_Size = blocks[blck_count_simult].second;
1028 blck_count_simult++;
1029 }
1030 else if (i < static_cast<int>(prologue+blocks.size()+epilogue))
1031 {
1032 Blck_Size = 1;
1033 MFS_Size = 1;
1034 }
1035
1036 Lag = Lead = 0;
1037 set<pair<int, int>> endo;
1038 for (count_equ = first_count_equ; count_equ < Blck_Size+first_count_equ; count_equ++)
1039 {
1040 endo.clear();
1041 equations[equation_reordered[count_equ]]->collectEndogenous(endo);
1042 for (const auto &it : endo)
1043 {
1044 int curr_variable = it.first;
1045 int curr_lag = it.second;
1046 if (linear_decomposition)
1047 {
1048 if (dynamic_jacobian.find({ curr_lag, equation_reordered[count_equ], curr_variable }) != dynamic_jacobian.end())
1049 {
1050 if (curr_lag > Lead)
1051 Lead = curr_lag;
1052 else if (-curr_lag > Lag)
1053 Lag = -curr_lag;
1054 }
1055 }
1056 else
1057 {
1058 if (find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable)
1059 != variable_reordered.begin()+(first_count_equ+Blck_Size)
1060 && dynamic_jacobian.find({ curr_lag, equation_reordered[count_equ], curr_variable }) != dynamic_jacobian.end())
1061 {
1062 if (curr_lag > Lead)
1063 Lead = curr_lag;
1064 else if (-curr_lag > Lag)
1065 Lag = -curr_lag;
1066 }
1067 }
1068 }
1069 }
1070 if (Lag > 0 && Lead > 0)
1071 {
1072 if (Blck_Size == 1)
1073 Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
1074 else
1075 Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE;
1076 }
1077 else if (Blck_Size > 1)
1078 {
1079 if (Lead > 0)
1080 Simulation_Type = SOLVE_BACKWARD_COMPLETE;
1081 else
1082 Simulation_Type = SOLVE_FORWARD_COMPLETE;
1083 }
1084 else
1085 {
1086 if (Lead > 0)
1087 Simulation_Type = SOLVE_BACKWARD_SIMPLE;
1088 else
1089 Simulation_Type = SOLVE_FORWARD_SIMPLE;
1090 }
1091 l_n_static = n_static[i];
1092 l_n_forward = n_forward[i];
1093 l_n_backward = n_backward[i];
1094 l_n_mixed = n_mixed[i];
1095 if (Blck_Size == 1)
1096 {
1097 if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
1098 {
1099 if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
1100 Simulation_Type = EVALUATE_BACKWARD;
1101 else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
1102 Simulation_Type = EVALUATE_FORWARD;
1103 }
1104 if (i > 0)
1105 {
1106 bool is_lead = false, is_lag = false;
1107 int c_Size = get<2>(block_type_size_mfs[block_type_size_mfs.size()-1]);
1108 int first_equation = get<1>(block_type_size_mfs[block_type_size_mfs.size()-1]);
1109 if (c_Size > 0 && ((prev_Type == EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
1110 || (prev_Type == EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag)))
1111 {
1112 for (int j = first_equation; j < first_equation+c_Size; j++)
1113 {
1114 auto it = dynamic_jacobian.find({ -1, equation_reordered[eq], variable_reordered[j] });
1115 if (it != dynamic_jacobian.end())
1116 is_lag = true;
1117 it = dynamic_jacobian.find({ +1, equation_reordered[eq], variable_reordered[j] });
1118 if (it != dynamic_jacobian.end())
1119 is_lead = true;
1120 }
1121 }
1122 if ((prev_Type == EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
1123 || (prev_Type == EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag))
1124 {
1125 //merge the current block with the previous one
1126 BlockSimulationType c_Type = get<0>(block_type_size_mfs[block_type_size_mfs.size()-1]);
1127 c_Size++;
1128 block_type_size_mfs[block_type_size_mfs.size()-1] = { c_Type, first_equation, c_Size, c_Size };
1129 if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
1130 Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
1131 if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
1132 Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
1133 block_lag_lead[block_type_size_mfs.size()-1] = { Lag, Lead };
1134 auto tmp = block_col_type[block_col_type.size()-1];
1135 block_col_type[block_col_type.size()-1] = { get<0>(tmp)+l_n_static, get<1>(tmp)+l_n_forward, get<2>(tmp)+l_n_backward, get<3>(tmp)+l_n_mixed };
1136 }
1137 else
1138 {
1139 block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
1140 block_lag_lead.emplace_back(Lag, Lead);
1141 block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
1142 }
1143 }
1144 else
1145 {
1146 block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
1147 block_lag_lead.emplace_back(Lag, Lead);
1148 block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
1149 }
1150 }
1151 else
1152 {
1153 block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
1154 block_lag_lead.emplace_back(Lag, Lead);
1155 block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
1156 }
1157 prev_Type = Simulation_Type;
1158 eq += Blck_Size;
1159 }
1160 return block_type_size_mfs;
1161 }
1162
1163 vector<bool>
equationLinear(map<tuple<int,int,int>,expr_t> first_order_endo_derivatives) const1164 ModelTree::equationLinear(map<tuple<int, int, int>, expr_t> first_order_endo_derivatives) const
1165 {
1166 vector<bool> is_linear(symbol_table.endo_nbr(), true);
1167 for (const auto &it : first_order_endo_derivatives)
1168 {
1169 expr_t Id = it.second;
1170 set<pair<int, int>> endogenous;
1171 Id->collectEndogenous(endogenous);
1172 if (endogenous.size() > 0)
1173 {
1174 int eq = get<0>(it.first);
1175 is_linear[eq] = false;
1176 }
1177 }
1178 return is_linear;
1179 }
1180
1181 vector<bool>
BlockLinear(const blocks_derivatives_t & blocks_derivatives,const vector<int> & variable_reordered) const1182 ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vector<int> &variable_reordered) const
1183 {
1184 unsigned int nb_blocks = getNbBlocks();
1185 vector<bool> blocks_linear(nb_blocks, true);
1186 for (unsigned int block = 0; block < nb_blocks; block++)
1187 {
1188 BlockSimulationType simulation_type = getBlockSimulationType(block);
1189 int block_size = getBlockSize(block);
1190 block_derivatives_equation_variable_laglead_nodeid_t derivatives_block = blocks_derivatives[block];
1191 int first_variable_position = getBlockFirstEquation(block);
1192 if (simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
1193 for (const auto &[ignore, ignore2, lag, d1] : derivatives_block)
1194 {
1195 if (lag == 0)
1196 {
1197 set<pair<int, int>> endogenous;
1198 d1->collectEndogenous(endogenous);
1199 if (endogenous.size() > 0)
1200 for (int l = 0; l < block_size; l++)
1201 if (endogenous.find({ variable_reordered[first_variable_position+l], 0 }) != endogenous.end())
1202 {
1203 blocks_linear[block] = false;
1204 goto the_end;
1205 }
1206 }
1207 }
1208 else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
1209 for (const auto &[ignore, ignore2, lag, d1] : derivatives_block)
1210 {
1211 set<pair<int, int>> endogenous;
1212 d1->collectEndogenous(endogenous);
1213 if (endogenous.size() > 0)
1214 for (int l = 0; l < block_size; l++)
1215 if (endogenous.find({ variable_reordered[first_variable_position+l], lag }) != endogenous.end())
1216 {
1217 blocks_linear[block] = false;
1218 goto the_end;
1219 }
1220 }
1221 the_end:
1222 ;
1223 }
1224 return blocks_linear;
1225 }
1226
1227 int
equation_number() const1228 ModelTree::equation_number() const
1229 {
1230 return (equations.size());
1231 }
1232
1233 void
writeDerivative(ostream & output,int eq,int symb_id,int lag,ExprNodeOutputType output_type,const temporary_terms_t & temporary_terms) const1234 ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
1235 ExprNodeOutputType output_type,
1236 const temporary_terms_t &temporary_terms) const
1237 {
1238 if (auto it = derivatives[1].find({ eq, getDerivID(symb_id, lag) });
1239 it != derivatives[1].end())
1240 it->second->writeOutput(output, output_type, temporary_terms, {});
1241 else
1242 output << 0;
1243 }
1244
1245 void
computeDerivatives(int order,const set<int> & vars)1246 ModelTree::computeDerivatives(int order, const set<int> &vars)
1247 {
1248 assert(order >= 1);
1249
1250 computed_derivs_order = order;
1251
1252 // Do not shrink the vectors, since they have a minimal size of 4 (see constructor)
1253 derivatives.resize(max(static_cast<size_t>(order+1), derivatives.size()));
1254 NNZDerivatives.resize(max(static_cast<size_t>(order+1), NNZDerivatives.size()), 0);
1255
1256 // First-order derivatives
1257 for (int var : vars)
1258 for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
1259 {
1260 expr_t d1 = equations[eq]->getDerivative(var);
1261 if (d1 == Zero)
1262 continue;
1263 derivatives[1][{ eq, var }] = d1;
1264 ++NNZDerivatives[1];
1265 }
1266
1267 // Higher-order derivatives
1268 for (int o = 2; o <= order; o++)
1269 for (const auto &it : derivatives[o-1])
1270 for (int var : vars)
1271 {
1272 if (it.first.back() > var)
1273 continue;
1274
1275 expr_t d = it.second->getDerivative(var);
1276 if (d == Zero)
1277 continue;
1278
1279 vector<int> indices{it.first};
1280 indices.push_back(var);
1281 // At this point, indices of endogenous variables are sorted in non-decreasing order
1282 derivatives[o][indices] = d;
1283 // We output symmetric elements at order = 2
1284 if (o == 2 && indices[1] != indices[2])
1285 NNZDerivatives[o] += 2;
1286 else
1287 NNZDerivatives[o]++;
1288 }
1289 }
1290
1291 void
computeTemporaryTerms(bool is_matlab,bool no_tmp_terms)1292 ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
1293 {
1294 /* Collect all model local variables appearing in equations (and only those,
1295 because printing unused model local variables can lead to a crash,
1296 see Dynare/dynare#101).
1297 Then store them in a dedicated structure (temporary_terms_mlv), that will
1298 be treated as the rest of temporary terms. */
1299 temporary_terms_mlv.clear();
1300 set<int> used_local_vars;
1301 for (auto &equation : equations)
1302 equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
1303 for (int used_local_var : used_local_vars)
1304 {
1305 VariableNode *v = AddVariable(used_local_var);
1306 temporary_terms_mlv[v] = local_variables_table.find(used_local_var)->second;
1307 }
1308
1309 // Compute the temporary terms in equations and derivatives
1310 map<pair<int, int>, temporary_terms_t> temp_terms_map;
1311 map<expr_t, pair<int, pair<int, int>>> reference_count;
1312
1313 for (auto &equation : equations)
1314 equation->computeTemporaryTerms({ 0, 0 },
1315 temp_terms_map,
1316 reference_count,
1317 is_matlab);
1318
1319 for (int order = 1; order < static_cast<int>(derivatives.size()); order++)
1320 for (const auto &it : derivatives[order])
1321 it.second->computeTemporaryTerms({ 0, order },
1322 temp_terms_map,
1323 reference_count,
1324 is_matlab);
1325
1326 /* If the user has specified the notmpterms option, clear all temporary
1327 terms, except those that correspond to external functions (since they are
1328 not optional) */
1329 if (no_tmp_terms)
1330 for (auto &it : temp_terms_map)
1331 // The following loop can be simplified with std::erase_if() in C++20
1332 for (auto it2 = it.second.begin(); it2 != it.second.end();)
1333 if (!dynamic_cast<AbstractExternalFunctionNode *>(*it2))
1334 it2 = it.second.erase(it2);
1335 else
1336 ++it2;
1337
1338 // Fill the (now obsolete) temporary_terms structure
1339 temporary_terms.clear();
1340 for (const auto &it : temp_terms_map)
1341 temporary_terms.insert(it.second.begin(), it.second.end());
1342
1343 // Fill the new structure
1344 temporary_terms_derivatives.clear();
1345 temporary_terms_derivatives.resize(derivatives.size());
1346 for (int order = 0; order < static_cast<int>(derivatives.size()); order++)
1347 temporary_terms_derivatives[order] = move(temp_terms_map[{ 0, order }]);
1348
1349 // Compute indices in MATLAB/Julia vector
1350 int idx = 0;
1351 for (auto &it : temporary_terms_mlv)
1352 temporary_terms_idxs[it.first] = idx++;
1353 for (int order = 0; order < static_cast<int>(derivatives.size()); order++)
1354 for (const auto &it : temporary_terms_derivatives[order])
1355 temporary_terms_idxs[it] = idx++;
1356 }
1357
1358 void
writeModelLocalVariableTemporaryTerms(temporary_terms_t & temp_term_union,const temporary_terms_idxs_t & tt_idxs,ostream & output,ExprNodeOutputType output_type,deriv_node_temp_terms_t & tef_terms) const1359 ModelTree::writeModelLocalVariableTemporaryTerms(temporary_terms_t &temp_term_union,
1360 const temporary_terms_idxs_t &tt_idxs,
1361 ostream &output, ExprNodeOutputType output_type,
1362 deriv_node_temp_terms_t &tef_terms) const
1363 {
1364 temporary_terms_t tto;
1365 for (auto it : temporary_terms_mlv)
1366 tto.insert(it.first);
1367
1368 for (auto &it : temporary_terms_mlv)
1369 {
1370 it.second->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
1371
1372 if (isJuliaOutput(output_type))
1373 output << " @inbounds const ";
1374
1375 it.first->writeOutput(output, output_type, tto, tt_idxs, tef_terms);
1376 output << " = ";
1377 it.second->writeOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
1378
1379 if (isCOutput(output_type) || isMatlabOutput(output_type))
1380 output << ";";
1381 output << endl;
1382
1383 /* We put in temp_term_union the VariableNode corresponding to the MLV,
1384 not its definition, so that when equations use the MLV,
1385 T(XXX) is printed instead of the MLV name */
1386 temp_term_union.insert(it.first);
1387 }
1388 }
1389
1390 void
writeTemporaryTerms(const temporary_terms_t & tt,temporary_terms_t & temp_term_union,const temporary_terms_idxs_t & tt_idxs,ostream & output,ExprNodeOutputType output_type,deriv_node_temp_terms_t & tef_terms) const1391 ModelTree::writeTemporaryTerms(const temporary_terms_t &tt,
1392 temporary_terms_t &temp_term_union,
1393 const temporary_terms_idxs_t &tt_idxs,
1394 ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const
1395 {
1396 for (auto it : tt)
1397 {
1398 if (dynamic_cast<AbstractExternalFunctionNode *>(it))
1399 it->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
1400
1401 if (isJuliaOutput(output_type))
1402 output << " @inbounds ";
1403
1404 it->writeOutput(output, output_type, tt, tt_idxs, tef_terms);
1405 output << " = ";
1406 it->writeOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
1407
1408 if (isCOutput(output_type) || isMatlabOutput(output_type))
1409 output << ";";
1410 output << endl;
1411
1412 temp_term_union.insert(it);
1413 }
1414 }
1415
1416 void
writeJsonTemporaryTerms(const temporary_terms_t & tt,temporary_terms_t & temp_term_union,ostream & output,deriv_node_temp_terms_t & tef_terms,const string & concat) const1417 ModelTree::writeJsonTemporaryTerms(const temporary_terms_t &tt,
1418 temporary_terms_t &temp_term_union,
1419 ostream &output,
1420 deriv_node_temp_terms_t &tef_terms, const string &concat) const
1421 {
1422 // Local var used to keep track of temp nodes already written
1423 bool wrote_term = false;
1424 temporary_terms_t tt2 = temp_term_union;
1425
1426 output << R"("external_functions_temporary_terms_)" << concat << R"(": [)";
1427 for (auto it : tt)
1428 {
1429 if (dynamic_cast<AbstractExternalFunctionNode *>(it))
1430 {
1431 if (wrote_term)
1432 output << ", ";
1433 vector<string> efout;
1434 it->writeJsonExternalFunctionOutput(efout, tt2, tef_terms);
1435 for (auto it1 = efout.begin(); it1 != efout.end(); ++it1)
1436 {
1437 if (it1 != efout.begin())
1438 output << ", ";
1439 output << *it1;
1440 }
1441 wrote_term = true;
1442 }
1443 tt2.insert(it);
1444 }
1445
1446 wrote_term = false;
1447 output << "]"
1448 << R"(, "temporary_terms_)" << concat << R"(": [)";
1449 for (const auto &it : tt)
1450 {
1451 if (wrote_term)
1452 output << ", ";
1453 output << R"({"temporary_term": ")";
1454 it->writeJsonOutput(output, tt, tef_terms);
1455 output << R"(")"
1456 << R"(, "value": ")";
1457 it->writeJsonOutput(output, temp_term_union, tef_terms);
1458 output << R"("})" << endl;
1459 wrote_term = true;
1460
1461 temp_term_union.insert(it);
1462 }
1463 output << "]";
1464 }
1465
1466 void
fixNestedParenthesis(ostringstream & output,map<string,string> & tmp_paren_vars,bool & message_printed) const1467 ModelTree::fixNestedParenthesis(ostringstream &output, map<string, string> &tmp_paren_vars, bool &message_printed) const
1468 {
1469 string str = output.str();
1470 if (!testNestedParenthesis(str))
1471 return;
1472 int open = 0;
1473 int first_open_paren = 0;
1474 int matching_paren = 0;
1475 bool hit_limit = false;
1476 int i1 = 0;
1477 for (size_t i = 0; i < str.length(); i++)
1478 {
1479 if (str.at(i) == '(')
1480 {
1481 if (open == 0)
1482 first_open_paren = i;
1483 open++;
1484 }
1485 else if (str.at(i) == ')')
1486 {
1487 open--;
1488 if (open == 0)
1489 matching_paren = i;
1490 }
1491 if (open > 32)
1492 hit_limit = true;
1493
1494 if (hit_limit && open == 0)
1495 {
1496 if (!message_printed)
1497 {
1498 cerr << "Warning: A .m file created by Dynare will have more than 32 nested parenthesis. MATLAB cannot support this. " << endl
1499 << " We are going to modify, albeit inefficiently, this output to have fewer than 32 nested parenthesis. " << endl
1500 << " It would hence behoove you to use the use_dll option of the model block to circumnavigate this problem." << endl
1501 << " If you have not yet set up a compiler on your system, see the MATLAB documentation for doing so." << endl
1502 << " For Windows, see: https://www.mathworks.com/help/matlab/matlab_external/install-mingw-support-package.html" << endl << endl;
1503 message_printed = true;
1504 }
1505 string str1 = str.substr(first_open_paren, matching_paren - first_open_paren + 1);
1506 string repstr, varname;
1507 while (testNestedParenthesis(str1))
1508 {
1509 size_t open_paren_idx = string::npos;
1510 size_t match_paren_idx = string::npos;
1511 size_t last_open_paren = string::npos;
1512 for (size_t j = 0; j < str1.length(); j++)
1513 {
1514 if (str1.at(j) == '(')
1515 {
1516 // don't match, e.g. y(1)
1517 if (size_t idx = str1.find_last_of("*/-+", j - 1);
1518 j == 0 || (idx != string::npos && idx == j - 1))
1519 open_paren_idx = j;
1520 last_open_paren = j;
1521 }
1522 else if (str1.at(j) == ')')
1523 {
1524 // don't match, e.g. y(1)
1525 if (size_t idx = str1.find_last_not_of("0123456789", j - 1);
1526 idx != string::npos && idx != last_open_paren)
1527 match_paren_idx = j;
1528 }
1529
1530 if (open_paren_idx != string::npos && match_paren_idx != string::npos)
1531 {
1532 string val = str1.substr(open_paren_idx, match_paren_idx - open_paren_idx + 1);
1533 if (auto it = tmp_paren_vars.find(val);
1534 it == tmp_paren_vars.end())
1535 {
1536 ostringstream ptvstr;
1537 ptvstr << i1++;
1538 varname = "paren32_tmp_var_" + ptvstr.str();
1539 repstr = repstr + varname + " = " + val + ";\n";
1540 tmp_paren_vars[val] = varname;
1541 }
1542 else
1543 varname = it->second;
1544 str1.replace(open_paren_idx, match_paren_idx - open_paren_idx + 1, varname);
1545 break;
1546 }
1547 }
1548 }
1549 if (auto it = tmp_paren_vars.find(str1);
1550 it == tmp_paren_vars.end())
1551 {
1552 ostringstream ptvstr;
1553 ptvstr << i1++;
1554 varname = "paren32_tmp_var_" + ptvstr.str();
1555 repstr = repstr + varname + " = " + str1 + ";\n";
1556 }
1557 else
1558 varname = it->second;
1559 str.replace(first_open_paren, matching_paren - first_open_paren + 1, varname);
1560 size_t insertLoc = str.find_last_of("\n", first_open_paren);
1561 str.insert(insertLoc + 1, repstr);
1562 hit_limit = false;
1563 i = -1;
1564 first_open_paren = matching_paren = open = 0;
1565 }
1566 }
1567 output.str(str);
1568 }
1569
1570 bool
testNestedParenthesis(const string & str) const1571 ModelTree::testNestedParenthesis(const string &str) const
1572 {
1573 int open = 0;
1574 for (char i : str)
1575 {
1576 if (i == '(')
1577 open++;
1578 else if (i == ')')
1579 open--;
1580 if (open > 32)
1581 return true;
1582 }
1583 return false;
1584 }
1585
1586 void
compileTemporaryTerms(ostream & code_file,unsigned int & instruction_number,const temporary_terms_t & tt,map_idx_t map_idx,bool dynamic,bool steady_dynamic) const1587 ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const
1588 {
1589 // Local var used to keep track of temp nodes already written
1590 temporary_terms_t tt2;
1591 // To store the functions that have already been written in the form TEF* = ext_fun();
1592 deriv_node_temp_terms_t tef_terms;
1593 for (auto it : tt)
1594 {
1595 if (dynamic_cast<AbstractExternalFunctionNode *>(it))
1596 {
1597 it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
1598 }
1599
1600 FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(map_idx.find(it->idx)->second));
1601 fnumexpr.write(code_file, instruction_number);
1602 it->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
1603 if (dynamic)
1604 {
1605 FSTPT_ fstpt(static_cast<int>(map_idx.find(it->idx)->second));
1606 fstpt.write(code_file, instruction_number);
1607 }
1608 else
1609 {
1610 FSTPST_ fstpst(static_cast<int>(map_idx.find(it->idx)->second));
1611 fstpst.write(code_file, instruction_number);
1612 }
1613 // Insert current node into tt2
1614 tt2.insert(it);
1615 }
1616 }
1617
1618 void
writeJsonModelLocalVariables(ostream & output,bool write_tef_terms,deriv_node_temp_terms_t & tef_terms) const1619 ModelTree::writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, deriv_node_temp_terms_t &tef_terms) const
1620 {
1621 /* Collect all model local variables appearing in equations, and print only
1622 them. Printing unused model local variables can lead to a crash (see
1623 ticket #101). */
1624 set<int> used_local_vars;
1625
1626 // Use an empty set for the temporary terms
1627 const temporary_terms_t tt;
1628
1629 for (auto equation : equations)
1630 equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
1631
1632 output << R"("model_local_variables": [)";
1633 bool printed = false;
1634 for (int it : local_variables_vector)
1635 if (used_local_vars.find(it) != used_local_vars.end())
1636 {
1637 if (printed)
1638 output << ", ";
1639 else
1640 printed = true;
1641
1642 int id = it;
1643 expr_t value = local_variables_table.find(id)->second;
1644 if (write_tef_terms)
1645 {
1646 vector<string> efout;
1647 value->writeJsonExternalFunctionOutput(efout, tt, tef_terms);
1648 for (auto it1 = efout.begin(); it1 != efout.end(); ++it1)
1649 {
1650 if (it1 != efout.begin())
1651 output << ", ";
1652 output << *it1;
1653 }
1654
1655 if (!efout.empty())
1656 output << ", ";
1657 }
1658
1659 output << R"({"variable": ")" << symbol_table.getName(id)
1660 << R"(", "value": ")";
1661 value->writeJsonOutput(output, tt, tef_terms);
1662 output << R"("})" << endl;
1663 }
1664 output << "]";
1665 }
1666
1667 void
writeModelEquations(ostream & output,ExprNodeOutputType output_type) const1668 ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const
1669 {
1670 temporary_terms_t tt;
1671 temporary_terms_idxs_t ttidxs;
1672 writeModelEquations(output, output_type, tt);
1673 }
1674
1675 void
writeModelEquations(ostream & output,ExprNodeOutputType output_type,const temporary_terms_t & temporary_terms) const1676 ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type,
1677 const temporary_terms_t &temporary_terms) const
1678 {
1679 for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
1680 {
1681 BinaryOpNode *eq_node = equations[eq];
1682 expr_t lhs = eq_node->arg1;
1683 expr_t rhs = eq_node->arg2;
1684
1685 // Test if the right hand side of the equation is empty.
1686 double vrhs = 1.0;
1687 try
1688 {
1689 vrhs = rhs->eval(eval_context_t());
1690 }
1691 catch (ExprNode::EvalException &e)
1692 {
1693 }
1694
1695 if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
1696 if (isJuliaOutput(output_type))
1697 {
1698 output << " @inbounds residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
1699 << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
1700 << RIGHT_ARRAY_SUBSCRIPT(output_type)
1701 << " = (";
1702 lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
1703 output << ") - (";
1704 rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
1705 output << ")" << endl;
1706 }
1707 else
1708 {
1709 output << "lhs = ";
1710 lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
1711 output << ";" << endl
1712 << "rhs = ";
1713 rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
1714 output << ";" << endl
1715 << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
1716 << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
1717 << RIGHT_ARRAY_SUBSCRIPT(output_type)
1718 << " = lhs - rhs;" << endl;
1719 }
1720 else // The right hand side of the equation is empty ==> residual=lhs;
1721 {
1722 if (isJuliaOutput(output_type))
1723 output << " @inbounds ";
1724 output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
1725 << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
1726 << RIGHT_ARRAY_SUBSCRIPT(output_type)
1727 << " = ";
1728 lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
1729 output << ";" << endl;
1730 }
1731 }
1732 }
1733
1734 void
compileModelEquations(ostream & code_file,unsigned int & instruction_number,const temporary_terms_t & tt,const map_idx_t & map_idx,bool dynamic,bool steady_dynamic) const1735 ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
1736 {
1737 for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
1738 {
1739 BinaryOpNode *eq_node = equations[eq];
1740 expr_t lhs = eq_node->arg1;
1741 expr_t rhs = eq_node->arg2;
1742 FNUMEXPR_ fnumexpr(ModelEquation, eq);
1743 fnumexpr.write(code_file, instruction_number);
1744 // Test if the right hand side of the equation is empty.
1745 double vrhs = 1.0;
1746 try
1747 {
1748 vrhs = rhs->eval(eval_context_t());
1749 }
1750 catch (ExprNode::EvalException &e)
1751 {
1752 }
1753
1754 if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
1755 {
1756 lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
1757 rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
1758
1759 FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
1760 fbinary.write(code_file, instruction_number);
1761
1762 FSTPR_ fstpr(eq);
1763 fstpr.write(code_file, instruction_number);
1764 }
1765 else // The right hand side of the equation is empty ==> residual=lhs;
1766 {
1767 lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
1768 FSTPR_ fstpr(eq);
1769 fstpr.write(code_file, instruction_number);
1770 }
1771 }
1772 }
1773
1774 void
Write_Inf_To_Bin_File(const string & filename,int & u_count_int,bool & file_open,bool is_two_boundaries,int block_mfs) const1775 ModelTree::Write_Inf_To_Bin_File(const string &filename,
1776 int &u_count_int, bool &file_open, bool is_two_boundaries, int block_mfs) const
1777 {
1778 int j;
1779 std::ofstream SaveCode;
1780 if (file_open)
1781 SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate);
1782 else
1783 SaveCode.open(filename, ios::out | ios::binary);
1784 if (!SaveCode.is_open())
1785 {
1786 cerr << R"(Error : Can't open file ")" << filename << R"(" for writing)" << endl;
1787 exit(EXIT_FAILURE);
1788 }
1789 u_count_int = 0;
1790 for (const auto & [indices, d1] : derivatives[1])
1791 {
1792 int deriv_id = indices[1];
1793 if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
1794 {
1795 int eq = indices[0];
1796 int symb = getSymbIDByDerivID(deriv_id);
1797 int var = symbol_table.getTypeSpecificID(symb);
1798 int lag = getLagByDerivID(deriv_id);
1799 SaveCode.write(reinterpret_cast<char *>(&eq), sizeof(eq));
1800 int varr = var + lag * block_mfs;
1801 SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
1802 SaveCode.write(reinterpret_cast<char *>(&lag), sizeof(lag));
1803 int u = u_count_int + block_mfs;
1804 SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
1805 u_count_int++;
1806 }
1807 }
1808 if (is_two_boundaries)
1809 u_count_int += symbol_table.endo_nbr();
1810 for (j = 0; j < static_cast<int>(symbol_table.endo_nbr()); j++)
1811 SaveCode.write(reinterpret_cast<char *>(&j), sizeof(j));
1812 for (j = 0; j < static_cast<int>(symbol_table.endo_nbr()); j++)
1813 SaveCode.write(reinterpret_cast<char *>(&j), sizeof(j));
1814 SaveCode.close();
1815 }
1816
1817 void
writeLatexModelFile(const string & mod_basename,const string & latex_basename,ExprNodeOutputType output_type,bool write_equation_tags) const1818 ModelTree::writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const
1819 {
1820 filesystem::create_directories(mod_basename + "/latex");
1821
1822 ofstream output, content_output;
1823 string filename = mod_basename + "/latex/" + latex_basename + ".tex";
1824 string content_filename = mod_basename + "/latex/" + latex_basename + "_content" + ".tex";
1825 output.open(filename, ios::out | ios::binary);
1826 if (!output.is_open())
1827 {
1828 cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
1829 exit(EXIT_FAILURE);
1830 }
1831
1832 content_output.open(content_filename, ios::out | ios::binary);
1833 if (!content_output.is_open())
1834 {
1835 cerr << "ERROR: Can't open file " << content_filename << " for writing" << endl;
1836 exit(EXIT_FAILURE);
1837 }
1838
1839 output << R"(\documentclass[10pt,a4paper]{article})" << endl
1840 << R"(\usepackage[landscape]{geometry})" << endl
1841 << R"(\usepackage{fullpage})" << endl
1842 << R"(\usepackage{amsfonts})" << endl
1843 << R"(\usepackage{breqn})" << endl
1844 << R"(\begin{document})" << endl
1845 << R"(\footnotesize)" << endl;
1846
1847 // Write model local variables
1848 for (int id : local_variables_vector)
1849 {
1850 expr_t value = local_variables_table.find(id)->second;
1851
1852 content_output << R"(\begin{dmath*})" << endl
1853 << symbol_table.getTeXName(id) << " = ";
1854 // Use an empty set for the temporary terms
1855 value->writeOutput(content_output, output_type);
1856 content_output << endl << R"(\end{dmath*})" << endl;
1857 }
1858
1859 for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
1860 {
1861 content_output << "% Equation " << eq + 1 << endl;
1862 if (write_equation_tags)
1863 {
1864 auto escape_special_latex_symbols
1865 = [](string str)
1866 {
1867 const regex special_latex_chars (R"([&%$#_{}])");
1868 const regex backslash (R"(\\)");
1869 const regex tilde (R"(~)");
1870 const regex carrot (R"(\^)");
1871 const regex textbackslash (R"(\\textbackslash)");
1872 str = regex_replace(str, backslash, R"(\textbackslash)");
1873 str = regex_replace(str, special_latex_chars, R"(\$&)");
1874 str = regex_replace(str, carrot, R"(\^{})");
1875 str = regex_replace(str, tilde, R"(\textasciitilde{})");
1876 return regex_replace(str, textbackslash, R"(\textbackslash{})");
1877 };
1878 bool wrote_eq_tag = false;
1879 for (const auto & [tagged_eq, tag_pair] : equation_tags)
1880 if (tagged_eq == eq)
1881 {
1882 if (!wrote_eq_tag)
1883 content_output << R"(\noindent[)";
1884 else
1885 content_output << ", ";
1886
1887 content_output << escape_special_latex_symbols(tag_pair.first);
1888
1889 if (!(tag_pair.second.empty()))
1890 content_output << "= `" << escape_special_latex_symbols(tag_pair.second) << "'";
1891
1892 wrote_eq_tag = true;
1893 }
1894
1895 if (wrote_eq_tag)
1896 content_output << "]" << endl;
1897 }
1898
1899 content_output << R"(\begin{dmath})" << endl;
1900 // Here it is necessary to cast to superclass ExprNode, otherwise the overloaded writeOutput() method is not found
1901 dynamic_cast<ExprNode *>(equations[eq])->writeOutput(content_output, output_type);
1902 content_output << endl << R"(\end{dmath})" << endl;
1903 }
1904
1905 output << R"(\include{)" << latex_basename + "_content" << "}" << endl
1906 << R"(\end{document})" << endl;
1907
1908 output.close();
1909 content_output.close();
1910 }
1911
1912 void
addEquation(expr_t eq,int lineno)1913 ModelTree::addEquation(expr_t eq, int lineno)
1914 {
1915 auto beq = dynamic_cast<BinaryOpNode *>(eq);
1916 assert(beq && beq->op_code == BinaryOpcode::equal);
1917
1918 equations.push_back(beq);
1919 equations_lineno.push_back(lineno);
1920 }
1921
1922 vector<int>
includeExcludeEquations(set<pair<string,string>> & eqs,bool exclude_eqs,vector<BinaryOpNode * > & equations,vector<int> & equations_lineno,vector<pair<int,pair<string,string>>> & equation_tags,multimap<pair<string,string>,int> & equation_tags_xref,bool static_equations) const1923 ModelTree::includeExcludeEquations(set<pair<string, string>> &eqs, bool exclude_eqs,
1924 vector<BinaryOpNode *> &equations, vector<int> &equations_lineno,
1925 vector<pair<int, pair<string, string>>> &equation_tags,
1926 multimap<pair<string, string>, int> &equation_tags_xref, bool static_equations) const
1927 {
1928 vector<int> excluded_vars;
1929 if (equations.empty())
1930 return excluded_vars;
1931
1932 // Get equation numbers of tags
1933 set<int> tag_eqns;
1934 for (auto &it : eqs)
1935 if (equation_tags_xref.find(it) != equation_tags_xref.end())
1936 {
1937 auto range = equation_tags_xref.equal_range(it);
1938 for_each(range.first, range.second, [&tag_eqns](auto &x) { tag_eqns.insert(x.second); });
1939 eqs.erase(it);
1940 }
1941 if (tag_eqns.empty())
1942 return excluded_vars;
1943
1944 set<int> eqns;
1945 if (exclude_eqs)
1946 eqns = tag_eqns;
1947 else
1948 for (size_t i = 0; i < equations.size(); i++)
1949 if (tag_eqns.find(i) == tag_eqns.end())
1950 eqns.insert(i);
1951
1952 // remove from equations, equations_lineno, equation_tags, equation_tags_xref
1953 vector<BinaryOpNode *> new_eqns;
1954 vector<int> new_equations_lineno;
1955 map<int, int> old_eqn_num_2_new;
1956 for (size_t i = 0; i < equations.size(); i++)
1957 if (eqns.find(i) != eqns.end())
1958 {
1959 bool found = false;
1960 for (const auto & [tagged_eq, tag_pair] : equation_tags)
1961 if (tagged_eq == static_cast<int>(i) && tag_pair.first == "endogenous")
1962 {
1963 found = true;
1964 excluded_vars.push_back(symbol_table.getID(tag_pair.second));
1965 break;
1966 }
1967 if (!found)
1968 {
1969 set<pair<int, int>> result;
1970 equations[i]->arg1->collectDynamicVariables(SymbolType::endogenous, result);
1971 if (result.size() == 1)
1972 excluded_vars.push_back(result.begin()->first);
1973 else
1974 {
1975 cerr << "ERROR: Equation " << i
1976 << " has been excluded but does not have a single variable on LHS or `endogenous` tag" << endl;
1977 exit(EXIT_FAILURE);
1978 }
1979 }
1980 }
1981 else
1982 {
1983 new_eqns.emplace_back(equations[i]);
1984 old_eqn_num_2_new[i] = new_eqns.size() - 1;
1985 new_equations_lineno.emplace_back(equations_lineno[i]);
1986 }
1987 int n_excl = equations.size() - new_eqns.size();
1988
1989 equations = new_eqns;
1990 equations_lineno = new_equations_lineno;
1991
1992 equation_tags.erase(remove_if(equation_tags.begin(), equation_tags.end(),
1993 [&](const auto &it) { return eqns.find(it.first) != eqns.end(); }),
1994 equation_tags.end());
1995 for (auto &it : old_eqn_num_2_new)
1996 for (auto &it1 : equation_tags)
1997 if (it1.first == it.first)
1998 it1.first = it.second;
1999
2000 equation_tags_xref.clear();
2001 for (const auto &it : equation_tags)
2002 equation_tags_xref.emplace(it.second, it.first);
2003
2004 if (!static_equations)
2005 for (size_t i = 0; i < excluded_vars.size(); i++)
2006 for (size_t j = i+1; j < excluded_vars.size(); j++)
2007 if (excluded_vars[i] == excluded_vars[j])
2008 {
2009 cerr << "Error: Variable " << symbol_table.getName(i) << " was excluded twice"
2010 << " via in/exclude_eqs option" << endl;
2011 exit(EXIT_FAILURE);
2012 }
2013
2014 cout << "Excluded " << n_excl << (static_equations ? " static " : " dynamic ")
2015 << "equation" << (n_excl > 1 ? "s" : "") << " via in/exclude_eqs option" << endl;
2016
2017 return excluded_vars;
2018 }
2019
2020 void
simplifyEquations()2021 ModelTree::simplifyEquations()
2022 {
2023 size_t last_subst_table_size = 0;
2024 map<VariableNode *, NumConstNode *> subst_table;
2025 // Equations with tags are excluded, in particular because of MCPs, see dynare#1697
2026 findConstantEquationsWithoutTags(subst_table);
2027 while (subst_table.size() != last_subst_table_size)
2028 {
2029 last_subst_table_size = subst_table.size();
2030 for (auto &[id, definition] : local_variables_table)
2031 definition = definition->replaceVarsInEquation(subst_table);
2032 for (auto &equation : equations)
2033 equation = dynamic_cast<BinaryOpNode *>(equation->replaceVarsInEquation(subst_table));
2034 subst_table.clear();
2035 findConstantEquationsWithoutTags(subst_table);
2036 }
2037 }
2038
2039 void
findConstantEquationsWithoutTags(map<VariableNode *,NumConstNode * > & subst_table) const2040 ModelTree::findConstantEquationsWithoutTags(map<VariableNode *, NumConstNode *> &subst_table) const
2041 {
2042 for (size_t i = 0; i < equations.size(); i++)
2043 if (getEquationTags(i).empty())
2044 equations[i]->findConstantEquations(subst_table);
2045 }
2046
2047 void
addEquation(expr_t eq,int lineno,const vector<pair<string,string>> & eq_tags)2048 ModelTree::addEquation(expr_t eq, int lineno, const vector<pair<string, string>> &eq_tags)
2049 {
2050 int n = equations.size();
2051 for (const auto &eq_tag : eq_tags)
2052 {
2053 equation_tags.emplace_back(n, eq_tag);
2054 equation_tags_xref.emplace(eq_tag, n);
2055 }
2056 addEquation(eq, lineno);
2057 }
2058
2059 void
addAuxEquation(expr_t eq)2060 ModelTree::addAuxEquation(expr_t eq)
2061 {
2062 auto beq = dynamic_cast<BinaryOpNode *>(eq);
2063 assert(beq && beq->op_code == BinaryOpcode::equal);
2064
2065 aux_equations.push_back(beq);
2066 }
2067
2068 void
addTrendVariables(const vector<int> & trend_vars,expr_t growth_factor)2069 ModelTree::addTrendVariables(const vector<int> &trend_vars, expr_t growth_factor) noexcept(false)
2070 {
2071 for (int id : trend_vars)
2072 if (trend_symbols_map.find(id) != trend_symbols_map.end())
2073 throw TrendException(symbol_table.getName(id));
2074 else
2075 trend_symbols_map[id] = growth_factor;
2076 }
2077
2078 void
addNonstationaryVariables(const vector<int> & nonstationary_vars,bool log_deflator,expr_t deflator)2079 ModelTree::addNonstationaryVariables(const vector<int> &nonstationary_vars, bool log_deflator, expr_t deflator) noexcept(false)
2080 {
2081 for (int id : nonstationary_vars)
2082 if (nonstationary_symbols_map.find(id) != nonstationary_symbols_map.end())
2083 throw TrendException(symbol_table.getName(id));
2084 else
2085 nonstationary_symbols_map[id] = { log_deflator, deflator };
2086 }
2087
2088 void
initializeVariablesAndEquations()2089 ModelTree::initializeVariablesAndEquations()
2090 {
2091 for (size_t j = 0; j < equations.size(); j++)
2092 equation_reordered.push_back(j);
2093
2094 for (int j = 0; j < symbol_table.endo_nbr(); j++)
2095 variable_reordered.push_back(j);
2096 }
2097
2098 void
set_cutoff_to_zero()2099 ModelTree::set_cutoff_to_zero()
2100 {
2101 cutoff = 0;
2102 }
2103
2104 void
jacobianHelper(ostream & output,int eq_nb,int col_nb,ExprNodeOutputType output_type) const2105 ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
2106 {
2107 if (isJuliaOutput(output_type))
2108 output << " @inbounds ";
2109 output << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
2110 if (isMatlabOutput(output_type) || isJuliaOutput(output_type))
2111 output << eq_nb + 1 << "," << col_nb + 1;
2112 else
2113 output << eq_nb + col_nb *equations.size();
2114 output << RIGHT_ARRAY_SUBSCRIPT(output_type);
2115 }
2116
2117 void
sparseHelper(int order,ostream & output,int row_nb,int col_nb,ExprNodeOutputType output_type) const2118 ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
2119 {
2120 output << "v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
2121 if (isMatlabOutput(output_type) || isJuliaOutput(output_type))
2122 output << row_nb + 1 << "," << col_nb + 1;
2123 else
2124 output << row_nb + col_nb * NNZDerivatives[order];
2125 output << RIGHT_ARRAY_SUBSCRIPT(output_type);
2126 }
2127
2128 void
computeParamsDerivatives(int paramsDerivsOrder)2129 ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
2130 {
2131 assert(paramsDerivsOrder >= 1);
2132
2133 set<int> deriv_id_set;
2134 addAllParamDerivId(deriv_id_set);
2135
2136 // First-order derivatives w.r.t. params
2137 for (int param : deriv_id_set)
2138 {
2139 for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
2140 {
2141 expr_t d = equations[eq]->getDerivative(param);
2142 if (d == Zero)
2143 continue;
2144 params_derivatives[{ 0, 1 }][{ eq, param }] = d;
2145 }
2146
2147 for (int endoOrd = 1; endoOrd < static_cast<int>(derivatives.size()); endoOrd++)
2148 for (const auto &[indices, dprev] : derivatives[endoOrd])
2149 {
2150 expr_t d = dprev->getDerivative(param);
2151 if (d == Zero)
2152 continue;
2153 vector<int> new_indices = indices;
2154 new_indices.push_back(param);
2155 params_derivatives[{ endoOrd, 1 }][new_indices] = d;
2156 }
2157 }
2158
2159 // Higher-order derivatives w.r.t. parameters
2160 for (int endoOrd = 0; endoOrd < static_cast<int>(derivatives.size()); endoOrd++)
2161 for (int paramOrd = 2; paramOrd <= paramsDerivsOrder; paramOrd++)
2162 for (const auto &[indices, dprev] : params_derivatives[{ endoOrd, paramOrd-1 }])
2163 for (int param : deriv_id_set)
2164 {
2165 if (indices.back() > param)
2166 continue;
2167
2168 expr_t d = dprev->getDerivative(param);
2169 if (d == Zero)
2170 continue;
2171 vector<int> new_indices = indices;
2172 new_indices.push_back(param);
2173 // At this point, indices of both endogenous and parameters are sorted in non-decreasing order
2174 params_derivatives[{ endoOrd, paramOrd }][new_indices] = d;
2175 }
2176 }
2177
2178 void
computeParamsDerivativesTemporaryTerms()2179 ModelTree::computeParamsDerivativesTemporaryTerms()
2180 {
2181 map<expr_t, pair<int, pair<int, int>>> reference_count;
2182
2183 /* The temp terms should be constructed in the same order as the for loops in
2184 {Static,Dynamic}Model::write{Json,}ParamsDerivativesFile() */
2185 params_derivs_temporary_terms.clear();
2186 for (const auto &[order, derivs] : params_derivatives)
2187 for (const auto &[indices, d] : derivs)
2188 d->computeTemporaryTerms(order, params_derivs_temporary_terms,
2189 reference_count, true);
2190
2191 int idx = 0;
2192 for (auto &[mlv, value] : temporary_terms_mlv)
2193 params_derivs_temporary_terms_idxs[mlv] = idx++;
2194 for (const auto &[order, tts] : params_derivs_temporary_terms)
2195 for (const auto &tt : tts)
2196 params_derivs_temporary_terms_idxs[tt] = idx++;
2197 }
2198
2199 bool
isNonstationary(int symb_id) const2200 ModelTree::isNonstationary(int symb_id) const
2201 {
2202 return nonstationary_symbols_map.find(symb_id) != nonstationary_symbols_map.end();
2203 }
2204
2205 void
writeJsonModelEquations(ostream & output,bool residuals) const2206 ModelTree::writeJsonModelEquations(ostream &output, bool residuals) const
2207 {
2208 if (residuals)
2209 output << endl << R"("residuals":[)" << endl;
2210 else
2211 output << endl << R"("model":[)" << endl;
2212 for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
2213 {
2214 if (eq > 0)
2215 output << ", ";
2216
2217 BinaryOpNode *eq_node = equations[eq];
2218 expr_t lhs = eq_node->arg1;
2219 expr_t rhs = eq_node->arg2;
2220
2221 if (residuals)
2222 {
2223 output << R"({"residual": {)"
2224 << R"("lhs": ")";
2225 lhs->writeJsonOutput(output, temporary_terms, {});
2226 output << R"(")";
2227
2228 output << R"(, "rhs": ")";
2229 rhs->writeJsonOutput(output, temporary_terms, {});
2230 output << R"(")";
2231 try
2232 {
2233 // Test if the right hand side of the equation is empty.
2234 if (rhs->eval(eval_context_t()) != 0)
2235 {
2236 output << R"(, "rhs": ")";
2237 rhs->writeJsonOutput(output, temporary_terms, {});
2238 output << R"(")";
2239 }
2240 }
2241 catch (ExprNode::EvalException &e)
2242 {
2243 }
2244 output << "}";
2245 }
2246 else
2247 {
2248 output << R"({"lhs": ")";
2249 lhs->writeJsonOutput(output, {}, {});
2250 output << R"(", "rhs": ")";
2251 rhs->writeJsonOutput(output, {}, {});
2252 output << R"(")"
2253 << R"(, "line": )" << equations_lineno[eq];
2254
2255 if (auto eqtags = getEquationTags(eq);
2256 !eqtags.empty())
2257 {
2258 output << R"(, "tags": {)";
2259 int i = 0;
2260 for (const auto &[name, value] : eqtags)
2261 {
2262 if (i != 0)
2263 output << ", ";
2264 output << R"(")" << name << R"(": ")" << value << R"(")";
2265 i++;
2266 }
2267 output << "}";
2268 eqtags.clear();
2269 }
2270 }
2271 output << "}" << endl;
2272 }
2273 output << endl << "]" << endl;
2274 }
2275
2276 string
matlab_arch(const string & mexext)2277 ModelTree::matlab_arch(const string &mexext)
2278 {
2279 if (mexext == "mexglx")
2280 return "glnx86";
2281 else if (mexext == "mexa64")
2282 return "glnxa64";
2283 if (mexext == "mexw32")
2284 return "win32";
2285 else if (mexext == "mexw64")
2286 return "win64";
2287 else if (mexext == "mexmaci")
2288 {
2289 cerr << "32-bit MATLAB not supported on macOS" << endl;
2290 exit(EXIT_FAILURE);
2291 }
2292 else if (mexext == "mexmaci64")
2293 return "maci64";
2294 else
2295 {
2296 cerr << "ERROR: 'mexext' option to preprocessor incorrectly set, needed with 'use_dll'" << endl;
2297 exit(EXIT_FAILURE);
2298 }
2299 }
2300
2301 void
compileDll(const string & basename,const string & static_or_dynamic,const string & mexext,const filesystem::path & matlabroot,const filesystem::path & dynareroot) const2302 ModelTree::compileDll(const string &basename, const string &static_or_dynamic, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const
2303 {
2304 const string opt_flags = "-O3 -g0 --param ira-max-conflict-table-size=1 -fno-forward-propagate -fno-gcse -fno-dce -fno-dse -fno-tree-fre -fno-tree-pre -fno-tree-cselim -fno-tree-dse -fno-tree-dce -fno-tree-pta -fno-gcse-after-reload";
2305
2306 filesystem::path compiler;
2307 ostringstream flags;
2308 string libs;
2309
2310 if (mexext == "mex")
2311 {
2312 // Octave
2313 compiler = matlabroot / "bin" / "mkoctfile";
2314 flags << "--mex";
2315 }
2316 else
2317 {
2318 // MATLAB
2319 compiler = "gcc";
2320 string arch = matlab_arch(mexext);
2321 auto include_dir = matlabroot / "extern" / "include";
2322 flags << "-I " << include_dir;
2323 auto bin_dir = matlabroot / "bin" / arch;
2324 flags << " -L " << bin_dir;
2325 flags << " -fexceptions -DNDEBUG";
2326 libs = "-lmex -lmx";
2327 if (mexext == "mexglx" || mexext == "mexa64")
2328 {
2329 // GNU/Linux
2330 flags << " -D_GNU_SOURCE -fPIC -pthread"
2331 << " -shared -Wl,--no-undefined -Wl,-rpath-link," << bin_dir;
2332 libs += " -lm -lstdc++";
2333
2334 if (mexext == "mexglx")
2335 flags << " -D_FILE_OFFSET_BITS=64 -m32";
2336 else
2337 flags << " -fno-omit-frame-pointer";
2338 }
2339 else if (mexext == "mexw32" || mexext == "mexw64")
2340 {
2341 // Windows
2342 flags << " -static-libgcc -static-libstdc++ -shared";
2343 // Put the MinGW environment shipped with Dynare in the path
2344 auto mingwpath = dynareroot / (string{"mingw"} + (mexext == "mexw32" ? "32" : "64")) / "bin";
2345 string newpath = "PATH=" + mingwpath.string() + ';' + string{getenv("PATH")};
2346 if (putenv(const_cast<char *>(newpath.c_str())) != 0)
2347 {
2348 cerr << "Can't set PATH" << endl;
2349 exit(EXIT_FAILURE);
2350 }
2351 }
2352 else
2353 {
2354 // macOS
2355 #ifdef __APPLE__
2356 char dynare_m_path[PATH_MAX];
2357 uint32_t size = PATH_MAX;
2358 string gcc_relative_path = "";
2359 if (_NSGetExecutablePath(dynare_m_path, &size) == 0)
2360 {
2361 string str = dynare_m_path;
2362 gcc_relative_path = str.substr(0, str.find_last_of("/")) + "/../../.brew/bin/gcc-9";
2363 }
2364
2365 if (filesystem::exists(gcc_relative_path))
2366 compiler = gcc_relative_path;
2367 else if (filesystem::exists("/usr/local/bin/gcc-9"))
2368 compiler = "/usr/local/bin/gcc-9";
2369 else
2370 {
2371 cerr << "ERROR: You must install gcc-9 on your system before using the `use_dll` option of Dynare. "
2372 << "You can do this via the Dynare installation package." << endl;
2373 exit(EXIT_FAILURE);
2374 }
2375 #endif
2376 flags << " -fno-common -arch x86_64 -mmacosx-version-min=10.9 -Wl,-twolevel_namespace -undefined error -bundle";
2377 libs += " -lm -lstdc++";
2378 }
2379 }
2380
2381 auto model_dir = filesystem::path{basename} / "model" / "src";
2382 filesystem::path main_src{model_dir / (static_or_dynamic + ".c")},
2383 mex_src{model_dir / (static_or_dynamic + "_mex.c")};
2384
2385 filesystem::path mex_dir{"+" + basename};
2386 filesystem::path binary{mex_dir / (static_or_dynamic + "." + mexext)};
2387
2388 ostringstream cmd;
2389
2390 #ifdef _WIN32
2391 /* On Windows, system() hands the command over to "cmd.exe /C". We need to
2392 enclose the whole command line within double quotes if we want the inner
2393 quotes to be correctly handled. See "cmd /?" for more details. */
2394 cmd << '"';
2395 #endif
2396
2397 if (user_set_compiler.empty())
2398 cmd << compiler << " ";
2399 else
2400 if (!filesystem::exists(user_set_compiler))
2401 {
2402 cerr << "Error: The specified compiler '" << user_set_compiler << "' cannot be found on your system" << endl;
2403 exit(EXIT_FAILURE);
2404 }
2405 else
2406 cmd << user_set_compiler << " ";
2407
2408 if (user_set_subst_flags.empty())
2409 cmd << opt_flags << " " << flags.str() << " ";
2410 else
2411 cmd << user_set_subst_flags << " ";
2412
2413 if (!user_set_add_flags.empty())
2414 cmd << user_set_add_flags << " ";
2415
2416 cmd << main_src << " " << mex_src << " -o " << binary << " ";
2417
2418 if (user_set_subst_libs.empty())
2419 cmd << libs;
2420 else
2421 cmd << user_set_subst_libs;
2422
2423 if (!user_set_add_libs.empty())
2424 cmd << " " << user_set_add_libs;
2425
2426 #ifdef _WIN32
2427 cmd << '"';
2428 #endif
2429
2430 cout << "Compiling " << static_or_dynamic << " MEX..." << endl << cmd.str() << endl;
2431
2432 if (system(cmd.str().c_str()))
2433 {
2434 cerr << "Compilation failed" << endl;
2435 exit(EXIT_FAILURE);
2436 }
2437 }
2438
2439 void
reorderAuxiliaryEquations()2440 ModelTree::reorderAuxiliaryEquations()
2441 {
2442 using namespace boost;
2443
2444 // Create the mapping between auxiliary variables and auxiliary equations
2445 int n = static_cast<int>(aux_equations.size());
2446 map<int, int> auxEndoToEq;
2447 for (int i = 0; i < n; i++)
2448 {
2449 auto varexpr = dynamic_cast<VariableNode *>(aux_equations[i]->arg1);
2450 assert(varexpr && symbol_table.getType(varexpr->symb_id) == SymbolType::endogenous);
2451 auxEndoToEq[varexpr->symb_id] = i;
2452 }
2453 assert(static_cast<int>(auxEndoToEq.size()) == n);
2454
2455 /* Construct the directed acyclic graph where auxiliary equations are
2456 vertices and edges represent dependency relationships. */
2457 using Graph = adjacency_list<vecS, vecS, directedS>;
2458 Graph g(n);
2459 for (int i = 0; i < n; i++)
2460 {
2461 set<int> endos;
2462 aux_equations[i]->collectVariables(SymbolType::endogenous, endos);
2463 for (int endo : endos)
2464 if (auto it = auxEndoToEq.find(endo);
2465 it != auxEndoToEq.end() && it->second != i)
2466 add_edge(i, it->second, g);
2467 }
2468
2469 // Topological sort of the graph
2470 using Vertex = graph_traits<Graph>::vertex_descriptor;
2471 vector<Vertex> ordered;
2472 topological_sort(g, back_inserter(ordered));
2473
2474 // Reorder auxiliary equations accordingly
2475 auto aux_equations_old = aux_equations;
2476 auto index = get(vertex_index, g); // Maps vertex descriptors to their index
2477 for (int i = 0; i < n; i++)
2478 aux_equations[i] = aux_equations_old[index[ordered[i]]];
2479 }
2480