...
 
/*
* Copyright © 2007-2019 Dynare Team
* Copyright © 2007-2020 Dynare Team
*
* This file is part of Dynare.
*
......@@ -105,33 +105,33 @@ enum Tags
};
enum BlockType
enum class BlockType
{
SIMULTANS, //!< Simultaneous time separable block
PROLOGUE, //!< Prologue block (one equation at the beginning, later merged)
EPILOGUE, //!< Epilogue block (one equation at the beginning, later merged)
SIMULTAN //!< Simultaneous time unseparable block
simultans, //!< Simultaneous time separable block
prologue, //!< Prologue block (one equation at the beginning, later merged)
epilogue, //!< Epilogue block (one equation at the beginning, later merged)
simultan //!< Simultaneous time unseparable block
};
enum EquationType
enum class EquationType
{
E_UNKNOWN, //!< Unknown equation type
E_EVALUATE, //!< Simple evaluation, normalized variable on left-hand side
E_EVALUATE_S, //!< Simple evaluation, normalize using the first order derivative
E_SOLVE //!< No simple evaluation of the equation, it has to be solved
unknown, //!< Unknown equation type
evaluate, //!< Simple evaluation, normalized variable on left-hand side
evaluate_s, //!< Simple evaluation, normalize using the first order derivative
solve //!< No simple evaluation of the equation, it has to be solved
};
enum BlockSimulationType
enum class BlockSimulationType
{
UNKNOWN, //!< Unknown simulation type
EVALUATE_FORWARD, //!< Simple evaluation, normalized variable on left-hand side, forward
EVALUATE_BACKWARD, //!< Simple evaluation, normalized variable on left-hand side, backward
SOLVE_FORWARD_SIMPLE, //!< Block of one equation, newton solver needed, forward
SOLVE_BACKWARD_SIMPLE, //!< Block of one equation, newton solver needed, backward
SOLVE_TWO_BOUNDARIES_SIMPLE, //!< Block of one equation, newton solver needed, forward & ackward
SOLVE_FORWARD_COMPLETE, //!< Block of several equations, newton solver needed, forward
SOLVE_BACKWARD_COMPLETE, //!< Block of several equations, newton solver needed, backward
SOLVE_TWO_BOUNDARIES_COMPLETE //!< Block of several equations, newton solver needed, forward and backwar
unknown, //!< Unknown simulation type
evaluateForward, //!< Simple evaluation, normalized variable on left-hand side, forward
evaluateBackward, //!< Simple evaluation, normalized variable on left-hand side, backward
solveForwardSimple, //!< Block of one equation, newton solver needed, forward
solveBackwardSimple, //!< Block of one equation, newton solver needed, backward
solveTwoBoundariesSimple, //!< Block of one equation, newton solver needed, forward & ackward
solveForwardComplete, //!< Block of several equations, newton solver needed, forward
solveBackwardComplete, //!< Block of several equations, newton solver needed, backward
solveTwoBoundariesComplete //!< Block of several equations, newton solver needed, forward and backwar
};
//! Enumeration of possible symbol types
......@@ -1428,7 +1428,7 @@ private:
unsigned int nb_col_det_exo_jacob, nb_col_exo_jacob, nb_col_other_endo_jacob;
public:
inline
FBEGINBLOCK_() : size{0}, type{UNKNOWN},
FBEGINBLOCK_() : size{0}, type{static_cast<uint8_t>(BlockSimulationType::unknown)},
is_linear{false}, endo_nbr{0}, Max_Lag{0}, Max_Lead{0}, u_count_int{0}, nb_col_jacob{0}
{
}
......@@ -1577,8 +1577,10 @@ public:
CompileCode.write(reinterpret_cast<char *>(&variable[i]), sizeof(variable[0]));
CompileCode.write(reinterpret_cast<char *>(&equation[i]), sizeof(equation[0]));
}
if (type == SOLVE_TWO_BOUNDARIES_SIMPLE || type == SOLVE_TWO_BOUNDARIES_COMPLETE
|| type == SOLVE_BACKWARD_COMPLETE || type == SOLVE_FORWARD_COMPLETE)
if (type == static_cast<uint8_t>(BlockSimulationType::solveTwoBoundariesSimple)
|| type == static_cast<uint8_t>(BlockSimulationType::solveTwoBoundariesComplete)
|| type == static_cast<uint8_t>(BlockSimulationType::solveBackwardComplete)
|| type == static_cast<uint8_t>(BlockSimulationType::solveForwardComplete))
{
CompileCode.write(reinterpret_cast<char *>(&is_linear), sizeof(is_linear));
CompileCode.write(reinterpret_cast<char *>(&endo_nbr), sizeof(endo_nbr));
......@@ -1617,8 +1619,10 @@ public:
memcpy(&bc.Equation, code, sizeof(bc.Equation)); code += sizeof(bc.Equation);
Block_Contain_.push_back(bc);
}
if (type == SOLVE_TWO_BOUNDARIES_SIMPLE || type == SOLVE_TWO_BOUNDARIES_COMPLETE
|| type == SOLVE_BACKWARD_COMPLETE || type == SOLVE_FORWARD_COMPLETE)
if (type == static_cast<uint8_t>(BlockSimulationType::solveTwoBoundariesSimple)
|| type == static_cast<uint8_t>(BlockSimulationType::solveTwoBoundariesComplete)
|| type == static_cast<uint8_t>(BlockSimulationType::solveBackwardComplete)
|| type == static_cast<uint8_t>(BlockSimulationType::solveForwardComplete))
{
memcpy(&is_linear, code, sizeof(is_linear)); code += sizeof(is_linear);
memcpy(&endo_nbr, code, sizeof(endo_nbr)); code += sizeof(endo_nbr);
......
......@@ -419,31 +419,40 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
<< "% Warning : this file is generated automatically by Dynare" << endl
<< "% from model file (.mod)" << endl << endl
<< "%/" << endl;
if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
if (simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
{
output << "function [y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)" << endl;
}
else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
else if (simulation_type == BlockSimulationType::solveForwardComplete
|| simulation_type == BlockSimulationType::solveBackwardComplete)
output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE)
else if (simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple)
output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
else
output << "function [residual, y, g1, g2, g3, b, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)" << endl;
BlockType block_type;
if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
block_type = SIMULTAN;
else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
block_type = SIMULTANS;
else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
|| simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
block_type = BlockType::simultan;
else if (simulation_type == BlockSimulationType::solveForwardComplete
|| simulation_type == BlockSimulationType::solveBackwardComplete)
block_type = BlockType::simultans;
else if ((simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
&& getBlockFirstEquation(block) < prologue)
block_type = PROLOGUE;
else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
|| simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
block_type = BlockType::prologue;
else if ((simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
&& getBlockFirstEquation(block) >= equations.size() - epilogue)
block_type = EPILOGUE;
block_type = BlockType::epilogue;
else
block_type = SIMULTANS;
block_type = BlockType::simultans;
output << " % ////////////////////////////////////////////////////////////////////////" << endl
<< " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
<< " //" << endl
......@@ -451,7 +460,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
<< BlockSim(simulation_type) << " //" << endl
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
//The Temporary terms
if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
if (simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
{
output << " if(jacobian_eval)" << endl
<< " g1 = spalloc(" << block_mfs << ", " << count_col_endo << ", " << nze << ");" << endl
......@@ -468,7 +478,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
<< " g1_xd=spalloc(" << block_size << ", " << count_col_exo_det << ", " << nze_exo_det << ");" << endl
<< " g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");" << endl
<< " else" << endl;
if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
output << " g1 = spalloc(" << block_mfs << "*Periods, "
<< block_mfs << "*(Periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
<< ", " << nze << "*Periods);" << endl;
......@@ -486,7 +497,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
tmp_output << " T" << it;
output << " global" << tmp_output.str() << ";" << endl;
}
if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
{
temporary_terms_t tt2;
for (int i = 0; i < static_cast<int>(block_size); i++)
......@@ -505,16 +517,21 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
}
}
}
if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
if (simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete)
output << " residual=zeros(" << block_mfs << ",1);" << endl;
else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
else if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
output << " residual=zeros(" << block_mfs << ",y_kmin+periods);" << endl;
if (simulation_type == EVALUATE_BACKWARD)
if (simulation_type == BlockSimulationType::evaluateBackward)
output << " for it_ = (y_kmin+periods):y_kmin+1" << endl;
if (simulation_type == EVALUATE_FORWARD)
if (simulation_type == BlockSimulationType::evaluateForward)
output << " for it_ = y_kmin+1:(y_kmin+periods)" << endl;
if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
{
output << " b = zeros(periods*y_size,1);" << endl
<< " for it_ = y_kmin+1:(periods+y_kmin)" << endl
......@@ -524,7 +541,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
sps = " ";
}
else
if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
if (simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
sps = " ";
else
sps = "";
......@@ -562,19 +580,20 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
switch (simulation_type)
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
evaluation: if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
case BlockSimulationType::evaluateBackward:
case BlockSimulationType::evaluateForward:
evaluation: if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
output << " % equation " << getBlockEquationID(block, i)+1 << " variable : " << sModel
<< " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
output << " ";
if (equ_type == E_EVALUATE)
if (equ_type == EquationType::evaluate)
{
output << tmp_output.str();
output << " = ";
rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
}
else if (equ_type == E_EVALUATE_S)
else if (equ_type == EquationType::evaluate_s)
{
output << "%" << tmp_output.str();
output << " = ";
......@@ -598,10 +617,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
}
output << ";" << endl;
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
case BlockSimulationType::solveBackwardSimple:
case BlockSimulationType::solveForwardSimple:
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
if (i < block_recursive)
goto evaluation;
feedback_variables.push_back(variable_ID);
......@@ -609,8 +628,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
<< " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(SymbolType::endogenous, variable_ID) << endl;
output << " " << "residual(" << i+1-block_recursive << ") = (";
goto end;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case BlockSimulationType::solveTwoBoundariesComplete:
case BlockSimulationType::solveTwoBoundariesSimple:
if (i < block_recursive)
goto evaluation;
feedback_variables.push_back(variable_ID);
......@@ -628,17 +647,21 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
output << ");" << endl;
#ifdef CONDITION
if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
output << " condition(" << i+1 << ")=0;" << endl;
#endif
}
}
// The Jacobian if we have to solve the block
if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple
|| simulation_type == BlockSimulationType::solveTwoBoundariesComplete)
output << " " << sps << "% Jacobian " << endl << " if jacobian_eval" << endl;
else
if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE
|| simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
if (simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete)
output << " % Jacobian " << endl << " if jacobian_eval" << endl;
else
output << " % Jacobian " << endl << " if jacobian_eval" << endl;
......@@ -744,15 +767,15 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
switch (simulation_type)
{
case EVALUATE_FORWARD:
case EVALUATE_BACKWARD:
case BlockSimulationType::evaluateForward:
case BlockSimulationType::evaluateBackward:
output << " end;" << endl
<< " end;" << endl;
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
case BlockSimulationType::solveBackwardSimple:
case BlockSimulationType::solveForwardSimple:
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
output << " else" << endl;
for (const auto &it : blocks_derivatives[block])
{
......@@ -775,8 +798,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
}
output << " end;" << endl;
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case BlockSimulationType::solveTwoBoundariesSimple:
case BlockSimulationType::solveTwoBoundariesComplete:
output << " else" << endl;
for (const auto &it : blocks_derivatives[block])
{
......@@ -892,13 +915,13 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
int u_count_int = 0;
BlockSimulationType simulation_type;
if ((max_endo_lag > 0) && (max_endo_lead > 0))
simulation_type = SOLVE_TWO_BOUNDARIES_COMPLETE;
simulation_type = BlockSimulationType::solveTwoBoundariesComplete;
else if ((max_endo_lag >= 0) && (max_endo_lead == 0))
simulation_type = SOLVE_FORWARD_COMPLETE;
simulation_type = BlockSimulationType::solveForwardComplete;
else
simulation_type = SOLVE_BACKWARD_COMPLETE;
simulation_type = BlockSimulationType::solveBackwardComplete;
Write_Inf_To_Bin_File(basename + "/model/bytecode/dynamic.bin", u_count_int, file_open, simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE, symbol_table.endo_nbr());
Write_Inf_To_Bin_File(basename + "/model/bytecode/dynamic.bin", u_count_int, file_open, simulation_type == BlockSimulationType::solveTwoBoundariesComplete, symbol_table.endo_nbr());
file_open = true;
//Temporary variables declaration
......@@ -1178,11 +1201,13 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
int block_max_lag = max_leadlag_block[block].first;
int block_max_lead = max_leadlag_block[block].second;
if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE
|| simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple
|| simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete)
{
Write_Inf_To_Bin_File_Block(basename, block, u_count_int, file_open,
simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE, linear_decomposition);
simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple, linear_decomposition);
file_open = true;
}
map<tuple<int, int, int>, expr_t> tmp_block_endo_derivative;
......@@ -1311,14 +1336,14 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
switch (simulation_type)
{
evaluation:
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
case BlockSimulationType::evaluateBackward:
case BlockSimulationType::evaluateForward:
equ_type = getBlockEquationType(block, i);
{
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
fnumexpr.write(code_file, instruction_number);
}
if (equ_type == E_EVALUATE)
if (equ_type == EquationType::evaluate)
{
eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
lhs = eq_node->arg1;
......@@ -1326,7 +1351,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, true, false);
}
else if (equ_type == E_EVALUATE_S)
else if (equ_type == EquationType::evaluate_s)
{
eq_node = static_cast<BinaryOpNode *>(getBlockEquationRenormalizedExpr(block, i));
lhs = eq_node->arg1;
......@@ -1335,10 +1360,10 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, true, false);
}
break;
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
case BlockSimulationType::solveTwoBoundariesComplete:
case BlockSimulationType::solveTwoBoundariesSimple:
if (i < static_cast<int>(block_recursive))
goto evaluation;
variable_ID = getBlockVariableID(block, i);
......@@ -1371,13 +1396,13 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
fjmp_if_eval.write(code_file, instruction_number);
int prev_instruction_number = instruction_number;
// The Jacobian if we have to solve the block determinsitic block
if (simulation_type != EVALUATE_BACKWARD
&& simulation_type != EVALUATE_FORWARD)
if (simulation_type != BlockSimulationType::evaluateBackward
&& simulation_type != BlockSimulationType::evaluateForward)
{
switch (simulation_type)
{
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case BlockSimulationType::solveBackwardSimple:
case BlockSimulationType::solveForwardSimple:
{
FNUMEXPR_ fnumexpr(FirstEndoDerivative, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0);
fnumexpr.write(code_file, instruction_number);
......@@ -1389,10 +1414,10 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
}
break;
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
case BlockSimulationType::solveTwoBoundariesComplete:
case BlockSimulationType::solveTwoBoundariesSimple:
count_u = feedback_variables.size();
for (const auto &it : blocks_derivatives[block])
{
......@@ -1403,7 +1428,9 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
unsigned int varr = getBlockVariableID(block, var);
if (eq >= block_recursive and var >= block_recursive)
{
if (lag != 0 && (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE))
if (lag != 0
&& (simulation_type == BlockSimulationType::solveForwardComplete
|| simulation_type == BlockSimulationType::solveBackwardComplete))
continue;
if (!Uf[eqr].Ufl)
{
......@@ -1927,7 +1954,8 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
block_recursive = block_size - block_mfs;
BlockSimulationType simulation_type = getBlockSimulationType(block);
if (simulation_type == EVALUATE_FORWARD || simulation_type == EVALUATE_BACKWARD)
if (simulation_type == BlockSimulationType::evaluateForward
|| simulation_type == BlockSimulationType::evaluateBackward)
{
for (unsigned int ik = 0; ik < block_size; ik++)
{
......@@ -1948,23 +1976,23 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
switch (simulation_type)
{
case EVALUATE_FORWARD:
case EVALUATE_BACKWARD:
case BlockSimulationType::evaluateForward:
case BlockSimulationType::evaluateBackward:
mDynamicModelFile << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 1, it_-1, 1);" << endl
<< " residual(y_index_eq)=ys(y_index)-y(it_, y_index);" << endl;
break;
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_SIMPLE:
case BlockSimulationType::solveForwardSimple:
case BlockSimulationType::solveBackwardSimple:
mDynamicModelFile << " [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_, 1);" << endl
<< " residual(y_index_eq)=r;" << endl;
break;
case SOLVE_FORWARD_COMPLETE:
case SOLVE_BACKWARD_COMPLETE:
case BlockSimulationType::solveForwardComplete:
case BlockSimulationType::solveBackwardComplete:
mDynamicModelFile << " [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_, 1);" << endl
<< " residual(y_index_eq)=r;" << endl;
break;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case BlockSimulationType::solveTwoBoundariesComplete:
case BlockSimulationType::solveTwoBoundariesSimple:
mDynamicModelFile << " [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_-" << max_lag << ", 1, " << max_lag << ", " << block_recursive << "," << "options_.periods" << ");" << endl
<< " residual(y_index_eq)=r(:,M_.maximum_lag+1);" << endl;
break;
......@@ -2017,7 +2045,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
unsigned int block_recursive = block_size - block_mfs;
BlockSimulationType simulation_type = getBlockSimulationType(block);
if (simulation_type == EVALUATE_FORWARD && block_size)
if (simulation_type == BlockSimulationType::evaluateForward && block_size)
{
if (open_par)
mDynamicModelFile << " end" << endl;
......@@ -2043,7 +2071,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
<< " return;" << endl
<< " end;" << endl;
}
else if (simulation_type == EVALUATE_BACKWARD && block_size)
else if (simulation_type == BlockSimulationType::evaluateBackward && block_size)
{
if (open_par)
mDynamicModelFile << " end" << endl;
......@@ -2069,7 +2097,8 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
<< " return;" << endl
<< " end;" << endl;
}
else if ((simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_FORWARD_SIMPLE) && block_size)
else if ((simulation_type == BlockSimulationType::solveForwardComplete
|| simulation_type == BlockSimulationType::solveForwardSimple) && block_size)
{
if (open_par)
mDynamicModelFile << " end" << endl;
......@@ -2099,7 +2128,8 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
<< " return;" << endl
<< " end;" << endl;
}
else if ((simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_SIMPLE) && block_size)
else if ((simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveBackwardSimple) && block_size)
{
if (open_par)
mDynamicModelFile << " end" << endl;
......@@ -2130,7 +2160,8 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
<< " return;" << endl
<< " end;" << endl;
}
else if ((simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE) && block_size)
else if ((simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple) && block_size)
{
if (open_par)
mDynamicModelFile << " end" << endl;
......@@ -3173,7 +3204,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
for (const auto &it : other_endo_block[block])
for (int it1 : it.second)
other_endogenous.insert(it1);
output << "block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";" << endl
output << "block_structure.block(" << block+1 << ").Simulation_Type = " << static_cast<int>(simulation_type) << ";" << endl
<< "block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";" << endl
<< "block_structure.block(" << block+1 << ").maximum_lead = " << max_lead << ";" << endl
<< "block_structure.block(" << block+1 << ").maximum_endo_lag = " << max_lag_endo << ";" << endl
......@@ -4832,15 +4863,15 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
tie(blocks, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed)
= select_non_linear_equations_and_variables(is_equation_linear);
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, 0);
equationTypeDetermination(first_order_endo_derivatives, 0);
prologue = 0;
epilogue = 0;
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
computeChainRuleJacobian();
blocks_linear = BlockLinear();
determineLinearBlocks();
collect_block_first_order_derivatives();
......@@ -4861,7 +4892,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, mfs);
equationTypeDetermination(first_order_endo_derivatives, mfs);
cout << "Finding the optimal block decomposition of the model ..." << endl;
......@@ -4869,13 +4900,13 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
tie(blocks, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false, true);
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
printBlockDecomposition(blocks);
computeChainRuleJacobian();
blocks_linear = BlockLinear();
determineLinearBlocks();
collect_block_first_order_derivatives();
......@@ -5015,7 +5046,8 @@ DynamicModel::get_Derivatives(int block)
int max_lag, max_lead;
map<tuple<int, int, int, int, int>, int> Derivatives;
BlockSimulationType simulation_type = getBlockSimulationType(block);
if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
if (simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
{
max_lag = 1;
max_lead = 1;
......@@ -5045,7 +5077,8 @@ DynamicModel::get_Derivatives(int block)
if (OK)
{
if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursive)
if (getBlockEquationType(block, eq) == EquationType::evaluate_s
&& eq < block_nb_recursive)
//It's a normalized equation, we have to recompute the derivative using chain rule derivative function
Derivatives[{ lag, eq, var, eqr, varr }] = 1;
else
......@@ -5085,7 +5118,7 @@ DynamicModel::computeChainRuleJacobian()
int block_nb_recursives = block_size - block_nb_mfs;
for (int i = 0; i < block_nb_recursives; i++)
{
if (getBlockEquationType(block, i) == E_EVALUATE_S)
if (getBlockEquationType(block, i) == EquationType::evaluate_s)
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
else
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
......@@ -5101,7 +5134,8 @@ DynamicModel::computeChainRuleJacobian()
first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else if (Deriv_type == 2)
{
if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursives)
if (getBlockEquationType(block, eq) == EquationType::evaluate_s
&& eq < block_nb_recursives)
first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else
first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
......
......@@ -619,7 +619,7 @@ public:
bool
isBlockEquationRenormalized(int block_number, int equation_number) const override
{
return equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first == E_EVALUATE_S;
return equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first == EquationType::evaluate_s;
};
//! Return the expr_t of the equation equation_number belonging to the block block_number
expr_t
......
......@@ -136,19 +136,19 @@ ExprNode::collectVariables(SymbolType type, set<int> &result) const
void
ExprNode::collectEndogenous(set<pair<int, int>> &result) const
{
set<pair<int, int>> symb_ids;
collectDynamicVariables(SymbolType::endogenous, symb_ids);
for (const auto &symb_id : symb_ids)
result.emplace(datatree.symbol_table.getTypeSpecificID(symb_id.first), symb_id.second);
set<pair<int, int>> symb_ids_and_lags;
collectDynamicVariables(SymbolType::endogenous, symb_ids_and_lags);
for (const auto &[symb_id, lag] : symb_ids_and_lags)
result.emplace(datatree.symbol_table.getTypeSpecificID(symb_id), lag);
}
void
ExprNode::collectExogenous(set<pair<int, int>> &result) const
{
set<pair<int, int>> symb_ids;
collectDynamicVariables(SymbolType::exogenous, symb_ids);
for (const auto &symb_id : symb_ids)
result.emplace(datatree.symbol_table.getTypeSpecificID(symb_id.first), symb_id.second);
set<pair<int, int>> symb_ids_and_lags;
collectDynamicVariables(SymbolType::exogenous, symb_ids_and_lags);
for (const auto &[symb_id, lag] : symb_ids_and_lags)
result.emplace(datatree.symbol_table.getTypeSpecificID(symb_id), lag);
}
void
......
......@@ -264,41 +264,13 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
BipartiteGraph g(2 * n);
// Fill in the graph
set<pair<int, int>> endo;
for (const auto &it : contemporaneous_jacobian)
add_edge(it.first.first + n, it.first.second, g);
for (const auto &[eq_and_endo, val] : contemporaneous_jacobian)
add_edge(eq_and_endo.first + n, eq_and_endo.second, g);
// Compute maximum cardinality matching
vector<int> mate_map(2*n);
#if 1
bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
#else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
fill(mate_map.begin(), mate_map.end(), boost::graph_traits<BipartiteGraph>::null_vertex());
auto natural_endo2eqs = computeNormalizedEquations();
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
if (natural_endo2eqs.count(i) == 0)
continue;
int j = natural_endo2eqs.find(i)->second;
put(&mate_map[0], i, n+j);
put(&mate_map[0], n+j, i);
}
boost::edmonds_augmenting_path_finder<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type> augmentor(g, &mate_map[0], get(boost::vertex_index, g));
while (augmentor.augment_matching())
{
};
augmentor.get_current_matching(&mate_map[0]);
bool check = boost::maximum_cardinality_matching_verifier<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(boost::vertex_index, g));
#endif
assert(check);
......@@ -312,29 +284,6 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
endo2eq.resize(equations.size());
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), [=](int i) { return i-n; });
#ifdef DEBUG
auto natural_endo2eqs = computeNormalizedEquations(natural_endo2eqs);
int n1 = 0, n2 = 0;
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
if (natural_endo2eqs.count(i) == 0)
continue;
n1++;
auto x = natural_endo2eqs.equal_range(i);
if (find_if(x.first, x.second, [=](auto y) { return y.second == endo2eq[i]; }) == x.second)
cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, i))
<< " not used." << endl;
else
n2++;
}
cout << "Used " << n2 << " natural normalizations out of " << n1 << ", for a total of " << n << " equations." << endl;
#endif
// Check if all variables are normalized
if (auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
it != mate_map.begin() + n)
......@@ -357,28 +306,27 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
int n = equations.size();
// compute the maximum value of each row of the contemporaneous Jacobian matrix
//jacob_map normalized_contemporaneous_jacobian;
// Compute the maximum value of each row of the contemporaneous Jacobian matrix
jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
vector<double> max_val(n, 0.0);
for (const auto &it : contemporaneous_jacobian)
if (fabs(it.second) > max_val[it.first.first])
max_val[it.first.first] = fabs(it.second);
for (const auto &[eq_and_endo, val] : contemporaneous_jacobian)
max_val[eq_and_endo.first] = max(max_val[eq_and_endo.first], fabs(val));
for (auto &iter : normalized_contemporaneous_jacobian)
iter.second /= max_val[iter.first.first];
for (auto &[eq_and_endo, val] : normalized_contemporaneous_jacobian)
val /= max_val[eq_and_endo.first];
//We start with the highest value of the cutoff and try to normalize the model
// We start with the highest value of the cutoff and try to normalize the model
double current_cutoff = 0.99999999;
const double cutoff_lower_limit = 1e-19;
int suppressed = 0;
while (!check && current_cutoff > 1e-19)
while (!check && current_cutoff > cutoff_lower_limit)
{
jacob_map_t tmp_normalized_contemporaneous_jacobian;
int suppress = 0;
for (auto &iter : normalized_contemporaneous_jacobian)
if (fabs(iter.second) > max(current_cutoff, cutoff))
tmp_normalized_contemporaneous_jacobian[{ iter.first.first, iter.first.second }] = iter.second;
for (auto &[eq_and_endo, val] : normalized_contemporaneous_jacobian)
if (fabs(val) > max(current_cutoff, cutoff))
tmp_normalized_contemporaneous_jacobian[eq_and_endo] = val;
else
suppress++;
......@@ -389,7 +337,7 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
{
current_cutoff /= 2;
// In this last case try to normalize with the complete jacobian
if (current_cutoff <= 1e-19)
if (current_cutoff <= cutoff_lower_limit)
check = computeNormalization(normalized_contemporaneous_jacobian, false);
}
}
......@@ -399,34 +347,33 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
//if no non-singular normalization can be found, try to find a normalization even with a potential singularity
jacob_map_t tmp_normalized_contemporaneous_jacobian;
set<pair<int, int>> endo;
for (int i = 0; i < n; i++)
{
endo.clear();
equations[i]->collectEndogenous(endo);
for (const auto &it : endo)
tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1;
set<pair<int, int>> endos_and_lags;
equations[i]->collectEndogenous(endos_and_lags);
for (const auto &[endo, lag] : endos_and_lags)
tmp_normalized_contemporaneous_jacobian[{ i, endo }] = 1;
}
check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
if (check)
{
// Update the jacobian matrix
for (const auto &[key, ignore] : tmp_normalized_contemporaneous_jacobian)
for (const auto &[eq_and_endo, ignore] : tmp_normalized_contemporaneous_jacobian)
{
if (static_jacobian.find({ key.first, key.second }) == static_jacobian.end())
static_jacobian[{ key.first, key.second }] = 0;
if (dynamic_jacobian.find({ 0, key.first, key.second }) == dynamic_jacobian.end())
dynamic_jacobian[{ 0, key.first, key.second }] = nullptr;
if (contemporaneous_jacobian.find({ key.first, key.second }) == contemporaneous_jacobian.end())
contemporaneous_jacobian[{ key.first, key.second }] = 0;
if (static_jacobian.find(eq_and_endo) == static_jacobian.end())
static_jacobian[eq_and_endo] = 0;
if (dynamic_jacobian.find({ 0, eq_and_endo.first, eq_and_endo.second }) == dynamic_jacobian.end())
dynamic_jacobian[{ 0, eq_and_endo.first, eq_and_endo.second }] = nullptr;
if (contemporaneous_jacobian.find(eq_and_endo) == contemporaneous_jacobian.end())
contemporaneous_jacobian[eq_and_endo] = 0;
try
{
if (derivatives[1].find({ key.first, getDerivID(symbol_table.getID(SymbolType::endogenous, key.second), 0) }) == derivatives[1].end())
derivatives[1][{ key.first, getDerivID(symbol_table.getID(SymbolType::endogenous, key.second), 0) }] = Zero;
if (derivatives[1].find({ eq_and_endo.first, getDerivID(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second), 0) }) == derivatives[1].end())
derivatives[1][{ eq_and_endo.first, getDerivID(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second), 0) }] = Zero;
}
catch (DataTree::UnknownDerivIDException &e)
{
cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, key.second))
cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, eq_and_endo.second))
<< " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block." << endl;
exit(EXIT_FAILURE);
}
......@@ -441,36 +388,11 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
}
}
multimap<int, int>
ModelTree::computeNormalizedEquations() const
{
multimap<int, int> endo2eqs;
for (size_t i = 0; i < equations.size(); i++)
{
auto lhs = dynamic_cast<VariableNode *>(equations[i]->arg1);
if (!lhs)
continue;
int symb_id = lhs->symb_id;
if (symbol_table.getType(symb_id) != SymbolType::endogenous)
continue;
set<pair<int, int>> endo;
equations[i]->arg2->collectEndogenous(endo);
if (endo.find({ symbol_table.getTypeSpecificID(symb_id), 0 }) != endo.end())
continue;
endo2eqs.emplace(symbol_table.getTypeSpecificID(symb_id), static_cast<int>(i));
cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << i+1 << endl;
}
return endo2eqs;
}
pair<ModelTree::jacob_map_t, ModelTree::jacob_map_t>
ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose)
{
jacob_map_t contemporaneous_jacobian, static_jacobian;
int nb_elements_contemparenous_Jacobian = 0;
int nb_elements_contemporaneous_jacobian = 0;
set<vector<int>> jacobian_elements_to_delete;
for (const auto &[indices, d1] : derivatives[1])
{
......@@ -500,20 +422,22 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double
if (fabs(val) < cutoff)
{
if (verbose)
cout << "the coefficient related to variable " << var << " with lag " << lag << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")" << endl;
cout << "The coefficient related to variable " << var << " with lag " << lag << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")." << endl;
jacobian_elements_to_delete.insert({ eq, deriv_id });
}
else
{
if (lag == 0)
{
nb_elements_contemparenous_Jacobian++;
nb_elements_contemporaneous_jacobian++;
contemporaneous_jacobian[{ eq, var }] = val;
}
if (static_jacobian.find({ eq, var }) != static_jacobian.end())
static_jacobian[{ eq, var }] += val;
else
static_jacobian[{ eq, var }] = val;
dynamic_jacobian[{ lag, eq, var }] = d1;
}
}
......@@ -525,8 +449,8 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double
if (jacobian_elements_to_delete.size() > 0)
{
cout << jacobian_elements_to_delete.size() << " elements among " << derivatives[1].size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded" << endl
<< "The contemporaneous incidence matrix has " << nb_elements_contemparenous_Jacobian << " elements" << endl;
cout << jacobian_elements_to_delete.size() << " elements among " << derivatives[1].size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded." << endl
<< "The contemporaneous incidence matrix has " << nb_elements_contemporaneous_jacobian << " elements." << endl;
}
return { contemporaneous_jacobian, static_jacobian };
......@@ -728,20 +652,21 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg)
}
}
equation_type_and_normalized_equation_t
ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs) const
void
ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs)
{
expr_t lhs;
BinaryOpNode *eq_node;
EquationType Equation_Simulation_Type;
equation_type_and_normalized_equation_t V_Equation_Simulation_Type(equations.size());
equation_type_and_normalized_equation.clear();
equation_type_and_normalized_equation.resize(equations.size());
for (unsigned int i = 0; i < equations.size(); i++)
{
int eq = equation_reordered[i];
int var = variable_reordered[i];
eq_node = equations[eq];
lhs = eq_node->arg1;
Equation_Simulation_Type = E_SOLVE;
Equation_Simulation_Type = EquationType::solve;
pair<bool, expr_t> res;
if (auto derivative = first_order_endo_derivatives.find({ eq, var, 0 });
derivative != first_order_endo_derivatives.end())
......@@ -751,7 +676,7 @@ ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &fi
auto d_endo_variable = result.find({ var, 0 });
//Determine whether the equation could be evaluated rather than to be solved
if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, variable_reordered[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
Equation_Simulation_Type = E_EVALUATE;
Equation_Simulation_Type = EquationType::evaluate;
else
{
vector<tuple<int, expr_t, expr_t>> List_of_Op_RHS;
......@@ -759,18 +684,17 @@ ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &fi
if (mfs == 2)
{
if (d_endo_variable == result.end() && res.second)
Equation_Simulation_Type = E_EVALUATE_S;
Equation_Simulation_Type = EquationType::evaluate_s;
}
else if (mfs == 3)
{
if (res.second) // The equation could be solved analytically
Equation_Simulation_Type = E_EVALUATE_S;
Equation_Simulation_Type = EquationType::evaluate_s;
}
}
}
V_Equation_Simulation_Type[eq] = { Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second) };
equation_type_and_normalized_equation[eq] = { Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second) };
}
return V_Equation_Simulation_Type;
}
pair<lag_lead_vector_t, lag_lead_vector_t>
......@@ -912,7 +836,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
if (select_feedback_variable)
{
for (int i = 0; i < n; i++)
if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
if (Equation_Type[equation_reordered[i+prologue]].first == EquationType::solve
|| variable_lag_lead[variable_reordered[i+prologue]].second > 0
|| variable_lag_lead[variable_reordered[i+prologue]].first > 0
|| equation_lag_lead[equation_reordered[i+prologue]].second > 0
......@@ -922,7 +846,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
}
else
for (int i = 0; i < n; i++)
if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
if (Equation_Type[equation_reordered[i+prologue]].first == EquationType::solve || mfs == 0)
add_edge(vertex(i, G2), vertex(i, G2), G2);
//Determines the dynamic structure of each equation
......@@ -1053,7 +977,9 @@ ModelTree::printBlockDecomposition(const vector<pair<int, int>> &blocks) const
for (unsigned int block = 0; block < Nb_TotalBlocks; block++)
{
BlockSimulationType simulation_type = getBlockSimulationType(block);
if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
if (simulation_type == BlockSimulationType::solveForwardComplete
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesComplete)
{
Nb_SimulBlocks++;
int size = getBlockSize(block);
......@@ -1072,15 +998,15 @@ ModelTree::printBlockDecomposition(const vector<pair<int, int>> &blocks) const
<< " and " << Nb_feedback_variable << " feedback variable(s)." << endl;
}
block_type_firstequation_size_mfs_t
void
ModelTree::reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<unsigned int> &n_static, const vector<unsigned int> &n_forward, const vector<unsigned int> &n_backward, const vector<unsigned int> &n_mixed, bool linear_decomposition)
{
int i = 0;
int count_equ = 0, blck_count_simult = 0;
int Blck_Size, MFS_Size;
int Lead, Lag;
block_type_firstequation_size_mfs_t block_type_size_mfs;
BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
block_type_firstequation_size_mfs.clear();
BlockSimulationType Simulation_Type, prev_Type = BlockSimulationType::unknown;
int eq = 0;
unsigned int l_n_static = 0, l_n_forward = 0, l_n_backward = 0, l_n_mixed = 0;
for (i = 0; i < static_cast<int>(prologue+blocks.size()+epilogue); i++)
......@@ -1140,23 +1066,23 @@ ModelTree::reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &blocks
if (Lag > 0 && Lead > 0)
{
if (Blck_Size == 1)
Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
Simulation_Type = BlockSimulationType::solveTwoBoundariesSimple;
else
Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE;
Simulation_Type = BlockSimulationType::solveTwoBoundariesComplete;
}
else if (Blck_Size > 1)
{
if (Lead > 0)
Simulation_Type = SOLVE_BACKWARD_COMPLETE;
Simulation_Type = BlockSimulationType::solveBackwardComplete;
else
Simulation_Type = SOLVE_FORWARD_COMPLETE;
Simulation_Type = BlockSimulationType::solveForwardComplete;
}
else
{
if (Lead > 0)
Simulation_Type = SOLVE_BACKWARD_SIMPLE;
Simulation_Type = BlockSimulationType::solveBackwardSimple;
else
Simulation_Type = SOLVE_FORWARD_SIMPLE;
Simulation_Type = BlockSimulationType::solveForwardSimple;
}
l_n_static = n_static[i];
l_n_forward = n_forward[i];
......@@ -1164,20 +1090,22 @@ ModelTree::reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &blocks
l_n_mixed = n_mixed[i];
if (Blck_Size == 1)
{
if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
if (Equation_Type[equation_reordered[eq]].first == EquationType::evaluate
|| Equation_Type[equation_reordered[eq]].first == EquationType::evaluate_s)
{
if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
Simulation_Type = EVALUATE_BACKWARD;
else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
Simulation_Type = EVALUATE_FORWARD;
if (Simulation_Type == BlockSimulationType::solveBackwardSimple)
Simulation_Type = BlockSimulationType::evaluateBackward;
else if (Simulation_Type == BlockSimulationType::solveForwardSimple)
Simulation_Type = BlockSimulationType::evaluateForward;
}
if (i > 0)
{
bool is_lead = false, is_lag = false;
int c_Size = get<2>(block_type_size_mfs[block_type_size_mfs.size()-1]);
int first_equation = get<1>(block_type_size_mfs[block_type_size_mfs.size()-1]);
if (c_Size > 0 && ((prev_Type == EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
|| (prev_Type == EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag)))
int c_Size = get<2>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
int first_equation = get<1>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
if (c_Size > 0
&& ((prev_Type == BlockSimulationType::evaluateForward && Simulation_Type == BlockSimulationType::evaluateForward && !is_lead)
|| (prev_Type == BlockSimulationType::evaluateBackward && Simulation_Type == BlockSimulationType::evaluateBackward && !is_lag)))
{
for (int j = first_equation; j < first_equation+c_Size; j++)
{
......@@ -1189,45 +1117,44 @@ ModelTree::reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &blocks
is_lead = true;
}
}
if ((prev_Type == EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
|| (prev_Type == EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag))
if ((prev_Type == BlockSimulationType::evaluateForward && Simulation_Type == BlockSimulationType::evaluateForward && !is_lead)
|| (prev_Type == BlockSimulationType::evaluateBackward && Simulation_Type == BlockSimulationType::evaluateBackward && !is_lag))
{
//merge the current block with the previous one
BlockSimulationType c_Type = get<0>(block_type_size_mfs[block_type_size_mfs.size()-1]);
BlockSimulationType c_Type = get<0>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
c_Size++;
block_type_size_mfs[block_type_size_mfs.size()-1] = { c_Type, first_equation, c_Size, c_Size };
if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
block_lag_lead[block_type_size_mfs.size()-1] = { Lag, Lead };
block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1] = { c_Type, first_equation, c_Size, c_Size };
if (block_lag_lead[block_type_firstequation_size_mfs.size()-1].first > Lag)
Lag = block_lag_lead[block_type_firstequation_size_mfs.size()-1].first;
if (block_lag_lead[block_type_firstequation_size_mfs.size()-1].second > Lead)
Lead = block_lag_lead[block_type_firstequation_size_mfs.size()-1].second;
block_lag_lead[block_type_firstequation_size_mfs.size()-1] = { Lag, Lead };
auto tmp = block_col_type[block_col_type.size()-1];
block_col_type[block_col_type.size()-1] = { get<0>(tmp)+l_n_static, get<1>(tmp)+l_n_forward, get<2>(tmp)+l_n_backward, get<3>(tmp)+l_n_mixed };
}
else
{
block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
}
}
else
{
block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
}
}
else
{
block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
}
prev_Type = Simulation_Type;
eq += Blck_Size;
}
return block_type_size_mfs;
}
vector<bool>
......@@ -1248,18 +1175,20 @@ ModelTree::equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_e
return is_linear;
}
vector<bool>
ModelTree::BlockLinear() const
void
ModelTree::determineLinearBlocks()
{
unsigned int nb_blocks = getNbBlocks();
vector<bool> blocks_linear(nb_blocks, true);
blocks_linear.clear();
blocks_linear.resize(nb_blocks, true);
for (unsigned int block = 0; block < nb_blocks; block++)
{
BlockSimulationType simulation_type = getBlockSimulationType(block);
int block_size = getBlockSize(block);
block_derivatives_equation_variable_laglead_nodeid_t derivatives_block = blocks_derivatives[block];
int first_variable_position = getBlockFirstEquation(block);
if (simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
if (simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete)
for (const auto &[ignore, ignore2, lag, d1] : derivatives_block)
{
if (lag == 0)
......@@ -1275,7 +1204,8 @@ ModelTree::BlockLinear() const
}
}
}
else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
else if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
for (const auto &[ignore, ignore2, lag, d1] : derivatives_block)
{
set<pair<int, int>> endogenous;
......@@ -1291,7 +1221,6 @@ ModelTree::BlockLinear() const
the_end:
;
}
return blocks_linear;
}
int
......
......@@ -282,11 +282,9 @@ protected:
void computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian);
//! Try to find a natural normalization if all equations are matched to an endogenous variable on the LHS
bool computeNaturalNormalization();
//! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS)
multimap<int, int> computeNormalizedEquations() const;
//! Evaluate the jacobian and suppress all the elements below the cutoff
//! Evaluate the jacobian (w.r.t. endogenous) and suppress all the elements below the cutoff
/*! Returns a pair (contemporaneous_jacobian, static_jacobian). Also fills
dynamic_jacobian. */
dynamic_jacobian. External functions are evaluated to 1. */
pair<jacob_map_t, jacob_map_t> evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose);
//! Select and reorder the non linear equations of the model
/*! Returns a tuple (blocks, equation_lag_lead, variable_lag_lead, n_static,
......@@ -294,14 +292,14 @@ protected:
tuple<vector<pair<int, int>>, lag_lead_vector_t, lag_lead_vector_t, vector<unsigned int>, vector<unsigned int>, vector<unsigned int>, vector<unsigned int>> select_non_linear_equations_and_variables(const vector<bool> &is_equation_linear);
//! Search the equations and variables belonging to the prologue and the epilogue of the model
void computePrologueAndEpilogue(const jacob_map_t &static_jacobian);
//! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
equation_type_and_normalized_equation_t equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs) const;
//! Determine the type of each equation of model and try to normalize the unnormalized equation
void equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs);
//! Compute the block decomposition and for a non-recusive block find the minimum feedback set
/*! Returns a tuple (blocks, equation_lag_lead, variable_lag_lead, n_static,
n_forward, n_backward, n_mixed) */
tuple<vector<pair<int, int>>, lag_lead_vector_t, lag_lead_vector_t, vector<unsigned int>, vector<unsigned int>, vector<unsigned int>, vector<unsigned int>> computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable);
//! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<unsigned int> &n_static, const vector<unsigned int> &n_forward, const vector<unsigned int> &n_backward, const vector<unsigned int> &n_mixed, bool linear_decomposition);
void reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<unsigned int> &n_static, const vector<unsigned int> &n_forward, const vector<unsigned int> &n_backward, const vector<unsigned int> &n_mixed, bool linear_decomposition);
//! Determine the maximum number of lead and lag for the endogenous variable in a bloc
/*! Returns a pair { equation_lead,lag, variable_lead_lag } */
pair<lag_lead_vector_t, lag_lead_vector_t> getVariableLeadLagByBlock(const vector<int> &components_set, int nb_blck_sim) const;
......@@ -310,7 +308,7 @@ protected:
//! Print an abstract of the block structure of the model
void printBlockDecomposition(const vector<pair<int, int>> &blocks) const;
//! Determine for each block if it is linear or not
vector<bool> BlockLinear() const;
void determineLinearBlocks();
//! Remove equations specified by exclude_eqs
vector<int> includeExcludeEquations(set<pair<string, string>> &eqs, bool exclude_eqs,
vector<BinaryOpNode *> &equations, vector<int> &equations_lineno,
......@@ -437,61 +435,64 @@ public:
}
inline static string
c_Equation_Type(int type)
c_Equation_Type(EquationType type)
{
vector<string> c_Equation_Type =
switch (type)
{
"E_UNKNOWN ",
"E_EVALUATE ",
"E_EVALUATE_S",
"E_SOLVE "
};
return c_Equation_Type[type];
};
case EquationType::evaluate:
return "EVALUATE ";
case EquationType::evaluate_s:
return "EVALUATE_S";
case EquationType::solve:
return "SOLVE ";
default:
return "UNKNOWN ";
}
}
inline static string
BlockType0(BlockType type)
{
switch (type)
{
case SIMULTANS:
case BlockType::simultans:
return "SIMULTANEOUS TIME SEPARABLE ";
case PROLOGUE:
case BlockType::prologue:
return "PROLOGUE ";
case EPILOGUE:
case BlockType::epilogue:
return "EPILOGUE ";
case SIMULTAN:
case BlockType::simultan:
return "SIMULTANEOUS TIME UNSEPARABLE";
default:
return "UNKNOWN ";
}
};
}
inline static string
BlockSim(int type)
BlockSim(BlockSimulationType type)
{
switch (type)
{
case EVALUATE_FORWARD:
case BlockSimulationType::evaluateForward:
return "EVALUATE FORWARD ";
case EVALUATE_BACKWARD:
case BlockSimulationType::evaluateBackward:
return "EVALUATE BACKWARD ";
case SOLVE_FORWARD_SIMPLE:
case BlockSimulationType::solveForwardSimple:
return "SOLVE FORWARD SIMPLE ";
case SOLVE_BACKWARD_SIMPLE:
case BlockSimulationType::solveBackwardSimple:
return "SOLVE BACKWARD SIMPLE ";
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case BlockSimulationType::solveTwoBoundariesSimple:
return "SOLVE TWO BOUNDARIES SIMPLE ";
case SOLVE_FORWARD_COMPLETE:
case BlockSimulationType::solveForwardComplete:
return "SOLVE FORWARD COMPLETE ";
case SOLVE_BACKWARD_COMPLETE:
case BlockSimulationType::solveBackwardComplete:
return "SOLVE BACKWARD COMPLETE ";
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case BlockSimulationType::solveTwoBoundariesComplete:
return "SOLVE TWO BOUNDARIES COMPLETE";
default:
return "UNKNOWN ";
}
};
}
};
#endif
......@@ -305,24 +305,30 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
<< "% Warning : this file is generated automatically by Dynare" << endl
<< "% from model file (.mod)" << endl << endl
<< "%/" << endl;
if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
if (simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
output << "function y = static_" << block+1 << "(y, x, params)" << endl;
else
output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)" << endl;
BlockType block_type;
if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
block_type = SIMULTANS;
else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
|| simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
if (simulation_type == BlockSimulationType::solveForwardComplete
|| simulation_type == BlockSimulationType::solveBackwardComplete)
block_type = BlockType::simultans;
else if ((simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
&& getBlockFirstEquation(block) < prologue)
block_type = PROLOGUE;
else if ((simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_SIMPLE
|| simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
block_type = BlockType::prologue;
else if ((simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::evaluateBackward
|| simulation_type == BlockSimulationType::evaluateForward)
&& getBlockFirstEquation(block) >= equations.size() - epilogue)
block_type = EPILOGUE;
block_type = BlockType::epilogue;
else
block_type = SIMULTANS;
block_type = BlockType::simultans;
output << " % ////////////////////////////////////////////////////////////////////////" << endl
<< " % //" << string(" Block ").substr(int (log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
<< " //" << endl
......@@ -331,7 +337,8 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
output << " global options_;" << endl;
//The Temporary terms
if (simulation_type != EVALUATE_BACKWARD && simulation_type != EVALUATE_FORWARD)
if (simulation_type != BlockSimulationType::evaluateBackward
&& simulation_type != BlockSimulationType::evaluateForward)
output << " g1 = spalloc(" << block_mfs << ", " << block_mfs << ", " << derivative_endo[block].size() << ");" << endl;
if (v_temporary_terms_inuse[block].size())
......@@ -342,7 +349,8 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
output << " global" << tmp_output.str() << ";" << endl;
}
if (simulation_type != EVALUATE_BACKWARD && simulation_type != EVALUATE_FORWARD)
if (simulation_type != BlockSimulationType::evaluateBackward
&& simulation_type != BlockSimulationType::evaluateForward)
output << " residual=zeros(" << block_mfs << ",1);" << endl;
// The equations
......@@ -381,19 +389,19 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
switch (simulation_type)
{
case EVALUATE_BACKWARD: