2009-04-14 16:39:53 +02:00
|
|
|
/*
|
2022-04-19 17:06:37 +02:00
|
|
|
* Copyright © 2003-2022 Dynare Team
|
2009-04-14 16:39:53 +02:00
|
|
|
*
|
|
|
|
* This file is part of Dynare.
|
|
|
|
*
|
|
|
|
* Dynare is free software: you can redistribute it and/or modify
|
|
|
|
* it under the terms of the GNU General Public License as published by
|
|
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
|
|
* (at your option) any later version.
|
|
|
|
*
|
|
|
|
* Dynare is distributed in the hope that it will be useful,
|
|
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
* GNU General Public License for more details.
|
|
|
|
*
|
|
|
|
* You should have received a copy of the GNU General Public License
|
2021-06-09 16:52:20 +02:00
|
|
|
* along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
2009-04-14 16:39:53 +02:00
|
|
|
*/
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
#include <iostream>
|
|
|
|
#include <cmath>
|
2009-04-14 16:47:57 +02:00
|
|
|
#include <cstdlib>
|
2009-04-27 19:15:14 +02:00
|
|
|
#include <cassert>
|
2009-04-28 19:11:48 +02:00
|
|
|
#include <algorithm>
|
2021-04-23 17:30:38 +02:00
|
|
|
#include <sstream>
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2018-06-27 15:01:31 +02:00
|
|
|
#include "StaticModel.hh"
|
2018-10-10 17:08:54 +02:00
|
|
|
#include "DynamicModel.hh"
|
2009-04-28 19:11:48 +02:00
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
StaticModel::StaticModel(SymbolTable &symbol_table_arg,
|
2010-02-22 17:33:38 +01:00
|
|
|
NumericalConstants &num_constants_arg,
|
2018-09-14 17:04:06 +02:00
|
|
|
ExternalFunctionsTable &external_functions_table_arg) :
|
2018-10-04 17:18:27 +02:00
|
|
|
ModelTree{symbol_table_arg, num_constants_arg, external_functions_table_arg}
|
2009-12-16 14:21:31 +01:00
|
|
|
{
|
|
|
|
}
|
2009-04-14 16:47:57 +02:00
|
|
|
|
2018-10-09 18:27:19 +02:00
|
|
|
StaticModel::StaticModel(const StaticModel &m) :
|
2020-05-06 17:13:47 +02:00
|
|
|
ModelTree{m}
|
2018-10-09 18:27:19 +02:00
|
|
|
{
|
|
|
|
}
|
|
|
|
|
|
|
|
StaticModel &
|
|
|
|
StaticModel::operator=(const StaticModel &m)
|
|
|
|
{
|
|
|
|
ModelTree::operator=(m);
|
|
|
|
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
2018-10-10 17:08:54 +02:00
|
|
|
StaticModel::StaticModel(const DynamicModel &m) :
|
2019-12-20 16:59:30 +01:00
|
|
|
ModelTree{m.symbol_table, m.num_constants, m.external_functions_table}
|
2018-10-10 17:08:54 +02:00
|
|
|
{
|
|
|
|
// Convert model local variables (need to be done first)
|
|
|
|
for (int it : m.local_variables_vector)
|
|
|
|
AddLocalVariable(it, m.local_variables_table.find(it)->second->toStatic(*this));
|
|
|
|
|
|
|
|
// Convert equations
|
|
|
|
int static_only_index = 0;
|
2020-02-20 15:29:10 +01:00
|
|
|
set<int> dynamic_equations = m.equation_tags.getDynamicEqns();
|
2018-10-10 17:08:54 +02:00
|
|
|
for (int i = 0; i < static_cast<int>(m.equations.size()); i++)
|
2020-02-20 15:29:10 +01:00
|
|
|
try
|
|
|
|
{
|
|
|
|
// If equation is dynamic, replace it by an equation marked [static]
|
2022-05-04 16:01:34 +02:00
|
|
|
if (dynamic_equations.contains(i))
|
2018-10-10 17:08:54 +02:00
|
|
|
{
|
2020-02-20 15:29:10 +01:00
|
|
|
auto [static_only_equations,
|
|
|
|
static_only_equations_lineno,
|
|
|
|
static_only_equations_equation_tags] = m.getStaticOnlyEquationsInfo();
|
|
|
|
|
|
|
|
addEquation(static_only_equations[static_only_index]->toStatic(*this),
|
|
|
|
static_only_equations_lineno[static_only_index],
|
|
|
|
static_only_equations_equation_tags.getTagsByEqn(static_only_index));
|
|
|
|
static_only_index++;
|
2018-10-10 17:08:54 +02:00
|
|
|
}
|
2020-02-20 15:29:10 +01:00
|
|
|
else
|
|
|
|
addEquation(m.equations[i]->toStatic(*this),
|
|
|
|
m.equations_lineno[i],
|
|
|
|
m.equation_tags.getTagsByEqn(i));
|
|
|
|
}
|
|
|
|
catch (DataTree::DivisionByZeroException)
|
|
|
|
{
|
2021-12-10 12:44:42 +01:00
|
|
|
cerr << "...division by zero error encountered when converting equation " << i << " to static" << endl;
|
2020-02-20 15:29:10 +01:00
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
2018-10-10 17:08:54 +02:00
|
|
|
|
|
|
|
// Convert auxiliary equations
|
|
|
|
for (auto aux_eq : m.aux_equations)
|
|
|
|
addAuxEquation(aux_eq->toStatic(*this));
|
2019-12-02 19:21:14 +01:00
|
|
|
|
|
|
|
user_set_add_flags = m.user_set_add_flags;
|
|
|
|
user_set_subst_flags = m.user_set_subst_flags;
|
|
|
|
user_set_add_libs = m.user_set_add_libs;
|
|
|
|
user_set_subst_libs = m.user_set_subst_libs;
|
|
|
|
user_set_compiler = m.user_set_compiler;
|
2018-10-10 17:08:54 +02:00
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
void
|
2022-06-23 14:28:13 +02:00
|
|
|
StaticModel::writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
|
2009-12-16 14:21:31 +01:00
|
|
|
{
|
2019-12-16 19:42:59 +01:00
|
|
|
if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) });
|
|
|
|
it != derivatives[1].end())
|
2022-07-06 16:44:51 +02:00
|
|
|
it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms, temporary_terms_idxs, tef_terms);
|
2009-12-16 14:21:31 +01:00
|
|
|
else
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FLDZ_{};
|
2009-12-16 14:21:31 +01:00
|
|
|
}
|
2009-04-14 16:39:53 +02:00
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
void
|
2022-06-23 14:28:13 +02:00
|
|
|
StaticModel::writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
|
2009-12-16 14:21:31 +01:00
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
if (auto it = blocks_derivatives[blk].find({ eq, var, lag });
|
|
|
|
it != blocks_derivatives[blk].end())
|
2022-07-06 16:44:51 +02:00
|
|
|
it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms, temporary_terms_idxs, tef_terms);
|
2009-12-16 14:21:31 +01:00
|
|
|
else
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FLDZ_{};
|
2009-12-16 14:21:31 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2020-05-19 16:56:53 +02:00
|
|
|
StaticModel::writeStaticPerBlockMFiles(const string &basename) const
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-05-06 17:13:47 +02:00
|
|
|
temporary_terms_t temporary_terms; // Temp terms written so far
|
2020-05-05 17:49:58 +02:00
|
|
|
|
2020-05-19 16:56:53 +02:00
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-05-19 16:56:53 +02:00
|
|
|
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
2009-12-16 18:13:23 +01:00
|
|
|
|
2020-05-19 16:56:53 +02:00
|
|
|
string filename = packageDir(basename + ".block") + "/static_" + to_string(blk+1) + ".m";
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream output{filename, ios::out | ios::binary};
|
2020-06-23 15:13:04 +02:00
|
|
|
if (!output.is_open())
|
|
|
|
{
|
|
|
|
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
2018-11-23 17:19:59 +01:00
|
|
|
output << "%" << endl
|
2020-05-19 16:56:53 +02:00
|
|
|
<< "% " << filename << " : Computes static version of one block" << endl
|
2018-11-23 17:19:59 +01:00
|
|
|
<< "%" << endl
|
|
|
|
<< "% Warning : this file is generated automatically by Dynare" << endl
|
|
|
|
<< "% from model file (.mod)" << endl << endl
|
2020-05-19 16:56:53 +02:00
|
|
|
<< "%" << endl;
|
2020-03-20 17:31:14 +01:00
|
|
|
if (simulation_type == BlockSimulationType::evaluateBackward
|
|
|
|
|| simulation_type == BlockSimulationType::evaluateForward)
|
2020-05-25 18:35:36 +02:00
|
|
|
output << "function [y, T] = static_" << blk+1 << "(y, x, params, T)" << endl;
|
2009-12-16 18:13:23 +01:00
|
|
|
else
|
2020-10-06 17:16:50 +02:00
|
|
|
output << "function [residual, y, T, g1] = static_" << blk+1 << "(y, x, params, T)" << endl;
|
2009-12-16 18:13:23 +01:00
|
|
|
|
|
|
|
output << " % ////////////////////////////////////////////////////////////////////////" << endl
|
2022-06-22 12:47:11 +02:00
|
|
|
<< " % //" << " Block "s.substr(static_cast<int>(log10(blk + 1))) << blk+1
|
2020-04-30 12:48:16 +02:00
|
|
|
<< " //" << endl
|
2009-12-16 18:13:23 +01:00
|
|
|
<< " % // Simulation type "
|
|
|
|
<< BlockSim(simulation_type) << " //" << endl
|
2020-05-19 16:56:53 +02:00
|
|
|
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
|
|
|
|
|
2020-03-20 17:31:14 +01:00
|
|
|
if (simulation_type != BlockSimulationType::evaluateBackward
|
|
|
|
&& simulation_type != BlockSimulationType::evaluateForward)
|
2020-06-22 12:46:33 +02:00
|
|
|
output << " residual=zeros(" << blocks[blk].mfs_size << ",1);" << endl
|
|
|
|
<< " g1_i=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
|
|
|
|
<< " g1_j=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
|
|
|
|
<< " g1_v=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
|
|
|
|
<< endl;
|
2020-05-13 16:58:19 +02:00
|
|
|
|
2022-07-12 14:13:27 +02:00
|
|
|
writeStaticPerBlockHelper<ExprNodeOutputType::matlabStaticModel>(blk, output, temporary_terms);
|
2020-05-13 16:58:19 +02:00
|
|
|
|
2020-06-22 12:46:33 +02:00
|
|
|
if (simulation_type != BlockSimulationType::evaluateBackward
|
|
|
|
&& simulation_type != BlockSimulationType::evaluateForward)
|
|
|
|
output << endl
|
|
|
|
<< " g1=sparse(g1_i, g1_j, g1_v, " << blocks[blk].mfs_size << "," << blocks[blk].mfs_size << ");" << endl;
|
2020-06-19 16:09:38 +02:00
|
|
|
|
2010-12-17 14:21:35 +01:00
|
|
|
output << "end" << endl;
|
2009-12-16 18:13:23 +01:00
|
|
|
output.close();
|
|
|
|
}
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-06-23 15:13:04 +02:00
|
|
|
void
|
|
|
|
StaticModel::writeStaticPerBlockCFiles(const string &basename) const
|
|
|
|
{
|
|
|
|
temporary_terms_t temporary_terms; // Temp terms written so far
|
|
|
|
|
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
|
|
|
{
|
|
|
|
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
|
|
|
|
|
|
|
string filename = basename + "/model/src/static_" + to_string(blk+1) + ".c";
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream output{filename, ios::out | ios::binary};
|
2020-06-23 15:13:04 +02:00
|
|
|
if (!output.is_open())
|
|
|
|
{
|
|
|
|
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
output << "/* Block " << blk+1 << endl
|
|
|
|
<< " " << BlockSim(simulation_type) << " */" << endl
|
|
|
|
<< endl
|
|
|
|
<< "#include <math.h>" << endl
|
|
|
|
<< "#include <stdlib.h>" << endl
|
|
|
|
<< R"(#include "mex.h")" << endl
|
|
|
|
<< endl;
|
|
|
|
|
|
|
|
// Write function definition if BinaryOpcode::powerDeriv is used
|
|
|
|
writePowerDerivHeader(output);
|
|
|
|
|
|
|
|
output << endl;
|
|
|
|
|
|
|
|
if (simulation_type == BlockSimulationType::evaluateBackward
|
|
|
|
|| simulation_type == BlockSimulationType::evaluateForward)
|
2020-09-03 17:52:12 +02:00
|
|
|
output << "void static_" << blk+1 << "(double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl;
|
2020-06-23 15:13:04 +02:00
|
|
|
else
|
2020-10-06 17:16:50 +02:00
|
|
|
output << "void static_" << blk+1 << "(double *restrict y, const double *restrict x, const double *restrict params, double *restrict T, double *restrict residual, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v)" << endl;
|
2020-06-23 15:13:04 +02:00
|
|
|
output << '{' << endl;
|
|
|
|
|
2022-07-12 14:13:27 +02:00
|
|
|
writeStaticPerBlockHelper<ExprNodeOutputType::CStaticModel>(blk, output, temporary_terms);
|
2020-06-23 15:13:04 +02:00
|
|
|
|
|
|
|
output << '}' << endl
|
|
|
|
<< endl;
|
|
|
|
|
|
|
|
ostringstream header;
|
|
|
|
if (simulation_type == BlockSimulationType::evaluateBackward
|
|
|
|
|| simulation_type == BlockSimulationType::evaluateForward)
|
|
|
|
{
|
|
|
|
header << "void static_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, mxArray *T)";
|
|
|
|
output << header.str() << endl
|
|
|
|
<< '{' << endl
|
|
|
|
<< " static_" << blk+1 << "(mxGetPr(y), mxGetPr(x), mxGetPr(params), mxGetPr(T));" << endl
|
|
|
|
<< '}' << endl;
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
2020-10-06 17:16:50 +02:00
|
|
|
header << "void static_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, mxArray *T, mxArray **residual, mxArray **g1)";
|
2020-06-23 15:13:04 +02:00
|
|
|
output << header.str() << endl
|
|
|
|
<< '{' << endl
|
|
|
|
<< " *residual = mxCreateDoubleMatrix(" << blocks[blk].mfs_size << ",1,mxREAL);" << endl
|
|
|
|
<< " mxArray *g1_i = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
|
|
|
|
<< " mxArray *g1_j = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
|
|
|
|
<< " mxArray *g1_v = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
|
|
|
|
<< " static_" << blk+1 << "(mxGetPr(y), mxGetPr(x), mxGetPr(params), mxGetPr(T), mxGetPr(*residual), mxGetPr(g1_i), mxGetPr(g1_j), mxGetPr(g1_v));" << endl
|
|
|
|
<< " mxArray *plhs[1];" << endl
|
|
|
|
<< " mxArray *m = mxCreateDoubleScalar(" << blocks[blk].mfs_size << ");" << endl
|
|
|
|
<< " mxArray *n = mxCreateDoubleScalar(" << blocks[blk].mfs_size << ");" << endl
|
|
|
|
<< " mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
|
|
|
|
<< R"( mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
|
|
|
|
<< " *g1 = plhs[0];" << endl
|
|
|
|
<< " mxDestroyArray(g1_i);" << endl
|
|
|
|
<< " mxDestroyArray(g1_j);" << endl
|
|
|
|
<< " mxDestroyArray(g1_v);" << endl
|
|
|
|
<< " mxDestroyArray(m);" << endl
|
|
|
|
<< " mxDestroyArray(n);" << endl
|
|
|
|
<< '}' << endl;
|
|
|
|
}
|
|
|
|
|
|
|
|
output.close();
|
|
|
|
|
|
|
|
filename = basename + "/model/src/static_" + to_string(blk+1) + ".h";
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream header_output{filename, ios::out | ios::binary};
|
2020-06-23 15:13:04 +02:00
|
|
|
if (!header_output.is_open())
|
|
|
|
{
|
|
|
|
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
header_output << header.str() << ';' << endl;
|
|
|
|
header_output.close();
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
void
|
2020-05-19 17:43:35 +02:00
|
|
|
StaticModel::writeStaticBytecode(const string &basename) const
|
2010-01-22 11:03:29 +01:00
|
|
|
{
|
|
|
|
ostringstream tmp_output;
|
|
|
|
bool file_open = false;
|
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
BytecodeWriter code_file{basename + "/model/bytecode/static.cod"};
|
2010-01-22 11:03:29 +01:00
|
|
|
int count_u;
|
|
|
|
int u_count_int = 0;
|
|
|
|
|
2020-05-19 17:43:35 +02:00
|
|
|
writeBytecodeBinFile(basename + "/model/bytecode/static.bin", u_count_int, file_open, false);
|
2010-01-22 11:03:29 +01:00
|
|
|
file_open = true;
|
|
|
|
|
2020-05-06 17:13:47 +02:00
|
|
|
// Compute the union of temporary terms from residuals and 1st derivatives
|
|
|
|
temporary_terms_t temporary_terms = temporary_terms_derivatives[0];
|
2022-06-08 13:59:10 +02:00
|
|
|
temporary_terms.insert(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end());
|
2020-05-06 17:13:47 +02:00
|
|
|
|
2010-01-22 11:03:29 +01:00
|
|
|
//Temporary variables declaration
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FDIMST_{static_cast<int>(temporary_terms.size())}
|
|
|
|
<< FBEGINBLOCK_{symbol_table.endo_nbr(),
|
|
|
|
BlockSimulationType::solveForwardComplete,
|
|
|
|
0,
|
|
|
|
symbol_table.endo_nbr(),
|
|
|
|
endo_idx_block2orig,
|
|
|
|
eq_idx_block2orig,
|
|
|
|
false,
|
|
|
|
symbol_table.endo_nbr(),
|
|
|
|
0,
|
|
|
|
0,
|
|
|
|
u_count_int,
|
|
|
|
symbol_table.endo_nbr()};
|
2010-01-22 11:03:29 +01:00
|
|
|
|
2020-05-13 12:53:47 +02:00
|
|
|
temporary_terms_t temporary_terms_union;
|
2022-05-20 11:43:02 +02:00
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
2022-07-06 16:44:51 +02:00
|
|
|
writeBytecodeTemporaryTerms(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
|
2010-01-22 11:03:29 +01:00
|
|
|
|
2022-07-06 16:44:51 +02:00
|
|
|
writeBytecodeModelEquations(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
|
2010-01-22 11:03:29 +01:00
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FENDEQU_{};
|
2010-01-22 11:03:29 +01:00
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
// Get the current code_file position and jump if evaluating
|
|
|
|
int pos_jmpifeval = code_file.getInstructionCounter();
|
|
|
|
code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being
|
2010-10-18 17:28:21 +02:00
|
|
|
|
2018-11-23 17:19:59 +01:00
|
|
|
vector<vector<pair<int, int>>> my_derivatives(symbol_table.endo_nbr());
|
2010-01-22 11:03:29 +01:00
|
|
|
count_u = symbol_table.endo_nbr();
|
2019-12-16 19:42:59 +01:00
|
|
|
for (const auto & [indices, d1] : derivatives[1])
|
2010-01-22 11:03:29 +01:00
|
|
|
{
|
2019-12-16 19:42:59 +01:00
|
|
|
int deriv_id = indices[1];
|
2018-07-17 18:34:07 +02:00
|
|
|
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
|
2010-01-22 11:03:29 +01:00
|
|
|
{
|
2020-03-24 18:26:06 +01:00
|
|
|
int eq = indices[0];
|
2022-07-12 15:28:52 +02:00
|
|
|
int var { getTypeSpecificIDByDerivID(deriv_id) };
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, var};
|
2018-11-15 16:39:53 +01:00
|
|
|
if (!my_derivatives[eq].size())
|
|
|
|
my_derivatives[eq].clear();
|
|
|
|
my_derivatives[eq].emplace_back(var, count_u);
|
2010-01-22 11:03:29 +01:00
|
|
|
|
2022-07-06 16:44:51 +02:00
|
|
|
d1->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
|
2010-01-22 11:03:29 +01:00
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FSTPSU_{count_u};
|
2010-01-22 11:03:29 +01:00
|
|
|
count_u++;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
for (int i = 0; i < symbol_table.endo_nbr(); i++)
|
|
|
|
{
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FLDR_{i};
|
2018-11-15 16:39:53 +01:00
|
|
|
if (my_derivatives[i].size())
|
2010-01-22 11:03:29 +01:00
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
for (bool printed_something{false};
|
|
|
|
const auto &it : my_derivatives[i])
|
2010-01-22 11:03:29 +01:00
|
|
|
{
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FLDSU_{it.second}
|
|
|
|
<< FLDSV_{SymbolType::endogenous, it.first}
|
|
|
|
<< FBINARY_{BinaryOpcode::times};
|
2022-06-03 16:24:26 +02:00
|
|
|
if (exchange(printed_something, true))
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FBINARY_{BinaryOpcode::plus};
|
2010-01-22 11:03:29 +01:00
|
|
|
}
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FBINARY_{BinaryOpcode::minus};
|
2010-01-22 11:03:29 +01:00
|
|
|
}
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FSTPSU_{i};
|
2010-01-22 11:03:29 +01:00
|
|
|
}
|
2022-06-23 14:28:13 +02:00
|
|
|
|
|
|
|
// Jump unconditionally after the block
|
|
|
|
int pos_jmp = code_file.getInstructionCounter();
|
|
|
|
code_file << FJMP_{0}; // Use 0 as jump offset for the time being
|
|
|
|
// Update jump offset for previous JMPIFEVAL
|
|
|
|
code_file.overwriteInstruction(pos_jmpifeval, FJMPIFEVAL_{pos_jmp-pos_jmpifeval});
|
2010-10-18 17:28:21 +02:00
|
|
|
|
2018-11-23 17:19:59 +01:00
|
|
|
temporary_terms_t tt2, tt3;
|
2010-10-18 17:28:21 +02:00
|
|
|
|
|
|
|
// The Jacobian if we have to solve the block determinsitic bloc
|
2019-12-16 19:42:59 +01:00
|
|
|
for (const auto & [indices, d1] : derivatives[1])
|
2010-10-18 17:28:21 +02:00
|
|
|
{
|
2019-12-16 19:42:59 +01:00
|
|
|
int deriv_id = indices[1];
|
2018-07-17 18:34:07 +02:00
|
|
|
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
|
2010-10-18 17:28:21 +02:00
|
|
|
{
|
2020-03-24 18:26:06 +01:00
|
|
|
int eq = indices[0];
|
2022-07-12 15:28:52 +02:00
|
|
|
int var { getTypeSpecificIDByDerivID(deriv_id) };
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, var};
|
2018-11-15 16:39:53 +01:00
|
|
|
if (!my_derivatives[eq].size())
|
|
|
|
my_derivatives[eq].clear();
|
|
|
|
my_derivatives[eq].emplace_back(var, count_u);
|
2010-10-18 17:28:21 +02:00
|
|
|
|
2022-07-06 16:44:51 +02:00
|
|
|
d1->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FSTPG2_{eq, var};
|
2010-10-18 17:28:21 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
// Update jump offset for previous JMP
|
|
|
|
int pos_end_block = code_file.getInstructionCounter();
|
|
|
|
code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1});
|
|
|
|
|
|
|
|
code_file << FENDBLOCK_{} << FEND_{};
|
2010-01-22 11:03:29 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2020-05-19 17:43:35 +02:00
|
|
|
StaticModel::writeStaticBlockBytecode(const string &basename) const
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
|
|
|
struct Uff_l
|
2009-12-16 14:21:31 +01:00
|
|
|
{
|
2009-12-16 18:13:23 +01:00
|
|
|
int u, var, lag;
|
|
|
|
Uff_l *pNext;
|
|
|
|
};
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2009-12-16 18:13:23 +01:00
|
|
|
struct Uff
|
|
|
|
{
|
|
|
|
Uff_l *Ufl, *Ufl_First;
|
|
|
|
};
|
|
|
|
|
|
|
|
int i, v;
|
|
|
|
string tmp_s;
|
|
|
|
ostringstream tmp_output;
|
2018-06-04 12:52:14 +02:00
|
|
|
expr_t lhs = nullptr, rhs = nullptr;
|
2009-12-16 18:13:23 +01:00
|
|
|
BinaryOpNode *eq_node;
|
|
|
|
Uff Uf[symbol_table.endo_nbr()];
|
2010-09-16 19:18:45 +02:00
|
|
|
map<expr_t, int> reference_count;
|
2009-12-16 18:13:23 +01:00
|
|
|
vector<int> feedback_variables;
|
|
|
|
bool file_open = false;
|
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
BytecodeWriter code_file{basename + "/model/bytecode/static.cod"};
|
2009-12-16 18:13:23 +01:00
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
//Temporary variables declaration
|
|
|
|
code_file << FDIMST_{static_cast<int>(blocks_temporary_terms_idxs.size())};
|
2009-12-16 18:13:23 +01:00
|
|
|
|
2020-05-13 12:53:47 +02:00
|
|
|
temporary_terms_t temporary_terms_union;
|
|
|
|
|
2020-04-23 16:04:02 +02:00
|
|
|
for (int block = 0; block < static_cast<int>(blocks.size()); block++)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
|
|
|
feedback_variables.clear();
|
|
|
|
if (block > 0)
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FENDBLOCK_{};
|
2009-12-16 18:13:23 +01:00
|
|
|
int count_u;
|
|
|
|
int u_count_int = 0;
|
2020-04-23 16:04:02 +02:00
|
|
|
BlockSimulationType simulation_type = blocks[block].simulation_type;
|
|
|
|
int block_size = blocks[block].size;
|
|
|
|
int block_mfs = blocks[block].mfs_size;
|
|
|
|
int block_recursive = blocks[block].getRecursiveSize();
|
2009-12-16 18:13:23 +01:00
|
|
|
|
2020-03-20 17:31:14 +01:00
|
|
|
if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple
|
|
|
|
|| simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|
|
|
|
|| simulation_type == BlockSimulationType::solveBackwardComplete
|
|
|
|
|| simulation_type == BlockSimulationType::solveForwardComplete)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-05-19 17:43:35 +02:00
|
|
|
writeBlockBytecodeBinFile(basename, block, u_count_int, file_open);
|
2009-12-16 18:13:23 +01:00
|
|
|
file_open = true;
|
|
|
|
}
|
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FBEGINBLOCK_{block_mfs,
|
|
|
|
simulation_type,
|
|
|
|
blocks[block].first_equation,
|
|
|
|
block_size,
|
|
|
|
endo_idx_block2orig,
|
|
|
|
eq_idx_block2orig,
|
|
|
|
blocks[block].linear,
|
|
|
|
symbol_table.endo_nbr(),
|
|
|
|
0,
|
|
|
|
0,
|
|
|
|
u_count_int,
|
|
|
|
block_size};
|
|
|
|
|
|
|
|
// Get the current code_file position and jump if evaluating
|
|
|
|
int pos_jmpifeval = code_file.getInstructionCounter();
|
|
|
|
code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being
|
2010-12-31 16:46:09 +01:00
|
|
|
|
2020-05-06 17:13:47 +02:00
|
|
|
//The Temporary terms
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
2020-05-13 16:58:19 +02:00
|
|
|
/* Keep a backup of temporary_terms_union here, since temp. terms are
|
2020-05-13 12:53:47 +02:00
|
|
|
written a second time below. This is probably unwanted… */
|
2020-05-13 16:58:19 +02:00
|
|
|
temporary_terms_t ttu_old = temporary_terms_union;
|
|
|
|
|
|
|
|
auto write_eq_tt = [&](int eq)
|
|
|
|
{
|
|
|
|
for (auto it : blocks_temporary_terms[block][eq])
|
|
|
|
{
|
|
|
|
if (dynamic_cast<AbstractExternalFunctionNode *>(it))
|
2022-07-06 16:44:51 +02:00
|
|
|
it->writeBytecodeExternalFunctionOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
2020-05-13 16:58:19 +02:00
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)};
|
2022-07-06 16:44:51 +02:00
|
|
|
it->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FSTPST_{blocks_temporary_terms_idxs.at(it)};
|
2020-05-13 16:58:19 +02:00
|
|
|
temporary_terms_union.insert(it);
|
|
|
|
}
|
|
|
|
};
|
2020-05-06 17:13:47 +02:00
|
|
|
|
|
|
|
for (i = 0; i < block_size; i++)
|
|
|
|
{
|
2020-05-13 16:58:19 +02:00
|
|
|
write_eq_tt(i);
|
|
|
|
|
2010-10-11 19:21:32 +02:00
|
|
|
// The equations
|
2009-12-16 18:13:23 +01:00
|
|
|
int variable_ID, equation_ID;
|
|
|
|
EquationType equ_type;
|
|
|
|
switch (simulation_type)
|
|
|
|
{
|
|
|
|
evaluation:
|
2020-03-20 17:31:14 +01:00
|
|
|
case BlockSimulationType::evaluateBackward:
|
|
|
|
case BlockSimulationType::evaluateForward:
|
2009-12-16 18:13:23 +01:00
|
|
|
equ_type = getBlockEquationType(block, i);
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
|
2020-03-20 18:00:56 +01:00
|
|
|
if (equ_type == EquationType::evaluate)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-05-07 15:20:11 +02:00
|
|
|
eq_node = getBlockEquationExpr(block, i);
|
2018-11-28 14:32:26 +01:00
|
|
|
lhs = eq_node->arg1;
|
|
|
|
rhs = eq_node->arg2;
|
2022-07-06 16:44:51 +02:00
|
|
|
rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
2009-12-16 18:13:23 +01:00
|
|
|
}
|
2020-05-20 11:49:32 +02:00
|
|
|
else if (equ_type == EquationType::evaluateRenormalized)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-05-07 15:20:11 +02:00
|
|
|
eq_node = getBlockEquationRenormalizedExpr(block, i);
|
2018-11-28 14:32:26 +01:00
|
|
|
lhs = eq_node->arg1;
|
|
|
|
rhs = eq_node->arg2;
|
2022-07-06 16:44:51 +02:00
|
|
|
rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
2009-12-16 18:13:23 +01:00
|
|
|
}
|
|
|
|
break;
|
2020-03-20 17:31:14 +01:00
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
2020-03-24 18:26:06 +01:00
|
|
|
if (i < block_recursive)
|
2009-12-16 18:13:23 +01:00
|
|
|
goto evaluation;
|
|
|
|
variable_ID = getBlockVariableID(block, i);
|
|
|
|
equation_ID = getBlockEquationID(block, i);
|
|
|
|
feedback_variables.push_back(variable_ID);
|
2018-06-04 12:52:14 +02:00
|
|
|
Uf[equation_ID].Ufl = nullptr;
|
2009-12-16 18:13:23 +01:00
|
|
|
goto end;
|
|
|
|
default:
|
|
|
|
end:
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
|
2020-05-07 15:20:11 +02:00
|
|
|
eq_node = getBlockEquationExpr(block, i);
|
2018-11-28 14:32:26 +01:00
|
|
|
lhs = eq_node->arg1;
|
|
|
|
rhs = eq_node->arg2;
|
2022-07-06 16:44:51 +02:00
|
|
|
lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
2009-12-16 18:13:23 +01:00
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
|
2009-12-16 18:13:23 +01:00
|
|
|
}
|
|
|
|
}
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FENDEQU_{};
|
2010-12-10 11:50:27 +01:00
|
|
|
|
2009-12-16 18:13:23 +01:00
|
|
|
// The Jacobian if we have to solve the block
|
2020-03-20 17:31:14 +01:00
|
|
|
if (simulation_type != BlockSimulationType::evaluateBackward
|
|
|
|
&& simulation_type != BlockSimulationType::evaluateForward)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-05-13 16:58:19 +02:00
|
|
|
// Write temporary terms for derivatives
|
|
|
|
write_eq_tt(blocks[block].size);
|
|
|
|
|
2009-12-16 18:13:23 +01:00
|
|
|
switch (simulation_type)
|
|
|
|
{
|
2020-03-20 17:31:14 +01:00
|
|
|
case BlockSimulationType::solveBackwardSimple:
|
|
|
|
case BlockSimulationType::solveForwardSimple:
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, 0, 0};
|
|
|
|
writeBytecodeDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
code_file << FSTPG_{0};
|
2009-12-16 18:13:23 +01:00
|
|
|
break;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-03-20 17:31:14 +01:00
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
2009-12-16 18:13:23 +01:00
|
|
|
count_u = feedback_variables.size();
|
2020-04-24 12:29:02 +02:00
|
|
|
for (const auto &[indices, ignore2] : blocks_derivatives[block])
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
auto [eq, var, ignore] = indices;
|
2020-03-24 18:26:06 +01:00
|
|
|
int eqr = getBlockEquationID(block, eq);
|
|
|
|
int varr = getBlockVariableID(block, var);
|
2010-04-28 16:03:32 +02:00
|
|
|
if (eq >= block_recursive && var >= block_recursive)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
|
|
|
if (!Uf[eqr].Ufl)
|
|
|
|
{
|
2019-04-23 11:07:32 +02:00
|
|
|
Uf[eqr].Ufl = static_cast<Uff_l *>(malloc(sizeof(Uff_l)));
|
2009-12-16 18:13:23 +01:00
|
|
|
Uf[eqr].Ufl_First = Uf[eqr].Ufl;
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
2019-04-23 11:07:32 +02:00
|
|
|
Uf[eqr].Ufl->pNext = static_cast<Uff_l *>(malloc(sizeof(Uff_l)));
|
2009-12-16 18:13:23 +01:00
|
|
|
Uf[eqr].Ufl = Uf[eqr].Ufl->pNext;
|
|
|
|
}
|
2018-06-04 12:52:14 +02:00
|
|
|
Uf[eqr].Ufl->pNext = nullptr;
|
2009-12-16 18:13:23 +01:00
|
|
|
Uf[eqr].Ufl->u = count_u;
|
|
|
|
Uf[eqr].Ufl->var = varr;
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr};
|
|
|
|
writeBytecodeChainRuleDerivative(code_file, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
code_file << FSTPSU_{count_u};
|
2009-12-16 18:13:23 +01:00
|
|
|
count_u++;
|
|
|
|
}
|
|
|
|
}
|
2020-03-24 18:26:06 +01:00
|
|
|
for (i = 0; i < block_size; i++)
|
2022-06-02 10:50:21 +02:00
|
|
|
if (i >= block_recursive)
|
|
|
|
{
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FLDR_{i-block_recursive} << FLDZ_{};
|
2022-06-02 10:50:21 +02:00
|
|
|
|
|
|
|
v = getBlockEquationID(block, i);
|
|
|
|
for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FLDSU_{Uf[v].Ufl->u}
|
|
|
|
<< FLDSV_{SymbolType::endogenous, Uf[v].Ufl->var}
|
|
|
|
<< FBINARY_{BinaryOpcode::times}
|
|
|
|
<< FCUML_{};
|
2022-06-02 10:50:21 +02:00
|
|
|
Uf[v].Ufl = Uf[v].Ufl_First;
|
|
|
|
while (Uf[v].Ufl)
|
|
|
|
{
|
|
|
|
Uf[v].Ufl_First = Uf[v].Ufl->pNext;
|
|
|
|
free(Uf[v].Ufl);
|
|
|
|
Uf[v].Ufl = Uf[v].Ufl_First;
|
|
|
|
}
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FBINARY_{BinaryOpcode::minus}
|
|
|
|
<< FSTPSU_{i - block_recursive};
|
2022-06-02 10:50:21 +02:00
|
|
|
}
|
2009-12-16 18:13:23 +01:00
|
|
|
break;
|
|
|
|
default:
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
2010-10-11 19:21:32 +02:00
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
// Jump unconditionally after the block
|
|
|
|
int pos_jmp = code_file.getInstructionCounter();
|
|
|
|
code_file << FJMP_{0}; // Use 0 as jump offset for the time being
|
|
|
|
// Update jump offset for previous JMPIFEVAL
|
|
|
|
code_file.overwriteInstruction(pos_jmpifeval, FJMPIFEVAL_{pos_jmp-pos_jmpifeval});
|
2010-10-11 19:21:32 +02:00
|
|
|
|
2020-05-06 17:13:47 +02:00
|
|
|
tef_terms.clear();
|
2020-05-13 16:58:19 +02:00
|
|
|
temporary_terms_union = ttu_old;
|
2010-10-11 19:21:32 +02:00
|
|
|
|
2020-05-06 17:13:47 +02:00
|
|
|
for (i = 0; i < block_size; i++)
|
|
|
|
{
|
2020-05-13 16:58:19 +02:00
|
|
|
write_eq_tt(i);
|
|
|
|
|
2010-10-11 19:21:32 +02:00
|
|
|
// The equations
|
|
|
|
int variable_ID, equation_ID;
|
|
|
|
EquationType equ_type;
|
|
|
|
switch (simulation_type)
|
|
|
|
{
|
|
|
|
evaluation_l:
|
2020-03-20 17:31:14 +01:00
|
|
|
case BlockSimulationType::evaluateBackward:
|
|
|
|
case BlockSimulationType::evaluateForward:
|
2010-10-11 19:21:32 +02:00
|
|
|
equ_type = getBlockEquationType(block, i);
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
|
2020-03-20 18:00:56 +01:00
|
|
|
if (equ_type == EquationType::evaluate)
|
2010-10-11 19:21:32 +02:00
|
|
|
{
|
2020-05-07 15:20:11 +02:00
|
|
|
eq_node = getBlockEquationExpr(block, i);
|
2018-11-28 14:32:26 +01:00
|
|
|
lhs = eq_node->arg1;
|
|
|
|
rhs = eq_node->arg2;
|
2022-07-06 16:44:51 +02:00
|
|
|
rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
2010-10-11 19:21:32 +02:00
|
|
|
}
|
2020-05-20 11:49:32 +02:00
|
|
|
else if (equ_type == EquationType::evaluateRenormalized)
|
2010-10-11 19:21:32 +02:00
|
|
|
{
|
2020-05-07 15:20:11 +02:00
|
|
|
eq_node = getBlockEquationRenormalizedExpr(block, i);
|
2018-11-28 14:32:26 +01:00
|
|
|
lhs = eq_node->arg1;
|
|
|
|
rhs = eq_node->arg2;
|
2022-07-06 16:44:51 +02:00
|
|
|
rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
2010-10-11 19:21:32 +02:00
|
|
|
}
|
|
|
|
break;
|
2020-03-20 17:31:14 +01:00
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
2020-03-24 18:26:06 +01:00
|
|
|
if (i < block_recursive)
|
2010-10-11 19:21:32 +02:00
|
|
|
goto evaluation_l;
|
|
|
|
variable_ID = getBlockVariableID(block, i);
|
|
|
|
equation_ID = getBlockEquationID(block, i);
|
|
|
|
feedback_variables.push_back(variable_ID);
|
2018-06-04 12:52:14 +02:00
|
|
|
Uf[equation_ID].Ufl = nullptr;
|
2010-10-11 19:21:32 +02:00
|
|
|
goto end_l;
|
|
|
|
default:
|
|
|
|
end_l:
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
|
2020-05-07 15:20:11 +02:00
|
|
|
eq_node = getBlockEquationExpr(block, i);
|
2018-11-28 14:32:26 +01:00
|
|
|
lhs = eq_node->arg1;
|
|
|
|
rhs = eq_node->arg2;
|
2022-07-06 16:44:51 +02:00
|
|
|
lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
2010-10-11 19:21:32 +02:00
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
|
2010-10-11 19:21:32 +02:00
|
|
|
}
|
|
|
|
}
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FENDEQU_{};
|
2010-10-11 19:21:32 +02:00
|
|
|
|
|
|
|
// The Jacobian if we have to solve the block determinsitic bloc
|
2020-05-13 16:58:19 +02:00
|
|
|
|
|
|
|
// Write temporary terms for derivatives
|
|
|
|
write_eq_tt(blocks[block].size);
|
|
|
|
|
2010-10-11 19:21:32 +02:00
|
|
|
switch (simulation_type)
|
|
|
|
{
|
2020-03-20 17:31:14 +01:00
|
|
|
case BlockSimulationType::solveBackwardSimple:
|
|
|
|
case BlockSimulationType::solveForwardSimple:
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, 0, 0};
|
|
|
|
writeBytecodeDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
code_file << FSTPG2_{0, 0};
|
2010-10-11 19:21:32 +02:00
|
|
|
break;
|
2020-03-20 17:31:14 +01:00
|
|
|
case BlockSimulationType::evaluateBackward:
|
|
|
|
case BlockSimulationType::evaluateForward:
|
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
2010-10-11 19:21:32 +02:00
|
|
|
count_u = feedback_variables.size();
|
2020-04-24 12:29:02 +02:00
|
|
|
for (const auto &[indices, ignore2] : blocks_derivatives[block])
|
2010-10-11 19:21:32 +02:00
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
auto &[eq, var, ignore] = indices;
|
2020-03-24 18:26:06 +01:00
|
|
|
int eqr = getBlockEquationID(block, eq);
|
|
|
|
int varr = getBlockVariableID(block, var);
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, 0};
|
|
|
|
writeBytecodeChainRuleDerivative(code_file, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
code_file << FSTPG2_{eq, var};
|
2010-10-11 19:21:32 +02:00
|
|
|
}
|
|
|
|
break;
|
|
|
|
default:
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
// Set codefile position to previous JMP_ and set the number of instructions to jump
|
2022-06-23 14:28:13 +02:00
|
|
|
// Update jump offset for previous JMP
|
|
|
|
int pos_end_block = code_file.getInstructionCounter();
|
|
|
|
code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1});
|
2009-12-16 18:13:23 +01:00
|
|
|
}
|
2022-06-23 14:28:13 +02:00
|
|
|
code_file << FENDBLOCK_{} << FEND_{};
|
2009-12-16 18:13:23 +01:00
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
void
|
2020-05-19 17:43:35 +02:00
|
|
|
StaticModel::writeBlockBytecodeBinFile(const string &basename, int num,
|
|
|
|
int &u_count_int, bool &file_open) const
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
|
|
|
int j;
|
|
|
|
std::ofstream SaveCode;
|
2018-06-27 15:01:31 +02:00
|
|
|
string filename = basename + "/model/bytecode/static.bin";
|
2009-12-16 18:13:23 +01:00
|
|
|
if (file_open)
|
2018-06-27 15:01:31 +02:00
|
|
|
SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate);
|
2009-12-16 18:13:23 +01:00
|
|
|
else
|
2018-06-27 15:01:31 +02:00
|
|
|
SaveCode.open(filename, ios::out | ios::binary);
|
2009-12-16 18:13:23 +01:00
|
|
|
if (!SaveCode.is_open())
|
|
|
|
{
|
2018-06-27 15:01:31 +02:00
|
|
|
cerr << "Error : Can't open file " << filename << " for writing" << endl;
|
2009-12-16 18:13:23 +01:00
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
u_count_int = 0;
|
2020-04-23 16:04:02 +02:00
|
|
|
int block_size = blocks[num].size;
|
|
|
|
int block_mfs = blocks[num].mfs_size;
|
|
|
|
int block_recursive = blocks[num].getRecursiveSize();
|
2020-04-24 12:29:02 +02:00
|
|
|
for (const auto &[indices, ignore2] : blocks_derivatives[num])
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
auto [eq, var, ignore] = indices;
|
2009-12-16 18:13:23 +01:00
|
|
|
int lag = 0;
|
2010-04-28 16:03:32 +02:00
|
|
|
if (eq >= block_recursive && var >= block_recursive)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
|
|
|
int v = eq - block_recursive;
|
|
|
|
SaveCode.write(reinterpret_cast<char *>(&v), sizeof(v));
|
|
|
|
int varr = var - block_recursive;
|
|
|
|
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
|
|
|
|
SaveCode.write(reinterpret_cast<char *>(&lag), sizeof(lag));
|
|
|
|
int u = u_count_int + block_mfs;
|
|
|
|
SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
|
|
|
|
u_count_int++;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2020-03-24 18:26:06 +01:00
|
|
|
for (j = block_recursive; j < block_size; j++)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-03-24 18:26:06 +01:00
|
|
|
int varr = getBlockVariableID(num, j);
|
2009-12-16 18:13:23 +01:00
|
|
|
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
|
|
|
|
}
|
2020-03-24 18:26:06 +01:00
|
|
|
for (j = block_recursive; j < block_size; j++)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-03-24 18:26:06 +01:00
|
|
|
int eqr = getBlockEquationID(num, j);
|
2009-12-16 18:13:23 +01:00
|
|
|
SaveCode.write(reinterpret_cast<char *>(&eqr), sizeof(eqr));
|
|
|
|
}
|
|
|
|
SaveCode.close();
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
void
|
2022-05-19 14:10:22 +02:00
|
|
|
StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context, bool no_tmp_terms, bool block)
|
2009-12-16 14:21:31 +01:00
|
|
|
{
|
2011-06-22 11:34:25 +02:00
|
|
|
initializeVariablesAndEquations();
|
|
|
|
|
2016-03-18 14:55:14 +01:00
|
|
|
vector<BinaryOpNode *> neweqs;
|
2020-03-24 18:26:06 +01:00
|
|
|
for (int eq = 0; eq < static_cast<int>(equations.size() - aux_equations.size()); eq++)
|
2016-03-21 11:51:48 +01:00
|
|
|
{
|
|
|
|
expr_t eq_tmp = equations[eq]->substituteStaticAuxiliaryVariable();
|
|
|
|
neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
|
|
|
|
}
|
2016-03-18 14:55:14 +01:00
|
|
|
|
2019-12-20 16:59:30 +01:00
|
|
|
for (auto &aux_equation : aux_equations)
|
2016-03-21 20:42:35 +01:00
|
|
|
{
|
2018-06-04 12:26:16 +02:00
|
|
|
expr_t eq_tmp = aux_equation->substituteStaticAuxiliaryDefinition();
|
2016-03-21 20:42:35 +01:00
|
|
|
neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
|
|
|
|
}
|
2017-06-14 07:01:31 +02:00
|
|
|
|
2016-03-18 14:55:14 +01:00
|
|
|
equations.clear();
|
2017-06-14 07:01:31 +02:00
|
|
|
copy(neweqs.begin(), neweqs.end(), back_inserter(equations));
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2022-01-21 14:31:29 +01:00
|
|
|
/* In both MATLAB and Julia, tensors for higher-order derivatives are stored
|
|
|
|
in matrices whose columns correspond to variable multi-indices. Since we
|
|
|
|
currently are limited to 32-bit signed integers (hence 31 bits) for matrix
|
|
|
|
indices, check that we will not overflow (see #89). Note that such a check
|
|
|
|
is not needed for parameter derivatives, since tensors for those are not
|
|
|
|
stored as matrices. This check is implemented at this place for symmetry
|
|
|
|
with DynamicModel::computingPass(). */
|
|
|
|
if (log2(symbol_table.endo_nbr())*derivsOrder >= numeric_limits<int>::digits)
|
|
|
|
{
|
|
|
|
cerr << "ERROR: The static derivatives matrix is too large. Please decrease the approximation order." << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
|
2018-11-22 14:32:40 +01:00
|
|
|
// Compute derivatives w.r. to all endogenous
|
|
|
|
set<int> vars;
|
2009-12-16 18:13:23 +01:00
|
|
|
for (int i = 0; i < symbol_table.endo_nbr(); i++)
|
2016-03-21 11:51:48 +01:00
|
|
|
{
|
2018-07-17 18:34:07 +02:00
|
|
|
int id = symbol_table.getID(SymbolType::endogenous, i);
|
2016-03-21 20:42:35 +01:00
|
|
|
vars.insert(getDerivID(id, 0));
|
2017-06-14 07:01:31 +02:00
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
// Launch computations
|
2018-12-20 17:04:28 +01:00
|
|
|
cout << "Computing static model derivatives (order " << derivsOrder << ")." << endl;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2018-11-22 14:32:40 +01:00
|
|
|
computeDerivatives(derivsOrder, vars);
|
2013-11-22 21:09:04 +01:00
|
|
|
|
2016-05-18 12:26:19 +02:00
|
|
|
if (paramsDerivsOrder > 0)
|
2012-11-29 18:07:48 +01:00
|
|
|
{
|
2018-12-20 17:04:28 +01:00
|
|
|
cout << "Computing static model derivatives w.r.t. parameters (order " << paramsDerivsOrder << ")." << endl;
|
2016-05-18 12:26:19 +02:00
|
|
|
computeParamsDerivatives(paramsDerivsOrder);
|
2012-11-29 18:07:48 +01:00
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
if (block)
|
|
|
|
{
|
2020-05-06 18:09:44 +02:00
|
|
|
auto contemporaneous_jacobian = evaluateAndReduceJacobian(eval_context);
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2022-04-19 17:06:37 +02:00
|
|
|
computeNonSingularNormalization(contemporaneous_jacobian, false);
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-04-30 12:48:16 +02:00
|
|
|
auto [prologue, epilogue] = computePrologueAndEpilogue();
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-03-30 14:51:53 +02:00
|
|
|
auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-03-20 15:23:23 +01:00
|
|
|
equationTypeDetermination(first_order_endo_derivatives, mfs);
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2018-12-20 17:04:28 +01:00
|
|
|
cout << "Finding the optimal block decomposition of the model ..." << endl;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-04-30 12:48:16 +02:00
|
|
|
computeBlockDecomposition(prologue, epilogue);
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-04-30 12:48:16 +02:00
|
|
|
reduceBlockDecomposition();
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-04-15 17:56:28 +02:00
|
|
|
printBlockDecomposition();
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-03-19 17:46:10 +01:00
|
|
|
computeChainRuleJacobian();
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-03-20 15:23:23 +01:00
|
|
|
determineLinearBlocks();
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
if (!no_tmp_terms)
|
2020-05-06 17:13:47 +02:00
|
|
|
computeBlockTemporaryTerms();
|
2009-12-16 14:21:31 +01:00
|
|
|
}
|
|
|
|
else
|
2010-01-22 11:03:29 +01:00
|
|
|
{
|
2018-09-25 15:56:52 +02:00
|
|
|
computeTemporaryTerms(true, no_tmp_terms);
|
2018-11-16 18:13:34 +01:00
|
|
|
|
|
|
|
/* Must be called after computeTemporaryTerms(), because it depends on
|
|
|
|
temporary_terms_mlv to be filled */
|
|
|
|
if (paramsDerivsOrder > 0 && !no_tmp_terms)
|
|
|
|
computeParamsDerivativesTemporaryTerms();
|
2010-01-22 11:03:29 +01:00
|
|
|
}
|
2009-04-14 16:39:53 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2018-06-27 15:01:31 +02:00
|
|
|
StaticModel::writeStaticMFile(const string &basename) const
|
2009-04-14 16:39:53 +02:00
|
|
|
{
|
2022-07-11 17:33:09 +02:00
|
|
|
auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::matlabStaticModel>();
|
|
|
|
|
|
|
|
ostringstream init_output, end_output;
|
|
|
|
init_output << "residual = zeros(" << equations.size() << ", 1);";
|
|
|
|
end_output << "if ~isreal(residual)" << endl
|
|
|
|
<< " residual = real(residual)+imag(residual).^2;" << endl
|
|
|
|
<< "end";
|
2022-07-12 17:04:41 +02:00
|
|
|
writeStaticMFileHelper(basename, "static_resid", "residual", "static_resid_tt",
|
2022-07-11 17:33:09 +02:00
|
|
|
temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(),
|
|
|
|
"", init_output, end_output,
|
|
|
|
d_output[0], tt_output[0]);
|
|
|
|
|
|
|
|
init_output.str("");
|
|
|
|
end_output.str("");
|
|
|
|
init_output << "g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");";
|
|
|
|
end_output << "if ~isreal(g1)" << endl
|
|
|
|
<< " g1 = real(g1)+2*imag(g1);" << endl
|
|
|
|
<< "end";
|
2022-07-12 17:04:41 +02:00
|
|
|
writeStaticMFileHelper(basename, "static_g1", "g1", "static_g1_tt",
|
2022-07-11 17:33:09 +02:00
|
|
|
temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(),
|
|
|
|
"static_resid_tt",
|
|
|
|
init_output, end_output,
|
|
|
|
d_output[1], tt_output[1]);
|
2022-07-12 17:04:41 +02:00
|
|
|
writeStaticMWrapperFunction(basename, "g1");
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
|
|
// For order ≥ 2
|
|
|
|
int ncols{symbol_table.endo_nbr()};
|
|
|
|
int ntt{static_cast<int>(temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size())};
|
|
|
|
for (size_t i{2}; i < derivatives.size(); i++)
|
|
|
|
{
|
|
|
|
ncols *= symbol_table.endo_nbr();
|
|
|
|
ntt += temporary_terms_derivatives[i].size();
|
|
|
|
string gname{"g" + to_string(i)};
|
|
|
|
string gprevname{"g" + to_string(i-1)};
|
|
|
|
|
|
|
|
init_output.str("");
|
|
|
|
end_output.str("");
|
|
|
|
if (derivatives[i].size())
|
|
|
|
{
|
|
|
|
init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl
|
|
|
|
<< gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl
|
|
|
|
<< gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl;
|
|
|
|
end_output << gname << " = sparse("
|
|
|
|
<< gname << "_i," << gname << "_j," << gname << "_v,"
|
|
|
|
<< equations.size() << "," << ncols << ");";
|
|
|
|
}
|
|
|
|
else
|
|
|
|
init_output << gname << " = sparse([],[],[]," << equations.size() << "," << ncols << ");";
|
2022-07-12 17:04:41 +02:00
|
|
|
writeStaticMFileHelper(basename, "static_" + gname, gname,
|
2022-07-11 17:33:09 +02:00
|
|
|
"static_" + gname + "_tt",
|
|
|
|
ntt,
|
|
|
|
"static_" + gprevname + "_tt",
|
|
|
|
init_output, end_output,
|
|
|
|
d_output[i], tt_output[i]);
|
|
|
|
if (i <= 3)
|
2022-07-12 17:04:41 +02:00
|
|
|
writeStaticMWrapperFunction(basename, gname);
|
2022-07-11 17:33:09 +02:00
|
|
|
}
|
|
|
|
|
2022-07-12 17:04:41 +02:00
|
|
|
writeStaticMCompatFile(basename);
|
2018-03-27 17:14:30 +02:00
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2018-03-27 17:14:30 +02:00
|
|
|
void
|
2022-07-12 17:04:41 +02:00
|
|
|
StaticModel::writeStaticMWrapperFunction(const string &basename, const string &ending) const
|
2018-03-27 17:14:30 +02:00
|
|
|
{
|
|
|
|
string name;
|
|
|
|
if (ending == "g1")
|
2018-06-27 15:01:31 +02:00
|
|
|
name = "static_resid_g1";
|
2018-03-27 17:14:30 +02:00
|
|
|
else if (ending == "g2")
|
2018-06-27 15:01:31 +02:00
|
|
|
name = "static_resid_g1_g2";
|
2018-03-27 17:14:30 +02:00
|
|
|
else if (ending == "g3")
|
2018-06-27 15:01:31 +02:00
|
|
|
name = "static_resid_g1_g2_g3";
|
2018-03-27 17:14:30 +02:00
|
|
|
|
2018-06-27 15:01:31 +02:00
|
|
|
string filename = packageDir(basename) + "/" + name + ".m";
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream output{filename, ios::out | ios::binary};
|
2009-12-16 14:21:31 +01:00
|
|
|
if (!output.is_open())
|
|
|
|
{
|
2018-03-27 17:14:30 +02:00
|
|
|
cerr << "Error: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
|
|
|
|
if (ending == "g1")
|
2018-05-23 17:32:50 +02:00
|
|
|
output << "function [residual, g1] = " << name << "(T, y, x, params, T_flag)" << endl
|
|
|
|
<< "% function [residual, g1] = " << name << "(T, y, x, params, T_flag)" << endl;
|
2018-03-27 17:14:30 +02:00
|
|
|
else if (ending == "g2")
|
2018-05-23 17:32:50 +02:00
|
|
|
output << "function [residual, g1, g2] = " << name << "(T, y, x, params, T_flag)" << endl
|
|
|
|
<< "% function [residual, g1, g2] = " << name << "(T, y, x, params, T_flag)" << endl;
|
2018-03-27 17:14:30 +02:00
|
|
|
else if (ending == "g3")
|
2018-05-23 17:32:50 +02:00
|
|
|
output << "function [residual, g1, g2, g3] = " << name << "(T, y, x, params, T_flag)" << endl
|
|
|
|
<< "% function [residual, g1, g2, g3] = " << name << "(T, y, x, params, T_flag)" << endl;
|
2018-03-27 17:14:30 +02:00
|
|
|
|
|
|
|
output << "%" << endl
|
|
|
|
<< "% Wrapper function automatically created by Dynare" << endl
|
2018-05-23 17:32:50 +02:00
|
|
|
<< "%" << endl
|
|
|
|
<< endl
|
|
|
|
<< " if T_flag" << endl
|
2018-06-27 15:01:31 +02:00
|
|
|
<< " T = " << basename << ".static_" << ending << "_tt(T, y, x, params);" << endl
|
2018-05-23 17:32:50 +02:00
|
|
|
<< " end" << endl;
|
2018-03-27 17:14:30 +02:00
|
|
|
|
|
|
|
if (ending == "g1")
|
2018-06-27 15:01:31 +02:00
|
|
|
output << " residual = " << basename << ".static_resid(T, y, x, params, false);" << endl
|
|
|
|
<< " g1 = " << basename << ".static_g1(T, y, x, params, false);" << endl;
|
2018-03-27 17:14:30 +02:00
|
|
|
else if (ending == "g2")
|
2018-06-27 15:01:31 +02:00
|
|
|
output << " [residual, g1] = " << basename << ".static_resid_g1(T, y, x, params, false);" << endl
|
|
|
|
<< " g2 = " << basename << ".static_g2(T, y, x, params, false);" << endl;
|
2018-03-27 17:14:30 +02:00
|
|
|
else if (ending == "g3")
|
2018-06-27 15:01:31 +02:00
|
|
|
output << " [residual, g1, g2] = " << basename << ".static_resid_g1_g2(T, y, x, params, false);" << endl
|
|
|
|
<< " g3 = " << basename << ".static_g3(T, y, x, params, false);" << endl;
|
2018-03-27 17:14:30 +02:00
|
|
|
|
|
|
|
output << endl << "end" << endl;
|
|
|
|
output.close();
|
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2022-07-12 17:04:41 +02:00
|
|
|
StaticModel::writeStaticMFileHelper(const string &basename,
|
2018-06-27 15:01:31 +02:00
|
|
|
const string &name, const string &retvalname,
|
2018-03-27 17:14:30 +02:00
|
|
|
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) const
|
|
|
|
{
|
2018-06-27 15:01:31 +02:00
|
|
|
string filename = packageDir(basename) + "/" + name_tt + ".m";
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream output{filename, ios::out | ios::binary};
|
2018-03-27 17:14:30 +02:00
|
|
|
if (!output.is_open())
|
|
|
|
{
|
|
|
|
cerr << "Error: Can't open file " << filename << " for writing" << endl;
|
2009-12-16 14:21:31 +01:00
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
|
2018-03-27 17:14:30 +02:00
|
|
|
output << "function T = " << name_tt << "(T, y, x, params)" << endl
|
|
|
|
<< "% function T = " << name_tt << "(T, y, x, params)" << endl
|
2009-07-07 16:20:48 +02:00
|
|
|
<< "%" << endl
|
2018-03-27 17:14:30 +02:00
|
|
|
<< "% File created by Dynare Preprocessor from .mod file" << endl
|
2009-07-07 16:20:48 +02:00
|
|
|
<< "%" << endl
|
2018-03-27 17:14:30 +02:00
|
|
|
<< "% Inputs:" << endl
|
|
|
|
<< "% T [#temp variables by 1] double vector of temporary terms to be filled by function" << endl
|
|
|
|
<< "% y [M_.endo_nbr by 1] double vector of endogenous variables in declaration order" << endl
|
|
|
|
<< "% x [M_.exo_nbr by 1] double vector of exogenous variables in declaration order" << endl
|
|
|
|
<< "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
|
2013-07-28 11:32:14 +02:00
|
|
|
<< "%" << endl
|
2018-03-27 17:14:30 +02:00
|
|
|
<< "% Output:" << endl
|
|
|
|
<< "% T [#temp variables by 1] double vector of temporary terms" << endl
|
|
|
|
<< "%" << endl << endl
|
|
|
|
<< "assert(length(T) >= " << ttlen << ");" << endl
|
|
|
|
<< endl;
|
|
|
|
|
|
|
|
if (!previous_tt_name.empty())
|
2018-06-27 15:01:31 +02:00
|
|
|
output << "T = " << basename << "." << previous_tt_name << "(T, y, x, params);" << endl << endl;
|
2018-03-27 17:14:30 +02:00
|
|
|
|
|
|
|
output << s_tt.str() << endl
|
|
|
|
<< "end" << endl;
|
|
|
|
output.close();
|
|
|
|
|
2018-06-27 15:01:31 +02:00
|
|
|
filename = packageDir(basename) + "/" + name + ".m";
|
2018-06-27 15:12:12 +02:00
|
|
|
output.open(filename, ios::out | ios::binary);
|
2018-03-27 17:14:30 +02:00
|
|
|
if (!output.is_open())
|
|
|
|
{
|
|
|
|
cerr << "Error: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
|
|
|
|
output << "function " << retvalname << " = " << name << "(T, y, x, params, T_flag)" << endl
|
|
|
|
<< "% function " << retvalname << " = " << name << "(T, y, x, params, T_flag)" << endl
|
2013-07-28 11:32:14 +02:00
|
|
|
<< "%" << endl
|
2018-03-27 17:14:30 +02:00
|
|
|
<< "% File created by Dynare Preprocessor from .mod file" << endl
|
2017-06-14 07:01:31 +02:00
|
|
|
<< "%" << endl
|
2018-03-27 17:14:30 +02:00
|
|
|
<< "% Inputs:" << endl
|
|
|
|
<< "% T [#temp variables by 1] double vector of temporary terms to be filled by function" << endl
|
|
|
|
<< "% y [M_.endo_nbr by 1] double vector of endogenous variables in declaration order" << endl
|
|
|
|
<< "% x [M_.exo_nbr by 1] double vector of exogenous variables in declaration order" << endl
|
|
|
|
<< "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
|
|
|
|
<< "% to evaluate the model" << endl
|
|
|
|
<< "% T_flag boolean boolean flag saying whether or not to calculate temporary terms" << endl
|
|
|
|
<< "%" << endl
|
|
|
|
<< "% Output:" << endl
|
|
|
|
<< "% " << retvalname << endl
|
|
|
|
<< "%" << endl << endl;
|
|
|
|
|
|
|
|
if (!name_tt.empty())
|
|
|
|
output << "if T_flag" << endl
|
2018-06-27 15:01:31 +02:00
|
|
|
<< " T = " << basename << "." << name_tt << "(T, y, x, params);" << endl
|
2018-03-27 17:14:30 +02:00
|
|
|
<< "end" << endl;
|
|
|
|
|
|
|
|
output << init_s.str() << endl
|
|
|
|
<< s.str()
|
|
|
|
<< end_s.str() << endl
|
|
|
|
<< "end" << endl;
|
2011-08-18 12:44:11 +02:00
|
|
|
output.close();
|
|
|
|
}
|
|
|
|
|
2018-05-23 16:10:26 +02:00
|
|
|
void
|
2022-07-12 17:04:41 +02:00
|
|
|
StaticModel::writeStaticMCompatFile(const string &basename) const
|
2018-05-23 16:10:26 +02:00
|
|
|
{
|
2018-06-27 15:01:31 +02:00
|
|
|
string filename = packageDir(basename) + "/static.m";
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream output{filename, ios::out | ios::binary};
|
2018-05-23 16:10:26 +02:00
|
|
|
if (!output.is_open())
|
|
|
|
{
|
|
|
|
cerr << "Error: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
2018-11-15 16:39:53 +01:00
|
|
|
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();
|
2018-05-23 16:10:26 +02:00
|
|
|
|
2018-06-27 15:01:31 +02:00
|
|
|
output << "function [residual, g1, g2, g3] = static(y, x, params)" << endl
|
2018-05-23 16:10:26 +02:00
|
|
|
<< " T = NaN(" << ntt << ", 1);" << endl
|
|
|
|
<< " if nargout <= 1" << endl
|
2018-06-27 15:01:31 +02:00
|
|
|
<< " residual = " << basename << ".static_resid(T, y, x, params, true);" << endl
|
2018-05-23 16:10:26 +02:00
|
|
|
<< " elseif nargout == 2" << endl
|
2018-06-27 15:01:31 +02:00
|
|
|
<< " [residual, g1] = " << basename << ".static_resid_g1(T, y, x, params, true);" << endl
|
2018-05-23 16:10:26 +02:00
|
|
|
<< " elseif nargout == 3" << endl
|
2018-06-27 15:01:31 +02:00
|
|
|
<< " [residual, g1, g2] = " << basename << ".static_resid_g1_g2(T, y, x, params, true);" << endl
|
2018-05-23 16:10:26 +02:00
|
|
|
<< " else" << endl
|
2018-06-27 15:01:31 +02:00
|
|
|
<< " [residual, g1, g2, g3] = " << basename << ".static_resid_g1_g2_g3(T, y, x, params, true);" << endl
|
2018-05-23 16:10:26 +02:00
|
|
|
<< " end" << endl
|
|
|
|
<< "end" << endl;
|
|
|
|
|
|
|
|
output.close();
|
|
|
|
}
|
|
|
|
|
2011-08-18 12:44:11 +02:00
|
|
|
void
|
2018-06-27 15:01:31 +02:00
|
|
|
StaticModel::writeStaticCFile(const string &basename) const
|
2011-08-18 12:44:11 +02:00
|
|
|
{
|
|
|
|
// Writing comments and function definition command
|
2022-07-11 17:33:09 +02:00
|
|
|
string filename{basename + "/model/src/static.c"};
|
2011-08-18 12:44:11 +02:00
|
|
|
|
2022-07-11 17:33:09 +02:00
|
|
|
int ntt{static_cast<int>(temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size())};
|
2018-09-25 19:15:22 +02:00
|
|
|
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream output{filename, ios::out | ios::binary};
|
2011-08-18 12:44:11 +02:00
|
|
|
if (!output.is_open())
|
|
|
|
{
|
|
|
|
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
|
|
|
|
output << "/*" << endl
|
|
|
|
<< " * " << filename << " : Computes static model for Dynare" << endl
|
|
|
|
<< " *" << endl
|
|
|
|
<< " * Warning : this file is generated automatically by Dynare" << endl
|
|
|
|
<< " * from model file (.mod)" << endl << endl
|
|
|
|
<< " */" << endl
|
2020-05-29 16:12:01 +02:00
|
|
|
<< endl
|
|
|
|
<< "#include <math.h>" << endl
|
|
|
|
<< "#include <stdlib.h>" << endl
|
|
|
|
<< R"(#include "mex.h")" << endl
|
|
|
|
<< endl;
|
2011-08-18 12:44:11 +02:00
|
|
|
|
2018-07-18 16:18:26 +02:00
|
|
|
// Write function definition if BinaryOpcode::powerDeriv is used
|
2020-05-29 16:12:01 +02:00
|
|
|
writePowerDeriv(output);
|
2011-08-18 12:44:11 +02:00
|
|
|
|
2018-09-25 19:15:22 +02:00
|
|
|
output << endl;
|
|
|
|
|
2022-07-11 17:33:09 +02:00
|
|
|
auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::CStaticModel>();
|
|
|
|
|
|
|
|
for (size_t i = 0; i < d_output.size(); i++)
|
|
|
|
{
|
|
|
|
string funcname{i == 0 ? "resid" : "g" + to_string(i)};
|
|
|
|
output << "void static_" << funcname << "_tt(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl
|
|
|
|
<< "{" << endl
|
|
|
|
<< tt_output[i].str()
|
|
|
|
<< "}" << endl
|
|
|
|
<< endl
|
|
|
|
<< "void static_" << funcname << "(const double *restrict y, const double *restrict x, const double *restrict params, const double *restrict T, ";
|
|
|
|
if (i == 0)
|
|
|
|
output << "double *restrict residual";
|
|
|
|
else if (i == 1)
|
|
|
|
output << "double *restrict g1";
|
|
|
|
else
|
|
|
|
output << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v";
|
|
|
|
output << ")" << endl
|
|
|
|
<< "{" << endl;
|
|
|
|
if (i == 0)
|
|
|
|
output << " double lhs, rhs;" << endl;
|
|
|
|
output << d_output[i].str()
|
|
|
|
<< "}" << endl
|
|
|
|
<< endl;
|
|
|
|
}
|
2018-09-25 19:15:22 +02:00
|
|
|
|
2020-05-29 16:12:01 +02:00
|
|
|
output << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
|
2011-08-18 12:44:11 +02:00
|
|
|
<< "{" << endl
|
2020-06-23 16:30:21 +02:00
|
|
|
<< " if (nrhs > 3)" << endl
|
|
|
|
<< R"( mexErrMsgTxt("Accepts at most 3 output arguments");)" << endl
|
|
|
|
<< " if (nrhs != 3)" << endl
|
|
|
|
<< R"( mexErrMsgTxt("Requires exactly 3 input arguments");)" << endl
|
2018-09-25 19:15:22 +02:00
|
|
|
<< " double *y = mxGetPr(prhs[0]);" << endl
|
|
|
|
<< " double *x = mxGetPr(prhs[1]);" << endl
|
|
|
|
<< " double *params = mxGetPr(prhs[2]);" << endl
|
|
|
|
<< endl
|
2020-05-29 16:12:01 +02:00
|
|
|
<< " double *T = (double *) malloc(sizeof(double)*" << ntt << ");" << endl
|
2011-08-18 12:44:11 +02:00
|
|
|
<< endl
|
|
|
|
<< " if (nlhs >= 1)" << endl
|
|
|
|
<< " {" << endl
|
|
|
|
<< " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
|
2018-09-25 19:15:22 +02:00
|
|
|
<< " double *residual = mxGetPr(plhs[0]);" << endl
|
2020-06-23 18:19:27 +02:00
|
|
|
<< " static_resid_tt(y, x, params, T);" << endl
|
|
|
|
<< " static_resid(y, x, params, T, residual);" << endl
|
2011-08-18 12:44:11 +02:00
|
|
|
<< " }" << endl
|
|
|
|
<< endl
|
|
|
|
<< " if (nlhs >= 2)" << endl
|
|
|
|
<< " {" << endl
|
|
|
|
<< " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl
|
2018-09-25 19:15:22 +02:00
|
|
|
<< " double *g1 = mxGetPr(plhs[1]);" << endl
|
2020-06-23 18:19:27 +02:00
|
|
|
<< " static_g1_tt(y, x, params, T);" << endl
|
|
|
|
<< " static_g1(y, x, params, T, g1);" << endl
|
2011-08-18 12:44:11 +02:00
|
|
|
<< " }" << endl
|
|
|
|
<< endl
|
|
|
|
<< " if (nlhs >= 3)" << endl
|
|
|
|
<< " {" << endl
|
2020-06-23 17:50:50 +02:00
|
|
|
<< " mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
|
|
|
|
<< " mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
|
|
|
|
<< " mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
|
2020-06-23 18:19:27 +02:00
|
|
|
<< " static_g2_tt(y, x, params, T);" << endl
|
|
|
|
<< " static_g2(y, x, params, T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
|
2020-06-23 17:50:50 +02:00
|
|
|
<< " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
|
|
|
|
<< " mxArray *n = mxCreateDoubleScalar(" << symbol_table.endo_nbr()*symbol_table.endo_nbr() << ");" << endl
|
|
|
|
<< " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
|
|
|
|
<< R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
|
|
|
|
<< " plhs[2] = plhs_sparse[0];" << endl
|
|
|
|
<< " mxDestroyArray(g2_i);" << endl
|
|
|
|
<< " mxDestroyArray(g2_j);" << endl
|
|
|
|
<< " mxDestroyArray(g2_v);" << endl
|
|
|
|
<< " mxDestroyArray(m);" << endl
|
|
|
|
<< " mxDestroyArray(n);" << endl
|
2011-08-18 12:44:11 +02:00
|
|
|
<< " }" << endl
|
|
|
|
<< endl
|
2018-09-25 19:15:22 +02:00
|
|
|
<< " free(T);" << endl
|
|
|
|
<< "}" << endl;
|
2020-05-29 16:12:01 +02:00
|
|
|
|
2009-07-07 16:20:48 +02:00
|
|
|
output.close();
|
2009-04-14 16:39:53 +02:00
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
void
|
2015-07-27 15:59:13 +02:00
|
|
|
StaticModel::writeStaticJuliaFile(const string &basename) const
|
2015-07-27 15:33:38 +02:00
|
|
|
{
|
2022-07-11 17:33:09 +02:00
|
|
|
auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::juliaStaticModel>();
|
|
|
|
|
|
|
|
stringstream output;
|
|
|
|
output << "module " << basename << "Static" << endl
|
|
|
|
<< "#" << endl
|
|
|
|
<< "# NB: this file was automatically generated by Dynare" << endl
|
|
|
|
<< "# from " << basename << ".mod" << endl
|
|
|
|
<< "#" << endl
|
|
|
|
<< "using StatsFuns" << endl << endl
|
|
|
|
<< "export tmp_nbr, static!, staticResid!, staticG1!, staticG2!, staticG3!" << endl << endl
|
|
|
|
<< "#=" << endl
|
|
|
|
<< "# The comments below apply to all functions contained in this module #" << endl
|
|
|
|
<< " NB: The arguments contained on the first line of the function" << endl
|
|
|
|
<< " definition are those that are modified in place" << endl << endl
|
|
|
|
<< "## Exported Functions ##" << endl
|
|
|
|
<< " static! : Wrapper function; computes residuals, Jacobian, Hessian," << endl
|
|
|
|
<< " and third order derivatives matroces depending on the arguments provided" << endl
|
|
|
|
<< " staticResid! : Computes the static model residuals" << endl
|
|
|
|
<< " staticG1! : Computes the static model Jacobian" << endl
|
|
|
|
<< " staticG2! : Computes the static model Hessian" << endl
|
|
|
|
<< " staticG3! : Computes the static model third derivatives" << endl << endl
|
|
|
|
<< "## Exported Variables ##" << endl
|
|
|
|
<< " tmp_nbr : Vector{Int}(4) respectively the number of temporary variables" << endl
|
|
|
|
<< " for the residuals, g1, g2 and g3." << endl << endl
|
|
|
|
<< "## Local Functions ##" << endl
|
|
|
|
<< " staticResidTT! : Computes the static model temporary terms for the residuals" << endl
|
|
|
|
<< " staticG1TT! : Computes the static model temporary terms for the Jacobian" << endl
|
|
|
|
<< " staticG2TT! : Computes the static model temporary terms for the Hessian" << endl
|
|
|
|
<< " staticG3TT! : Computes the static model temporary terms for the third derivatives" << endl << endl
|
|
|
|
<< "## Function Arguments ##" << endl
|
2022-07-12 11:15:51 +02:00
|
|
|
<< " T : Vector{<: Real}(num_temp_terms) temporary terms" << endl
|
|
|
|
<< " y : Vector{<: Real}(model_.endo_nbr) endogenous variables in declaration order" << endl
|
|
|
|
<< " x : Vector{<: Real}(model_.exo_nbr) exogenous variables in declaration order" << endl
|
|
|
|
<< " params : Vector{<: Real}(model_.param) parameter values in declaration order" << endl
|
|
|
|
<< " residual : Vector{<: Real}(model_.eq_nbr) residuals of the static model equations" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " in order of declaration of the equations. Dynare may prepend auxiliary equations," << endl
|
|
|
|
<< " see model.aux_vars" << endl
|
2022-07-12 11:15:51 +02:00
|
|
|
<< " g1 : Matrix{<: Real}(model.eq_nbr,model_.endo_nbr) Jacobian matrix of the static model equations" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " The columns and rows respectively correspond to the variables in declaration order and the" << endl
|
|
|
|
<< " equations in order of declaration" << endl
|
|
|
|
<< " g2 : spzeros(model.eq_nbr, model_.endo^2) Hessian matrix of the static model equations" << endl
|
|
|
|
<< " The columns and rows respectively correspond to the variables in declaration order and the" << endl
|
|
|
|
<< " equations in order of declaration" << endl
|
|
|
|
<< " g3 : spzeros(model.eq_nbr, model_.endo^3) Third order derivatives matrix of the static model equations" << endl
|
|
|
|
<< " The columns and rows respectively correspond to the variables in declaration order and the" << endl
|
|
|
|
<< " equations in order of declaration" << endl << endl
|
|
|
|
<< "## Remarks ##" << endl
|
|
|
|
<< " [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
|
|
|
|
<< " used for the different returned objects (residuals, jacobian, hessian or third order derivatives) is given by the elements in `tmp_nbr`" << endl
|
|
|
|
<< " exported vector. The first element is the number of temporaries used for the computation of the residuals, the second element is the" << endl
|
|
|
|
<< " number of temporaries used for the evaluation of the jacobian matrix, etc. If one calls the version of the static model computing the" << endl
|
|
|
|
<< " residuals, and the jacobian and hessian matrices, then `T` must have at least `sum(tmp_nbr[1:3])` elements." << endl
|
|
|
|
<< "=#" << endl << endl;
|
|
|
|
|
|
|
|
// Write the number of temporary terms
|
|
|
|
output << "tmp_nbr = zeros(Int,4)" << endl
|
|
|
|
<< "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << "# Number of temporary terms for the residuals" << endl
|
|
|
|
<< "tmp_nbr[2] = " << temporary_terms_derivatives[1].size() << "# Number of temporary terms for g1 (jacobian)" << endl
|
|
|
|
<< "tmp_nbr[3] = " << temporary_terms_derivatives[2].size() << "# Number of temporary terms for g2 (hessian)" << endl
|
|
|
|
<< "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl;
|
|
|
|
|
|
|
|
// staticResidTT!
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function staticResidTT!(T::Vector{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " @inbounds begin" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< tt_output[0].str()
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " end" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl << endl;
|
|
|
|
|
|
|
|
// static!
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function staticResid!(T::Vector{<: Real}, residual::Vector{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T0_flag::Bool)" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " @assert length(y) == " << symbol_table.endo_nbr() << endl
|
|
|
|
<< " @assert length(x) == " << symbol_table.exo_nbr() << endl
|
|
|
|
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
|
|
|
|
<< " @assert length(residual) == " << equations.size() << endl
|
|
|
|
<< " if T0_flag" << endl
|
|
|
|
<< " staticResidTT!(T, y, x, params)" << endl
|
|
|
|
<< " end" << endl
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " @inbounds begin" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< d_output[0].str()
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " end" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " if ~isreal(residual)" << endl
|
|
|
|
<< " residual = real(residual)+imag(residual).^2;" << endl
|
|
|
|
<< " end" << endl
|
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl << endl;
|
|
|
|
|
|
|
|
// staticG1TT!
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function staticG1TT!(T::Vector{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T0_flag::Bool)" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " if T0_flag" << endl
|
|
|
|
<< " staticResidTT!(T, y, x, params)" << endl
|
|
|
|
<< " end" << endl
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " @inbounds begin" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< tt_output[1].str()
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " end" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl << endl;
|
|
|
|
|
|
|
|
// staticG1!
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function staticG1!(T::Vector{<: Real}, g1::Matrix{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T1_flag::Bool, T0_flag::Bool)" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " @assert length(T) >= "
|
|
|
|
<< temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl
|
|
|
|
<< " @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr() << ")" << endl
|
|
|
|
<< " @assert length(y) == " << symbol_table.endo_nbr() << endl
|
|
|
|
<< " @assert length(x) == " << symbol_table.exo_nbr() << endl
|
|
|
|
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
|
|
|
|
<< " if T1_flag" << endl
|
|
|
|
<< " staticG1TT!(T, y, x, params, T0_flag)" << endl
|
|
|
|
<< " end" << endl
|
|
|
|
<< " fill!(g1, 0.0)" << endl
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " @inbounds begin" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< d_output[1].str()
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " end" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " if ~isreal(g1)" << endl
|
|
|
|
<< " g1 = real(g1)+2*imag(g1);" << endl
|
|
|
|
<< " end" << endl
|
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl << endl;
|
|
|
|
|
|
|
|
// staticG2TT!
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function staticG2TT!(T::Vector{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T1_flag::Bool, T0_flag::Bool)" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " if T1_flag" << endl
|
|
|
|
<< " staticG1TT!(T, y, x, params, TO_flag)" << endl
|
|
|
|
<< " end" << endl
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " @inbounds begin" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< tt_output[2].str()
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " end" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl << endl;
|
|
|
|
|
|
|
|
// staticG2!
|
|
|
|
int hessianColsNbr{symbol_table.endo_nbr() * symbol_table.endo_nbr()};
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function staticG2!(T::Vector{<: Real}, g2::Matrix{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " @assert length(T) >= "
|
|
|
|
<< temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl
|
|
|
|
<< " @assert size(g2) == (" << equations.size() << ", " << hessianColsNbr << ")" << endl
|
|
|
|
<< " @assert length(y) == " << symbol_table.endo_nbr() << endl
|
|
|
|
<< " @assert length(x) == " << symbol_table.exo_nbr() << endl
|
|
|
|
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
|
|
|
|
<< " if T2_flag" << endl
|
|
|
|
<< " staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl
|
|
|
|
<< " end" << endl
|
|
|
|
<< " fill!(g2, 0.0)" << endl
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " @inbounds begin" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< d_output[2].str()
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " end" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl << endl;
|
|
|
|
|
|
|
|
// staticG3TT!
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function staticG3TT!(T::Vector{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " if T2_flag" << endl
|
|
|
|
<< " staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl
|
|
|
|
<< " end" << endl
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " @inbounds begin" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< tt_output[3].str()
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " end" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl << endl;
|
|
|
|
|
|
|
|
// staticG3!
|
|
|
|
int ncols{hessianColsNbr * symbol_table.endo_nbr()};
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function staticG3!(T::Vector{<: Real}, g3::Matrix{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T3_flag::Bool, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " @assert length(T) >= "
|
|
|
|
<< 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
|
|
|
|
<< " @assert size(g3) == (" << equations.size() << ", " << ncols << ")" << endl
|
|
|
|
<< " @assert length(y) == " << symbol_table.endo_nbr() << endl
|
|
|
|
<< " @assert length(x) == " << symbol_table.exo_nbr() << endl
|
|
|
|
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
|
|
|
|
<< " if T3_flag" << endl
|
|
|
|
<< " staticG3TT!(T, y, x, params, T2_flag, T1_flag, T0_flag)" << endl
|
|
|
|
<< " end" << endl
|
|
|
|
<< " fill!(g3, 0.0)" << endl
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " @inbounds begin" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< d_output[3].str()
|
2022-07-12 13:07:33 +02:00
|
|
|
<< " end" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl << endl;
|
|
|
|
|
|
|
|
// static!
|
2022-07-12 11:15:51 +02:00
|
|
|
output << "function static!(T::Vector{<: Real}, residual::Vector{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " staticResid!(T, residual, y, x, params, true)" << endl
|
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl
|
|
|
|
<< endl
|
2022-07-12 11:15:51 +02:00
|
|
|
<< "function static!(T::Vector{<: Real}, residual::Vector{<: Real}, g1::Matrix{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " staticG1!(T, g1, y, x, params, true, true)" << endl
|
|
|
|
<< " staticResid!(T, residual, y, x, params, false)" << endl
|
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl
|
|
|
|
<< endl
|
2022-07-12 11:15:51 +02:00
|
|
|
<< "function static!(T::Vector{<: Real}, g1::Matrix{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " staticG1!(T, g1, y, x, params, true, false)" << endl
|
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl
|
|
|
|
<< endl
|
2022-07-12 11:15:51 +02:00
|
|
|
<< "function static!(T::Vector{<: Real}, residual::Vector{<: Real}, g1::Matrix{<: Real}, g2::Matrix{<: Real}," << endl
|
|
|
|
<< " y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
|
2022-07-11 17:33:09 +02:00
|
|
|
<< " staticG2!(T, g2, y, x, params, true, true, true)" << endl
|
|
|
|
<< " staticG1!(T, g1, y, x, params, false, false)" << endl
|
|
|
|
<< " staticResid!(T, residual, y, x, params, false)" << endl
|
|
|
|
<< " return nothing" << endl
|
|
|
|
<< "end" << endl
|
|
|
|
<< endl;
|
|
|
|
|
|
|
|
// Write function definition if BinaryOpcode::powerDeriv is used
|
|
|
|
writePowerDerivJulia(output);
|
|
|
|
|
|
|
|
output << "end" << endl;
|
|
|
|
|
|
|
|
writeToFileIfModified(output, basename + "Static.jl");
|
2015-07-27 15:33:38 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2022-05-19 14:15:43 +02:00
|
|
|
StaticModel::writeStaticFile(const string &basename, bool block, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2020-06-23 15:13:04 +02:00
|
|
|
filesystem::path model_dir{basename};
|
|
|
|
model_dir /= "model";
|
|
|
|
if (use_dll)
|
|
|
|
filesystem::create_directories(model_dir / "src");
|
2022-05-19 14:15:43 +02:00
|
|
|
filesystem::create_directories(model_dir / "bytecode");
|
2020-06-23 15:13:04 +02:00
|
|
|
|
|
|
|
if (block)
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2022-05-19 14:15:43 +02:00
|
|
|
writeStaticBlockBytecode(basename);
|
|
|
|
|
|
|
|
if (use_dll)
|
2020-06-23 15:13:04 +02:00
|
|
|
{
|
|
|
|
writeStaticPerBlockCFiles(basename);
|
|
|
|
writeStaticBlockCFile(basename);
|
2022-06-02 10:50:21 +02:00
|
|
|
vector<filesystem::path> src_files(blocks.size() + 1);
|
2020-06-23 15:13:04 +02:00
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
|
|
|
src_files[blk] = model_dir / "src" / ("static_" + to_string(blk+1) + ".c");
|
|
|
|
src_files[blocks.size()] = model_dir / "src" / "static.c";
|
|
|
|
compileMEX(basename, "static", mexext, src_files, matlabroot, dynareroot);
|
|
|
|
}
|
|
|
|
else if (julia)
|
|
|
|
{
|
|
|
|
cerr << "'block' option is not available with Julia" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
else // M-files
|
|
|
|
{
|
|
|
|
writeStaticPerBlockMFiles(basename);
|
|
|
|
writeStaticBlockMFile(basename);
|
|
|
|
}
|
2009-12-16 18:13:23 +01:00
|
|
|
}
|
2020-06-23 15:13:04 +02:00
|
|
|
else
|
2018-10-26 11:44:26 +02:00
|
|
|
{
|
2022-05-19 14:15:43 +02:00
|
|
|
writeStaticBytecode(basename);
|
|
|
|
|
|
|
|
if (use_dll)
|
2020-06-23 15:13:04 +02:00
|
|
|
{
|
|
|
|
writeStaticCFile(basename);
|
|
|
|
compileMEX(basename, "static", mexext, { model_dir / "src" / "static.c" },
|
|
|
|
matlabroot, dynareroot);
|
|
|
|
}
|
|
|
|
else if (julia)
|
|
|
|
writeStaticJuliaFile(basename);
|
|
|
|
else // M-files
|
|
|
|
writeStaticMFile(basename);
|
2018-10-26 11:44:26 +02:00
|
|
|
}
|
2020-06-23 15:13:04 +02:00
|
|
|
|
2016-04-04 17:11:03 +02:00
|
|
|
writeSetAuxiliaryVariables(basename, julia);
|
2009-12-16 18:13:23 +01:00
|
|
|
}
|
2009-04-27 19:15:14 +02:00
|
|
|
|
2016-08-12 11:50:57 +02:00
|
|
|
bool
|
|
|
|
StaticModel::exoPresentInEqs() const
|
|
|
|
{
|
2018-06-04 12:26:16 +02:00
|
|
|
for (auto equation : equations)
|
2021-07-20 12:10:58 +02:00
|
|
|
if (equation->hasExogenous())
|
|
|
|
return true;
|
2016-08-12 11:50:57 +02:00
|
|
|
return false;
|
|
|
|
}
|
|
|
|
|
2009-04-27 19:15:14 +02:00
|
|
|
void
|
2020-05-19 16:56:53 +02:00
|
|
|
StaticModel::writeStaticBlockMFile(const string &basename) const
|
2009-04-27 19:15:14 +02:00
|
|
|
{
|
2018-06-27 15:01:31 +02:00
|
|
|
string filename = packageDir(basename) + "/static.m";
|
2009-04-27 19:15:14 +02:00
|
|
|
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream output{filename, ios::out | ios::binary};
|
2009-12-16 14:21:31 +01:00
|
|
|
if (!output.is_open())
|
2009-04-27 19:15:14 +02:00
|
|
|
{
|
2009-12-16 14:21:31 +01:00
|
|
|
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
2009-04-27 19:15:14 +02:00
|
|
|
}
|
|
|
|
|
2020-06-19 18:33:53 +02:00
|
|
|
output << "function [residual, y, T, g1] = static(nblock, y, x, params, T)" << endl
|
2009-12-16 14:21:31 +01:00
|
|
|
<< " switch nblock" << endl;
|
2009-06-30 17:07:09 +02:00
|
|
|
|
2020-05-19 16:56:53 +02:00
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
2009-06-30 17:07:09 +02:00
|
|
|
{
|
2020-05-19 16:56:53 +02:00
|
|
|
output << " case " << blk+1 << endl;
|
2009-04-27 19:15:14 +02:00
|
|
|
|
2020-05-19 16:56:53 +02:00
|
|
|
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
2020-03-20 17:31:14 +01:00
|
|
|
if (simulation_type == BlockSimulationType::evaluateBackward
|
|
|
|
|| simulation_type == BlockSimulationType::evaluateForward)
|
2020-06-19 18:33:53 +02:00
|
|
|
output << " [y, T] = " << basename << ".block.static_" << blk+1 << "(y, x, params, T);" << endl
|
|
|
|
<< " residual = [];" << endl
|
|
|
|
<< " g1 = [];" << endl;
|
2009-12-16 18:13:23 +01:00
|
|
|
else
|
2020-10-06 17:16:50 +02:00
|
|
|
output << " [residual, y, T, g1] = " << basename << ".block.static_" << blk+1 << "(y, x, params, T);" << endl;
|
2010-09-17 16:53:27 +02:00
|
|
|
|
2009-04-28 19:11:48 +02:00
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
output << " end" << endl
|
|
|
|
<< "end" << endl;
|
|
|
|
output.close();
|
|
|
|
}
|
2009-04-28 19:11:48 +02:00
|
|
|
|
2020-06-23 15:13:04 +02:00
|
|
|
void
|
|
|
|
StaticModel::writeStaticBlockCFile(const string &basename) const
|
|
|
|
{
|
|
|
|
string filename = basename + "/model/src/static.c";
|
|
|
|
|
2022-07-11 16:09:07 +02:00
|
|
|
ofstream output{filename, ios::out | ios::binary};
|
2020-06-23 15:13:04 +02:00
|
|
|
if (!output.is_open())
|
|
|
|
{
|
|
|
|
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
}
|
|
|
|
|
|
|
|
output << "#include <math.h>" << endl
|
|
|
|
<< R"(#include "mex.h")" << endl;
|
|
|
|
|
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
|
|
|
output << R"(#include "static_)" << blk+1 << R"(.h")" << endl;
|
|
|
|
|
|
|
|
output << endl;
|
|
|
|
writePowerDeriv(output);
|
|
|
|
|
|
|
|
output << endl
|
|
|
|
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
|
|
|
|
<< "{" << endl
|
|
|
|
<< " if (nrhs != 5)" << endl
|
|
|
|
<< R"( mexErrMsgTxt("Requires exactly 5 input arguments");)" << endl
|
|
|
|
<< " if (nlhs > 4)" << endl
|
|
|
|
<< R"( mexErrMsgTxt("Accepts at most 4 output arguments");)" << endl
|
|
|
|
<< " int nblock = (int) mxGetScalar(prhs[0]);" << endl
|
|
|
|
<< " const mxArray *y = prhs[1], *x = prhs[2], *params = prhs[3], *T = prhs[4];" << endl
|
|
|
|
<< " mxArray *T_new = mxDuplicateArray(T);" << endl
|
|
|
|
<< " mxArray *y_new = mxDuplicateArray(y);" << endl
|
|
|
|
<< " mxArray *residual, *g1;" << endl
|
|
|
|
<< " switch (nblock)" << endl
|
|
|
|
<< " {" << endl;
|
|
|
|
|
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
|
|
|
{
|
|
|
|
output << " case " << blk+1 << ':' << endl;
|
|
|
|
|
|
|
|
BlockSimulationType simulation_type = blocks[blk].simulation_type;
|
|
|
|
|
|
|
|
if (simulation_type == BlockSimulationType::evaluateBackward
|
|
|
|
|| simulation_type == BlockSimulationType::evaluateForward)
|
|
|
|
output << " static_" << blk+1 << "_mx(y_new, x, params, T_new);" << endl
|
|
|
|
<< " residual = mxCreateDoubleMatrix(0,0,mxREAL);" << endl
|
|
|
|
<< " g1 = mxCreateDoubleMatrix(0,0,mxREAL);" << endl;
|
|
|
|
else
|
2020-10-06 17:16:50 +02:00
|
|
|
output << " static_" << blk+1 << "_mx(y_new, x, params, T_new, &residual, &g1);" << endl;
|
2020-06-23 15:13:04 +02:00
|
|
|
output << " break;" << endl;
|
|
|
|
}
|
|
|
|
output << " }" << endl
|
|
|
|
<< endl
|
|
|
|
<< " if (nlhs >= 1)" << endl
|
|
|
|
<< " plhs[0] = residual;" << endl
|
|
|
|
<< " else" << endl
|
|
|
|
<< " mxDestroyArray(residual);" << endl
|
|
|
|
<< " if (nlhs >= 2)" << endl
|
|
|
|
<< " plhs[1] = y_new;" << endl
|
|
|
|
<< " else" << endl
|
|
|
|
<< " mxDestroyArray(y_new);" << endl
|
|
|
|
<< " if (nlhs >= 3)" << endl
|
|
|
|
<< " plhs[2] = T_new;" << endl
|
|
|
|
<< " else" << endl
|
|
|
|
<< " mxDestroyArray(T_new);" << endl
|
|
|
|
<< " if (nlhs >= 4)" << endl
|
|
|
|
<< " plhs[3] = g1;" << endl
|
|
|
|
<< " else" << endl
|
|
|
|
<< " mxDestroyArray(g1);" << endl
|
|
|
|
<< "}" << endl;
|
|
|
|
output.close();
|
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
void
|
2020-05-20 11:44:40 +02:00
|
|
|
StaticModel::writeDriverOutput(ostream &output, bool block) const
|
2009-12-16 14:21:31 +01:00
|
|
|
{
|
2019-04-12 15:41:52 +02:00
|
|
|
output << "M_.static_tmp_nbr = [";
|
2019-12-20 16:59:30 +01:00
|
|
|
for (const auto &temporary_terms_derivative : temporary_terms_derivatives)
|
2019-12-16 19:42:59 +01:00
|
|
|
output << temporary_terms_derivative.size() << "; ";
|
2019-04-12 15:41:52 +02:00
|
|
|
output << "];" << endl;
|
2018-03-28 18:53:02 +02:00
|
|
|
|
2020-06-05 17:11:29 +02:00
|
|
|
/* Write mapping between model local variables and indices in the temporary
|
|
|
|
terms vector (dynare#1722) */
|
|
|
|
output << "M_.model_local_variables_static_tt_idxs = {" << endl;
|
|
|
|
for (auto [mlv, value] : temporary_terms_mlv)
|
|
|
|
output << " '" << symbol_table.getName(mlv->symb_id) << "', "
|
|
|
|
<< temporary_terms_idxs.at(mlv)+1 << ';' << endl;
|
|
|
|
output << "};" << endl;
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
if (!block)
|
|
|
|
return;
|
2009-06-30 17:07:09 +02:00
|
|
|
|
2020-05-20 11:35:14 +02:00
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
2012-06-06 16:36:56 +02:00
|
|
|
{
|
2020-05-20 11:35:14 +02:00
|
|
|
output << "block_structure_stat.block(" << blk+1 << ").Simulation_Type = " << static_cast<int>(blocks[blk].simulation_type) << ";" << endl
|
|
|
|
<< "block_structure_stat.block(" << blk+1 << ").endo_nbr = " << blocks[blk].size << ";" << endl
|
|
|
|
<< "block_structure_stat.block(" << blk+1 << ").mfs = " << blocks[blk].mfs_size << ";" << endl
|
|
|
|
<< "block_structure_stat.block(" << blk+1 << ").equation = [";
|
|
|
|
for (int eq = 0; eq < blocks[blk].size; eq++)
|
|
|
|
output << " " << getBlockEquationID(blk, eq)+1;
|
2020-05-05 17:49:58 +02:00
|
|
|
output << "];" << endl
|
2020-05-20 11:35:14 +02:00
|
|
|
<< "block_structure_stat.block(" << blk+1 << ").variable = [";
|
|
|
|
for (int var = 0; var < blocks[blk].size; var++)
|
|
|
|
output << " " << getBlockVariableID(blk, var)+1;
|
2020-05-05 17:49:58 +02:00
|
|
|
output << "];" << endl;
|
2012-06-06 16:36:56 +02:00
|
|
|
}
|
2020-05-05 17:49:58 +02:00
|
|
|
output << "M_.block_structure_stat.block = block_structure_stat.block;" << endl
|
|
|
|
<< "M_.block_structure_stat.variable_reordered = [";
|
|
|
|
for (int i = 0; i < symbol_table.endo_nbr(); i++)
|
2020-04-17 14:55:55 +02:00
|
|
|
output << " " << endo_idx_block2orig[i]+1;
|
2018-11-23 17:19:59 +01:00
|
|
|
output << "];" << endl
|
|
|
|
<< "M_.block_structure_stat.equation_reordered = [";
|
2020-05-05 17:49:58 +02:00
|
|
|
for (int i = 0; i < symbol_table.endo_nbr(); i++)
|
2020-04-17 14:55:55 +02:00
|
|
|
output << " " << eq_idx_block2orig[i]+1;
|
2018-11-23 17:19:59 +01:00
|
|
|
output << "];" << endl;
|
2017-06-14 07:01:31 +02:00
|
|
|
|
2020-05-05 17:49:58 +02:00
|
|
|
set<pair<int, int>> row_incidence;
|
|
|
|
for (const auto &[indices, d1] : derivatives[1])
|
|
|
|
if (int deriv_id = indices[1];
|
|
|
|
getTypeByDerivID(deriv_id) == SymbolType::endogenous)
|
|
|
|
{
|
|
|
|
int eq = indices[0];
|
2022-07-12 15:28:52 +02:00
|
|
|
int var { getTypeSpecificIDByDerivID(deriv_id) };
|
2020-05-05 17:49:58 +02:00
|
|
|
row_incidence.emplace(eq, var);
|
|
|
|
}
|
2020-05-20 11:35:14 +02:00
|
|
|
output << "M_.block_structure_stat.incidence.sparse_IM = [" << endl;
|
2020-05-05 17:49:58 +02:00
|
|
|
for (auto [eq, var] : row_incidence)
|
2020-05-20 11:35:14 +02:00
|
|
|
output << " " << eq+1 << " " << var+1 << ";" << endl;
|
2020-05-25 18:35:36 +02:00
|
|
|
output << "];" << endl
|
|
|
|
<< "M_.block_structure_stat.tmp_nbr = " << blocks_temporary_terms_idxs.size()
|
|
|
|
<< ";" << endl;
|
2009-04-28 19:11:48 +02:00
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
SymbolType
|
2018-06-04 12:50:53 +02:00
|
|
|
StaticModel::getTypeByDerivID(int deriv_id) const noexcept(false)
|
2009-12-16 14:21:31 +01:00
|
|
|
{
|
2012-11-29 18:07:48 +01:00
|
|
|
if (deriv_id < symbol_table.endo_nbr())
|
2018-07-17 18:34:07 +02:00
|
|
|
return SymbolType::endogenous;
|
2012-11-29 18:07:48 +01:00
|
|
|
else if (deriv_id < symbol_table.endo_nbr() + symbol_table.param_nbr())
|
2018-07-17 18:34:07 +02:00
|
|
|
return SymbolType::parameter;
|
2012-11-29 18:07:48 +01:00
|
|
|
else
|
|
|
|
throw UnknownDerivIDException();
|
2009-04-27 19:15:14 +02:00
|
|
|
}
|
2009-04-30 15:14:33 +02:00
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
int
|
2022-06-24 17:10:12 +02:00
|
|
|
StaticModel::getLagByDerivID([[maybe_unused]] int deriv_id) const noexcept(false)
|
2009-04-30 15:14:33 +02:00
|
|
|
{
|
2009-12-16 14:21:31 +01:00
|
|
|
return 0;
|
2009-04-30 15:14:33 +02:00
|
|
|
}
|
2009-06-30 17:07:09 +02:00
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
int
|
2018-06-04 12:50:53 +02:00
|
|
|
StaticModel::getSymbIDByDerivID(int deriv_id) const noexcept(false)
|
2009-06-30 17:07:09 +02:00
|
|
|
{
|
2012-11-29 18:07:48 +01:00
|
|
|
if (deriv_id < symbol_table.endo_nbr())
|
2018-07-17 18:34:07 +02:00
|
|
|
return symbol_table.getID(SymbolType::endogenous, deriv_id);
|
2012-11-29 18:07:48 +01:00
|
|
|
else if (deriv_id < symbol_table.endo_nbr() + symbol_table.param_nbr())
|
2018-07-17 18:34:07 +02:00
|
|
|
return symbol_table.getID(SymbolType::parameter, deriv_id - symbol_table.endo_nbr());
|
2012-11-29 18:07:48 +01:00
|
|
|
else
|
|
|
|
throw UnknownDerivIDException();
|
2009-06-30 17:07:09 +02:00
|
|
|
}
|
|
|
|
|
2022-07-12 15:28:52 +02:00
|
|
|
int
|
|
|
|
StaticModel::getTypeSpecificIDByDerivID(int deriv_id) const
|
|
|
|
{
|
|
|
|
if (deriv_id < symbol_table.endo_nbr())
|
|
|
|
return deriv_id;
|
|
|
|
else if (deriv_id < symbol_table.endo_nbr() + symbol_table.param_nbr())
|
|
|
|
return deriv_id - symbol_table.endo_nbr();
|
|
|
|
else
|
|
|
|
throw UnknownDerivIDException();
|
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
int
|
2022-06-24 17:10:12 +02:00
|
|
|
StaticModel::getDerivID(int symb_id, [[maybe_unused]] int lag) const noexcept(false)
|
2009-06-30 17:07:09 +02:00
|
|
|
{
|
2018-07-17 18:34:07 +02:00
|
|
|
if (symbol_table.getType(symb_id) == SymbolType::endogenous)
|
2012-11-29 18:07:48 +01:00
|
|
|
return symbol_table.getTypeSpecificID(symb_id);
|
2018-07-17 18:34:07 +02:00
|
|
|
else if (symbol_table.getType(symb_id) == SymbolType::parameter)
|
2012-11-29 18:07:48 +01:00
|
|
|
return symbol_table.getTypeSpecificID(symb_id) + symbol_table.endo_nbr();
|
2009-12-16 14:21:31 +01:00
|
|
|
else
|
2022-06-08 12:45:33 +02:00
|
|
|
/* See the special treatment in VariableNode::prepareForDerivation(),
|
|
|
|
VariableNode::computeDerivative() and VariableNode::getChainRuleDerivative() */
|
|
|
|
throw UnknownDerivIDException{};
|
2009-12-16 14:21:31 +01:00
|
|
|
}
|
2009-06-30 17:07:09 +02:00
|
|
|
|
2012-11-29 18:07:48 +01:00
|
|
|
void
|
|
|
|
StaticModel::addAllParamDerivId(set<int> &deriv_id_set)
|
|
|
|
{
|
|
|
|
for (int i = 0; i < symbol_table.param_nbr(); i++)
|
|
|
|
deriv_id_set.insert(i + symbol_table.endo_nbr());
|
|
|
|
}
|
|
|
|
|
2009-07-07 16:20:48 +02:00
|
|
|
void
|
2020-03-19 17:46:10 +01:00
|
|
|
StaticModel::computeChainRuleJacobian()
|
2009-07-07 16:20:48 +02:00
|
|
|
{
|
2020-04-23 16:04:02 +02:00
|
|
|
int nb_blocks = blocks.size();
|
2020-03-19 17:46:10 +01:00
|
|
|
blocks_derivatives.resize(nb_blocks);
|
2020-04-24 12:29:02 +02:00
|
|
|
for (int blk = 0; blk < nb_blocks; blk++)
|
2009-07-07 16:20:48 +02:00
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
int nb_recursives = blocks[blk].getRecursiveSize();
|
|
|
|
|
2020-10-02 18:31:55 +02:00
|
|
|
map<int, BinaryOpNode *> recursive_vars;
|
2020-04-24 12:29:02 +02:00
|
|
|
for (int i = 0; i < nb_recursives; i++)
|
2009-07-10 16:22:40 +02:00
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
int deriv_id = getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(blk, i)), 0);
|
2020-05-20 11:49:32 +02:00
|
|
|
if (getBlockEquationType(blk, i) == EquationType::evaluateRenormalized)
|
2020-04-24 12:29:02 +02:00
|
|
|
recursive_vars[deriv_id] = getBlockEquationRenormalizedExpr(blk, i);
|
|
|
|
else
|
|
|
|
recursive_vars[deriv_id] = getBlockEquationExpr(blk, i);
|
2009-07-10 16:22:40 +02:00
|
|
|
}
|
2020-04-24 12:29:02 +02:00
|
|
|
|
|
|
|
assert(blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesSimple
|
|
|
|
&& blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesComplete);
|
|
|
|
|
|
|
|
int size = blocks[blk].size;
|
|
|
|
|
|
|
|
for (int eq = nb_recursives; eq < size; eq++)
|
2009-07-07 16:20:48 +02:00
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
int eq_orig = getBlockEquationID(blk, eq);
|
|
|
|
for (int var = nb_recursives; var < size; var++)
|
2009-12-16 14:21:31 +01:00
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
int var_orig = getBlockVariableID(blk, var);
|
|
|
|
expr_t d1 = equations[eq_orig]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, var_orig), 0), recursive_vars);
|
2020-05-06 14:21:33 +02:00
|
|
|
if (d1 != Zero)
|
|
|
|
blocks_derivatives[blk][{ eq, var, 0 }] = d1;
|
2009-07-07 17:58:17 +02:00
|
|
|
}
|
2009-07-07 16:20:48 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
void
|
2017-08-24 15:35:10 +02:00
|
|
|
StaticModel::writeLatexFile(const string &basename, bool write_equation_tags) const
|
2009-12-16 18:13:23 +01:00
|
|
|
{
|
2019-07-11 17:33:53 +02:00
|
|
|
writeLatexModelFile(basename, "static", ExprNodeOutputType::latexStaticModel, write_equation_tags);
|
2009-12-16 18:13:23 +01:00
|
|
|
}
|
2009-07-07 17:58:17 +02:00
|
|
|
|
2009-09-30 17:10:31 +02:00
|
|
|
void
|
2010-04-27 17:04:52 +02:00
|
|
|
StaticModel::writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type) const
|
2009-09-30 17:10:31 +02:00
|
|
|
{
|
2018-06-04 12:26:16 +02:00
|
|
|
for (auto aux_equation : aux_equations)
|
2009-09-30 17:10:31 +02:00
|
|
|
{
|
2018-06-04 12:26:16 +02:00
|
|
|
dynamic_cast<ExprNode *>(aux_equation)->writeOutput(output, output_type);
|
2009-09-30 17:10:31 +02:00
|
|
|
output << ";" << endl;
|
|
|
|
}
|
|
|
|
}
|
2011-07-24 20:52:03 +02:00
|
|
|
|
2017-06-14 07:01:31 +02:00
|
|
|
void
|
2019-12-16 19:42:59 +01:00
|
|
|
StaticModel::writeSetAuxiliaryVariables(const string &basename, bool julia) const
|
2011-07-24 20:52:03 +02:00
|
|
|
{
|
2017-08-29 10:58:39 +02:00
|
|
|
ostringstream output_func_body;
|
2021-01-09 14:45:31 +01:00
|
|
|
ExprNodeOutputType output_type = julia ? ExprNodeOutputType::juliaStaticModel : ExprNodeOutputType::matlabStaticModel;
|
|
|
|
writeAuxVarRecursiveDefinitions(output_func_body, output_type);
|
2017-08-29 10:58:39 +02:00
|
|
|
|
|
|
|
if (output_func_body.str().empty())
|
|
|
|
return;
|
2017-06-14 07:01:31 +02:00
|
|
|
|
2021-08-18 16:52:35 +02:00
|
|
|
string func_name = julia ? "set_auxiliary_variables!" : "set_auxiliary_variables";
|
2015-08-17 15:36:18 +02:00
|
|
|
string comment = julia ? "#" : "%";
|
2011-07-24 20:52:03 +02:00
|
|
|
|
2021-04-23 17:30:38 +02:00
|
|
|
stringstream output;
|
2021-04-23 17:55:28 +02:00
|
|
|
if (julia)
|
|
|
|
output << "module " << basename << "SetAuxiliaryVariables" << endl
|
|
|
|
<< "export " << func_name << endl;
|
2021-01-09 14:45:31 +01:00
|
|
|
output << "function ";
|
|
|
|
if (!julia)
|
|
|
|
output << "y = ";
|
2021-04-23 17:55:28 +02:00
|
|
|
output << func_name << "(y, x, params)" << endl
|
2015-08-17 15:36:18 +02:00
|
|
|
<< comment << endl
|
|
|
|
<< comment << " Status : Computes static model for Dynare" << endl
|
|
|
|
<< comment << endl
|
|
|
|
<< comment << " Warning : this file is generated automatically by Dynare" << endl
|
2022-07-12 14:31:30 +02:00
|
|
|
<< comment << " from model file (.mod)" << endl << endl;
|
|
|
|
if (julia)
|
|
|
|
output << "@inbounds begin" << endl;
|
|
|
|
output << output_func_body.str()
|
2021-04-23 17:55:28 +02:00
|
|
|
<< "end" << endl;
|
2021-01-09 14:45:31 +01:00
|
|
|
if (julia)
|
2022-07-12 14:31:30 +02:00
|
|
|
output << "end" << endl
|
|
|
|
<< "end" << endl;
|
2018-03-02 18:39:16 +01:00
|
|
|
|
2021-04-23 17:55:28 +02:00
|
|
|
writeToFileIfModified(output, julia ? basename + "SetAuxiliaryVariables.jl" : packageDir(basename) + "/" + func_name + ".m");
|
2016-04-04 17:11:03 +02:00
|
|
|
}
|
2015-05-12 16:36:03 +02:00
|
|
|
|
2016-04-04 17:11:03 +02:00
|
|
|
void
|
|
|
|
StaticModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const
|
|
|
|
{
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
2018-06-04 12:26:16 +02:00
|
|
|
for (auto aux_equation : aux_equations)
|
|
|
|
if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
|
2018-09-05 18:27:13 +02:00
|
|
|
dynamic_cast<ExprNode *>(aux_equation)->writeExternalFunctionOutput(output, ExprNodeOutputType::matlabStaticModel, {}, {}, tef_terms);
|
2018-06-04 12:26:16 +02:00
|
|
|
for (auto aux_equation : aux_equations)
|
2011-07-24 20:52:03 +02:00
|
|
|
{
|
2018-06-04 12:26:16 +02:00
|
|
|
dynamic_cast<ExprNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->writeOutput(output, output_type);
|
2011-07-24 20:52:03 +02:00
|
|
|
output << ";" << endl;
|
|
|
|
}
|
|
|
|
}
|
2012-11-29 18:07:48 +01:00
|
|
|
|
2017-08-30 11:32:01 +02:00
|
|
|
void
|
|
|
|
StaticModel::writeLatexAuxVarRecursiveDefinitions(ostream &output) const
|
|
|
|
{
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
temporary_terms_t temporary_terms;
|
2019-12-20 16:59:30 +01:00
|
|
|
temporary_terms_idxs_t temporary_terms_idxs;
|
2018-06-04 12:26:16 +02:00
|
|
|
for (auto aux_equation : aux_equations)
|
|
|
|
if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
|
2018-09-05 18:27:13 +02:00
|
|
|
dynamic_cast<ExprNode *>(aux_equation)->writeExternalFunctionOutput(output, ExprNodeOutputType::latexStaticModel,
|
2019-12-20 16:59:30 +01:00
|
|
|
temporary_terms, temporary_terms_idxs, tef_terms);
|
2018-06-04 12:26:16 +02:00
|
|
|
for (auto aux_equation : aux_equations)
|
2017-08-30 11:32:01 +02:00
|
|
|
{
|
2019-04-03 16:32:52 +02:00
|
|
|
output << R"(\begin{dmath})" << endl;
|
2018-09-05 18:27:13 +02:00
|
|
|
dynamic_cast<ExprNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->writeOutput(output, ExprNodeOutputType::latexStaticModel);
|
2019-04-03 16:32:52 +02:00
|
|
|
output << endl << R"(\end{dmath})" << endl;
|
2017-08-30 11:32:01 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2017-10-16 17:24:55 +02:00
|
|
|
void
|
|
|
|
StaticModel::writeJsonAuxVarRecursiveDefinitions(ostream &output) const
|
|
|
|
{
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
temporary_terms_t temporary_terms;
|
|
|
|
|
2018-06-04 12:26:16 +02:00
|
|
|
for (auto aux_equation : aux_equations)
|
|
|
|
if (dynamic_cast<ExprNode *>(aux_equation)->containsExternalFunction())
|
2017-10-16 17:24:55 +02:00
|
|
|
{
|
|
|
|
vector<string> efout;
|
2018-06-04 12:26:16 +02:00
|
|
|
dynamic_cast<ExprNode *>(aux_equation)->writeJsonExternalFunctionOutput(efout,
|
2019-12-16 19:42:59 +01:00
|
|
|
temporary_terms,
|
|
|
|
tef_terms,
|
|
|
|
false);
|
2022-06-03 16:24:26 +02:00
|
|
|
for (bool printed_something{false};
|
|
|
|
const auto &it : efout)
|
2017-10-16 17:24:55 +02:00
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
if (exchange(printed_something, true))
|
2017-10-16 17:24:55 +02:00
|
|
|
output << ", ";
|
2022-06-03 16:24:26 +02:00
|
|
|
output << it;
|
2017-10-16 17:24:55 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2018-06-04 12:26:16 +02:00
|
|
|
for (auto aux_equation : aux_equations)
|
2017-10-16 17:24:55 +02:00
|
|
|
{
|
2019-04-03 16:32:52 +02:00
|
|
|
output << R"(, {"lhs": ")";
|
2018-11-28 14:32:26 +01:00
|
|
|
aux_equation->arg1->writeJsonOutput(output, temporary_terms, tef_terms, false);
|
2019-04-03 16:32:52 +02:00
|
|
|
output << R"(", "rhs": ")";
|
2018-11-28 14:32:26 +01:00
|
|
|
dynamic_cast<BinaryOpNode *>(aux_equation->substituteStaticAuxiliaryDefinition())->arg2->writeJsonOutput(output, temporary_terms, tef_terms, false);
|
2019-04-03 16:32:52 +02:00
|
|
|
output << R"("})";
|
2017-10-16 17:24:55 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2017-02-27 15:40:34 +01:00
|
|
|
void
|
|
|
|
StaticModel::writeJsonOutput(ostream &output) const
|
|
|
|
{
|
2020-06-05 16:07:52 +02:00
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
writeJsonModelLocalVariables(output, false, tef_terms);
|
|
|
|
output << ", ";
|
2017-02-27 15:40:34 +01:00
|
|
|
writeJsonModelEquations(output, false);
|
|
|
|
}
|
|
|
|
|
2017-02-20 12:18:11 +01:00
|
|
|
void
|
2017-03-02 18:34:18 +01:00
|
|
|
StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) const
|
2017-02-20 12:18:11 +01:00
|
|
|
{
|
2019-12-20 16:59:30 +01:00
|
|
|
ostringstream model_local_vars_output; // Used for storing model local vars
|
|
|
|
vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
|
2017-02-20 12:18:11 +01:00
|
|
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
2019-04-18 14:34:48 +02:00
|
|
|
temporary_terms_t temp_term_union;
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2020-06-05 16:07:52 +02:00
|
|
|
writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms, "");
|
|
|
|
d_output[0] << ", ";
|
|
|
|
writeJsonModelEquations(d_output[0], true);
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
int ncols = symbol_table.endo_nbr();
|
|
|
|
for (size_t i = 1; i < derivatives.size(); i++)
|
2017-02-20 12:18:11 +01:00
|
|
|
{
|
2019-04-18 17:07:55 +02:00
|
|
|
string matrix_name = i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative";
|
|
|
|
writeJsonTemporaryTerms(temporary_terms_derivatives[i], temp_term_union, d_output[i], tef_terms, matrix_name);
|
|
|
|
temp_term_union.insert(temporary_terms_derivatives[i].begin(), temporary_terms_derivatives[i].end());
|
|
|
|
d_output[i] << R"(, ")" << matrix_name << R"(": {)"
|
|
|
|
<< R"( "nrows": )" << equations.size()
|
|
|
|
<< R"(, "ncols": )" << ncols
|
|
|
|
<< R"(, "entries": [)";
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2022-06-03 16:24:26 +02:00
|
|
|
for (bool printed_something{false};
|
|
|
|
const auto &[vidx, d] : derivatives[i])
|
2019-04-18 17:07:55 +02:00
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
if (exchange(printed_something, true))
|
2019-04-18 17:07:55 +02:00
|
|
|
d_output[i] << ", ";
|
2017-03-02 18:34:18 +01:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
int eq = vidx[0];
|
2017-03-02 18:34:18 +01:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
int col_idx = 0;
|
|
|
|
for (size_t j = 1; j < vidx.size(); j++)
|
|
|
|
{
|
|
|
|
col_idx *= symbol_table.endo_nbr();
|
2022-07-12 15:18:36 +02:00
|
|
|
col_idx += getJacobianCol(vidx[j]);
|
2019-04-18 17:07:55 +02:00
|
|
|
}
|
2017-03-02 18:34:18 +01:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
if (writeDetails)
|
|
|
|
d_output[i] << R"({"eq": )" << eq + 1;
|
|
|
|
else
|
|
|
|
d_output[i] << R"({"row": )" << eq + 1;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
d_output[i] << R"(, "col": )" << (i > 1 ? "[" : "") << col_idx + 1;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian
|
|
|
|
{
|
2022-07-12 15:18:36 +02:00
|
|
|
int col_idx_sym = getJacobianCol(vidx[2]) * symbol_table.endo_nbr() + getJacobianCol(vidx[1]);
|
2019-04-18 17:07:55 +02:00
|
|
|
d_output[i] << ", " << col_idx_sym + 1;
|
|
|
|
}
|
|
|
|
if (i > 1)
|
|
|
|
d_output[i] << "]";
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
if (writeDetails)
|
|
|
|
for (size_t j = 1; j < vidx.size(); j++)
|
2022-07-12 15:28:52 +02:00
|
|
|
d_output[i] << R"(, "var)" << (i > 1 ? to_string(j) : "") << R"(": ")" << getNameByDerivID(vidx[j]) << R"(")";
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2019-04-18 17:07:55 +02:00
|
|
|
d_output[i] << R"(, "val": ")";
|
|
|
|
d->writeJsonOutput(d_output[i], temp_term_union, tef_terms);
|
|
|
|
d_output[i] << R"("})" << endl;
|
2017-03-02 18:34:18 +01:00
|
|
|
}
|
2019-12-20 16:59:30 +01:00
|
|
|
d_output[i] << "]}";
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-12-20 16:59:30 +01:00
|
|
|
ncols *= symbol_table.endo_nbr();
|
2017-02-20 12:18:11 +01:00
|
|
|
}
|
|
|
|
|
2017-03-02 18:34:18 +01:00
|
|
|
if (writeDetails)
|
2019-04-03 16:32:52 +02:00
|
|
|
output << R"("static_model": {)";
|
2017-03-02 18:34:18 +01:00
|
|
|
else
|
2019-04-03 16:32:52 +02:00
|
|
|
output << R"("static_model_simple": {)";
|
2019-04-18 17:07:55 +02:00
|
|
|
output << model_local_vars_output.str();
|
|
|
|
for (const auto &it : d_output)
|
|
|
|
output << ", " << it.str();
|
|
|
|
output << "}";
|
2017-02-20 12:18:11 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2017-03-02 18:34:18 +01:00
|
|
|
StaticModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) const
|
2017-02-20 12:18:11 +01:00
|
|
|
{
|
2018-11-15 16:39:53 +01:00
|
|
|
if (!params_derivatives.size())
|
2017-02-20 12:18:11 +01:00
|
|
|
return;
|
|
|
|
|
2019-12-20 16:59:30 +01:00
|
|
|
ostringstream model_local_vars_output; // Used for storing model local vars
|
|
|
|
ostringstream model_output; // Used for storing model temp vars and equations
|
|
|
|
ostringstream jacobian_output; // Used for storing jacobian equations
|
|
|
|
ostringstream hessian_output; // Used for storing Hessian equations
|
|
|
|
ostringstream hessian1_output; // Used for storing Hessian equations
|
|
|
|
ostringstream third_derivs_output; // Used for storing third order derivatives equations
|
|
|
|
ostringstream third_derivs1_output; // Used for storing third order derivatives equations
|
2017-02-20 12:18:11 +01:00
|
|
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
2020-06-05 16:07:52 +02:00
|
|
|
writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2018-11-16 18:24:06 +01:00
|
|
|
temporary_terms_t temp_term_union;
|
2018-11-30 12:22:13 +01:00
|
|
|
for (const auto &it : params_derivs_temporary_terms)
|
2019-04-18 14:34:48 +02:00
|
|
|
writeJsonTemporaryTerms(it.second, temp_term_union, model_output, tef_terms, "all");
|
2018-11-16 18:24:06 +01:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
jacobian_output << R"("deriv_wrt_params": {)"
|
|
|
|
<< R"( "neqs": )" << equations.size()
|
|
|
|
<< R"(, "nparamcols": )" << symbol_table.param_nbr()
|
|
|
|
<< R"(, "entries": [)";
|
2022-06-03 16:24:26 +02:00
|
|
|
for (bool printed_something{false};
|
|
|
|
const auto &[vidx, d] : params_derivatives.find({ 0, 1 })->second)
|
2017-02-20 12:18:11 +01:00
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
if (exchange(printed_something, true))
|
2017-02-20 12:18:11 +01:00
|
|
|
jacobian_output << ", ";
|
|
|
|
|
2022-06-03 16:24:26 +02:00
|
|
|
auto [eq, param] = vectorToTuple<2>(vidx);
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2022-07-12 15:28:52 +02:00
|
|
|
int param_col { getTypeSpecificIDByDerivID(param) + 1 };
|
2017-03-02 18:34:18 +01:00
|
|
|
|
|
|
|
if (writeDetails)
|
2019-04-03 16:32:52 +02:00
|
|
|
jacobian_output << R"({"eq": )" << eq + 1;
|
2017-03-02 18:34:18 +01:00
|
|
|
else
|
2019-04-03 16:32:52 +02:00
|
|
|
jacobian_output << R"({"row": )" << eq + 1;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
|
|
|
if (writeDetails)
|
2019-04-03 16:32:52 +02:00
|
|
|
jacobian_output << R"(, "param_col": )" << param_col;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2022-07-12 15:28:52 +02:00
|
|
|
jacobian_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
jacobian_output << R"(, "val": ")";
|
2022-06-03 16:24:26 +02:00
|
|
|
d->writeJsonOutput(jacobian_output, temp_term_union, tef_terms);
|
2019-04-03 16:32:52 +02:00
|
|
|
jacobian_output << R"("})" << endl;
|
2017-02-20 12:18:11 +01:00
|
|
|
}
|
|
|
|
jacobian_output << "]}";
|
2018-11-15 16:39:53 +01:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian_output << R"("deriv_jacobian_wrt_params": {)"
|
|
|
|
<< R"( "neqs": )" << equations.size()
|
|
|
|
<< R"(, "nvarcols": )" << symbol_table.endo_nbr()
|
|
|
|
<< R"(, "nparamcols": )" << symbol_table.param_nbr()
|
|
|
|
<< R"(, "entries": [)";
|
2022-06-03 16:24:26 +02:00
|
|
|
for (bool printed_something{false};
|
|
|
|
const auto &[vidx, d] : params_derivatives.find({ 1, 1 })->second)
|
2017-02-20 12:18:11 +01:00
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
if (exchange(printed_something, true))
|
2017-02-20 12:18:11 +01:00
|
|
|
hessian_output << ", ";
|
|
|
|
|
2022-06-03 16:24:26 +02:00
|
|
|
auto [eq, var, param] = vectorToTuple<3>(vidx);
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2022-07-12 15:28:52 +02:00
|
|
|
int var_col { getTypeSpecificIDByDerivID(var) + 1 };
|
|
|
|
int param_col { getTypeSpecificIDByDerivID(param) + 1 };
|
2017-03-02 18:34:18 +01:00
|
|
|
|
|
|
|
if (writeDetails)
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian_output << R"({"eq": )" << eq + 1;
|
2017-03-02 18:34:18 +01:00
|
|
|
else
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian_output << R"({"row": )" << eq + 1;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
|
|
|
if (writeDetails)
|
2022-07-12 15:28:52 +02:00
|
|
|
hessian_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")"
|
|
|
|
<< R"(, "param": ")" << getNameByDerivID(param) << R"(")";
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian_output << R"(, "var_col": )" << var_col
|
|
|
|
<< R"(, "param_col": )" << param_col
|
|
|
|
<< R"(, "val": ")";
|
2022-06-03 16:24:26 +02:00
|
|
|
d->writeJsonOutput(hessian_output, temp_term_union, tef_terms);
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian_output << R"("})" << endl;
|
2017-02-20 12:18:11 +01:00
|
|
|
}
|
|
|
|
hessian_output << "]}";
|
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian1_output << R"("second_deriv_residuals_wrt_params": {)"
|
|
|
|
<< R"( "nrows": )" << equations.size()
|
|
|
|
<< R"(, "nparam1cols": )" << symbol_table.param_nbr()
|
|
|
|
<< R"(, "nparam2cols": )" << symbol_table.param_nbr()
|
|
|
|
<< R"(, "entries": [)";
|
2022-06-03 16:24:26 +02:00
|
|
|
for (bool printed_something{false};
|
|
|
|
const auto &[vidx, d] : params_derivatives.find({ 0, 2 })->second)
|
2017-02-20 12:18:11 +01:00
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
if (exchange(printed_something, true))
|
2017-02-20 12:18:11 +01:00
|
|
|
hessian1_output << ", ";
|
|
|
|
|
2022-06-03 16:24:26 +02:00
|
|
|
auto [eq, param1, param2] = vectorToTuple<3>(vidx);
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2022-07-12 15:28:52 +02:00
|
|
|
int param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
|
|
|
|
int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
|
2017-03-02 18:34:18 +01:00
|
|
|
|
|
|
|
if (writeDetails)
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian1_output << R"({"eq": )" << eq + 1;
|
2017-03-02 18:34:18 +01:00
|
|
|
else
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian1_output << R"({"row": )" << eq + 1;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian1_output << R"(, "param1_col": )" << param1_col
|
|
|
|
<< R"(, "param2_col": )" << param2_col;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
|
|
|
if (writeDetails)
|
2022-07-12 15:28:52 +02:00
|
|
|
hessian1_output << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
|
|
|
|
<< R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian1_output << R"(, "val": ")";
|
2022-06-03 16:24:26 +02:00
|
|
|
d->writeJsonOutput(hessian1_output, temp_term_union, tef_terms);
|
2019-04-03 16:32:52 +02:00
|
|
|
hessian1_output << R"("})" << endl;
|
2017-02-20 12:18:11 +01:00
|
|
|
}
|
|
|
|
hessian1_output << "]}";
|
2018-11-15 16:39:53 +01:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs_output << R"("second_deriv_jacobian_wrt_params": {)"
|
|
|
|
<< R"( "neqs": )" << equations.size()
|
|
|
|
<< R"(, "nvarcols": )" << symbol_table.endo_nbr()
|
|
|
|
<< R"(, "nparam1cols": )" << symbol_table.param_nbr()
|
|
|
|
<< R"(, "nparam2cols": )" << symbol_table.param_nbr()
|
|
|
|
<< R"(, "entries": [)";
|
2022-06-03 16:24:26 +02:00
|
|
|
for (bool printed_something{false};
|
|
|
|
const auto &[vidx, d] : params_derivatives.find({ 1, 2 })->second)
|
2017-02-20 12:18:11 +01:00
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
if (exchange(printed_something, true))
|
2017-02-20 12:18:11 +01:00
|
|
|
third_derivs_output << ", ";
|
|
|
|
|
2022-06-03 16:24:26 +02:00
|
|
|
auto [eq, var, param1, param2] = vectorToTuple<4>(vidx);
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2022-07-12 15:28:52 +02:00
|
|
|
int var_col { getTypeSpecificIDByDerivID(var) + 1 };
|
|
|
|
int param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
|
|
|
|
int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
|
2017-03-02 18:34:18 +01:00
|
|
|
|
|
|
|
if (writeDetails)
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs_output << R"({"eq": )" << eq + 1;
|
2017-03-02 18:34:18 +01:00
|
|
|
else
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs_output << R"({"row": )" << eq + 1;
|
|
|
|
third_derivs_output << R"(, "var_col": )" << var_col
|
|
|
|
<< R"(, "param1_col": )" << param1_col
|
|
|
|
<< R"(, "param2_col": )" << param2_col;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
|
|
|
if (writeDetails)
|
2022-07-12 15:28:52 +02:00
|
|
|
third_derivs_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")"
|
|
|
|
<< R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
|
|
|
|
<< R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs_output << R"(, "val": ")";
|
2022-06-03 16:24:26 +02:00
|
|
|
d->writeJsonOutput(third_derivs_output, temp_term_union, tef_terms);
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs_output << R"("})" << endl;
|
2017-02-20 12:18:11 +01:00
|
|
|
}
|
|
|
|
third_derivs_output << "]}" << endl;
|
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs1_output << R"("derivative_hessian_wrt_params": {)"
|
|
|
|
<< R"( "neqs": )" << equations.size()
|
|
|
|
<< R"(, "nvar1cols": )" << symbol_table.endo_nbr()
|
|
|
|
<< R"(, "nvar2cols": )" << symbol_table.endo_nbr()
|
|
|
|
<< R"(, "nparamcols": )" << symbol_table.param_nbr()
|
|
|
|
<< R"(, "entries": [)";
|
2022-06-03 16:24:26 +02:00
|
|
|
for (bool printed_something{false};
|
|
|
|
const auto &[vidx, d] : params_derivatives.find({ 2, 1 })->second)
|
2017-02-20 12:18:11 +01:00
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
if (exchange(printed_something, true))
|
2017-02-20 12:18:11 +01:00
|
|
|
third_derivs1_output << ", ";
|
|
|
|
|
2022-06-03 16:24:26 +02:00
|
|
|
auto [eq, var1, var2, param] = vectorToTuple<4>(vidx);
|
2017-02-20 12:18:11 +01:00
|
|
|
|
2022-07-12 15:28:52 +02:00
|
|
|
int var1_col { getTypeSpecificIDByDerivID(var1) + 1 };
|
|
|
|
int var2_col { getTypeSpecificIDByDerivID(var2) + 1 };
|
|
|
|
int param_col { getTypeSpecificIDByDerivID(param) + 1 };
|
2017-03-02 18:34:18 +01:00
|
|
|
|
|
|
|
if (writeDetails)
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs1_output << R"({"eq": )" << eq + 1;
|
2017-03-02 18:34:18 +01:00
|
|
|
else
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs1_output << R"({"row": )" << eq + 1;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs1_output << R"(, "var1_col": )" << var1_col
|
|
|
|
<< R"(, "var2_col": )" << var2_col
|
|
|
|
<< R"(, "param_col": )" << param_col;
|
2017-06-29 12:42:28 +02:00
|
|
|
|
|
|
|
if (writeDetails)
|
2022-07-12 15:28:52 +02:00
|
|
|
third_derivs1_output << R"(, "var1": ")" << getNameByDerivID(var1) << R"(")"
|
|
|
|
<< R"(, "var2": ")" << getNameByDerivID(var2) << R"(")"
|
|
|
|
<< R"(, "param1": ")" << getNameByDerivID(param) << R"(")";
|
2017-06-29 12:42:28 +02:00
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs1_output << R"(, "val": ")";
|
2022-06-03 16:24:26 +02:00
|
|
|
d->writeJsonOutput(third_derivs1_output, temp_term_union, tef_terms);
|
2019-04-03 16:32:52 +02:00
|
|
|
third_derivs1_output << R"("})" << endl;
|
2017-02-20 12:18:11 +01:00
|
|
|
}
|
|
|
|
third_derivs1_output << "]}" << endl;
|
|
|
|
|
2017-03-02 18:34:18 +01:00
|
|
|
if (writeDetails)
|
2019-04-03 16:32:52 +02:00
|
|
|
output << R"("static_model_params_derivative": {)";
|
2017-03-02 18:34:18 +01:00
|
|
|
else
|
2019-04-03 16:32:52 +02:00
|
|
|
output << R"("static_model_params_derivatives_simple": {)";
|
2017-03-02 18:34:18 +01:00
|
|
|
output << model_local_vars_output.str()
|
2017-02-20 12:18:11 +01:00
|
|
|
<< ", " << model_output.str()
|
|
|
|
<< ", " << jacobian_output.str()
|
|
|
|
<< ", " << hessian_output.str()
|
|
|
|
<< ", " << hessian1_output.str()
|
|
|
|
<< ", " << third_derivs_output.str()
|
|
|
|
<< ", " << third_derivs1_output.str()
|
|
|
|
<< "}";
|
|
|
|
}
|