1 /*
2 * Copyright © 2003-2019 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 <iostream>
21 #include <cmath>
22 #include <cstdlib>
23 #include <cassert>
24 #include <algorithm>
25
26 #include "StaticModel.hh"
27 #include "DynamicModel.hh"
28
29 void
copyHelper(const StaticModel & m)30 StaticModel::copyHelper(const StaticModel &m)
31 {
32 auto f = [this](expr_t e) { return e->clone(*this); };
33
34 auto convert_vector_tt = [f](vector<temporary_terms_t> vtt)
35 {
36 vector<temporary_terms_t> vtt2;
37 for (const auto &tt : vtt)
38 {
39 temporary_terms_t tt2;
40 for (const auto &it : tt)
41 tt2.insert(f(it));
42 vtt2.push_back(tt2);
43 }
44 return vtt2;
45 };
46
47 for (const auto &it : m.v_temporary_terms)
48 v_temporary_terms.push_back(convert_vector_tt(it));
49 for (const auto &it : m.v_temporary_terms_local)
50 v_temporary_terms_local.push_back(convert_vector_tt(it));
51
52 for (const auto &it : m.first_chain_rule_derivatives)
53 first_chain_rule_derivatives[it.first] = f(it.second);
54
55 for (const auto &it : m.equation_type_and_normalized_equation)
56 equation_type_and_normalized_equation.emplace_back(it.first, f(it.second));
57
58 for (const auto &it : m.blocks_derivatives)
59 {
60 block_derivatives_equation_variable_laglead_nodeid_t v;
61 for (const auto &it2 : it)
62 v.emplace_back(get<0>(it2), get<1>(it2), get<2>(it2), f(get<3>(it2)));
63 blocks_derivatives.push_back(v);
64 }
65
66 for (const auto &it : m.dynamic_jacobian)
67 dynamic_jacobian[it.first] = f(it.second);
68
69 auto convert_derivative_t = [f](derivative_t dt)
70 {
71 derivative_t dt2;
72 for (const auto &it : dt)
73 dt2[it.first] = f(it.second);
74 return dt2;
75 };
76 for (const auto &it : m.derivative_endo)
77 derivative_endo.push_back(convert_derivative_t(it));
78 for (const auto &it : m.derivative_other_endo)
79 derivative_other_endo.push_back(convert_derivative_t(it));
80 for (const auto &it : m.derivative_exo)
81 derivative_exo.push_back(convert_derivative_t(it));
82 for (const auto &it : m.derivative_exo_det)
83 derivative_exo_det.push_back(convert_derivative_t(it));
84 }
85
StaticModel(SymbolTable & symbol_table_arg,NumericalConstants & num_constants_arg,ExternalFunctionsTable & external_functions_table_arg)86 StaticModel::StaticModel(SymbolTable &symbol_table_arg,
87 NumericalConstants &num_constants_arg,
88 ExternalFunctionsTable &external_functions_table_arg) :
89 ModelTree{symbol_table_arg, num_constants_arg, external_functions_table_arg}
90 {
91 }
92
StaticModel(const StaticModel & m)93 StaticModel::StaticModel(const StaticModel &m) :
94 ModelTree{m},
95 v_temporary_terms_inuse{m.v_temporary_terms_inuse},
96 map_idx{m.map_idx},
97 map_idx2{m.map_idx2},
98 global_temporary_terms{m.global_temporary_terms},
99 block_type_firstequation_size_mfs{m.block_type_firstequation_size_mfs},
100 blocks_linear{m.blocks_linear},
101 block_col_type{m.block_col_type},
102 endo_max_leadlag_block{m.endo_max_leadlag_block},
103 other_endo_max_leadlag_block{m.other_endo_max_leadlag_block},
104 exo_max_leadlag_block{m.exo_max_leadlag_block},
105 exo_det_max_leadlag_block{m.exo_det_max_leadlag_block},
106 max_leadlag_block{m.max_leadlag_block}
107 {
108 copyHelper(m);
109 }
110
111 StaticModel &
operator =(const StaticModel & m)112 StaticModel::operator=(const StaticModel &m)
113 {
114 ModelTree::operator=(m);
115
116 v_temporary_terms.clear();
117 v_temporary_terms_local.clear();
118
119 v_temporary_terms_inuse = m.v_temporary_terms_inuse;
120
121 first_chain_rule_derivatives.clear();
122
123 map_idx = m.map_idx;
124 map_idx2 = m.map_idx2;
125 global_temporary_terms = m.global_temporary_terms;
126
127 equation_type_and_normalized_equation.clear();
128
129 block_type_firstequation_size_mfs = m.block_type_firstequation_size_mfs;
130
131 blocks_derivatives.clear();
132 dynamic_jacobian.clear();
133
134 blocks_linear = m.blocks_linear;
135
136 derivative_endo.clear();
137 derivative_other_endo.clear();
138 derivative_exo.clear();
139 derivative_exo_det.clear();
140
141 block_col_type = m.block_col_type;
142 endo_max_leadlag_block = m.endo_max_leadlag_block;
143 other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
144 exo_max_leadlag_block = m.exo_max_leadlag_block;
145 exo_det_max_leadlag_block = m.exo_det_max_leadlag_block;
146 max_leadlag_block = m.max_leadlag_block;
147
148 copyHelper(m);
149
150 return *this;
151 }
152
StaticModel(const DynamicModel & m)153 StaticModel::StaticModel(const DynamicModel &m) :
154 ModelTree{m.symbol_table, m.num_constants, m.external_functions_table}
155 {
156 // Convert model local variables (need to be done first)
157 for (int it : m.local_variables_vector)
158 AddLocalVariable(it, m.local_variables_table.find(it)->second->toStatic(*this));
159
160 // Convert equations
161 int static_only_index = 0;
162 for (int i = 0; i < static_cast<int>(m.equations.size()); i++)
163 {
164 // Detect if equation is marked [dynamic]
165 bool is_dynamic_only = false;
166 vector<pair<string, string>> eq_tags;
167 for (const auto & [tagged_eq, tag_pair] : m.equation_tags)
168 if (tagged_eq == i)
169 {
170 eq_tags.push_back(tag_pair);
171 if (tag_pair.first == "dynamic")
172 is_dynamic_only = true;
173 }
174
175 try
176 {
177 // If yes, replace it by an equation marked [static]
178 if (is_dynamic_only)
179 {
180 auto [static_only_equations, static_only_equations_lineno, static_only_equations_equation_tags] = m.getStaticOnlyEquationsInfo();
181
182 addEquation(static_only_equations[static_only_index]->toStatic(*this), static_only_equations_lineno[static_only_index], static_only_equations_equation_tags[static_only_index]);
183 static_only_index++;
184 }
185 else
186 addEquation(m.equations[i]->toStatic(*this), m.equations_lineno[i], eq_tags);
187 }
188 catch (DataTree::DivisionByZeroException)
189 {
190 cerr << "...division by zero error encountred when converting equation " << i << " to static" << endl;
191 exit(EXIT_FAILURE);
192 }
193 }
194
195 // Convert auxiliary equations
196 for (auto aux_eq : m.aux_equations)
197 addAuxEquation(aux_eq->toStatic(*this));
198
199 user_set_add_flags = m.user_set_add_flags;
200 user_set_subst_flags = m.user_set_subst_flags;
201 user_set_add_libs = m.user_set_add_libs;
202 user_set_subst_libs = m.user_set_subst_libs;
203 user_set_compiler = m.user_set_compiler;
204 }
205
206 void
compileDerivative(ofstream & code_file,unsigned int & instruction_number,int eq,int symb_id,map_idx_t & map_idx,temporary_terms_t temporary_terms) const207 StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
208 {
209 if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) });
210 it != derivatives[1].end())
211 it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
212 else
213 {
214 FLDZ_ fldz;
215 fldz.write(code_file, instruction_number);
216 }
217 }
218
219 void
compileChainRuleDerivative(ofstream & code_file,unsigned int & instruction_number,int eqr,int varr,int lag,map_idx_t & map_idx,temporary_terms_t temporary_terms) const220 StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
221 {
222 if (auto it = first_chain_rule_derivatives.find({ eqr, varr, lag });
223 it != first_chain_rule_derivatives.end())
224 it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
225 else
226 {
227 FLDZ_ fldz;
228 fldz.write(code_file, instruction_number);
229 }
230 }
231
232 void
computeTemporaryTermsOrdered()233 StaticModel::computeTemporaryTermsOrdered()
234 {
235 map<expr_t, pair<int, int>> first_occurence;
236 map<expr_t, int> reference_count;
237 BinaryOpNode *eq_node;
238 first_chain_rule_derivatives_t::const_iterator it_chr;
239 ostringstream tmp_s;
240 map_idx.clear();
241
242 unsigned int nb_blocks = getNbBlocks();
243 v_temporary_terms = vector< vector<temporary_terms_t>>(nb_blocks);
244 v_temporary_terms_local = vector< vector<temporary_terms_t>>(nb_blocks);
245
246 v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
247
248 map_idx2 = vector<map_idx_t>(nb_blocks);
249
250 temporary_terms.clear();
251
252 //local temporay terms
253 for (unsigned int block = 0; block < nb_blocks; block++)
254 {
255 map<expr_t, int> reference_count_local;
256 map<expr_t, pair<int, int>> first_occurence_local;
257 temporary_terms_t temporary_terms_l;
258
259 unsigned int block_size = getBlockSize(block);
260 unsigned int block_nb_mfs = getBlockMfs(block);
261 unsigned int block_nb_recursives = block_size - block_nb_mfs;
262 v_temporary_terms_local[block] = vector<temporary_terms_t>(block_size);
263
264 for (unsigned int i = 0; i < block_size; i++)
265 {
266 if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
267 getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, i);
268 else
269 {
270 eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
271 eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, i);
272 }
273 }
274 for (const auto &it : blocks_derivatives[block])
275 {
276 expr_t id = get<3>(it);
277 id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, block_size-1);
278 }
279 v_temporary_terms_inuse[block] = {};
280 computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]);
281 }
282
283 // global temporay terms
284 for (unsigned int block = 0; block < nb_blocks; block++)
285 {
286 // Compute the temporary terms reordered
287 unsigned int block_size = getBlockSize(block);
288 unsigned int block_nb_mfs = getBlockMfs(block);
289 unsigned int block_nb_recursives = block_size - block_nb_mfs;
290 v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
291 for (unsigned int i = 0; i < block_size; i++)
292 {
293 if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
294 getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
295 else
296 {
297 eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
298 eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
299 }
300 }
301 for (const auto &it : blocks_derivatives[block])
302 {
303 expr_t id = get<3>(it);
304 id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
305 }
306 }
307
308 for (unsigned int block = 0; block < nb_blocks; block++)
309 {
310 // Collecte the temporary terms reordered
311 unsigned int block_size = getBlockSize(block);
312 unsigned int block_nb_mfs = getBlockMfs(block);
313 unsigned int block_nb_recursives = block_size - block_nb_mfs;
314 set<int> temporary_terms_in_use;
315 for (unsigned int i = 0; i < block_size; i++)
316 {
317 if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
318 getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
319 else
320 {
321 eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
322 eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
323 }
324 }
325 for (const auto &it : blocks_derivatives[block])
326 {
327 expr_t id = get<3>(it);
328 id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
329 }
330 for (int i = 0; i < static_cast<int>(getBlockSize(block)); i++)
331 for (const auto &it : v_temporary_terms[block][i])
332 it->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
333 v_temporary_terms_inuse[block] = temporary_terms_in_use;
334 }
335 computeTemporaryTermsMapping(temporary_terms, map_idx);
336 }
337
338 void
computeTemporaryTermsMapping(temporary_terms_t & temporary_terms,map_idx_t & map_idx)339 StaticModel::computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, map_idx_t &map_idx)
340 {
341 // Add a mapping form node ID to temporary terms order
342 int j = 0;
343 for (auto temporary_term : temporary_terms)
344 map_idx[temporary_term->idx] = j++;
345 }
346
347 void
writeModelEquationsOrdered_M(const string & basename) const348 StaticModel::writeModelEquationsOrdered_M(const string &basename) const
349 {
350 string tmp_s, sps;
351 ostringstream tmp_output, tmp1_output, global_output;
352 expr_t lhs = nullptr, rhs = nullptr;
353 BinaryOpNode *eq_node;
354 map<expr_t, int> reference_count;
355 temporary_terms_t local_temporary_terms;
356 ofstream output;
357 vector<int> feedback_variables;
358 deriv_node_temp_terms_t tef_terms;
359 ExprNodeOutputType local_output_type;
360
361 local_output_type = ExprNodeOutputType::matlabStaticModelSparse;
362 if (global_temporary_terms)
363 local_temporary_terms = temporary_terms;
364
365 //----------------------------------------------------------------------
366 //For each block
367 for (unsigned int block = 0; block < getNbBlocks(); block++)
368 {
369 //recursive_variables.clear();
370 feedback_variables.clear();
371 //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
372 BlockSimulationType simulation_type = getBlockSimulationType(block);
373 unsigned int block_size = getBlockSize(block);
374 unsigned int block_mfs = getBlockMfs(block);
375 unsigned int block_recursive = block_size - block_mfs;
376
377 tmp1_output.str("");
378 tmp1_output << packageDir(basename + ".block") << "/static_" << block+1 << ".m";
379 output.open(tmp1_output.str(), ios::out | ios::binary);
380 output << "%" << endl
381 << "% " << tmp1_output.str() << " : Computes static model for Dynare" << endl
382 << "%" << endl
383 << "% Warning : this file is generated automatically by Dynare" << endl
384 << "% from model file (.mod)" << endl << endl
385 << "%/" << endl;
386 if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
387 output << "function y = static_" << block+1 << "(y, x, params)" << endl;
388 else
389 output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)" << endl;
390
391 BlockType block_type;
392 if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
393 block_type = SIMULTANS;
394 else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
395 || simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
396 && getBlockFirstEquation(block) < prologue)
397 block_type = PROLOGUE;
398 else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
399 || simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
400 && getBlockFirstEquation(block) >= equations.size() - epilogue)
401 block_type = EPILOGUE;
402 else
403 block_type = SIMULTANS;
404 output << " % ////////////////////////////////////////////////////////////////////////" << endl
405 << " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
406 << " //" << endl
407 << " % // Simulation type "
408 << BlockSim(simulation_type) << " //" << endl
409 << " % ////////////////////////////////////////////////////////////////////////" << endl;
410 output << " global options_;" << endl;
411 //The Temporary terms
412 if (simulation_type != EVALUATE_BACKWARD && simulation_type != EVALUATE_FORWARD)
413 output << " g1 = spalloc(" << block_mfs << ", " << block_mfs << ", " << blocks_derivatives[block].size() << ");" << endl;
414
415 if (v_temporary_terms_inuse[block].size())
416 {
417 tmp_output.str("");
418 for (int it : v_temporary_terms_inuse[block])
419 tmp_output << " T" << it;
420 output << " global" << tmp_output.str() << ";" << endl;
421 }
422
423 if (simulation_type != EVALUATE_BACKWARD && simulation_type != EVALUATE_FORWARD)
424 output << " residual=zeros(" << block_mfs << ",1);" << endl;
425
426 // The equations
427 temporary_terms_idxs_t temporary_terms_idxs;
428 for (unsigned int i = 0; i < block_size; i++)
429 {
430 if (!global_temporary_terms)
431 local_temporary_terms = v_temporary_terms[block][i];
432 temporary_terms_t tt2;
433 if (v_temporary_terms[block].size())
434 {
435 output << " " << "% //Temporary variables" << endl;
436 for (auto it : v_temporary_terms[block][i])
437 {
438 if (dynamic_cast<AbstractExternalFunctionNode *>(it))
439 it->writeExternalFunctionOutput(output, local_output_type, tt2, {}, tef_terms);
440
441 output << " " << sps;
442 it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
443 output << " = ";
444 it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
445 // Insert current node into tt2
446 tt2.insert(it);
447 output << ";" << endl;
448 }
449 }
450
451 int variable_ID = getBlockVariableID(block, i);
452 int equation_ID = getBlockEquationID(block, i);
453 EquationType equ_type = getBlockEquationType(block, i);
454 string sModel = symbol_table.getName(symbol_table.getID(SymbolType::endogenous, variable_ID));
455 eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
456 lhs = eq_node->arg1;
457 rhs = eq_node->arg2;
458 tmp_output.str("");
459 lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
460 switch (simulation_type)
461 {
462 case EVALUATE_BACKWARD:
463 case EVALUATE_FORWARD:
464 evaluation:
465 output << " % equation " << getBlockEquationID(block, i)+1 << " variable : " << sModel
466 << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
467 output << " ";
468 if (equ_type == E_EVALUATE)
469 {
470 output << tmp_output.str();
471 output << " = ";
472 rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
473 }
474 else if (equ_type == E_EVALUATE_S)
475 {
476 output << "%" << tmp_output.str();
477 output << " = ";
478 if (isBlockEquationRenormalized(block, i))
479 {
480 rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
481 output << endl << " ";
482 tmp_output.str("");
483 eq_node = static_cast<BinaryOpNode *>(getBlockEquationRenormalizedExpr(block, i));
484 lhs = eq_node->arg1;
485 rhs = eq_node->arg2;
486 lhs->writeOutput(output, local_output_type, local_temporary_terms, {});
487 output << " = ";
488 rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
489 }
490 }
491 else
492 {
493 cerr << "Type mismatch for equation " << equation_ID+1 << endl;
494 exit(EXIT_FAILURE);
495 }
496 output << ";" << endl;
497 break;
498 case SOLVE_BACKWARD_SIMPLE:
499 case SOLVE_FORWARD_SIMPLE:
500 case SOLVE_BACKWARD_COMPLETE:
501 case SOLVE_FORWARD_COMPLETE:
502 if (i < block_recursive)
503 goto evaluation;
504 feedback_variables.push_back(variable_ID);
505 output << " % equation " << equation_ID+1 << " variable : " << sModel
506 << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
507 output << " " << "residual(" << i+1-block_recursive << ") = (";
508 goto end;
509 default:
510 end:
511 output << tmp_output.str();
512 output << ") - (";
513 rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
514 output << ");" << endl;
515 }
516 }
517 // The Jacobian if we have to solve the block
518 if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE
519 || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
520 output << " " << sps << "% Jacobian " << endl;
521 switch (simulation_type)
522 {
523 case SOLVE_BACKWARD_SIMPLE:
524 case SOLVE_FORWARD_SIMPLE:
525 case SOLVE_BACKWARD_COMPLETE:
526 case SOLVE_FORWARD_COMPLETE:
527 for (const auto &it : blocks_derivatives[block])
528 {
529 unsigned int eq, var;
530 expr_t id;
531 tie(eq, var, ignore, id) = it;
532 unsigned int eqr = getBlockEquationID(block, eq);
533 unsigned int varr = getBlockVariableID(block, var);
534 output << " g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = ";
535 id->writeOutput(output, local_output_type, local_temporary_terms, {});
536 output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
537 << "(" << 0
538 << ") " << varr+1
539 << ", equation=" << eqr+1 << endl;
540 }
541 break;
542 default:
543 break;
544 }
545 output << "end" << endl;
546 output.close();
547 }
548 }
549
550 void
writeModelEquationsCode(const string & basename,map_idx_t map_idx) const551 StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx) const
552 {
553
554 ostringstream tmp_output;
555 ofstream code_file;
556 unsigned int instruction_number = 0;
557 bool file_open = false;
558
559 filesystem::create_directories(basename + "/model/bytecode");
560
561 string main_name = basename + "/model/bytecode/static.cod";
562 code_file.open(main_name, ios::out | ios::binary | ios::ate);
563 if (!code_file.is_open())
564 {
565 cerr << R"(Error : Can't open file ")" << main_name << R"(" for writing)" << endl;
566 exit(EXIT_FAILURE);
567 }
568 int count_u;
569 int u_count_int = 0;
570
571 Write_Inf_To_Bin_File(basename + "/model/bytecode/static.bin", u_count_int, file_open, false, symbol_table.endo_nbr());
572 file_open = true;
573
574 //Temporary variables declaration
575 FDIMST_ fdimst(temporary_terms.size());
576 fdimst.write(code_file, instruction_number);
577 FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
578 SOLVE_FORWARD_COMPLETE,
579 0,
580 symbol_table.endo_nbr(),
581 variable_reordered,
582 equation_reordered,
583 false,
584 symbol_table.endo_nbr(),
585 0,
586 0,
587 u_count_int,
588 symbol_table.endo_nbr());
589 fbeginblock.write(code_file, instruction_number);
590
591 // Add a mapping form node ID to temporary terms order
592 int j = 0;
593 for (auto temporary_term : temporary_terms)
594 map_idx[temporary_term->idx] = j++;
595 compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, false, false);
596
597 compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, false, false);
598
599 FENDEQU_ fendequ;
600 fendequ.write(code_file, instruction_number);
601
602 // Get the current code_file position and jump if eval = true
603 streampos pos1 = code_file.tellp();
604 FJMPIFEVAL_ fjmp_if_eval(0);
605 fjmp_if_eval.write(code_file, instruction_number);
606 int prev_instruction_number = instruction_number;
607
608 vector<vector<pair<int, int>>> my_derivatives(symbol_table.endo_nbr());
609 count_u = symbol_table.endo_nbr();
610 for (const auto & [indices, d1] : derivatives[1])
611 {
612 int deriv_id = indices[1];
613 if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
614 {
615 unsigned int eq = indices[0];
616 int symb = getSymbIDByDerivID(deriv_id);
617 unsigned int var = symbol_table.getTypeSpecificID(symb);
618 FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
619 fnumexpr.write(code_file, instruction_number);
620 if (!my_derivatives[eq].size())
621 my_derivatives[eq].clear();
622 my_derivatives[eq].emplace_back(var, count_u);
623
624 d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
625
626 FSTPSU_ fstpsu(count_u);
627 fstpsu.write(code_file, instruction_number);
628 count_u++;
629 }
630 }
631 for (int i = 0; i < symbol_table.endo_nbr(); i++)
632 {
633 FLDR_ fldr(i);
634 fldr.write(code_file, instruction_number);
635 if (my_derivatives[i].size())
636 {
637 for (auto it = my_derivatives[i].begin(); it != my_derivatives[i].end(); ++it)
638 {
639 FLDSU_ fldsu(it->second);
640 fldsu.write(code_file, instruction_number);
641 FLDSV_ fldsv{static_cast<int>(SymbolType::endogenous), static_cast<unsigned int>(it->first)};
642 fldsv.write(code_file, instruction_number);
643 FBINARY_ fbinary{static_cast<int>(BinaryOpcode::times)};
644 fbinary.write(code_file, instruction_number);
645 if (it != my_derivatives[i].begin())
646 {
647 FBINARY_ fbinary{static_cast<int>(BinaryOpcode::plus)};
648 fbinary.write(code_file, instruction_number);
649 }
650 }
651 FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
652 fbinary.write(code_file, instruction_number);
653 }
654 FSTPSU_ fstpsu(i);
655 fstpsu.write(code_file, instruction_number);
656 }
657 // Get the current code_file position and jump = true
658 streampos pos2 = code_file.tellp();
659 FJMP_ fjmp(0);
660 fjmp.write(code_file, instruction_number);
661 // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
662 streampos pos3 = code_file.tellp();
663 code_file.seekp(pos1);
664 FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
665 fjmp_if_eval1.write(code_file, instruction_number);
666 code_file.seekp(pos3);
667 prev_instruction_number = instruction_number;
668
669 temporary_terms_t tt2, tt3;
670
671 // The Jacobian if we have to solve the block determinsitic bloc
672 for (const auto & [indices, d1] : derivatives[1])
673 {
674 int deriv_id = indices[1];
675 if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
676 {
677 unsigned int eq = indices[0];
678 int symb = getSymbIDByDerivID(deriv_id);
679 unsigned int var = symbol_table.getTypeSpecificID(symb);
680 FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
681 fnumexpr.write(code_file, instruction_number);
682 if (!my_derivatives[eq].size())
683 my_derivatives[eq].clear();
684 my_derivatives[eq].emplace_back(var, count_u);
685
686 d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
687 FSTPG2_ fstpg2(eq, var);
688 fstpg2.write(code_file, instruction_number);
689 }
690 }
691
692 // Set codefile position to previous JMP_ and set the number of instructions to jump
693 pos1 = code_file.tellp();
694 code_file.seekp(pos2);
695 FJMP_ fjmp1(instruction_number - prev_instruction_number);
696 fjmp1.write(code_file, instruction_number);
697 code_file.seekp(pos1);
698
699 FENDBLOCK_ fendblock;
700 fendblock.write(code_file, instruction_number);
701 FEND_ fend;
702 fend.write(code_file, instruction_number);
703 code_file.close();
704 }
705
706 void
writeModelEquationsCode_Block(const string & basename,map_idx_t map_idx,vector<map_idx_t> map_idx2) const707 StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map_idx, vector<map_idx_t> map_idx2) const
708 {
709 struct Uff_l
710 {
711 int u, var, lag;
712 Uff_l *pNext;
713 };
714
715 struct Uff
716 {
717 Uff_l *Ufl, *Ufl_First;
718 };
719
720 int i, v;
721 string tmp_s;
722 ostringstream tmp_output;
723 ofstream code_file;
724 unsigned int instruction_number = 0;
725 expr_t lhs = nullptr, rhs = nullptr;
726 BinaryOpNode *eq_node;
727 Uff Uf[symbol_table.endo_nbr()];
728 map<expr_t, int> reference_count;
729 vector<int> feedback_variables;
730 deriv_node_temp_terms_t tef_terms;
731 bool file_open = false;
732
733 filesystem::create_directories(basename + "/model/bytecode");
734
735 string main_name = basename + "/model/bytecode/static.cod";
736 code_file.open(main_name, ios::out | ios::binary | ios::ate);
737 if (!code_file.is_open())
738 {
739 cerr << R"(Error : Can't open file ")" << main_name << R"(" for writing)" << endl;
740 exit(EXIT_FAILURE);
741 }
742 //Temporary variables declaration
743
744 FDIMST_ fdimst(temporary_terms.size());
745 fdimst.write(code_file, instruction_number);
746
747 for (unsigned int block = 0; block < getNbBlocks(); block++)
748 {
749 feedback_variables.clear();
750 if (block > 0)
751 {
752 FENDBLOCK_ fendblock;
753 fendblock.write(code_file, instruction_number);
754 }
755 int count_u;
756 int u_count_int = 0;
757 BlockSimulationType simulation_type = getBlockSimulationType(block);
758 unsigned int block_size = getBlockSize(block);
759 unsigned int block_mfs = getBlockMfs(block);
760 unsigned int block_recursive = block_size - block_mfs;
761
762 if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE
763 || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
764 {
765 Write_Inf_To_Bin_File_Block(basename, block, u_count_int, file_open);
766 file_open = true;
767 }
768
769 FBEGINBLOCK_ fbeginblock(block_mfs,
770 simulation_type,
771 getBlockFirstEquation(block),
772 block_size,
773 variable_reordered,
774 equation_reordered,
775 blocks_linear[block],
776 symbol_table.endo_nbr(),
777 0,
778 0,
779 u_count_int,
780 block_size);
781
782 fbeginblock.write(code_file, instruction_number);
783
784 // Get the current code_file position and jump if eval = true
785 streampos pos1 = code_file.tellp();
786 FJMPIFEVAL_ fjmp_if_eval(0);
787 fjmp_if_eval.write(code_file, instruction_number);
788 int prev_instruction_number = instruction_number;
789
790 for (i = 0; i < static_cast<int>(block_size); i++)
791 {
792 //The Temporary terms
793 temporary_terms_t tt2;
794 if (v_temporary_terms[block].size())
795 {
796 for (auto it : v_temporary_terms[block][i])
797 {
798 if (dynamic_cast<AbstractExternalFunctionNode *>(it))
799 it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
800
801 FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(map_idx.find(it->idx)->second));
802 fnumexpr.write(code_file, instruction_number);
803 it->compile(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
804 FSTPST_ fstpst(static_cast<int>(map_idx.find(it->idx)->second));
805 fstpst.write(code_file, instruction_number);
806 // Insert current node into tt2
807 tt2.insert(it);
808 }
809 }
810
811 // The equations
812 int variable_ID, equation_ID;
813 EquationType equ_type;
814 switch (simulation_type)
815 {
816 evaluation:
817 case EVALUATE_BACKWARD:
818 case EVALUATE_FORWARD:
819 equ_type = getBlockEquationType(block, i);
820 {
821 FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
822 fnumexpr.write(code_file, instruction_number);
823 }
824 if (equ_type == E_EVALUATE)
825 {
826 eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
827 lhs = eq_node->arg1;
828 rhs = eq_node->arg2;
829 rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
830 lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
831 }
832 else if (equ_type == E_EVALUATE_S)
833 {
834 eq_node = static_cast<BinaryOpNode *>(getBlockEquationRenormalizedExpr(block, i));
835 lhs = eq_node->arg1;
836 rhs = eq_node->arg2;
837 rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
838 lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
839 }
840 break;
841 case SOLVE_BACKWARD_COMPLETE:
842 case SOLVE_FORWARD_COMPLETE:
843 if (i < static_cast<int>(block_recursive))
844 goto evaluation;
845 variable_ID = getBlockVariableID(block, i);
846 equation_ID = getBlockEquationID(block, i);
847 feedback_variables.push_back(variable_ID);
848 Uf[equation_ID].Ufl = nullptr;
849 goto end;
850 default:
851 end:
852 FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
853 fnumexpr.write(code_file, instruction_number);
854 eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
855 lhs = eq_node->arg1;
856 rhs = eq_node->arg2;
857 lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
858 rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
859
860 FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
861 fbinary.write(code_file, instruction_number);
862
863 FSTPR_ fstpr(i - block_recursive);
864 fstpr.write(code_file, instruction_number);
865 }
866 }
867 FENDEQU_ fendequ;
868 fendequ.write(code_file, instruction_number);
869
870 // The Jacobian if we have to solve the block
871 if (simulation_type != EVALUATE_BACKWARD
872 && simulation_type != EVALUATE_FORWARD)
873 {
874 switch (simulation_type)
875 {
876 case SOLVE_BACKWARD_SIMPLE:
877 case SOLVE_FORWARD_SIMPLE:
878 {
879 FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
880 fnumexpr.write(code_file, instruction_number);
881 }
882 compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx, temporary_terms);
883 {
884 FSTPG_ fstpg(0);
885 fstpg.write(code_file, instruction_number);
886 }
887 break;
888
889 case SOLVE_BACKWARD_COMPLETE:
890 case SOLVE_FORWARD_COMPLETE:
891 count_u = feedback_variables.size();
892 for (const auto &it : blocks_derivatives[block])
893 {
894 unsigned int eq, var;
895 tie(eq, var, ignore, ignore) = it;
896 unsigned int eqr = getBlockEquationID(block, eq);
897 unsigned int varr = getBlockVariableID(block, var);
898 if (eq >= block_recursive && var >= block_recursive)
899 {
900 if (!Uf[eqr].Ufl)
901 {
902 Uf[eqr].Ufl = static_cast<Uff_l *>(malloc(sizeof(Uff_l)));
903 Uf[eqr].Ufl_First = Uf[eqr].Ufl;
904 }
905 else
906 {
907 Uf[eqr].Ufl->pNext = static_cast<Uff_l *>(malloc(sizeof(Uff_l)));
908 Uf[eqr].Ufl = Uf[eqr].Ufl->pNext;
909 }
910 Uf[eqr].Ufl->pNext = nullptr;
911 Uf[eqr].Ufl->u = count_u;
912 Uf[eqr].Ufl->var = varr;
913 FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
914 fnumexpr.write(code_file, instruction_number);
915 compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx, temporary_terms);
916 FSTPSU_ fstpsu(count_u);
917 fstpsu.write(code_file, instruction_number);
918 count_u++;
919 }
920 }
921 for (i = 0; i < static_cast<int>(block_size); i++)
922 {
923 if (i >= static_cast<int>(block_recursive))
924 {
925 FLDR_ fldr(i-block_recursive);
926 fldr.write(code_file, instruction_number);
927
928 FLDZ_ fldz;
929 fldz.write(code_file, instruction_number);
930
931 v = getBlockEquationID(block, i);
932 for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
933 {
934 FLDSU_ fldsu(Uf[v].Ufl->u);
935 fldsu.write(code_file, instruction_number);
936 FLDSV_ fldsv{static_cast<int>(SymbolType::endogenous), static_cast<unsigned int>(Uf[v].Ufl->var)};
937 fldsv.write(code_file, instruction_number);
938
939 FBINARY_ fbinary{static_cast<int>(BinaryOpcode::times)};
940 fbinary.write(code_file, instruction_number);
941
942 FCUML_ fcuml;
943 fcuml.write(code_file, instruction_number);
944 }
945 Uf[v].Ufl = Uf[v].Ufl_First;
946 while (Uf[v].Ufl)
947 {
948 Uf[v].Ufl_First = Uf[v].Ufl->pNext;
949 free(Uf[v].Ufl);
950 Uf[v].Ufl = Uf[v].Ufl_First;
951 }
952 FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
953 fbinary.write(code_file, instruction_number);
954
955 FSTPSU_ fstpsu(i - block_recursive);
956 fstpsu.write(code_file, instruction_number);
957
958 }
959 }
960 break;
961 default:
962 break;
963 }
964 }
965
966 // Get the current code_file position and jump = true
967 streampos pos2 = code_file.tellp();
968 FJMP_ fjmp(0);
969 fjmp.write(code_file, instruction_number);
970 // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
971 streampos pos3 = code_file.tellp();
972 code_file.seekp(pos1);
973 FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
974 fjmp_if_eval1.write(code_file, instruction_number);
975 code_file.seekp(pos3);
976 prev_instruction_number = instruction_number;
977
978 temporary_terms_t tt2, tt3;
979 deriv_node_temp_terms_t tef_terms2;
980
981 for (i = 0; i < static_cast<int>(block_size); i++)
982 {
983 if (v_temporary_terms_local[block].size())
984 {
985 for (const auto &it : v_temporary_terms_local[block][i])
986 {
987 if (dynamic_cast<AbstractExternalFunctionNode *>(it))
988 it->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
989
990 FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(map_idx2[block].find(it->idx)->second));
991 fnumexpr.write(code_file, instruction_number);
992
993 it->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);
994
995 FSTPST_ fstpst(static_cast<int>(map_idx2[block].find(it->idx)->second));
996 fstpst.write(code_file, instruction_number);
997 // Insert current node into tt2
998 tt3.insert(it);
999 tt2.insert(it);
1000 }
1001 }
1002
1003 // The equations
1004 int variable_ID, equation_ID;
1005 EquationType equ_type;
1006 switch (simulation_type)
1007 {
1008 evaluation_l:
1009 case EVALUATE_BACKWARD:
1010 case EVALUATE_FORWARD:
1011 equ_type = getBlockEquationType(block, i);
1012 {
1013 FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
1014 fnumexpr.write(code_file, instruction_number);
1015 }
1016 if (equ_type == E_EVALUATE)
1017 {
1018 eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
1019 lhs = eq_node->arg1;
1020 rhs = eq_node->arg2;
1021 rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
1022 lhs->compile(code_file, instruction_number, true, tt2, map_idx2[block], false, false);
1023 }
1024 else if (equ_type == E_EVALUATE_S)
1025 {
1026 eq_node = static_cast<BinaryOpNode *>(getBlockEquationRenormalizedExpr(block, i));
1027 lhs = eq_node->arg1;
1028 rhs = eq_node->arg2;
1029 rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
1030 lhs->compile(code_file, instruction_number, true, tt2, map_idx2[block], false, false);
1031 }
1032 break;
1033 case SOLVE_BACKWARD_COMPLETE:
1034 case SOLVE_FORWARD_COMPLETE:
1035 if (i < static_cast<int>(block_recursive))
1036 goto evaluation_l;
1037 variable_ID = getBlockVariableID(block, i);
1038 equation_ID = getBlockEquationID(block, i);
1039 feedback_variables.push_back(variable_ID);
1040 Uf[equation_ID].Ufl = nullptr;
1041 goto end_l;
1042 default:
1043 end_l:
1044 FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
1045 fnumexpr.write(code_file, instruction_number);
1046 eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
1047 lhs = eq_node->arg1;
1048 rhs = eq_node->arg2;
1049 lhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
1050 rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
1051
1052 FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
1053 fbinary.write(code_file, instruction_number);
1054
1055 FSTPR_ fstpr(i - block_recursive);
1056 fstpr.write(code_file, instruction_number);
1057 }
1058 }
1059 FENDEQU_ fendequ_l;
1060 fendequ_l.write(code_file, instruction_number);
1061
1062 // The Jacobian if we have to solve the block determinsitic bloc
1063 switch (simulation_type)
1064 {
1065 case SOLVE_BACKWARD_SIMPLE:
1066 case SOLVE_FORWARD_SIMPLE:
1067 {
1068 FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
1069 fnumexpr.write(code_file, instruction_number);
1070 }
1071 compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx2[block], tt2 /*temporary_terms*/);
1072 {
1073 FSTPG2_ fstpg2(0, 0);
1074 fstpg2.write(code_file, instruction_number);
1075 }
1076 break;
1077 case EVALUATE_BACKWARD:
1078 case EVALUATE_FORWARD:
1079 case SOLVE_BACKWARD_COMPLETE:
1080 case SOLVE_FORWARD_COMPLETE:
1081 count_u = feedback_variables.size();
1082 for (const auto &it : blocks_derivatives[block])
1083 {
1084 unsigned int eq, var;
1085 tie(eq, var, ignore, ignore) = it;
1086 unsigned int eqr = getBlockEquationID(block, eq);
1087 unsigned int varr = getBlockVariableID(block, var);
1088 FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0);
1089 fnumexpr.write(code_file, instruction_number);
1090
1091 compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx2[block], tt2 /*temporary_terms*/);
1092
1093 FSTPG2_ fstpg2(eq, var);
1094 fstpg2.write(code_file, instruction_number);
1095 }
1096 break;
1097 default:
1098 break;
1099 }
1100 // Set codefile position to previous JMP_ and set the number of instructions to jump
1101 pos1 = code_file.tellp();
1102 code_file.seekp(pos2);
1103 FJMP_ fjmp1(instruction_number - prev_instruction_number);
1104 fjmp1.write(code_file, instruction_number);
1105 code_file.seekp(pos1);
1106 }
1107 FENDBLOCK_ fendblock;
1108 fendblock.write(code_file, instruction_number);
1109 FEND_ fend;
1110 fend.write(code_file, instruction_number);
1111 code_file.close();
1112 }
1113
1114 void
Write_Inf_To_Bin_File_Block(const string & basename,int num,int & u_count_int,bool & file_open) const1115 StaticModel::Write_Inf_To_Bin_File_Block(const string &basename, int num,
1116 int &u_count_int, bool &file_open) const
1117 {
1118 int j;
1119 std::ofstream SaveCode;
1120 string filename = basename + "/model/bytecode/static.bin";
1121 if (file_open)
1122 SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate);
1123 else
1124 SaveCode.open(filename, ios::out | ios::binary);
1125 if (!SaveCode.is_open())
1126 {
1127 cerr << "Error : Can't open file " << filename << " for writing" << endl;
1128 exit(EXIT_FAILURE);
1129 }
1130 u_count_int = 0;
1131 unsigned int block_size = getBlockSize(num);
1132 unsigned int block_mfs = getBlockMfs(num);
1133 unsigned int block_recursive = block_size - block_mfs;
1134 for (const auto &it : blocks_derivatives[num])
1135 {
1136 unsigned int eq, var;
1137 tie(eq, var, ignore, ignore) = it;
1138 int lag = 0;
1139 if (eq >= block_recursive && var >= block_recursive)
1140 {
1141 int v = eq - block_recursive;
1142 SaveCode.write(reinterpret_cast<char *>(&v), sizeof(v));
1143 int varr = var - block_recursive;
1144 SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
1145 SaveCode.write(reinterpret_cast<char *>(&lag), sizeof(lag));
1146 int u = u_count_int + block_mfs;
1147 SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
1148 u_count_int++;
1149 }
1150 }
1151
1152 for (j = block_recursive; j < static_cast<int>(block_size); j++)
1153 {
1154 unsigned int varr = getBlockVariableID(num, j);
1155 SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
1156 }
1157 for (j = block_recursive; j < static_cast<int>(block_size); j++)
1158 {
1159 unsigned int eqr = getBlockEquationID(num, j);
1160 SaveCode.write(reinterpret_cast<char *>(&eqr), sizeof(eqr));
1161 }
1162 SaveCode.close();
1163 }
1164
1165 map<tuple<int, int, int>, expr_t>
collect_first_order_derivatives_endogenous()1166 StaticModel::collect_first_order_derivatives_endogenous()
1167 {
1168 map<tuple<int, int, int>, expr_t> endo_derivatives;
1169 for (auto & [indices, d1] : derivatives[1])
1170 {
1171 if (getTypeByDerivID(indices[1]) == SymbolType::endogenous)
1172 {
1173 int eq = indices[0];
1174 int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(indices[1]));
1175 int lag = 0;
1176 endo_derivatives[{ eq, var, lag }] = d1;
1177 }
1178 }
1179 return endo_derivatives;
1180 }
1181
1182 void
computingPass(int derivsOrder,int paramsDerivsOrder,const eval_context_t & eval_context,bool no_tmp_terms,bool block,bool bytecode)1183 StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool bytecode)
1184 {
1185 initializeVariablesAndEquations();
1186
1187 vector<BinaryOpNode *> neweqs;
1188 for (unsigned int eq = 0; eq < equations.size() - aux_equations.size(); eq++)
1189 {
1190 expr_t eq_tmp = equations[eq]->substituteStaticAuxiliaryVariable();
1191 neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
1192 }
1193
1194 for (auto &aux_equation : aux_equations)
1195 {
1196 expr_t eq_tmp = aux_equation->substituteStaticAuxiliaryDefinition();
1197 neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
1198 }
1199
1200 equations.clear();
1201 copy(neweqs.begin(), neweqs.end(), back_inserter(equations));
1202
1203 // Compute derivatives w.r. to all endogenous
1204 set<int> vars;
1205 for (int i = 0; i < symbol_table.endo_nbr(); i++)
1206 {
1207 int id = symbol_table.getID(SymbolType::endogenous, i);
1208 // if (!symbol_table.isAuxiliaryVariableButNotMultiplier(id))
1209 vars.insert(getDerivID(id, 0));
1210 }
1211
1212 // Launch computations
1213 cout << "Computing static model derivatives (order " << derivsOrder << ")." << endl;
1214
1215 computeDerivatives(derivsOrder, vars);
1216
1217 if (paramsDerivsOrder > 0)
1218 {
1219 cout << "Computing static model derivatives w.r.t. parameters (order " << paramsDerivsOrder << ")." << endl;
1220 computeParamsDerivatives(paramsDerivsOrder);
1221 }
1222
1223 if (block)
1224 {
1225 jacob_map_t contemporaneous_jacobian, static_jacobian;
1226 vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
1227
1228 // for each block contains pair<Size, Feddback_variable>
1229 vector<pair<int, int>> blocks;
1230
1231 evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
1232
1233 computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
1234
1235 computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
1236
1237 map<tuple<int, int, int>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
1238
1239 equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
1240
1241 cout << "Finding the optimal block decomposition of the model ..." << endl;
1242
1243 lag_lead_vector_t equation_lag_lead, variable_lag_lead;
1244
1245 computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, false, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed);
1246
1247 block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, false);
1248
1249 printBlockDecomposition(blocks);
1250
1251 computeChainRuleJacobian(blocks_derivatives);
1252
1253 blocks_linear = BlockLinear(blocks_derivatives, variable_reordered);
1254
1255 collect_block_first_order_derivatives();
1256
1257 global_temporary_terms = true;
1258 if (!no_tmp_terms)
1259 computeTemporaryTermsOrdered();
1260 }
1261 else
1262 {
1263 computeTemporaryTerms(true, no_tmp_terms);
1264 if (bytecode && !no_tmp_terms)
1265 computeTemporaryTermsMapping(temporary_terms, map_idx);
1266
1267 /* Must be called after computeTemporaryTerms(), because it depends on
1268 temporary_terms_mlv to be filled */
1269 if (paramsDerivsOrder > 0 && !no_tmp_terms)
1270 computeParamsDerivativesTemporaryTerms();
1271 }
1272 }
1273
1274 void
writeStaticMFile(const string & basename) const1275 StaticModel::writeStaticMFile(const string &basename) const
1276 {
1277 writeStaticModel(basename, false, false);
1278 }
1279
1280 void
writeWrapperFunctions(const string & basename,const string & ending) const1281 StaticModel::writeWrapperFunctions(const string &basename, const string &ending) const
1282 {
1283 string name;
1284 if (ending == "g1")
1285 name = "static_resid_g1";
1286 else if (ending == "g2")
1287 name = "static_resid_g1_g2";
1288 else if (ending == "g3")
1289 name = "static_resid_g1_g2_g3";
1290
1291 string filename = packageDir(basename) + "/" + name + ".m";
1292 ofstream output;
1293 output.open(filename, ios::out | ios::binary);
1294 if (!output.is_open())
1295 {
1296 cerr << "Error: Can't open file " << filename << " for writing" << endl;
1297 exit(EXIT_FAILURE);
1298 }
1299
1300 if (ending == "g1")
1301 output << "function [residual, g1] = " << name << "(T, y, x, params, T_flag)" << endl
1302 << "% function [residual, g1] = " << name << "(T, y, x, params, T_flag)" << endl;
1303 else if (ending == "g2")
1304 output << "function [residual, g1, g2] = " << name << "(T, y, x, params, T_flag)" << endl
1305 << "% function [residual, g1, g2] = " << name << "(T, y, x, params, T_flag)" << endl;
1306 else if (ending == "g3")
1307 output << "function [residual, g1, g2, g3] = " << name << "(T, y, x, params, T_flag)" << endl
1308 << "% function [residual, g1, g2, g3] = " << name << "(T, y, x, params, T_flag)" << endl;
1309
1310 output << "%" << endl
1311 << "% Wrapper function automatically created by Dynare" << endl
1312 << "%" << endl
1313 << endl
1314 << " if T_flag" << endl
1315 << " T = " << basename << ".static_" << ending << "_tt(T, y, x, params);" << endl
1316 << " end" << endl;
1317
1318 if (ending == "g1")
1319 output << " residual = " << basename << ".static_resid(T, y, x, params, false);" << endl
1320 << " g1 = " << basename << ".static_g1(T, y, x, params, false);" << endl;
1321 else if (ending == "g2")
1322 output << " [residual, g1] = " << basename << ".static_resid_g1(T, y, x, params, false);" << endl
1323 << " g2 = " << basename << ".static_g2(T, y, x, params, false);" << endl;
1324 else if (ending == "g3")
1325 output << " [residual, g1, g2] = " << basename << ".static_resid_g1_g2(T, y, x, params, false);" << endl
1326 << " g3 = " << basename << ".static_g3(T, y, x, params, false);" << endl;
1327
1328 output << endl << "end" << endl;
1329 output.close();
1330 }
1331
1332 void
writeStaticModelHelper(const string & basename,const string & name,const string & retvalname,const string & name_tt,size_t ttlen,const string & previous_tt_name,const ostringstream & init_s,const ostringstream & end_s,const ostringstream & s,const ostringstream & s_tt) const1333 StaticModel::writeStaticModelHelper(const string &basename,
1334 const string &name, const string &retvalname,
1335 const string &name_tt, size_t ttlen,
1336 const string &previous_tt_name,
1337 const ostringstream &init_s, const ostringstream &end_s,
1338 const ostringstream &s, const ostringstream &s_tt) const
1339 {
1340 string filename = packageDir(basename) + "/" + name_tt + ".m";
1341 ofstream output;
1342 output.open(filename, ios::out | ios::binary);
1343 if (!output.is_open())
1344 {
1345 cerr << "Error: Can't open file " << filename << " for writing" << endl;
1346 exit(EXIT_FAILURE);
1347 }
1348
1349 output << "function T = " << name_tt << "(T, y, x, params)" << endl
1350 << "% function T = " << name_tt << "(T, y, x, params)" << endl
1351 << "%" << endl
1352 << "% File created by Dynare Preprocessor from .mod file" << endl
1353 << "%" << endl
1354 << "% Inputs:" << endl
1355 << "% T [#temp variables by 1] double vector of temporary terms to be filled by function" << endl
1356 << "% y [M_.endo_nbr by 1] double vector of endogenous variables in declaration order" << endl
1357 << "% x [M_.exo_nbr by 1] double vector of exogenous variables in declaration order" << endl
1358 << "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
1359 << "%" << endl
1360 << "% Output:" << endl
1361 << "% T [#temp variables by 1] double vector of temporary terms" << endl
1362 << "%" << endl << endl
1363 << "assert(length(T) >= " << ttlen << ");" << endl
1364 << endl;
1365
1366 if (!previous_tt_name.empty())
1367 output << "T = " << basename << "." << previous_tt_name << "(T, y, x, params);" << endl << endl;
1368
1369 output << s_tt.str() << endl
1370 << "end" << endl;
1371 output.close();
1372
1373 filename = packageDir(basename) + "/" + name + ".m";
1374 output.open(filename, ios::out | ios::binary);
1375 if (!output.is_open())
1376 {
1377 cerr << "Error: Can't open file " << filename << " for writing" << endl;
1378 exit(EXIT_FAILURE);
1379 }
1380
1381 output << "function " << retvalname << " = " << name << "(T, y, x, params, T_flag)" << endl
1382 << "% function " << retvalname << " = " << name << "(T, y, x, params, T_flag)" << endl
1383 << "%" << endl
1384 << "% File created by Dynare Preprocessor from .mod file" << endl
1385 << "%" << endl
1386 << "% Inputs:" << endl
1387 << "% T [#temp variables by 1] double vector of temporary terms to be filled by function" << endl
1388 << "% y [M_.endo_nbr by 1] double vector of endogenous variables in declaration order" << endl
1389 << "% x [M_.exo_nbr by 1] double vector of exogenous variables in declaration order" << endl
1390 << "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
1391 << "% to evaluate the model" << endl
1392 << "% T_flag boolean boolean flag saying whether or not to calculate temporary terms" << endl
1393 << "%" << endl
1394 << "% Output:" << endl
1395 << "% " << retvalname << endl
1396 << "%" << endl << endl;
1397
1398 if (!name_tt.empty())
1399 output << "if T_flag" << endl
1400 << " T = " << basename << "." << name_tt << "(T, y, x, params);" << endl
1401 << "end" << endl;
1402
1403 output << init_s.str() << endl
1404 << s.str()
1405 << end_s.str() << endl
1406 << "end" << endl;
1407 output.close();
1408 }
1409
1410 void
writeStaticMatlabCompatLayer(const string & basename) const1411 StaticModel::writeStaticMatlabCompatLayer(const string &basename) const
1412 {
1413 string filename = packageDir(basename) + "/static.m";
1414 ofstream output;
1415 output.open(filename, ios::out | ios::binary);
1416 if (!output.is_open())
1417 {
1418 cerr << "Error: Can't open file " << filename << " for writing" << endl;
1419 exit(EXIT_FAILURE);
1420 }
1421 int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size();
1422
1423 output << "function [residual, g1, g2, g3] = static(y, x, params)" << endl
1424 << " T = NaN(" << ntt << ", 1);" << endl
1425 << " if nargout <= 1" << endl
1426 << " residual = " << basename << ".static_resid(T, y, x, params, true);" << endl
1427 << " elseif nargout == 2" << endl
1428 << " [residual, g1] = " << basename << ".static_resid_g1(T, y, x, params, true);" << endl
1429 << " elseif nargout == 3" << endl
1430 << " [residual, g1, g2] = " << basename << ".static_resid_g1_g2(T, y, x, params, true);" << endl
1431 << " else" << endl
1432 << " [residual, g1, g2, g3] = " << basename << ".static_resid_g1_g2_g3(T, y, x, params, true);" << endl
1433 << " end" << endl
1434 << "end" << endl;
1435
1436 output.close();
1437 }
1438
1439 void
writeStaticModel(ostream & StaticOutput,bool use_dll,bool julia) const1440 StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const
1441 {
1442 writeStaticModel("", StaticOutput, use_dll, julia);
1443 }
1444
1445 void
writeStaticModel(const string & basename,bool use_dll,bool julia) const1446 StaticModel::writeStaticModel(const string &basename, bool use_dll, bool julia) const
1447 {
1448 ofstream StaticOutput;
1449 writeStaticModel(basename, StaticOutput, use_dll, julia);
1450 }
1451
1452 void
writeStaticModel(const string & basename,ostream & StaticOutput,bool use_dll,bool julia) const1453 StaticModel::writeStaticModel(const string &basename,
1454 ostream &StaticOutput, bool use_dll, bool julia) const
1455 {
1456 vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
1457 vector<ostringstream> tt_output(derivatives.size()); // Temp terms output (at all orders)
1458
1459 ExprNodeOutputType output_type = (use_dll ? ExprNodeOutputType::CStaticModel :
1460 julia ? ExprNodeOutputType::juliaStaticModel : ExprNodeOutputType::matlabStaticModel);
1461
1462 deriv_node_temp_terms_t tef_terms;
1463 temporary_terms_t temp_term_union;
1464
1465 writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_idxs,
1466 tt_output[0], output_type, tef_terms);
1467
1468 writeTemporaryTerms(temporary_terms_derivatives[0],
1469 temp_term_union,
1470 temporary_terms_idxs,
1471 tt_output[0], output_type, tef_terms);
1472
1473 writeModelEquations(d_output[0], output_type, temp_term_union);
1474
1475 int nrows = equations.size();
1476 int JacobianColsNbr = symbol_table.endo_nbr();
1477 int hessianColsNbr = JacobianColsNbr*JacobianColsNbr;
1478
1479 auto getJacobCol = [this](int var) { return symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)); };
1480
1481 // Write Jacobian w.r. to endogenous only
1482 if (!derivatives[1].empty())
1483 {
1484 writeTemporaryTerms(temporary_terms_derivatives[1],
1485 temp_term_union,
1486 temporary_terms_idxs,
1487 tt_output[1], output_type, tef_terms);
1488
1489 for (const auto & [indices, d1] : derivatives[1])
1490 {
1491 auto [eq, var] = vectorToTuple<2>(indices);
1492
1493 jacobianHelper(d_output[1], eq, getJacobCol(var), output_type);
1494 d_output[1] << "=";
1495 d1->writeOutput(d_output[1], output_type,
1496 temp_term_union, temporary_terms_idxs, tef_terms);
1497 d_output[1] << ";" << endl;
1498 }
1499 }
1500
1501 // Write derivatives for order ≥ 2
1502 for (size_t i = 2; i < derivatives.size(); i++)
1503 if (!derivatives[i].empty())
1504 {
1505 writeTemporaryTerms(temporary_terms_derivatives[i],
1506 temp_term_union,
1507 temporary_terms_idxs,
1508 tt_output[i], output_type, tef_terms);
1509
1510 /* When creating the sparse matrix (in MATLAB or C mode), since storage
1511 is in column-major order, output the first column, then the second,
1512 then the third. This gives a significant performance boost in use_dll
1513 mode (at both compilation and runtime), because it facilitates memory
1514 accesses and expression reusage. */
1515 ostringstream col0_output, col1_output, col2_output;
1516
1517 int k = 0; // Current line index in the 3-column matrix
1518 for (const auto &[vidx, d] : derivatives[i])
1519 {
1520 int eq = vidx[0];
1521
1522 int col_idx = 0;
1523 for (size_t j = 1; j < vidx.size(); j++)
1524 {
1525 col_idx *= JacobianColsNbr;
1526 col_idx += getJacobCol(vidx[j]);
1527 }
1528
1529 if (output_type == ExprNodeOutputType::juliaStaticModel)
1530 {
1531 d_output[i] << " @inbounds " << "g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = ";
1532 d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms);
1533 d_output[i] << endl;
1534 }
1535 else
1536 {
1537 sparseHelper(i, col0_output, k, 0, output_type);
1538 col0_output << "=" << eq + 1 << ";" << endl;
1539
1540 sparseHelper(i, col1_output, k, 1, output_type);
1541 col1_output << "=" << col_idx + 1 << ";" << endl;
1542
1543 sparseHelper(i, col2_output, k, 2, output_type);
1544 col2_output << "=";
1545 d->writeOutput(col2_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms);
1546 col2_output << ";" << endl;
1547
1548 k++;
1549 }
1550
1551 // Output symetric elements at order 2
1552 if (i == 2 && vidx[1] != vidx[2])
1553 {
1554 int col_idx_sym = getJacobCol(vidx[2]) * JacobianColsNbr + getJacobCol(vidx[1]);
1555
1556 if (output_type == ExprNodeOutputType::juliaStaticModel)
1557 d_output[2] << " @inbounds g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = "
1558 << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
1559 else
1560 {
1561 sparseHelper(2, col0_output, k, 0, output_type);
1562 col0_output << "=" << eq + 1 << ";" << endl;
1563
1564 sparseHelper(2, col1_output, k, 1, output_type);
1565 col1_output << "=" << col_idx_sym + 1 << ";" << endl;
1566
1567 sparseHelper(2, col2_output, k, 2, output_type);
1568 col2_output << "=";
1569 sparseHelper(2, col2_output, k-1, 2, output_type);
1570 col2_output << ";" << endl;
1571
1572 k++;
1573 }
1574 }
1575 }
1576 if (output_type != ExprNodeOutputType::juliaStaticModel)
1577 d_output[i] << col0_output.str() << col1_output.str() << col2_output.str();
1578 }
1579
1580 if (output_type == ExprNodeOutputType::matlabStaticModel)
1581 {
1582 // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201
1583 map<string, string> tmp_paren_vars;
1584 bool message_printed = false;
1585 for (auto &it : tt_output)
1586 fixNestedParenthesis(it, tmp_paren_vars, message_printed);
1587 for (auto &it : d_output)
1588 fixNestedParenthesis(it, tmp_paren_vars, message_printed);
1589
1590 ostringstream init_output, end_output;
1591 init_output << "residual = zeros(" << equations.size() << ", 1);";
1592 end_output << "if ~isreal(residual)" << endl
1593 << " residual = real(residual)+imag(residual).^2;" << endl
1594 << "end";
1595 writeStaticModelHelper(basename, "static_resid", "residual", "static_resid_tt",
1596 temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(),
1597 "", init_output, end_output,
1598 d_output[0], tt_output[0]);
1599
1600 init_output.str("");
1601 end_output.str("");
1602 init_output << "g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");";
1603 end_output << "if ~isreal(g1)" << endl
1604 << " g1 = real(g1)+2*imag(g1);" << endl
1605 << "end";
1606 writeStaticModelHelper(basename, "static_g1", "g1", "static_g1_tt",
1607 temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(),
1608 "static_resid_tt",
1609 init_output, end_output,
1610 d_output[1], tt_output[1]);
1611 writeWrapperFunctions(basename, "g1");
1612
1613 // For order ≥ 2
1614 int ncols = JacobianColsNbr;
1615 int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size();
1616 for (size_t i = 2; i < derivatives.size(); i++)
1617 {
1618 ncols *= JacobianColsNbr;
1619 ntt += temporary_terms_derivatives[i].size();
1620 string gname = "g" + to_string(i);
1621 string vname = "v" + to_string(i);
1622 string gprevname = "g" + to_string(i-1);
1623
1624 init_output.str("");
1625 end_output.str("");
1626 if (derivatives[i].size())
1627 {
1628 init_output << vname << " = zeros(" << NNZDerivatives[i] << ",3);";
1629 end_output << gname << " = sparse("
1630 << vname << "(:,1),"
1631 << vname << "(:,2),"
1632 << vname << "(:,3),"
1633 << nrows << "," << ncols << ");";
1634 }
1635 else
1636 init_output << gname << " = sparse([],[],[]," << nrows << "," << ncols << ");";
1637 writeStaticModelHelper(basename, "static_" + gname, gname,
1638 "static_" + gname + "_tt",
1639 ntt,
1640 "static_" + gprevname + "_tt",
1641 init_output, end_output,
1642 d_output[i], tt_output[i]);
1643 if (i <= 3)
1644 writeWrapperFunctions(basename, gname);
1645 }
1646
1647 writeStaticMatlabCompatLayer(basename);
1648 }
1649 else if (output_type == ExprNodeOutputType::CStaticModel)
1650 {
1651 for (size_t i = 0; i < d_output.size(); i++)
1652 {
1653 string funcname = i == 0 ? "resid" : "g" + to_string(i);
1654 string argname = i == 0 ? "residual" : i == 1 ? "g1" : "v" + to_string(i);
1655 StaticOutput << "void static_" << funcname << "_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T)" << endl
1656 << "{" << endl
1657 << tt_output[i].str()
1658 << "}" << endl
1659 << endl
1660 << "void static_" << funcname << "(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *" << argname << ")" << endl
1661 << "{" << endl;
1662 if (i == 0)
1663 StaticOutput << " double lhs, rhs;" << endl;
1664 StaticOutput << d_output[i].str()
1665 << "}" << endl
1666 << endl;
1667 }
1668 }
1669 else
1670 {
1671 string filename = basename + "Static.jl";
1672 ofstream output;
1673 output.open(filename, ios::out | ios::binary);
1674 if (!output.is_open())
1675 {
1676 cerr << "Error: Can't open file " << filename << " for writing" << endl;
1677 exit(EXIT_FAILURE);
1678 }
1679
1680 output << "module " << basename << "Static" << endl
1681 << "#" << endl
1682 << "# NB: this file was automatically generated by Dynare" << endl
1683 << "# from " << basename << ".mod" << endl
1684 << "#" << endl
1685 << "using Utils" << endl << endl
1686 << "export tmp_nbr, static!, staticResid!, staticG1!, staticG2!, staticG3!" << endl << endl
1687 << "#=" << endl
1688 << "# The comments below apply to all functions contained in this module #" << endl
1689 << " NB: The arguments contained on the first line of the function" << endl
1690 << " definition are those that are modified in place" << endl << endl
1691 << "## Exported Functions ##" << endl
1692 << " static! : Wrapper function; computes residuals, Jacobian, Hessian," << endl
1693 << " and third order derivatives matroces depending on the arguments provided" << endl
1694 << " staticResid! : Computes the static model residuals" << endl
1695 << " staticG1! : Computes the static model Jacobian" << endl
1696 << " staticG2! : Computes the static model Hessian" << endl
1697 << " staticG3! : Computes the static model third derivatives" << endl << endl
1698 << "## Exported Variables ##" << endl
1699 << " tmp_nbr : Vector{Int}(4) respectively the number of temporary variables" << endl
1700 << " for the residuals, g1, g2 and g3." << endl << endl
1701 << "## Local Functions ##" << endl
1702 << " staticResidTT! : Computes the static model temporary terms for the residuals" << endl
1703 << " staticG1TT! : Computes the static model temporary terms for the Jacobian" << endl
1704 << " staticG2TT! : Computes the static model temporary terms for the Hessian" << endl
1705 << " staticG3TT! : Computes the static model temporary terms for the third derivatives" << endl << endl
1706 << "## Function Arguments ##" << endl
1707 << " T : Vector{Float64}(num_temp_terms) temporary terms" << endl
1708 << " y : Vector{Float64}(model_.endo_nbr) endogenous variables in declaration order" << endl
1709 << " x : Vector{Float64}(model_.exo_nbr) exogenous variables in declaration order" << endl
1710 << " params : Vector{Float64}(model_.param) parameter values in declaration order" << endl
1711 << " residual : Vector{Float64}(model_.eq_nbr) residuals of the static model equations" << endl
1712 << " in order of declaration of the equations. Dynare may prepend auxiliary equations," << endl
1713 << " see model.aux_vars" << endl
1714 << " g1 : Matrix{Float64}(model.eq_nbr,model_.endo_nbr) Jacobian matrix of the static model equations" << endl
1715 << " The columns and rows respectively correspond to the variables in declaration order and the" << endl
1716 << " equations in order of declaration" << endl
1717 << " g2 : spzeros(model.eq_nbr, model_.endo^2) Hessian matrix of the static model equations" << endl
1718 << " The columns and rows respectively correspond to the variables in declaration order and the" << endl
1719 << " equations in order of declaration" << endl
1720 << " g3 : spzeros(model.eq_nbr, model_.endo^3) Third order derivatives matrix of the static model equations" << endl
1721 << " The columns and rows respectively correspond to the variables in declaration order and the" << endl
1722 << " equations in order of declaration" << endl << endl
1723 << "## Remarks ##" << endl
1724 << " [1] The size of `T`, ie the value of `num_temp_terms`, depends on the version of the static model called. The number of temporary variables" << endl
1725 << " used for the different returned objects (residuals, jacobian, hessian or third order derivatives) is given by the elements in `tmp_nbr`" << endl
1726 << " exported vector. The first element is the number of temporaries used for the computation of the residuals, the second element is the" << endl
1727 << " number of temporaries used for the evaluation of the jacobian matrix, etc. If one calls the version of the static model computing the" << endl
1728 << " residuals, and the jacobian and hessian matrices, then `T` must have at least `sum(tmp_nbr[1:3])` elements." << endl
1729 << "=#" << endl << endl;
1730
1731 // Write the number of temporary terms
1732 output << "tmp_nbr = zeros(Int,4)" << endl
1733 << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << "# Number of temporary terms for the residuals" << endl
1734 << "tmp_nbr[2] = " << temporary_terms_derivatives[1].size() << "# Number of temporary terms for g1 (jacobian)" << endl
1735 << "tmp_nbr[3] = " << temporary_terms_derivatives[2].size() << "# Number of temporary terms for g2 (hessian)" << endl
1736 << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl;
1737
1738 // staticResidTT!
1739 output << "function staticResidTT!(T::Vector{Float64}," << endl
1740 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
1741 << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl
1742 << tt_output[0].str()
1743 << " return nothing" << endl
1744 << "end" << endl << endl;
1745
1746 // static!
1747 output << "function staticResid!(T::Vector{Float64}, residual::Vector{Float64}," << endl
1748 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T0_flag::Bool)" << endl
1749 << " @assert length(y) == " << symbol_table.endo_nbr() << endl
1750 << " @assert length(x) == " << symbol_table.exo_nbr() << endl
1751 << " @assert length(params) == " << symbol_table.param_nbr() << endl
1752 << " @assert length(residual) == " << equations.size() << endl
1753 << " if T0_flag" << endl
1754 << " staticResidTT!(T, y, x, params)" << endl
1755 << " end" << endl
1756 << d_output[0].str()
1757 << " if ~isreal(residual)" << endl
1758 << " residual = real(residual)+imag(residual).^2;" << endl
1759 << " end" << endl
1760 << " return nothing" << endl
1761 << "end" << endl << endl;
1762
1763 // staticG1TT!
1764 output << "function staticG1TT!(T::Vector{Float64}," << endl
1765 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T0_flag::Bool)" << endl
1766 << " if T0_flag" << endl
1767 << " staticResidTT!(T, y, x, params)" << endl
1768 << " end" << endl
1769 << tt_output[1].str()
1770 << " return nothing" << endl
1771 << "end" << endl << endl;
1772
1773 // staticG1!
1774 output << "function staticG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl
1775 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl
1776 << " @assert length(T) >= "
1777 << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl
1778 << " @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr() << ")" << endl
1779 << " @assert length(y) == " << symbol_table.endo_nbr() << endl
1780 << " @assert length(x) == " << symbol_table.exo_nbr() << endl
1781 << " @assert length(params) == " << symbol_table.param_nbr() << endl
1782 << " if T1_flag" << endl
1783 << " staticG1TT!(T, y, x, params, T0_flag)" << endl
1784 << " end" << endl
1785 << " fill!(g1, 0.0)" << endl
1786 << d_output[1].str()
1787 << " if ~isreal(g1)" << endl
1788 << " g1 = real(g1)+2*imag(g1);" << endl
1789 << " end" << endl
1790 << " return nothing" << endl
1791 << "end" << endl << endl;
1792
1793 // staticG2TT!
1794 output << "function staticG2TT!(T::Vector{Float64}," << endl
1795 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl
1796 << " if T1_flag" << endl
1797 << " staticG1TT!(T, y, x, params, TO_flag)" << endl
1798 << " end" << endl
1799 << tt_output[2].str()
1800 << " return nothing" << endl
1801 << "end" << endl << endl;
1802
1803 // staticG2!
1804 output << "function staticG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl
1805 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
1806 << " @assert length(T) >= "
1807 << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl
1808 << " @assert size(g2) == (" << equations.size() << ", " << hessianColsNbr << ")" << endl
1809 << " @assert length(y) == " << symbol_table.endo_nbr() << endl
1810 << " @assert length(x) == " << symbol_table.exo_nbr() << endl
1811 << " @assert length(params) == " << symbol_table.param_nbr() << endl
1812 << " if T2_flag" << endl
1813 << " staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl
1814 << " end" << endl
1815 << " fill!(g2, 0.0)" << endl
1816 << d_output[2].str()
1817 << " return nothing" << endl
1818 << "end" << endl << endl;
1819
1820 // staticG3TT!
1821 output << "function staticG3TT!(T::Vector{Float64}," << endl
1822 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
1823 << " if T2_flag" << endl
1824 << " staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl
1825 << " end" << endl
1826 << tt_output[3].str()
1827 << " return nothing" << endl
1828 << "end" << endl << endl;
1829
1830 // staticG3!
1831 int ncols = hessianColsNbr * JacobianColsNbr;
1832 output << "function staticG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl
1833 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T3_flag::Bool, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
1834 << " @assert length(T) >= "
1835 << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl
1836 << " @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
1837 << " @assert length(y) == " << symbol_table.endo_nbr() << endl
1838 << " @assert length(x) == " << symbol_table.exo_nbr() << endl
1839 << " @assert length(params) == " << symbol_table.param_nbr() << endl
1840 << " if T3_flag" << endl
1841 << " staticG3TT!(T, y, x, params, T2_flag, T1_flag, T0_flag)" << endl
1842 << " end" << endl
1843 << " fill!(g3, 0.0)" << endl
1844 << d_output[3].str()
1845 << " return nothing" << endl
1846 << "end" << endl << endl;
1847
1848 // static!
1849 output << "function static!(T::Vector{Float64}, residual::Vector{Float64}," << endl
1850 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
1851 << " staticResid!(T, residual, y, x, params, true)" << endl
1852 << " return nothing" << endl
1853 << "end" << endl
1854 << endl
1855 << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}," << endl
1856 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
1857 << " staticG1!(T, g1, y, x, params, true)" << endl
1858 << " staticResid!(T, residual, y, x, params, false)" << endl
1859 << " return nothing" << endl
1860 << "end" << endl
1861 << endl
1862 << "function static!(T::Vector{Float64}, g1::Matrix{Float64}," << endl
1863 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
1864 << " staticG1!(T, g1, y, x, params, true, false)" << endl
1865 << " return nothing" << endl
1866 << "end" << endl
1867 << endl
1868 << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl
1869 << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
1870 << " staticG2!(T, g2, y, x, params, true)" << endl
1871 << " staticG1!(T, g1, y, x, params, false)" << endl
1872 << " staticResid!(T, residual, y, x, params, false)" << endl
1873 << " return nothing" << endl
1874 << "end" << endl << endl
1875 << "end" << endl;
1876 output.close();
1877 }
1878 }
1879
1880 void
writeStaticCFile(const string & basename) const1881 StaticModel::writeStaticCFile(const string &basename) const
1882 {
1883 // Writing comments and function definition command
1884 filesystem::create_directories(basename + "/model/src");
1885 string filename = basename + "/model/src/static.c";
1886 string filename_mex = basename + "/model/src/static_mex.c";
1887
1888 int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size();
1889
1890 ofstream output;
1891 output.open(filename, ios::out | ios::binary);
1892 if (!output.is_open())
1893 {
1894 cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
1895 exit(EXIT_FAILURE);
1896 }
1897
1898 output << "/*" << endl
1899 << " * " << filename << " : Computes static model for Dynare" << endl
1900 << " *" << endl
1901 << " * Warning : this file is generated automatically by Dynare" << endl
1902 << " * from model file (.mod)" << endl << endl
1903 << " */" << endl
1904 << "#include <math.h>" << endl;
1905
1906 output << "#include <stdlib.h>" << endl;
1907
1908 if (external_functions_table.get_total_number_of_unique_model_block_external_functions())
1909 // External Matlab function, implies Static function will call mex
1910 output
1911 #ifndef __APPLE__
1912 << "#include <uchar.h>" << endl // For MATLAB ≤ R2011a
1913 #else
1914 << "typedef uint_least16_t char16_t;" << endl
1915 << "typedef uint_least32_t char32_t;" << endl // uchar.h does not exist on macOS
1916 #endif
1917 << R"(#include "mex.h")" << endl;
1918
1919 output << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
1920 << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
1921
1922 // Write function definition if BinaryOpcode::powerDeriv is used
1923 writePowerDerivCHeader(output);
1924
1925 output << endl;
1926
1927 // Writing the function body
1928 writeStaticModel(output, true, false);
1929
1930 output << endl;
1931
1932 writePowerDeriv(output);
1933 output.close();
1934
1935 output.open(filename_mex, ios::out | ios::binary);
1936 if (!output.is_open())
1937 {
1938 cerr << "ERROR: Can't open file " << filename_mex << " for writing" << endl;
1939 exit(EXIT_FAILURE);
1940 }
1941
1942 // Writing the gateway routine
1943 output << "/*" << endl
1944 << " * " << filename_mex << " : The gateway routine used to call the Static function "
1945 << "located in " << filename << endl
1946 << " *" << endl
1947 << " * Warning : this file is generated automatically by Dynare" << endl
1948 << " * from model file (.mod)" << endl << endl
1949 << " */" << endl
1950 << endl
1951 << "#include <stdlib.h>" << endl
1952 #ifndef __APPLE__
1953 << "#include <uchar.h>" << endl // For MATLAB ≤ R2011a
1954 #else
1955 << "typedef uint_least16_t char16_t;" << endl
1956 << "typedef uint_least32_t char32_t;" << endl // uchar.h does not exist on macOS
1957 #endif
1958 << R"(#include "mex.h")" << endl
1959 << endl
1960 << "void static_resid_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T);" << endl
1961 << "void static_resid(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *residual);" << endl
1962 << "void static_g1_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T);" << endl
1963 << "void static_g1(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *g1);" << endl
1964 << "void static_g2_tt(const double *y, const double *x, int nb_row_x, const double *params, double *T);" << endl
1965 << "void static_g2(const double *y, const double *x, int nb_row_x, const double *params, const double *T, double *v2);" << endl
1966 << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
1967 << "{" << endl
1968 << " /* Create a pointer to the input matrix y. */" << endl
1969 << " double *y = mxGetPr(prhs[0]);" << endl
1970 << endl
1971 << " /* Create a pointer to the input matrix x. */" << endl
1972 << " double *x = mxGetPr(prhs[1]);" << endl
1973 << endl
1974 << " /* Create a pointer to the input matrix params. */" << endl
1975 << " double *params = mxGetPr(prhs[2]);" << endl
1976 << endl
1977 << " /* Gets number of rows of matrix x. */" << endl
1978 << " int nb_row_x = mxGetM(prhs[1]);" << endl
1979 << endl
1980 << " double *T = (double *) malloc(sizeof(double)*" << ntt << ");"
1981 << endl
1982 << " if (nlhs >= 1)" << endl
1983 << " {" << endl
1984 << " /* Set the output pointer to the output matrix residual. */" << endl
1985 << " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
1986 << " double *residual = mxGetPr(plhs[0]);" << endl
1987 << " static_resid_tt(y, x, nb_row_x, params, T);" << endl
1988 << " static_resid(y, x, nb_row_x, params, T, residual);" << endl
1989 << " }" << endl
1990 << endl
1991 << " if (nlhs >= 2)" << endl
1992 << " {" << endl
1993 << " /* Set the output pointer to the output matrix g1. */" << endl
1994 << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl
1995 << " double *g1 = mxGetPr(plhs[1]);" << endl
1996 << " static_g1_tt(y, x, nb_row_x, params, T);" << endl
1997 << " static_g1(y, x, nb_row_x, params, T, g1);" << endl
1998 << " }" << endl
1999 << endl
2000 << " if (nlhs >= 3)" << endl
2001 << " {" << endl
2002 << " /* Set the output pointer to the output matrix v2. */" << endl
2003 << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3
2004 << ", mxREAL);" << endl
2005 << " double *v2 = mxGetPr(plhs[2]);" << endl
2006 << " static_g2_tt(y, x, nb_row_x, params, T);" << endl
2007 << " static_g2(y, x, nb_row_x, params, T, v2);" << endl
2008 << " }" << endl
2009 << endl
2010 << " free(T);" << endl
2011 << "}" << endl;
2012 output.close();
2013 }
2014
2015 void
writeStaticJuliaFile(const string & basename) const2016 StaticModel::writeStaticJuliaFile(const string &basename) const
2017 {
2018 writeStaticModel(basename, false, true);
2019 }
2020
2021 void
writeStaticFile(const string & basename,bool block,bool bytecode,bool use_dll,const string & mexext,const filesystem::path & matlabroot,const filesystem::path & dynareroot,bool julia) const2022 StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const
2023 {
2024 if (block && bytecode)
2025 writeModelEquationsCode_Block(basename, map_idx, map_idx2);
2026 else if (!block && bytecode)
2027 writeModelEquationsCode(basename, map_idx);
2028 else if (block && !bytecode)
2029 {
2030 writeModelEquationsOrdered_M(basename);
2031 writeStaticBlockMFSFile(basename);
2032 }
2033 else if (use_dll)
2034 {
2035 writeStaticCFile(basename);
2036 compileDll(basename, "static", mexext, matlabroot, dynareroot);
2037 }
2038 else if (julia)
2039 writeStaticJuliaFile(basename);
2040 else
2041 writeStaticMFile(basename);
2042 writeSetAuxiliaryVariables(basename, julia);
2043 }
2044
2045 bool
exoPresentInEqs() const2046 StaticModel::exoPresentInEqs() const
2047 {
2048 for (auto equation : equations)
2049 if (equation->containsExogenous())
2050 return true;
2051 return false;
2052 }
2053
2054 void
writeStaticBlockMFSFile(const string & basename) const2055 StaticModel::writeStaticBlockMFSFile(const string &basename) const
2056 {
2057 string filename = packageDir(basename) + "/static.m";
2058
2059 ofstream output;
2060 output.open(filename, ios::out | ios::binary);
2061 if (!output.is_open())
2062 {
2063 cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
2064 exit(EXIT_FAILURE);
2065 }
2066
2067 output << "function [residual, g1, y, var_index] = static(nblock, y, x, params)" << endl
2068 << " residual = [];" << endl
2069 << " g1 = [];" << endl
2070 << " var_index = [];" << endl << endl
2071 << " switch nblock" << endl;
2072
2073 unsigned int nb_blocks = getNbBlocks();
2074
2075 for (int b = 0; b < static_cast<int>(nb_blocks); b++)
2076 {
2077 set<int> local_var;
2078
2079 output << " case " << b+1 << endl;
2080
2081 BlockSimulationType simulation_type = getBlockSimulationType(b);
2082
2083 if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
2084 {
2085 output << " y_tmp = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl;
2086 ostringstream tmp;
2087 for (int i = 0; i < static_cast<int>(getBlockSize(b)); i++)
2088 tmp << " " << getBlockVariableID(b, i)+1;
2089 output << " var_index = [" << tmp.str() << "];" << endl
2090 << " residual = y(var_index) - y_tmp(var_index);" << endl
2091 << " y = y_tmp;" << endl;
2092 }
2093 else
2094 output << " [residual, y, g1] = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl;
2095
2096 }
2097 output << " end" << endl
2098 << "end" << endl;
2099 output.close();
2100 }
2101
2102 void
writeOutput(ostream & output,bool block) const2103 StaticModel::writeOutput(ostream &output, bool block) const
2104 {
2105 output << "M_.static_tmp_nbr = [";
2106 for (const auto &temporary_terms_derivative : temporary_terms_derivatives)
2107 output << temporary_terms_derivative.size() << "; ";
2108 output << "];" << endl;
2109
2110 /* Write mapping between model local variables and indices in the temporary
2111 terms vector (dynare#1722) */
2112 output << "M_.model_local_variables_static_tt_idxs = {" << endl;
2113 for (auto [mlv, value] : temporary_terms_mlv)
2114 output << " '" << symbol_table.getName(mlv->symb_id) << "', "
2115 << temporary_terms_idxs.at(mlv)+1 << ';' << endl;
2116 output << "};" << endl;
2117
2118 if (!block)
2119 return;
2120
2121 unsigned int nb_blocks = getNbBlocks();
2122 for (int b = 0; b < static_cast<int>(nb_blocks); b++)
2123 {
2124 BlockSimulationType simulation_type = getBlockSimulationType(b);
2125 unsigned int block_size = getBlockSize(b);
2126 ostringstream tmp_s, tmp_s_eq;
2127 tmp_s.str("");
2128 tmp_s_eq.str("");
2129 for (unsigned int i = 0; i < block_size; i++)
2130 {
2131 tmp_s << " " << getBlockVariableID(b, i)+1;
2132 tmp_s_eq << " " << getBlockEquationID(b, i)+1;
2133 }
2134 output << "block_structure_stat.block(" << b+1 << ").Simulation_Type = " << simulation_type << ";" << endl
2135 << "block_structure_stat.block(" << b+1 << ").endo_nbr = " << block_size << ";" << endl
2136 << "block_structure_stat.block(" << b+1 << ").mfs = " << getBlockMfs(b) << ";" << endl
2137 << "block_structure_stat.block(" << b+1 << ").equation = [" << tmp_s_eq.str() << "];" << endl
2138 << "block_structure_stat.block(" << b+1 << ").variable = [" << tmp_s.str() << "];" << endl;
2139 }
2140 output << "M_.block_structure_stat.block = block_structure_stat.block;" << endl;
2141 string cst_s;
2142 int nb_endo = symbol_table.endo_nbr();
2143 output << "M_.block_structure_stat.variable_reordered = [";
2144 for (int i = 0; i < nb_endo; i++)
2145 output << " " << variable_reordered[i]+1;
2146 output << "];" << endl
2147 << "M_.block_structure_stat.equation_reordered = [";
2148 for (int i = 0; i < nb_endo; i++)
2149 output << " " << equation_reordered[i]+1;
2150 output << "];" << endl;
2151
2152 map<pair<int, int>, int> row_incidence;
2153 for (const auto & [indices, d1] : derivatives[1])
2154 {
2155 int deriv_id = indices[1];
2156 if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
2157 {
2158 int eq = indices[0];
2159 int symb = getSymbIDByDerivID(deriv_id);
2160 int var = symbol_table.getTypeSpecificID(symb);
2161 //int lag = getLagByDerivID(deriv_id);
2162 row_incidence[{ eq, var }] = 1;
2163 }
2164 }
2165 output << "M_.block_structure_stat.incidence.sparse_IM = [";
2166 for (const auto &it : row_incidence)
2167 output << it.first.first+1 << " " << it.first.second+1 << ";" << endl;
2168 output << "];" << endl;
2169 }
2170
2171 SymbolType
getTypeByDerivID(int deriv_id) const2172 StaticModel::getTypeByDerivID(int deriv_id) const noexcept(false)
2173 {
2174 if (deriv_id < symbol_table.endo_nbr())
2175 return SymbolType::endogenous;
2176 else if (deriv_id < symbol_table.endo_nbr() + symbol_table.param_nbr())
2177 return SymbolType::parameter;
2178 else
2179 throw UnknownDerivIDException();
2180 }
2181
2182 int
getLagByDerivID(int deriv_id) const2183 StaticModel::getLagByDerivID(int deriv_id) const noexcept(false)
2184 {
2185 return 0;
2186 }
2187
2188 int
getSymbIDByDerivID(int deriv_id) const2189 StaticModel::getSymbIDByDerivID(int deriv_id) const noexcept(false)
2190 {
2191 if (deriv_id < symbol_table.endo_nbr())
2192 return symbol_table.getID(SymbolType::endogenous, deriv_id);
2193 else if (deriv_id < symbol_table.endo_nbr() + symbol_table.param_nbr())
2194 return symbol_table.getID(SymbolType::parameter, deriv_id - symbol_table.endo_nbr());
2195 else
2196 throw UnknownDerivIDException();
2197 }
2198
2199 int
getDerivID(int symb_id,int lag) const2200 StaticModel::getDerivID(int symb_id, int lag) const noexcept(false)
2201 {
2202 if (symbol_table.getType(symb_id) == SymbolType::endogenous)
2203 return symbol_table.getTypeSpecificID(symb_id);
2204 else if (symbol_table.getType(symb_id) == SymbolType::parameter)
2205 return symbol_table.getTypeSpecificID(symb_id) + symbol_table.endo_nbr();
2206 else
2207 return -1;
2208 }
2209
2210 void
addAllParamDerivId(set<int> & deriv_id_set)2211 StaticModel::addAllParamDerivId(set<int> &deriv_id_set)
2212 {
2213 for (int i = 0; i < symbol_table.param_nbr(); i++)
2214 deriv_id_set.insert(i + symbol_table.endo_nbr());
2215 }
2216
2217 map<tuple<int, int, int, int, int>, int>
get_Derivatives(int block)2218 StaticModel::get_Derivatives(int block)
2219 {
2220 map<tuple<int, int, int, int, int>, int> Derivatives;
2221 int block_size = getBlockSize(block);
2222 int block_nb_recursive = block_size - getBlockMfs(block);
2223 int lag = 0;
2224 for (int eq = 0; eq < block_size; eq++)
2225 {
2226 int eqr = getBlockEquationID(block, eq);
2227 for (int var = 0; var < block_size; var++)
2228 {
2229 int varr = getBlockVariableID(block, var);
2230 if (dynamic_jacobian.find({ lag, eqr, varr }) != dynamic_jacobian.end())
2231 {
2232 bool OK = true;
2233 if (auto its = Derivatives.find({ lag, eq, var, eqr, varr });
2234 its != Derivatives.end() && its->second == 2)
2235 OK = false;
2236
2237 if (OK)
2238 {
2239 if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursive)
2240 //It's a normalized equation, we have to recompute the derivative using chain rule derivative function
2241 Derivatives[{ lag, eq, var, eqr, varr }] = 1;
2242 else
2243 //It's a feedback equation we can use the derivatives
2244 Derivatives[{ lag, eq, var, eqr, varr }] = 0;
2245 }
2246 if (var < block_nb_recursive)
2247 {
2248 int eqs = getBlockEquationID(block, var);
2249 for (int vars = block_nb_recursive; vars < block_size; vars++)
2250 {
2251 int varrs = getBlockVariableID(block, vars);
2252 //A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation)
2253 if (Derivatives.find({ lag, var, vars, eqs, varrs }) != Derivatives.end())
2254 Derivatives[{ lag, eq, vars, eqr, varrs }] = 2;
2255 }
2256 }
2257 }
2258 }
2259 }
2260
2261 return Derivatives;
2262 }
2263
2264 void
computeChainRuleJacobian(blocks_derivatives_t & blocks_derivatives)2265 StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
2266 {
2267 map<int, expr_t> recursive_variables;
2268 unsigned int nb_blocks = getNbBlocks();
2269 blocks_derivatives = blocks_derivatives_t(nb_blocks);
2270 for (unsigned int block = 0; block < nb_blocks; block++)
2271 {
2272 block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
2273 recursive_variables.clear();
2274 BlockSimulationType simulation_type = getBlockSimulationType(block);
2275 int block_size = getBlockSize(block);
2276 int block_nb_mfs = getBlockMfs(block);
2277 int block_nb_recursives = block_size - block_nb_mfs;
2278 if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
2279 {
2280 blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
2281 for (int i = 0; i < block_nb_recursives; i++)
2282 {
2283 if (getBlockEquationType(block, i) == E_EVALUATE_S)
2284 recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
2285 else
2286 recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
2287 }
2288 auto Derivatives = get_Derivatives(block);
2289 for (const auto &it : Derivatives)
2290 {
2291 int Deriv_type = it.second;
2292 auto [lag, eq, var, eqr, varr] = it.first;
2293 if (Deriv_type == 0)
2294 first_chain_rule_derivatives[{ eqr, varr, lag }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
2295 else if (Deriv_type == 1)
2296 first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
2297 else if (Deriv_type == 2)
2298 {
2299 if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursives)
2300 first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
2301 else
2302 first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
2303 }
2304 tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]);
2305 }
2306 }
2307 else
2308 {
2309 blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
2310 for (int i = 0; i < block_nb_recursives; i++)
2311 {
2312 if (getBlockEquationType(block, i) == E_EVALUATE_S)
2313 recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
2314 else
2315 recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
2316 }
2317 for (int eq = block_nb_recursives; eq < block_size; eq++)
2318 {
2319 int eqr = getBlockEquationID(block, eq);
2320 for (int var = block_nb_recursives; var < block_size; var++)
2321 {
2322 int varr = getBlockVariableID(block, var);
2323 expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), 0), recursive_variables);
2324 if (d1 == Zero)
2325 continue;
2326 first_chain_rule_derivatives[{ eqr, varr, 0 }] = d1;
2327 tmp_derivatives.emplace_back(eq, var, 0, first_chain_rule_derivatives[{ eqr, varr, 0 }]);
2328 }
2329 }
2330 }
2331 blocks_derivatives[block] = tmp_derivatives;
2332 }
2333 }
2334
2335 void
collect_block_first_order_derivatives()2336 StaticModel::collect_block_first_order_derivatives()
2337 {
2338 //! vector for an equation or a variable indicates the block number
2339 vector<int> equation_2_block(equation_reordered.size()), variable_2_block(variable_reordered.size());
2340 unsigned int nb_blocks = getNbBlocks();
2341 for (unsigned int block = 0; block < nb_blocks; block++)
2342 {
2343 unsigned int block_size = getBlockSize(block);
2344 for (unsigned int i = 0; i < block_size; i++)
2345 {
2346 equation_2_block[getBlockEquationID(block, i)] = block;
2347 variable_2_block[getBlockVariableID(block, i)] = block;
2348 }
2349 }
2350 derivative_endo = vector<derivative_t>(nb_blocks);
2351 endo_max_leadlag_block = vector<pair<int, int>>(nb_blocks, { 0, 0 });
2352 max_leadlag_block = vector<pair<int, int>>(nb_blocks, { 0, 0 });
2353 for (auto &first_derivative : derivatives[1])
2354 {
2355 int eq = first_derivative.first[0];
2356 int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1]));
2357 int lag = 0;
2358 int block_eq = equation_2_block[eq];
2359 int block_var = variable_2_block[var];
2360 max_leadlag_block[block_eq] = { 0, 0 };
2361 max_leadlag_block[block_eq] = { 0, 0 };
2362 endo_max_leadlag_block[block_eq] = { 0, 0 };
2363 endo_max_leadlag_block[block_eq] = { 0, 0 };
2364 derivative_t tmp_derivative;
2365 lag_var_t lag_var;
2366 if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous && block_eq == block_var)
2367 {
2368 tmp_derivative = derivative_endo[block_eq];
2369 tmp_derivative[{ lag, eq, var }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
2370 derivative_endo[block_eq] = tmp_derivative;
2371 }
2372 }
2373 }
2374
2375 void
writeLatexFile(const string & basename,bool write_equation_tags) const2376 StaticModel::writeLatexFile(const string &basename, bool write_equation_tags) const
2377 {
2378 writeLatexModelFile(basename, "static", ExprNodeOutputType::latexStaticModel, write_equation_tags);
2379 }
2380
2381 void
writeAuxVarInitval(ostream & output,ExprNodeOutputType output_type) const2382 StaticModel::writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type) const
2383 {
2384 for (auto aux_equation : aux_equations)
2385 {
2386 dynamic_cast<ExprNode *>(aux_equation)->writeOutput(output, output_type);
2387 output << ";" << endl;
2388 }
2389 }
2390
2391 void
writeSetAuxiliaryVariables(const string & basename,bool julia) const2392 StaticModel::writeSetAuxiliaryVariables(const string &basename, bool julia) const
2393 {
2394 ostringstream output_func_body;
2395 writeAuxVarRecursiveDefinitions(output_func_body, ExprNodeOutputType::matlabStaticModel);
2396
2397 if (output_func_body.str().empty())
2398 return;
2399
2400 string func_name = julia ? basename + "_set_auxiliary_variables" : "set_auxiliary_variables";
2401 string filename = julia ? func_name + ".jl" : packageDir(basename) + "/" + func_name + ".m";
2402 string comment = julia ? "#" : "%";
2403
2404 ofstream output;
2405 output.open(filename, ios::out | ios::binary);
2406 if (!output.is_open())
2407 {
2408 cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
2409 exit(EXIT_FAILURE);
2410 }
2411
2412 output << "function y = " << func_name + "(y, x, params)" << endl
2413 << comment << endl
2414 << comment << " Status : Computes static model for Dynare" << endl
2415 << comment << endl
2416 << comment << " Warning : this file is generated automatically by Dynare" << endl
2417 << comment << " from model file (.mod)" << endl << endl
2418 << output_func_body.str();
2419
2420 output.close();
2421 }
2422
2423 void
writeAuxVarRecursiveDefinitions(ostream & output,ExprNodeOutputType output_type) const2424 StaticModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const
2425 {
2426 deriv_node_temp_terms_t tef_terms;
2427 for (auto aux_equation : aux_equations)
2428 if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
2429 dynamic_cast<ExprNode *>(aux_equation)->writeExternalFunctionOutput(output, ExprNodeOutputType::matlabStaticModel, {}, {}, tef_terms);
2430 for (auto aux_equation : aux_equations)
2431 {
2432 dynamic_cast<ExprNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->writeOutput(output, output_type);
2433 output << ";" << endl;
2434 }
2435 }
2436
2437 void
writeLatexAuxVarRecursiveDefinitions(ostream & output) const2438 StaticModel::writeLatexAuxVarRecursiveDefinitions(ostream &output) const
2439 {
2440 deriv_node_temp_terms_t tef_terms;
2441 temporary_terms_t temporary_terms;
2442 temporary_terms_idxs_t temporary_terms_idxs;
2443 for (auto aux_equation : aux_equations)
2444 if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
2445 dynamic_cast<ExprNode *>(aux_equation)->writeExternalFunctionOutput(output, ExprNodeOutputType::latexStaticModel,
2446 temporary_terms, temporary_terms_idxs, tef_terms);
2447 for (auto aux_equation : aux_equations)
2448 {
2449 output << R"(\begin{dmath})" << endl;
2450 dynamic_cast<ExprNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->writeOutput(output, ExprNodeOutputType::latexStaticModel);
2451 output << endl << R"(\end{dmath})" << endl;
2452 }
2453 }
2454
2455 void
writeJsonAuxVarRecursiveDefinitions(ostream & output) const2456 StaticModel::writeJsonAuxVarRecursiveDefinitions(ostream &output) const
2457 {
2458 deriv_node_temp_terms_t tef_terms;
2459 temporary_terms_t temporary_terms;
2460
2461 for (auto aux_equation : aux_equations)
2462 if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
2463 {
2464 vector<string> efout;
2465 dynamic_cast<ExprNode *>(aux_equation)->writeJsonExternalFunctionOutput(efout,
2466 temporary_terms,
2467 tef_terms,
2468 false);
2469 for (auto it = efout.begin(); it != efout.end(); ++it)
2470 {
2471 if (it != efout.begin())
2472 output << ", ";
2473 output << *it;
2474 }
2475 }
2476
2477 for (auto aux_equation : aux_equations)
2478 {
2479 output << R"(, {"lhs": ")";
2480 aux_equation->arg1->writeJsonOutput(output, temporary_terms, tef_terms, false);
2481 output << R"(", "rhs": ")";
2482 dynamic_cast<BinaryOpNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->arg2->writeJsonOutput(output, temporary_terms, tef_terms, false);
2483 output << R"("})";
2484 }
2485 }
2486
2487 void
writeParamsDerivativesFile(const string & basename,bool julia) const2488 StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) const
2489 {
2490 if (!params_derivatives.size())
2491 return;
2492
2493 ExprNodeOutputType output_type = (julia ? ExprNodeOutputType::juliaStaticModel : ExprNodeOutputType::matlabStaticModel);
2494
2495 ostringstream tt_output; // Used for storing temporary terms
2496 ostringstream jacobian_output; // Used for storing jacobian equations
2497 ostringstream hessian_output; // Used for storing Hessian equations
2498 ostringstream hessian1_output; // Used for storing Hessian equations
2499 ostringstream third_derivs_output; // Used for storing third order derivatives equations
2500 ostringstream third_derivs1_output; // Used for storing third order derivatives equations
2501
2502 temporary_terms_t temp_term_union;
2503 deriv_node_temp_terms_t tef_terms;
2504
2505 writeModelLocalVariableTemporaryTerms(temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
2506 for (const auto &it : params_derivs_temporary_terms)
2507 writeTemporaryTerms(it.second, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
2508
2509 for (const auto & [indices, d1] : params_derivatives.find({ 0, 1 })->second)
2510 {
2511 auto [eq, param] = vectorToTuple<2>(indices);
2512
2513 int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
2514
2515 jacobian_output << "rp" << LEFT_ARRAY_SUBSCRIPT(output_type)
2516 << eq+1 << ", " << param_col
2517 << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
2518 d1->writeOutput(jacobian_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
2519 jacobian_output << ";" << endl;
2520 }
2521
2522 for (const auto & [indices, d2] : params_derivatives.find({ 1, 1 })->second)
2523 {
2524 auto [eq, var, param] = vectorToTuple<3>(indices);
2525
2526 int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
2527 int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
2528
2529 hessian_output << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type)
2530 << eq+1 << ", " << var_col << ", " << param_col
2531 << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
2532 d2->writeOutput(hessian_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
2533 hessian_output << ";" << endl;
2534 }
2535
2536 int i = 1;
2537 for (const auto &[indices, d2] : params_derivatives.find({ 0, 2 })->second)
2538 {
2539 auto [eq, param1, param2] = vectorToTuple<3>(indices);
2540
2541 int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
2542 int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
2543
2544 hessian1_output << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
2545 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
2546 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
2547 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
2548 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
2549 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
2550 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
2551 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
2552 d2->writeOutput(hessian1_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
2553 hessian1_output << ";" << endl;
2554
2555 i++;
2556
2557 if (param1 != param2)
2558 {
2559 // Treat symmetric elements
2560 hessian1_output << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
2561 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
2562 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
2563 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
2564 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
2565 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
2566 << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
2567 << RIGHT_ARRAY_SUBSCRIPT(output_type)
2568 << "=rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",4"
2569 << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
2570 i++;
2571 }
2572 }
2573
2574 i = 1;
2575 for (const auto &[indices, d2] : params_derivatives.find({ 1, 2 })->second)
2576 {
2577 auto [eq, var, param1, param2] = vectorToTuple<4>(indices);
2578
2579 int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
2580 int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
2581 int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
2582
2583 third_derivs_output << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
2584 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
2585 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
2586 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
2587 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
2588 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
2589 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
2590 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
2591 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
2592 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
2593 d2->writeOutput(third_derivs_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
2594 third_derivs_output << ";" << endl;
2595
2596 i++;
2597
2598 if (param1 != param2)
2599 {
2600 // Treat symmetric elements
2601 third_derivs_output << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
2602 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
2603 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
2604 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
2605 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
2606 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
2607 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
2608 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
2609 << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
2610 << RIGHT_ARRAY_SUBSCRIPT(output_type)
2611 << "=gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
2612 << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
2613 i++;
2614 }
2615 }
2616
2617 i = 1;
2618 for (const auto &[indices, d2] : params_derivatives.find({ 2, 1 })->second)
2619 {
2620 auto [eq, var1, var2, param] = vectorToTuple<4>(indices);
2621
2622 int var1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var1)) + 1;
2623 int var2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var2)) + 1;
2624 int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
2625
2626 third_derivs1_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
2627 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
2628 << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
2629 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
2630 << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
2631 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
2632 << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
2633 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
2634 << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
2635 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
2636 d2->writeOutput(third_derivs1_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs, tef_terms);
2637 third_derivs1_output << ";" << endl;
2638
2639 i++;
2640
2641 if (var1 != var2)
2642 {
2643 // Treat symmetric elements
2644 third_derivs1_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
2645 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
2646 << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
2647 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
2648 << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
2649 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
2650 << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
2651 << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
2652 << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
2653 << RIGHT_ARRAY_SUBSCRIPT(output_type)
2654 << "=hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i-1 << ",5"
2655 << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
2656 i++;
2657 }
2658 }
2659
2660 ofstream paramsDerivsFile;
2661 string filename = julia ? basename + "StaticParamsDerivs.jl" : packageDir(basename) + "/static_params_derivs.m";
2662 paramsDerivsFile.open(filename, ios::out | ios::binary);
2663 if (!paramsDerivsFile.is_open())
2664 {
2665 cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
2666 exit(EXIT_FAILURE);
2667 }
2668
2669 if (!julia)
2670 {
2671 // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201
2672 map<string, string> tmp_paren_vars;
2673 bool message_printed = false;
2674 fixNestedParenthesis(tt_output, tmp_paren_vars, message_printed);
2675 fixNestedParenthesis(jacobian_output, tmp_paren_vars, message_printed);
2676 fixNestedParenthesis(hessian_output, tmp_paren_vars, message_printed);
2677 fixNestedParenthesis(hessian1_output, tmp_paren_vars, message_printed);
2678 fixNestedParenthesis(third_derivs_output, tmp_paren_vars, message_printed);
2679 fixNestedParenthesis(third_derivs1_output, tmp_paren_vars, message_printed);
2680
2681 paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = static_params_derivs(y, x, params)" << endl
2682 << "%" << endl
2683 << "% Status : Computes derivatives of the static model with respect to the parameters" << endl
2684 << "%" << endl
2685 << "% Inputs : " << endl
2686 << "% y [M_.endo_nbr by 1] double vector of endogenous variables in declaration order" << endl
2687 << "% x [M_.exo_nbr by 1] double vector of exogenous variables in declaration order" << endl
2688 << "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
2689 << "%" << endl
2690 << "% Outputs:" << endl
2691 << "% rp [M_.eq_nbr by #params] double Jacobian matrix of static model equations with respect to parameters " << endl
2692 << "% Dynare may prepend or append auxiliary equations, see M_.aux_vars" << endl
2693 << "% gp [M_.endo_nbr by M_.endo_nbr by #params] double Derivative of the Jacobian matrix of the static model equations with respect to the parameters" << endl
2694 << "% rows: variables in declaration order" << endl
2695 << "% rows: equations in order of declaration" << endl
2696 << "% rpp [#second_order_residual_terms by 4] double Hessian matrix of second derivatives of residuals with respect to parameters;" << endl
2697 << "% rows: respective derivative term" << endl
2698 << "% 1st column: equation number of the term appearing" << endl
2699 << "% 2nd column: number of the first parameter in derivative" << endl
2700 << "% 3rd column: number of the second parameter in derivative" << endl
2701 << "% 4th column: value of the Hessian term" << endl
2702 << "% gpp [#second_order_Jacobian_terms by 5] double Hessian matrix of second derivatives of the Jacobian with respect to the parameters;" << endl
2703 << "% rows: respective derivative term" << endl
2704 << "% 1st column: equation number of the term appearing" << endl
2705 << "% 2nd column: column number of variable in Jacobian of the static model" << endl
2706 << "% 3rd column: number of the first parameter in derivative" << endl
2707 << "% 4th column: number of the second parameter in derivative" << endl
2708 << "% 5th column: value of the Hessian term" << endl
2709 << "%" << endl
2710 << "%" << endl
2711 << "% Warning : this file is generated automatically by Dynare" << endl
2712 << "% from model file (.mod)" << endl << endl
2713 << "T = NaN(" << params_derivs_temporary_terms_idxs.size() << ",1);" << endl
2714 << tt_output.str()
2715 << "rp = zeros(" << equations.size() << ", "
2716 << symbol_table.param_nbr() << ");" << endl
2717 << jacobian_output.str()
2718 << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
2719 << symbol_table.param_nbr() << ");" << endl
2720 << hessian_output.str()
2721 << "if nargout >= 3" << endl
2722 << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
2723 << hessian1_output.str()
2724 << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
2725 << third_derivs_output.str()
2726 << "end" << endl
2727 << "if nargout >= 5" << endl
2728 << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
2729 << third_derivs1_output.str()
2730 << "end" << endl
2731 << "end" << endl;
2732 }
2733 else
2734 paramsDerivsFile << "module " << basename << "StaticParamsDerivs" << endl
2735 << "#" << endl
2736 << "# NB: this file was automatically generated by Dynare" << endl
2737 << "# from " << basename << ".mod" << endl
2738 << "#" << endl
2739 << "export params_derivs" << endl << endl
2740 << "function params_derivs(y, x, params)" << endl
2741 << tt_output.str()
2742 << "rp = zeros(" << equations.size() << ", "
2743 << symbol_table.param_nbr() << ");" << endl
2744 << jacobian_output.str()
2745 << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
2746 << symbol_table.param_nbr() << ");" << endl
2747 << hessian_output.str()
2748 << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
2749 << hessian1_output.str()
2750 << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
2751 << third_derivs_output.str()
2752 << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
2753 << third_derivs1_output.str()
2754 << "(rp, gp, rpp, gpp, hp)" << endl
2755 << "end" << endl
2756 << "end" << endl;
2757
2758 paramsDerivsFile.close();
2759 }
2760
2761 void
writeJsonOutput(ostream & output) const2762 StaticModel::writeJsonOutput(ostream &output) const
2763 {
2764 deriv_node_temp_terms_t tef_terms;
2765 writeJsonModelLocalVariables(output, false, tef_terms);
2766 output << ", ";
2767 writeJsonModelEquations(output, false);
2768 }
2769
2770 void
writeJsonComputingPassOutput(ostream & output,bool writeDetails) const2771 StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) const
2772 {
2773 ostringstream model_local_vars_output; // Used for storing model local vars
2774 vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
2775
2776 deriv_node_temp_terms_t tef_terms;
2777 temporary_terms_t temp_term_union;
2778
2779 writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
2780
2781 writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms, "");
2782 d_output[0] << ", ";
2783 writeJsonModelEquations(d_output[0], true);
2784
2785 auto getJacobCol = [this](int var) { return symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)); };
2786 int ncols = symbol_table.endo_nbr();
2787 for (size_t i = 1; i < derivatives.size(); i++)
2788 {
2789 string matrix_name = i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative";
2790 writeJsonTemporaryTerms(temporary_terms_derivatives[i], temp_term_union, d_output[i], tef_terms, matrix_name);
2791 temp_term_union.insert(temporary_terms_derivatives[i].begin(), temporary_terms_derivatives[i].end());
2792 d_output[i] << R"(, ")" << matrix_name << R"(": {)"
2793 << R"( "nrows": )" << equations.size()
2794 << R"(, "ncols": )" << ncols
2795 << R"(, "entries": [)";
2796
2797 for (auto it = derivatives[i].begin(); it != derivatives[i].end(); ++it)
2798 {
2799 if (it != derivatives[i].begin())
2800 d_output[i] << ", ";
2801
2802 const vector<int> &vidx = it->first;
2803 expr_t d = it->second;
2804 int eq = vidx[0];
2805
2806 int col_idx = 0;
2807 for (size_t j = 1; j < vidx.size(); j++)
2808 {
2809 col_idx *= symbol_table.endo_nbr();
2810 col_idx += getJacobCol(vidx[j]);
2811 }
2812
2813 if (writeDetails)
2814 d_output[i] << R"({"eq": )" << eq + 1;
2815 else
2816 d_output[i] << R"({"row": )" << eq + 1;
2817
2818 d_output[i] << R"(, "col": )" << (i > 1 ? "[" : "") << col_idx + 1;
2819
2820 if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian
2821 {
2822 int col_idx_sym = getJacobCol(vidx[2]) * symbol_table.endo_nbr() + getJacobCol(vidx[1]);
2823 d_output[i] << ", " << col_idx_sym + 1;
2824 }
2825 if (i > 1)
2826 d_output[i] << "]";
2827
2828 if (writeDetails)
2829 for (size_t j = 1; j < vidx.size(); j++)
2830 d_output[i] << R"(, "var)" << (i > 1 ? to_string(j) : "") << R"(": ")" << symbol_table.getName(getSymbIDByDerivID(vidx[j])) << R"(")";
2831
2832 d_output[i] << R"(, "val": ")";
2833 d->writeJsonOutput(d_output[i], temp_term_union, tef_terms);
2834 d_output[i] << R"("})" << endl;
2835 }
2836 d_output[i] << "]}";
2837
2838 ncols *= symbol_table.endo_nbr();
2839 }
2840
2841 if (writeDetails)
2842 output << R"("static_model": {)";
2843 else
2844 output << R"("static_model_simple": {)";
2845 output << model_local_vars_output.str();
2846 for (const auto &it : d_output)
2847 output << ", " << it.str();
2848 output << "}";
2849 }
2850
2851 void
writeJsonParamsDerivativesFile(ostream & output,bool writeDetails) const2852 StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const
2853 {
2854 if (!params_derivatives.size())
2855 return;
2856
2857 ostringstream model_local_vars_output; // Used for storing model local vars
2858 ostringstream model_output; // Used for storing model temp vars and equations
2859 ostringstream jacobian_output; // Used for storing jacobian equations
2860 ostringstream hessian_output; // Used for storing Hessian equations
2861 ostringstream hessian1_output; // Used for storing Hessian equations
2862 ostringstream third_derivs_output; // Used for storing third order derivatives equations
2863 ostringstream third_derivs1_output; // Used for storing third order derivatives equations
2864
2865 deriv_node_temp_terms_t tef_terms;
2866 writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
2867
2868 temporary_terms_t temp_term_union;
2869 for (const auto &it : params_derivs_temporary_terms)
2870 writeJsonTemporaryTerms(it.second, temp_term_union, model_output, tef_terms, "all");
2871
2872 jacobian_output << R"("deriv_wrt_params": {)"
2873 << R"( "neqs": )" << equations.size()
2874 << R"(, "nparamcols": )" << symbol_table.param_nbr()
2875 << R"(, "entries": [)";
2876 auto &rp = params_derivatives.find({ 0, 1 })->second;
2877 for (auto it = rp.begin(); it != rp.end(); ++it)
2878 {
2879 if (it != rp.begin())
2880 jacobian_output << ", ";
2881
2882 auto [eq, param] = vectorToTuple<2>(it->first);
2883 expr_t d1 = it->second;
2884
2885 int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
2886
2887 if (writeDetails)
2888 jacobian_output << R"({"eq": )" << eq + 1;
2889 else
2890 jacobian_output << R"({"row": )" << eq + 1;
2891
2892 if (writeDetails)
2893 jacobian_output << R"(, "param_col": )" << param_col;
2894
2895 jacobian_output << R"(, "param": ")" << symbol_table.getName(getSymbIDByDerivID(param)) << R"(")";
2896
2897 jacobian_output << R"(, "val": ")";
2898 d1->writeJsonOutput(jacobian_output, temp_term_union, tef_terms);
2899 jacobian_output << R"("})" << endl;
2900 }
2901 jacobian_output << "]}";
2902
2903 hessian_output << R"("deriv_jacobian_wrt_params": {)"
2904 << R"( "neqs": )" << equations.size()
2905 << R"(, "nvarcols": )" << symbol_table.endo_nbr()
2906 << R"(, "nparamcols": )" << symbol_table.param_nbr()
2907 << R"(, "entries": [)";
2908 auto &gp = params_derivatives.find({ 1, 1 })->second;
2909 for (auto it = gp.begin(); it != gp.end(); ++it)
2910 {
2911 if (it != gp.begin())
2912 hessian_output << ", ";
2913
2914 auto [eq, var, param] = vectorToTuple<3>(it->first);
2915 expr_t d2 = it->second;
2916
2917 int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
2918 int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
2919
2920 if (writeDetails)
2921 hessian_output << R"({"eq": )" << eq + 1;
2922 else
2923 hessian_output << R"({"row": )" << eq + 1;
2924
2925 if (writeDetails)
2926 hessian_output << R"(, "var": ")" << symbol_table.getName(getSymbIDByDerivID(var)) << R"(")"
2927 << R"(, "param": ")" << symbol_table.getName(getSymbIDByDerivID(param)) << R"(")";
2928
2929 hessian_output << R"(, "var_col": )" << var_col
2930 << R"(, "param_col": )" << param_col
2931 << R"(, "val": ")";
2932 d2->writeJsonOutput(hessian_output, temp_term_union, tef_terms);
2933 hessian_output << R"("})" << endl;
2934 }
2935 hessian_output << "]}";
2936
2937 hessian1_output << R"("second_deriv_residuals_wrt_params": {)"
2938 << R"( "nrows": )" << equations.size()
2939 << R"(, "nparam1cols": )" << symbol_table.param_nbr()
2940 << R"(, "nparam2cols": )" << symbol_table.param_nbr()
2941 << R"(, "entries": [)";
2942 auto &rpp = params_derivatives.find({ 0, 2 })->second;
2943 for (auto it = rpp.begin(); it != rpp.end(); ++it)
2944 {
2945 if (it != rpp.begin())
2946 hessian1_output << ", ";
2947
2948 auto [eq, param1, param2] = vectorToTuple<3>(it->first);
2949 expr_t d2 = it->second;
2950
2951 int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
2952 int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
2953
2954 if (writeDetails)
2955 hessian1_output << R"({"eq": )" << eq + 1;
2956 else
2957 hessian1_output << R"({"row": )" << eq + 1;
2958
2959 hessian1_output << R"(, "param1_col": )" << param1_col
2960 << R"(, "param2_col": )" << param2_col;
2961
2962 if (writeDetails)
2963 hessian1_output << R"(, "param1": ")" << symbol_table.getName(getSymbIDByDerivID(param1)) << R"(")"
2964 << R"(, "param2": ")" << symbol_table.getName(getSymbIDByDerivID(param2)) << R"(")";
2965
2966 hessian1_output << R"(, "val": ")";
2967 d2->writeJsonOutput(hessian1_output, temp_term_union, tef_terms);
2968 hessian1_output << R"("})" << endl;
2969 }
2970 hessian1_output << "]}";
2971
2972 third_derivs_output << R"("second_deriv_jacobian_wrt_params": {)"
2973 << R"( "neqs": )" << equations.size()
2974 << R"(, "nvarcols": )" << symbol_table.endo_nbr()
2975 << R"(, "nparam1cols": )" << symbol_table.param_nbr()
2976 << R"(, "nparam2cols": )" << symbol_table.param_nbr()
2977 << R"(, "entries": [)";
2978 auto &gpp = params_derivatives.find({ 1, 2 })->second;
2979 for (auto it = gpp.begin(); it != gpp.end(); ++it)
2980 {
2981 if (it != gpp.begin())
2982 third_derivs_output << ", ";
2983
2984 auto [eq, var, param1, param2] = vectorToTuple<4>(it->first);
2985 expr_t d2 = it->second;
2986
2987 int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
2988 int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
2989 int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
2990
2991 if (writeDetails)
2992 third_derivs_output << R"({"eq": )" << eq + 1;
2993 else
2994 third_derivs_output << R"({"row": )" << eq + 1;
2995 third_derivs_output << R"(, "var_col": )" << var_col
2996 << R"(, "param1_col": )" << param1_col
2997 << R"(, "param2_col": )" << param2_col;
2998
2999 if (writeDetails)
3000 third_derivs_output << R"(, "var": ")" << symbol_table.getName(var) << R"(")"
3001 << R"(, "param1": ")" << symbol_table.getName(getSymbIDByDerivID(param1)) << R"(")"
3002 << R"(, "param2": ")" << symbol_table.getName(getSymbIDByDerivID(param2)) << R"(")";
3003
3004 third_derivs_output << R"(, "val": ")";
3005 d2->writeJsonOutput(third_derivs_output, temp_term_union, tef_terms);
3006 third_derivs_output << R"("})" << endl;
3007 }
3008 third_derivs_output << "]}" << endl;
3009
3010 third_derivs1_output << R"("derivative_hessian_wrt_params": {)"
3011 << R"( "neqs": )" << equations.size()
3012 << R"(, "nvar1cols": )" << symbol_table.endo_nbr()
3013 << R"(, "nvar2cols": )" << symbol_table.endo_nbr()
3014 << R"(, "nparamcols": )" << symbol_table.param_nbr()
3015 << R"(, "entries": [)";
3016 auto &hp = params_derivatives.find({ 2, 1 })->second;
3017 for (auto it = hp.begin(); it != hp.end(); ++it)
3018 {
3019 if (it != hp.begin())
3020 third_derivs1_output << ", ";
3021
3022 auto [eq, var1, var2, param] = vectorToTuple<4>(it->first);
3023 expr_t d2 = it->second;
3024
3025 int var1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var1)) + 1;
3026 int var2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var2)) + 1;
3027 int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
3028
3029 if (writeDetails)
3030 third_derivs1_output << R"({"eq": )" << eq + 1;
3031 else
3032 third_derivs1_output << R"({"row": )" << eq + 1;
3033
3034 third_derivs1_output << R"(, "var1_col": )" << var1_col
3035 << R"(, "var2_col": )" << var2_col
3036 << R"(, "param_col": )" << param_col;
3037
3038 if (writeDetails)
3039 third_derivs1_output << R"(, "var1": ")" << symbol_table.getName(getSymbIDByDerivID(var1)) << R"(")"
3040 << R"(, "var2": ")" << symbol_table.getName(getSymbIDByDerivID(var2)) << R"(")"
3041 << R"(, "param1": ")" << symbol_table.getName(getSymbIDByDerivID(param)) << R"(")";
3042
3043 third_derivs1_output << R"(, "val": ")";
3044 d2->writeJsonOutput(third_derivs1_output, temp_term_union, tef_terms);
3045 third_derivs1_output << R"("})" << endl;
3046 }
3047 third_derivs1_output << "]}" << endl;
3048
3049 if (writeDetails)
3050 output << R"("static_model_params_derivative": {)";
3051 else
3052 output << R"("static_model_params_derivatives_simple": {)";
3053 output << model_local_vars_output.str()
3054 << ", " << model_output.str()
3055 << ", " << jacobian_output.str()
3056 << ", " << hessian_output.str()
3057 << ", " << hessian1_output.str()
3058 << ", " << third_derivs_output.str()
3059 << ", " << third_derivs1_output.str()
3060 << "}";
3061 }
3062