stepan
/
preprocessor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

preprocessor + bytecode DLL: various enhancements to block and bytecode options (changes by Ferhat) 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3244 ac1d8469-bf42-47a9-8791-bf33cf982152
issue#70
sebastien 2009-12-16 13:21:31 +00:00
parent 04aa1dbdb3
commit 160ec5d7ca
20 changed files with 3446 additions and 4887 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1166
BlockTriangular.cc
View File

File diff suppressed because it is too large Load Diff

200
BlockTriangular.hh
View File

@ -1,200 +0,0 @@
/*
* Copyright (C) 2007-2008 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _BLOCKTRIANGULAR_HH
#define _BLOCKTRIANGULAR_HH
#include <string>
#include "CodeInterpreter.hh"
#include "ExprNode.hh"
#include "SymbolTable.hh"
//#include "ModelNormalization.hh"
//#include "ModelBlocks.hh"
#include "IncidenceMatrix.hh"
#include "ModelTree.hh"
//! Sparse matrix of double to store the values of the Jacobian
typedef map<pair<int ,int >,double> jacob_map;
typedef vector<pair<BlockSimulationType, pair<int, int> > > t_type;
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
typedef vector<pair<EquationType, NodeID > > t_etype;
//! Vector describing variables: max_lag in the block, max_lead in the block
typedef vector<pair< int, int> > t_vtype;
typedef set<int> temporary_terms_inuse_type;
typedef vector<pair< pair<int, pair<int, int> >, pair<int, int> > > chain_rule_derivatives_type;
//! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
struct IM_compact
{
int size, u_init, u_finish, nb_endo, nb_other_endo, size_exo, size_other_endo;
int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Exogenous, *Exogenous_Index, *Equ_X, *Equ_X_Index;
int *u_other_endo, *Var_other_endo, *Equ_other_endo, *Var_Index_other_endo, *Equ_Index_other_endo;
};
//! One block of the model
struct Block
{
int Size, Sized, nb_exo, nb_exo_det, nb_other_endo, Nb_Recursives;
BlockType Type;
BlockSimulationType Simulation_Type;
int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo;
int Max_Lag_Endo, Max_Lead_Endo;
int Max_Lag_Other_Endo, Max_Lead_Other_Endo;
int Max_Lag_Exo, Max_Lead_Exo;
bool is_linear;
int *Equation, *Own_Derivative;
EquationType *Equation_Type;
NodeID *Equation_Normalized;
int *Variable, *Other_Endogenous, *Exogenous;
temporary_terms_type **Temporary_Terms_in_Equation;
//temporary_terms_type *Temporary_terms;
temporary_terms_inuse_type *Temporary_InUse;
IM_compact *IM_lead_lag;
chain_rule_derivatives_type *Chain_Rule_Derivatives;
int Code_Start, Code_Length;
};
//! The set of all blocks of the model
struct Model_Block
{
int Size, Periods;
Block* Block_List;
//int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block;
};
//! Creates the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition
class BlockTriangular
{
private:
//! Find equations and endogenous variables belonging to the prologue and epilogue of the model
void Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM0);
//! Allocates and fills the Model structure describing the content of each block
void Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block *ModelBlock, t_etype &Equation_Type, int recurs_Size, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! Finds a matching between equations and endogenous variables
bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, int verbose, bool *IM0, vector<int> &Index_Var_IM) const;
//! Decomposes into recurive blocks the non purely recursive equations and determines for each block the minimum feedback variables
void Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, t_etype &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs) const;
//! determines the type of each equation of the model (could be evaluated or need to be solved)
t_etype Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int mfs);
//! Tries to merge the consecutive blocks in a single block and determine the type of each block: recursive, simultaneous, ...
t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! Compute for each variable its maximum lead and lag in its block
t_vtype Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) const;
public:
SymbolTable &symbol_table;
/*Blocks blocks;
Normalization normalization;*/
IncidenceMatrix incidencematrix;
NumericalConstants &num_const;
DataTree *Normalized_Equation;
BlockTriangular(SymbolTable &symbol_table_arg, NumericalConstants &num_const_arg);
~BlockTriangular();
//! Frees the Model structure describing the content of each block
void Free_Block(Model_Block* ModelBlock) const;
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> get_Derivatives(Model_Block *ModelBlock, int Blck);
void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, int mfs, double cutoff);
void Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM_0 , jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &equation_simulation_type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, bool dynamic, int mfs, double cutoff);
vector<int> Index_Equ_IM;
vector<int> Index_Var_IM;
int prologue, epilogue;
bool bt_verbose;
Model_Block* ModelBlock;
int periods;
inline static std::string BlockType0(int type)
{
switch (type)
{
case 0:
return ("SIMULTANEOUS TIME SEPARABLE ");
break;
case 1:
return ("PROLOGUE ");
break;
case 2:
return ("EPILOGUE ");
break;
case 3:
return ("SIMULTANEOUS TIME UNSEPARABLE");
break;
default:
return ("UNKNOWN ");
break;
}
};
inline static std::string BlockSim(int type)
{
switch (type)
{
case EVALUATE_FORWARD:
//case EVALUATE_FORWARD_R:
return ("EVALUATE FORWARD ");
break;
case EVALUATE_BACKWARD:
//case EVALUATE_BACKWARD_R:
return ("EVALUATE BACKWARD ");
break;
case SOLVE_FORWARD_SIMPLE:
return ("SOLVE FORWARD SIMPLE ");
break;
case SOLVE_BACKWARD_SIMPLE:
return ("SOLVE BACKWARD SIMPLE ");
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
return ("SOLVE TWO BOUNDARIES SIMPLE ");
break;
case SOLVE_FORWARD_COMPLETE:
return ("SOLVE FORWARD COMPLETE ");
break;
case SOLVE_BACKWARD_COMPLETE:
return ("SOLVE BACKWARD COMPLETE ");
break;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
return ("SOLVE TWO BOUNDARIES COMPLETE");
break;
default:
return ("UNKNOWN ");
break;
}
};
inline static std::string c_Equation_Type(int type)
{
char c_Equation_Type[4][13]=
{
"E_UNKNOWN ",
"E_EVALUATE ",
"E_EVALUATE_S",
"E_SOLVE "
};
return(c_Equation_Type[type]);
};
};
#endif

28
CodeInterpreter.hh
View File

@ -100,10 +100,10 @@ enum Tags
enum BlockType
{
SIMULTANS = 0, //!< Simultaneous time separable block
PROLOGUE = 1, //!< Prologue block (one equation at the beginning, later merged)
EPILOGUE = 2, //!< Epilogue block (one equation at the beginning, later merged)
SIMULTAN = 3 //!< 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
@ -126,7 +126,7 @@ enum BlockSimulationType
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
SOLVE_TWO_BOUNDARIES_COMPLETE //!< Block of several equations, newton solver needed, forward and backwar
};
//! Enumeration of possible symbol types
@ -475,9 +475,8 @@ private:
uint8_t op_code;
int size;
uint8_t type;
int *variable;
int *equation;
int *own_derivatives;
vector<int> variable;
vector<int> equation;
bool is_linear;
vector<Block_contain_type> Block_Contain_;
int endo_nbr;
@ -485,11 +484,14 @@ private:
int Max_Lead;
int u_count_int;
public:
inline FBEGINBLOCK_(){ op_code = FBEGINBLOCK; size = 0; type = UNKNOWN; variable = NULL; equation = NULL; own_derivatives = NULL;
inline FBEGINBLOCK_(){ op_code = FBEGINBLOCK; size = 0; type = UNKNOWN; /*variable = NULL; equation = NULL;*/
is_linear = false; endo_nbr = 0; Max_Lag = 0; Max_Lead = 0; u_count_int = 0;};
inline FBEGINBLOCK_(const int size_arg, const BlockSimulationType type_arg, int *variable_arg, int *equation_arg, int *own_derivatives_arg,
bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int u_count_int_arg)
{ op_code = FBEGINBLOCK; size = size_arg; type = type_arg; variable = variable_arg; equation = equation_arg; own_derivatives = own_derivatives_arg;
inline FBEGINBLOCK_(unsigned int &size_arg, BlockSimulationType &type_arg, int unsigned first_element, int unsigned &block_size,
const vector<int> &variable_arg, const vector<int> &equation_arg,
bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg)
{ op_code = FBEGINBLOCK; size = size_arg; type = type_arg;
variable = vector<int>(variable_arg.begin()+first_element, variable_arg.begin()+(first_element+block_size));
equation = vector<int>(equation_arg.begin()+first_element, equation_arg.begin()+(first_element+block_size));
is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg;/*Block_Contain.clear();*/};
inline unsigned int get_size() { return size;};
inline uint8_t get_type() { return type;};
@ -508,7 +510,6 @@ public:
{
CompileCode.write(reinterpret_cast<char *>(&variable[i]), sizeof(variable[0]));
CompileCode.write(reinterpret_cast<char *>(&equation[i]), sizeof(equation[0]));
CompileCode.write(reinterpret_cast<char *>(&own_derivatives[i]), sizeof(own_derivatives[0]));
}
if (type==SOLVE_TWO_BOUNDARIES_SIMPLE || type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
type==SOLVE_BACKWARD_COMPLETE || type==SOLVE_FORWARD_COMPLETE)
@ -532,7 +533,6 @@ public:
Block_contain_type bc;
memcpy(&bc.Variable, code, sizeof(bc.Variable)); code += sizeof(bc.Variable);
memcpy(&bc.Equation, code, sizeof(bc.Equation)); code += sizeof(bc.Equation);
memcpy(&bc.Own_Derivative, code, sizeof(bc.Own_Derivative)); code += sizeof(bc.Own_Derivative);
Block_Contain_.push_back(bc);
}
if (type==SOLVE_TWO_BOUNDARIES_SIMPLE || type==SOLVE_TWO_BOUNDARIES_COMPLETE ||

4
ComputingTasks.cc
View File

@ -876,13 +876,13 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct)
void
PlannerObjectiveStatement::computingPass()
{
model_tree->computingPass(false, true, false);
model_tree->computingPass(eval_context_type(), false, true, false);
}
void
PlannerObjectiveStatement::writeOutput(ostream &output, const string &basename) const
{
model_tree->writeStaticFile(basename + "_objective", false);
model_tree->writeStaticFile(basename + "_objective", false, false);
}
BVARDensityStatement::BVARDensityStatement(int maxnlags_arg, const OptionsList &options_list_arg) :

2427
DynamicModel.cc
View File

File diff suppressed because it is too large Load Diff

114
DynamicModel.hh
View File

@ -25,8 +25,6 @@ using namespace std;
#include <fstream>
#include "StaticModel.hh"
#include "StaticDllModel.hh"
#include "BlockTriangular.hh"
//! Stores a dynamic model
class DynamicModel : public ModelTree
@ -76,10 +74,15 @@ private:
//! Temporary terms for the file containing parameters dervicatives
temporary_terms_type params_derivs_temporary_terms;
//! Temporary terms for block decomposed models
vector< vector<temporary_terms_type> > v_temporary_terms;
vector<temporary_terms_inuse_type> v_temporary_terms_inuse;
//! Store the derivatives or the chainrule derivatives:map<pair< equation, pair< variable, lead_lag >, NodeID>
typedef map< pair< int, pair< int, int> >, NodeID> first_chain_rule_derivatives_type;
first_chain_rule_derivatives_type first_chain_rule_derivatives;
//! Writes dynamic model file (Matlab version)
void writeDynamicMFile(const string &dynamic_basename) const;
//! Writes dynamic model file (C version)
@ -91,37 +94,38 @@ private:
/*! \todo add third derivatives handling in C output */
void writeDynamicModel(ostream &DynamicOutput, bool use_dll) const;
//! Writes the Block reordred structure of the model in M output
void writeModelEquationsOrdered_M(Model_Block *ModelBlock, const string &dynamic_basename) const;
void writeModelEquationsOrdered_M(const string &dynamic_basename) const;
//! Writes the code of the Block reordred structure of the model in virtual machine bytecode
void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, map_idx_type map_idx) const;
void writeModelEquationsCodeOrdered(const string file_name, const string bin_basename, map_idx_type map_idx) const;
//! Computes jacobian and prepares for equation normalization
/*! Using values from initval/endval blocks and parameter initializations:
- computes the jacobian for the model w.r. to contemporaneous variables
- removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff)
*/
void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m, bool dynamic);
void BlockLinear(Model_Block *ModelBlock);
//void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m, bool dynamic);
//! return a map on the block jacobian
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> get_Derivatives(int block);
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(t_blocks_derivatives &blocks_derivatives);
string reform(string name) const;
map_idx_type map_idx;
//! Build The incidence matrix form the modeltree
void BuildIncidenceMatrix();
void computeTemporaryTermsOrdered(Model_Block *ModelBlock);
void computeTemporaryTermsOrdered();
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, map_idx_type &map_idx) const;
//! Get the type corresponding to a derivation ID
SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the lag corresponding to a derivation ID
int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException);
virtual int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the symbol ID corresponding to a derivation ID
int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
virtual int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Compute the column indices of the dynamic Jacobian
void computeDynJacobianCols(bool jacobianExo);
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(Model_Block *ModelBlock);
//! Computes derivatives of the Jacobian w.r. to parameters
void computeParamsDerivatives();
//! Computes temporary terms for the file containing parameters derivatives
@ -145,10 +149,50 @@ private:
//! Write chain rule derivative of a recursive equation w.r. to a variable
void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
//! Collecte the derivatives w.r. to endogenous of the block, to endogenous of previouys blocks and to exogenous
void collect_block_first_order_derivatives();
//! Factorized code for substitutions of leads/lags
/*! \param[in] type determines which type of variables is concerned */
void substituteLeadLagInternal(aux_var_t type);
private:
//! Indicate if the temporary terms are computed for the overall model (true) or not (false). Default value true
bool global_temporary_terms;
//! vector of block reordered variables and equations
vector<int> equation_reordered, variable_reordered, inv_equation_reordered, inv_variable_reordered;
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
t_equation_type_and_normalized_equation equation_type_and_normalized_equation;
//! for each block contains pair< Simulation_Type, pair < Block_Size, Recursive_part_Size > >
t_block_type_firstequation_size_mfs block_type_firstequation_size_mfs;
//! for all blocks derivatives description
t_blocks_derivatives blocks_derivatives;
//! The jacobian without the elements below the cutoff
dynamic_jacob_map dynamic_jacobian;
//! Vector indicating if the block is linear in endogenous variable (true) or not (false)
vector<bool> blocks_linear;
//! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), NodeID>
typedef map<pair< int, pair<int, int> >, NodeID> t_derivative;
//! Vector of derivative for each blocks
vector<t_derivative> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
//!List for each block and for each lag-leag all the other endogenous variables and exogenous variables
typedef set<int> t_var;
typedef map<int, t_var> t_lag_var;
vector<t_lag_var> other_endo_block, exo_block, exo_det_block;
//!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
vector<pair<int, int> > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
public:
DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
//! Adds a variable node
@ -179,8 +223,6 @@ public:
//! Writes model initialization and lead/lag incidence matrix to output
void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll) const;
//! Complete set to block decompose the model
BlockTriangular block_triangular;
//! Adds informations for simulation in a binary file
void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename,
const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
@ -192,14 +234,12 @@ public:
/*! It assumes that the static model given in argument has just been allocated */
void toStatic(StaticModel &static_model) const;
void toStaticDll(StaticDllModel &static_model) const;
//! Writes LaTeX file with the equations of the dynamic model
void writeLatexFile(const string &basename) const;
virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException);
//! Returns true indicating that this is a dynamic model
virtual bool isDynamic() const { return true; };
@ -225,6 +265,40 @@ public:
//! Fills eval context with values of model local variables and auxiliary variables
void fillEvalContext(eval_context_type &eval_context) const;
//! Return the number of blocks
virtual unsigned int getNbBlocks() const {return(block_type_firstequation_size_mfs.size());};
//! Determine the simulation type of each block
virtual BlockSimulationType getBlockSimulationType(int block_number) const {return(block_type_firstequation_size_mfs[block_number].first.first);};
//! Return the first equation number of a block
virtual unsigned int getBlockFirstEquation(int block_number) const {return(block_type_firstequation_size_mfs[block_number].first.second);};
//! Return the size of the block block_number
virtual unsigned int getBlockSize(int block_number) const {return(block_type_firstequation_size_mfs[block_number].second.first);};
//! Return the number of feedback variable of the block block_number
virtual unsigned int getBlockMfs(int block_number) const {return(block_type_firstequation_size_mfs[block_number].second.second);};
//! Return the maximum lag in a block
virtual unsigned int getBlockMaxLag(int block_number) const {return(block_lag_lead[block_number].first);};
//! Return the maximum lead in a block
virtual unsigned int getBlockMaxLead(int block_number) const {return(block_lag_lead[block_number].second);};
//! Return the type of equation (equation_number) belonging to the block block_number
virtual EquationType getBlockEquationType(int block_number, int equation_number) const {return( equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first);};
//! Return true if the equation has been normalized
virtual bool isBlockEquationRenormalized(int block_number, int equation_number) const {return( equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);};
//! Return the NodeID of the equation equation_number belonging to the block block_number
virtual NodeID getBlockEquationNodeID(int block_number, int equation_number) const {return( equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);};
//! Return the NodeID of the renormalized equation equation_number belonging to the block block_number
virtual NodeID getBlockEquationRenormalizedNodeID(int block_number, int equation_number) const {return( equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);};
//! Return the original number of equation equation_number belonging to the block block_number
virtual int getBlockEquationID(int block_number, int equation_number) const {return( equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]);};
//! Return the original number of variable variable_number belonging to the block block_number
virtual int getBlockVariableID(int block_number, int variable_number) const {return( variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);};
//! Return the position of equation_number in the block number belonging to the block block_number
virtual int getBlockInitialEquationID(int block_number, int equation_number) const {return((int)inv_equation_reordered[equation_number] - (int)block_type_firstequation_size_mfs[block_number].first.second);};
//! Return the position of variable_number in the block number belonging to the block block_number
virtual int getBlockInitialVariableID(int block_number, int variable_number) const {return((int)inv_variable_reordered[variable_number] - (int)block_type_firstequation_size_mfs[block_number].first.second);};
};
#endif

97
ExprNode.cc
View File

@ -32,7 +32,6 @@ using namespace __gnu_cxx;
#include "ExprNode.hh"
#include "DataTree.hh"
#include "BlockTriangular.hh"
#include "ModFile.hh"
ExprNode::ExprNode(DataTree &datatree_arg) : datatree(datatree_arg), preparedForDerivation(false)
@ -127,9 +126,8 @@ ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const
vector<vector<temporary_terms_type> > &v_temporary_terms,
int equation) const
{
// Nothing to do for a terminal node
}
@ -242,18 +240,17 @@ NumConstNode::computeDerivative(int deriv_id)
}
void
NumConstNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const
NumConstNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const
{
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<NumConstNode *>(this));
if (it != temporary_terms.end())
ModelBlock->Block_List[Curr_Block].Temporary_InUse->insert(idx);
temporary_terms_inuse.insert(idx);
}
void
NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
//cout << "writeOutput constante\n";
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<NumConstNode *>(this));
if (it != temporary_terms.end())
if (output_type == oMatlabDynamicModelSparse)
@ -427,13 +424,13 @@ VariableNode::computeDerivative(int deriv_id)
}
void
VariableNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const
VariableNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const
{
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<VariableNode *>(this));
if (it != temporary_terms.end())
ModelBlock->Block_List[Curr_Block].Temporary_InUse->insert(idx);
temporary_terms_inuse.insert(idx);
if (type== eModelLocalVariable)
datatree.local_variables_table[symb_id]->collectTemporary_terms(temporary_terms, ModelBlock, Curr_Block);
datatree.local_variables_table[symb_id]->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
}
void
@ -486,7 +483,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case eModelLocalVariable:
case eModFileLocalVariable:
if (output_type==oMatlabDynamicModelSparse || output_type==oMatlabStaticModelSparse)
if (output_type==oMatlabDynamicModelSparse || output_type==oMatlabStaticModelSparse || output_type == oMatlabDynamicModelSparseLocalTemporaryTerms)
{
output << "(";
datatree.local_variables_table[symb_id]->writeOutput(output, output_type,temporary_terms);
@ -510,6 +507,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case oMatlabDynamicModelSparse:
case oMatlabDynamicModelSparseLocalTemporaryTerms:
i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
if (lag > 0)
output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
@ -535,6 +533,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
{
case oMatlabDynamicModel:
case oMatlabDynamicModelSparse:
case oMatlabDynamicModelSparseLocalTemporaryTerms:
if (lag > 0)
output << "x(it_+" << lag << ", " << i << ")";
else if (lag < 0)
@ -572,6 +571,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
{
case oMatlabDynamicModel:
case oMatlabDynamicModelSparse:
case oMatlabDynamicModelSparseLocalTemporaryTerms:
if (lag > 0)
output << "x(it_+" << lag << ", " << i << ")";
else if (lag < 0)
@ -695,12 +695,11 @@ VariableNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const
vector<vector<temporary_terms_type> > &v_temporary_terms,
int equation) const
{
if (type== eModelLocalVariable)
datatree.local_variables_table[symb_id]->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx);
datatree.local_variables_table[symb_id]->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
}
void
@ -1204,9 +1203,8 @@ UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const
vector< vector<temporary_terms_type> > &v_temporary_terms,
int equation) const
{
NodeID this2 = const_cast<UnaryOpNode *>(this);
map<NodeID, int>::iterator it = reference_count.find(this2);
@ -1214,7 +1212,7 @@ UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
{
reference_count[this2] = 1;
first_occurence[this2] = make_pair(Curr_block,equation);
arg->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx);
arg->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
}
else
{
@ -1222,19 +1220,19 @@ UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C)
{
temporary_terms.insert(this2);
ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2);
v_temporary_terms[first_occurence[this2].first][first_occurence[this2].second].insert(this2);
}
}
}
void
UnaryOpNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const
UnaryOpNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const
{
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode*>(this));
if (it != temporary_terms.end())
ModelBlock->Block_List[Curr_Block].Temporary_InUse->insert(idx);
temporary_terms_inuse.insert(idx);
else
arg->collectTemporary_terms(temporary_terms, ModelBlock, Curr_Block);
arg->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
}
void
@ -1326,6 +1324,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
cerr << "Steady State Operator not implemented for oCDynamicModel." << endl;
exit(EXIT_FAILURE);
case oMatlabDynamicModelSparse:
case oMatlabDynamicModelSparseLocalTemporaryTerms:
cerr << "Steady State Operator not implemented for oMatlabDynamicModelSparse." << endl;
exit(EXIT_FAILURE);
default:
@ -1989,9 +1988,8 @@ BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const
vector<vector<temporary_terms_type> > &v_temporary_terms,
int equation) const
{
NodeID this2 = const_cast<BinaryOpNode *>(this);
map<NodeID, int>::iterator it = reference_count.find(this2);
@ -1999,8 +1997,8 @@ BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
{
reference_count[this2] = 1;
first_occurence[this2] = make_pair(Curr_block, equation);
arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx);
arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx);
arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
}
else
{
@ -2008,7 +2006,7 @@ BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C)
{
temporary_terms.insert(this2);
ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2);
v_temporary_terms[first_occurence[this2].first][first_occurence[this2].second].insert(this2);
}
}
}
@ -2092,15 +2090,15 @@ BinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
}
void
BinaryOpNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const
BinaryOpNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const
{
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
if (it != temporary_terms.end())
ModelBlock->Block_List[Curr_Block].Temporary_InUse->insert(idx);
temporary_terms_inuse.insert(idx);
else
{
arg1->collectTemporary_terms(temporary_terms, ModelBlock, Curr_Block);
arg2->collectTemporary_terms(temporary_terms, ModelBlock, Curr_Block);
arg1->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
arg2->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
}
}
@ -2109,7 +2107,6 @@ void
BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
//cout << "writeOutput binary\n";
// If current node is a temporary term
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
if (it != temporary_terms.end())
@ -2897,9 +2894,8 @@ TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const
vector<vector<temporary_terms_type> > &v_temporary_terms,
int equation) const
{
NodeID this2 = const_cast<TrinaryOpNode *>(this);
map<NodeID, int>::iterator it = reference_count.find(this2);
@ -2907,9 +2903,9 @@ TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
{
reference_count[this2] = 1;
first_occurence[this2] = make_pair(Curr_block,equation);
arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx);
arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx);
arg3->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, ModelBlock, equation, map_idx);
arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
arg3->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
}
else
{
@ -2917,7 +2913,7 @@ TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
if (reference_count[this2] * cost(temporary_terms, false) > MIN_COST_C)
{
temporary_terms.insert(this2);
ModelBlock->Block_List[first_occurence[this2].first].Temporary_Terms_in_Equation[first_occurence[this2].second]->insert(this2);
v_temporary_terms[first_occurence[this2].first][first_occurence[this2].second].insert(this2);
}
}
}
@ -2972,16 +2968,16 @@ TrinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms
}
void
TrinaryOpNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const
TrinaryOpNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const
{
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<TrinaryOpNode *>(this));
if (it != temporary_terms.end())
ModelBlock->Block_List[Curr_Block].Temporary_InUse->insert(idx);
temporary_terms_inuse.insert(idx);
else
{
arg1->collectTemporary_terms(temporary_terms, ModelBlock, Curr_Block);
arg2->collectTemporary_terms(temporary_terms, ModelBlock, Curr_Block);
arg3->collectTemporary_terms(temporary_terms, ModelBlock, Curr_Block);
arg1->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
arg2->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
arg3->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
}
}
@ -3206,9 +3202,8 @@ UnknownFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const
vector< vector<temporary_terms_type> > &v_temporary_terms,
int equation) const
{
cerr << "UnknownFunctionNode::computeTemporaryTerms: not implemented" << endl;
exit(EXIT_FAILURE);
@ -3223,11 +3218,11 @@ UnknownFunctionNode::collectVariables(SymbolType type_arg, set<pair<int, int> >
}
void
UnknownFunctionNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const
UnknownFunctionNode::collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const
{
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<UnknownFunctionNode *>(this));
if (it != temporary_terms.end())
ModelBlock->Block_List[Curr_Block].Temporary_InUse->insert(idx);
temporary_terms_inuse.insert(idx);
else
{
//arg->collectTemporary_terms(temporary_terms, result);

57
ExprNode.hh
View File

@ -44,6 +44,9 @@ struct ExprNodeLess;
/*! They are ordered by index number thanks to ExprNodeLess */
typedef set<NodeID, ExprNodeLess> temporary_terms_type;
//! set of temporary terms used in a block
typedef set<int> temporary_terms_inuse_type;
typedef map<int,int> map_idx_type;
//! Type for evaluation contexts
@ -63,7 +66,8 @@ enum ExprNodeOutputType
oLatexDynamicModel, //!< LaTeX code, dynamic model declarations
oLatexDynamicSteadyStateOperator, //!< LaTeX code, dynamic model steady state declarations
oMatlabDynamicSteadyStateOperator, //!< Matlab code, dynamic model steady state declarations
oMatlabDynamicModelSparseSteadyStateOperator //!< Matlab code, dynamic block decomposed mode steady state declarations
oMatlabDynamicModelSparseSteadyStateOperator, //!< Matlab code, dynamic block decomposed model steady state declarations
oMatlabDynamicModelSparseLocalTemporaryTerms //!< Matlab code, dynamic block decomposed model local temporary_terms
};
#define IS_MATLAB(output_type) ((output_type) == oMatlabStaticModel \
@ -71,6 +75,7 @@ enum ExprNodeOutputType
|| (output_type) == oMatlabOutsideModel \
|| (output_type) == oMatlabStaticModelSparse \
|| (output_type) == oMatlabDynamicModelSparse \
|| (output_type) == oMatlabDynamicModelSparseLocalTemporaryTerms \
|| (output_type) == oMatlabDynamicSteadyStateOperator \
|| (output_type) == oMatlabDynamicModelSparseSteadyStateOperator)
@ -102,7 +107,7 @@ class ExprNode
{
friend class DataTree;
friend class DynamicModel;
friend class StaticDllModel;
friend class StaticModel;
friend class ExprNodeLess;
friend class NumConstNode;
friend class VariableNode;
@ -200,14 +205,13 @@ public:
*/
virtual void collectModelLocalVariables(set<int> &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const = 0;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const = 0;
virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const;
vector< vector<temporary_terms_type> > &v_temporary_terms,
int equation) const;
class EvalException
@ -245,7 +249,7 @@ public:
typedef map<const ExprNode *, const VariableNode *> subst_table_t;
//! Creates auxiliary endo lead variables corresponding to this expression
/*!
/*!
If maximum endogenous lead >= 3, this method will also create intermediary auxiliary var, and will add the equations of the form aux1 = aux2(+1) to the substitution table.
\pre This expression is assumed to have maximum endogenous lead >= 2
\param[in,out] subst_table The table to which new auxiliary variables and their correspondance will be added
@ -255,7 +259,7 @@ public:
VariableNode *createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
//! Creates auxiliary exo lead variables corresponding to this expression
/*!
/*!
If maximum exogenous lead >= 2, this method will also create intermediary auxiliary var, and will add the equations of the form aux1 = aux2(+1) to the substitution table.
\pre This expression is assumed to have maximum exogenous lead >= 1
\param[in,out] subst_table The table to which new auxiliary variables and their correspondance will be added
@ -332,7 +336,7 @@ public:
virtual void prepareForDerivation();
virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic, bool steady_dynamic) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
@ -367,10 +371,9 @@ public:
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
vector< vector<temporary_terms_type> > &v_temporary_terms,
int equation) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic, bool steady_dynamic) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
@ -409,11 +412,10 @@ public:
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const;
vector< vector<temporary_terms_type> > &v_temporary_terms,
int equation) const;
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic, bool steady_dynamic) const;
@ -458,11 +460,10 @@ public:
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const;
vector< vector<temporary_terms_type> > &v_temporary_terms,
int equation) const;
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic, bool steady_dynamic) const;
@ -511,11 +512,10 @@ public:
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const;
vector< vector<temporary_terms_type> > &v_temporary_terms,
int equation) const;
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic, bool steady_dynamic) const;
@ -552,11 +552,10 @@ public:
temporary_terms_type &temporary_terms,
map<NodeID, pair<int, int> > &first_occurence,
int Curr_block,
Model_Block *ModelBlock,
int equation,
map_idx_type &map_idx) const;
vector< vector<temporary_terms_type> > &v_temporary_terms,
int equation) const;
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, temporary_terms_inuse_type &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic, bool steady_dynamic) const;
virtual NodeID toStatic(DataTree &static_datatree) const;

245
IncidenceMatrix.cc
View File

@ -1,245 +0,0 @@
/*
* Copyright (C) 2007-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <cstdlib>
#include <cstring>
#include "IncidenceMatrix.hh"
IncidenceMatrix::IncidenceMatrix(const SymbolTable &symbol_table_arg) :
symbol_table(symbol_table_arg)
{
Model_Max_Lead = Model_Max_Lead_Endo = Model_Max_Lead_Exo = 0;
Model_Max_Lag = Model_Max_Lag_Endo = Model_Max_Lag_Exo = 0;
}
//------------------------------------------------------------------------------
//For a lead or a lag build the Incidence Matrix structures
bool*
IncidenceMatrix::Build_IM(int lead_lag, SymbolType type)
{
int size;
bool *IM;
if(type==eEndogenous)
{
size = symbol_table.endo_nbr() * symbol_table.endo_nbr() * sizeof(IM[0]);
List_IM[lead_lag] = IM = (bool*)malloc(size);
for(int i = 0; i< symbol_table.endo_nbr() * symbol_table.endo_nbr(); i++) IM[i] = 0;
if(lead_lag > 0)
{
if(lead_lag > Model_Max_Lead_Endo)
{
Model_Max_Lead_Endo = lead_lag;
if(lead_lag > Model_Max_Lead)
Model_Max_Lead = lead_lag;
}
}
else
{
if( -lead_lag > Model_Max_Lag_Endo)
{
Model_Max_Lag_Endo = -lead_lag;
if(-lead_lag > Model_Max_Lag)
Model_Max_Lag = -lead_lag;
}
}
}
else
{ //eExogenous
size = symbol_table.endo_nbr() * symbol_table.exo_nbr() * sizeof(IM[0]);
List_IM_X[lead_lag] = IM = (bool*)malloc(size);
for(int i = 0; i< symbol_table.endo_nbr() * symbol_table.exo_nbr(); i++) IM[i] = 0;
if(lead_lag > 0)
{
if(lead_lag > Model_Max_Lead_Exo)
{
Model_Max_Lead_Exo = lead_lag;
if(lead_lag > Model_Max_Lead)
Model_Max_Lead = lead_lag;
}
}
else
{
if( -lead_lag > Model_Max_Lag_Exo)
{
Model_Max_Lag_Exo = -lead_lag;
if(-lead_lag > Model_Max_Lag)
Model_Max_Lag = -lead_lag;
}
}
}
return (IM);
}
void
IncidenceMatrix::Free_IM() const
{
IncidenceList::const_iterator it = List_IM.begin();
for(it = List_IM.begin(); it != List_IM.end(); it++)
free(it->second);
for(it = List_IM_X.begin(); it != List_IM_X.end(); it++)
free(it->second);
}
//------------------------------------------------------------------------------
// Return the incidence matrix related to a lead or a lag
bool*
IncidenceMatrix::Get_IM(int lead_lag, SymbolType type) const
{
IncidenceList::const_iterator it;
if(type==eEndogenous)
{
it = List_IM.find(lead_lag);
if(it!=List_IM.end())
return(it->second);
else
return(NULL);
}
else //eExogenous
{
it = List_IM_X.find(lead_lag);
if(it!=List_IM_X.end())
return(it->second);
else
return(NULL);
}
}
//------------------------------------------------------------------------------
// Fill the incidence matrix related to a lead or a lag
void
IncidenceMatrix::fill_IM(int equation, int variable, int lead_lag, SymbolType type)
{
bool* Cur_IM;
Cur_IM = Get_IM(lead_lag, type);
if(equation >= symbol_table.endo_nbr())
{
cout << "Error : The model has more equations (at least " << equation + 1 << ") than declared endogenous variables (" << symbol_table.endo_nbr() << ")\n";
exit(EXIT_FAILURE);
}
if (!Cur_IM)
Cur_IM = Build_IM(lead_lag, type);
if(type==eEndogenous)
Cur_IM[equation*symbol_table.endo_nbr() + variable] = 1;
else
Cur_IM[equation*symbol_table.exo_nbr() + variable] = 1;
}
//------------------------------------------------------------------------------
// unFill the incidence matrix related to a lead or a lag
void
IncidenceMatrix::unfill_IM(int equation, int variable, int lead_lag, SymbolType type)
{
bool* Cur_IM;
Cur_IM = Get_IM(lead_lag, type);
if (!Cur_IM)
Cur_IM = Build_IM(lead_lag, type);
if(type==eEndogenous)
Cur_IM[equation*symbol_table.endo_nbr() + variable] = 0;
else
Cur_IM[equation*symbol_table.exo_nbr() + variable] = 0;
}
//------------------------------------------------------------------------------
//Print azn incidence matrix
void
IncidenceMatrix::Print_SIM(bool* IM, SymbolType type) const
{
int i, j, n;
if(type == eEndogenous)
n = symbol_table.endo_nbr();
else
n = symbol_table.exo_nbr();
for(i = 0;i < symbol_table.endo_nbr();i++)
{
cout << " ";
for(j = 0;j < n;j++)
cout << IM[i*n + j] << " ";
cout << "\n";
}
}
//------------------------------------------------------------------------------
//Print all incidence matrix
void
IncidenceMatrix::Print_IM(SymbolType type) const
{
IncidenceList::const_iterator it;
cout << "-------------------------------------------------------------------\n";
if(type == eEndogenous)
for(int k=-Model_Max_Lag_Endo; k <= Model_Max_Lead_Endo; k++)
{
it = List_IM.find(k);
if(it!=List_IM.end())
{
cout << "Incidence matrix for lead_lag = " << k << "\n";
Print_SIM(it->second, type);
}
}
else // eExogenous
for(int k=-Model_Max_Lag_Exo; k <= Model_Max_Lead_Exo; k++)
{
it = List_IM_X.find(k);
if(it!=List_IM_X.end())
{
cout << "Incidence matrix for lead_lag = " << k << "\n";
Print_SIM(it->second, type);
}
}
}
//------------------------------------------------------------------------------
// Swap rows and columns of the incidence matrix
void
IncidenceMatrix::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int n) const
{
int tmp_i, j;
bool tmp_b;
/* We exchange equation (row)...*/
if(pos1 != pos2)
{
tmp_i = Index_Equ_IM[pos1];
Index_Equ_IM[pos1] = Index_Equ_IM[pos2];
Index_Equ_IM[pos2] = tmp_i;
for(j = 0;j < n;j++)
{
tmp_b = SIM[pos1 * n + j];
SIM[pos1*n + j] = SIM[pos2 * n + j];
SIM[pos2*n + j] = tmp_b;
}
}
/* ...and variables (column)*/
if(pos1 != pos3)
{
tmp_i = Index_Var_IM[pos1];
Index_Var_IM[pos1] = Index_Var_IM[pos3];
Index_Var_IM[pos3] = tmp_i;
for(j = 0;j < n;j++)
{
tmp_b = SIM[j * n + pos1];
SIM[j*n + pos1] = SIM[j * n + pos3];
SIM[j*n + pos3] = tmp_b;
}
}
}

57
IncidenceMatrix.hh
View File

@ -1,57 +0,0 @@
/*
* Copyright (C) 2007-2008 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _INCIDENCEMATRIX_HH
#define _INCIDENCEMATRIX_HH
#include <map>
#include "ExprNode.hh"
#include "SymbolTable.hh"
//! List of incidence matrix (one matrix per lead/lag)
typedef bool* pbool;
typedef map<int,pbool> IncidenceList;
//! create and manage the incidence matrix
class IncidenceMatrix
{
public:
const SymbolTable &symbol_table;
IncidenceMatrix(const SymbolTable &symbol_table_arg);
bool* Build_IM(int lead_lag, SymbolType type);
bool* Get_IM(int lead_lag, SymbolType type) const;
void fill_IM(int equation, int variable_endo, int lead_lag, SymbolType type);
void unfill_IM(int equation, int variable_endo, int lead_lag, SymbolType type);
void Free_IM() const;
void Print_IM(SymbolType type) const;
void Print_SIM(bool* IM, SymbolType type) const;
void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int n) const;
int Model_Max_Lead, Model_Max_Lag;
int Model_Max_Lead_Endo, Model_Max_Lag_Endo, Model_Max_Lead_Exo, Model_Max_Lag_Exo;
private:
IncidenceList List_IM, List_IM_X;
};
#endif

6
Makefile.am
View File

@ -16,8 +16,6 @@ dynare_m_SOURCES = \
ModelTree.hh \
StaticModel.cc \
StaticModel.hh \
StaticDllModel.cc \
StaticDllModel.hh \
DynamicModel.cc \
DynamicModel.hh \
NumericalConstants.cc \
@ -44,10 +42,6 @@ dynare_m_SOURCES = \
ExprNode.hh \
MinimumFeedbackSet.cc \
MinimumFeedbackSet.hh \
IncidenceMatrix.cc \
IncidenceMatrix.hh \
BlockTriangular.cc \
BlockTriangular.hh \
DynareMain.cc \
DynareMain2.cc \
CodeInterpreter.hh \

31
ModFile.cc
View File

@ -25,7 +25,6 @@
ModFile::ModFile() : expressions_tree(symbol_table, num_constants),
static_model(symbol_table, num_constants),
static_dll_model(symbol_table, num_constants),
dynamic_model(symbol_table, num_constants),
linear(false), block(false), byte_code(false),
use_dll(false)
@ -187,16 +186,8 @@ ModFile::computingPass(bool no_tmp_terms)
if (dynamic_model.equation_number() > 0)
{
// Compute static model and its derivatives
if(byte_code)
{
dynamic_model.toStaticDll(static_dll_model);
static_dll_model.computingPass(global_eval_context, no_tmp_terms, block);
}
else
{
dynamic_model.toStatic(static_model);
static_model.computingPass(block, false, no_tmp_terms);
}
dynamic_model.toStatic(static_model);
static_model.computingPass(global_eval_context, no_tmp_terms, false, block);
// Set things to compute for dynamic model
if (dynamic_model_needed)
{
@ -356,22 +347,14 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all
InitValStatement *ivs = dynamic_cast<InitValStatement *>(*it);
if (ivs != NULL)
{
if (!byte_code)
static_model.writeAuxVarInitval(mOutputFile);
else
static_dll_model.writeAuxVarInitval(mOutputFile);
static_model.writeAuxVarInitval(mOutputFile);
ivs->writeOutputPostInit(mOutputFile);
}
// Special treatment for load params and steady state statement: insert initial values for the auxiliary variables
LoadParamsAndSteadyStateStatement *lpass = dynamic_cast<LoadParamsAndSteadyStateStatement *>(*it);
if (lpass)
{
if (!byte_code)
static_model.writeAuxVarInitval(mOutputFile);
else
static_dll_model.writeAuxVarInitval(mOutputFile);
}
static_model.writeAuxVarInitval(mOutputFile);
}
// Remove path for block option with M-files
@ -387,10 +370,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all
// Create static and dynamic files
if (dynamic_model.equation_number() > 0)
{
if(byte_code)
static_dll_model.writeStaticFile(basename, block);
else
static_model.writeStaticFile(basename, block);
static_model.writeStaticFile(basename, block, byte_code);
if(dynamic_model_needed)
{
dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll);

5
ModFile.hh
View File

@ -29,7 +29,6 @@ using namespace std;
#include "NumericalConstants.hh"
#include "NumericalInitialization.hh"
#include "StaticModel.hh"
#include "StaticDllModel.hh"
#include "DynamicModel.hh"
#include "Statement.hh"
@ -45,10 +44,8 @@ public:
NumericalConstants num_constants;
//! Expressions outside model block
DataTree expressions_tree;
//! Static model
StaticModel static_model;
//! Static Dll model
StaticDllModel static_dll_model;
StaticModel static_model;
//! Dynamic model
DynamicModel dynamic_model;
//! Option linear

806
ModelTree.cc
View File

@ -23,6 +23,812 @@
#include <fstream>
#include "ModelTree.hh"
#include "MinimumFeedbackSet.hh"
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
#include <boost/graph/strong_components.hpp>
#include <boost/graph/topological_sort.hpp>
using namespace boost;
using namespace MFS;
void
ModelTree::computeNormalization(const set<pair<int, int> > &endo_eqs_incidence) throw (NormalizationException)
{
const int n = equation_number();
assert(n == symbol_table.endo_nbr());
typedef adjacency_list<vecS, vecS, undirectedS> BipartiteGraph;
/*
Vertices 0 to n-1 are for endogenous (using type specific ID)
Vertices n to 2*n-1 are for equations (using equation no.)
*/
BipartiteGraph g(2 * n);
// Fill in the graph
set<pair<int, int> > endo;
for(set<pair<int, int> >::const_iterator it = endo_eqs_incidence.begin(); it != endo_eqs_incidence.end(); it++)
add_edge(it->first + n, it->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(), graph_traits<BipartiteGraph>::null_vertex());
multimap<int, int> natural_endo2eqs;
computeNormalizedEquations(natural_endo2eqs);
for(int i = 0; i < symbol_table.endo_nbr(); i++)
{
if (natural_endo2eqs.count(i) == 0)
continue;
int j = natural_endo2eqs.find(i)->second;
put(&mate_map[0], i, n+j);
put(&mate_map[0], n+j, i);
}
edmonds_augmenting_path_finder<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type> augmentor(g, &mate_map[0], get(vertex_index, g));
bool not_maximum_yet = true;
while(not_maximum_yet)
{
not_maximum_yet = augmentor.augment_matching();
}
augmentor.get_current_matching(&mate_map[0]);
bool check = maximum_cardinality_matching_verifier<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(vertex_index, g));
#endif
assert(check);
#ifdef DEBUG
for(int i = 0; i < n; i++)
cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
<< " matched with equation " << (mate_map[i]-n+1) << endl;
#endif
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
endo2eq.resize(equation_number());
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));
#ifdef DEBUG
multimap<int, int> natural_endo2eqs;
computeNormalizedEquations(natural_endo2eqs);
int n1 = 0, n2 = 0;
for(int i = 0; i < symbol_table.endo_nbr(); i++)
{
if (natural_endo2eqs.count(i) == 0)
continue;
n1++;
pair<multimap<int, int>::const_iterator, multimap<int, int>::const_iterator> x = natural_endo2eqs.equal_range(i);
if (find_if(x.first, x.second, compose1(bind2nd(equal_to<int>(), endo2eq[i]), select2nd<multimap<int, int>::value_type>())) == x.second)
cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
<< " not used." << endl;
else
n2++;
}
cout << "Used " << n2 << " natural normalizations out of " << n1 << ", for a total of " << n << " equations." << endl;
#endif
// Check if all variables are normalized
vector<int>::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
if (it != mate_map.begin() + n)
throw NormalizationException(symbol_table.getID(eEndogenous, it - mate_map.begin()));
}
void
ModelTree::computePossiblySingularNormalization(const jacob_map &contemporaneous_jacobian, bool try_symbolic)
{
cout << "Normalizing the model..." << endl;
set<pair<int, int> > endo_eqs_incidence;
for(jacob_map::const_iterator it = contemporaneous_jacobian.begin();
it != contemporaneous_jacobian.end(); it++)
endo_eqs_incidence.insert(make_pair(it->first.first, it->first.second));
try
{
computeNormalization(endo_eqs_incidence);
return;
}
catch(NormalizationException &e)
{
if (try_symbolic)
cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
else
{
cerr << "ERROR: Could not normalize the model. Variable "
<< symbol_table.getName(e.symb_id)
<< " is not in the maximum cardinality matching. Try to decrease the cutoff." << endl;
exit(EXIT_FAILURE);
}
}
// If no non-singular normalization can be found, try to find a normalization even with a potential singularity
if (try_symbolic)
{
endo_eqs_incidence.clear();
set<pair<int, int> > endo;
for(int i = 0; i < equation_number(); i++)
{
endo.clear();
equations[i]->collectEndogenous(endo);
for(set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
endo_eqs_incidence.insert(make_pair(i, it->first));
}
try
{
computeNormalization(endo_eqs_incidence);
}
catch(NormalizationException &e)
{
cerr << "ERROR: Could not normalize the model even with zero cutoff. Variable "
<< symbol_table.getName(e.symb_id)
<< " is not in the maximum cardinality matching." << endl;
exit(EXIT_FAILURE);
}
}
}
void
ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
{
for(int i = 0; i < equation_number(); i++)
{
VariableNode *lhs = dynamic_cast<VariableNode *>(equations[i]->get_arg1());
if (lhs == NULL)
continue;
int symb_id = lhs->get_symb_id();
if (symbol_table.getType(symb_id) != eEndogenous)
continue;
set<pair<int, int> > endo;
equations[i]->get_arg2()->collectEndogenous(endo);
if (endo.find(make_pair(symbol_table.getTypeSpecificID(symb_id), 0)) != endo.end())
continue;
endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i));
cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl;
}
}
void
ModelTree::evaluateAndReduceJacobian(const eval_context_type &eval_context, jacob_map &contemporaneous_jacobian, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian, double cutoff, bool verbose)
{
int nb_elements_contemparenous_Jacobian = 0;
set<pair<int, int> > jacobian_elements_to_delete;
for (first_derivatives_type::iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
int deriv_id = it->first.second;
if (getTypeByDerivID(deriv_id) == eEndogenous)
{
NodeID Id = it->second;
int eq = it->first.first;
int symb = getSymbIDByDerivID(deriv_id);
int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
double val = 0;
try
{
val = Id->eval(eval_context);
}
catch (ExprNode::EvalException &e)
{
cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
Id->writeOutput(cerr, oMatlabDynamicModelSparse, temporary_terms);
cerr << endl;
exit(EXIT_FAILURE);
}
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;
jacobian_elements_to_delete.insert(make_pair(eq, deriv_id));
}
else
{
if (lag == 0)
{
nb_elements_contemparenous_Jacobian++;
contemporaneous_jacobian[make_pair(eq,var)] = val;
}
if (static_jacobian.find(make_pair(eq, var)) != static_jacobian.end())
static_jacobian[make_pair(eq, var)] += val;
else
static_jacobian[make_pair(eq, var)] = val;
dynamic_jacobian[make_pair(lag, make_pair(eq, var))] = Id;
}
}
}
// Get rid of the elements of the Jacobian matrix below the cutoff
for(set<pair<int, int> >::const_iterator it = jacobian_elements_to_delete.begin(); it != jacobian_elements_to_delete.end(); it++)
first_derivatives.erase(*it);
if (jacobian_elements_to_delete.size()>0)
{
cout << jacobian_elements_to_delete.size() << " elements among " << first_derivatives.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;
}
}
void
ModelTree::computePrologueAndEpilogue(jacob_map &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered, unsigned int &prologue, unsigned int &epilogue)
{
vector<int> eq2endo(equation_number(),0);
equation_reordered.resize(equation_number());
variable_reordered.resize(equation_number());
bool *IM;
int n = equation_number();
IM = (bool*)calloc(n*n,sizeof(bool));
int i = 0;
for(vector<int>::const_iterator it=endo2eq.begin(); it != endo2eq.end(); it++, i++)
{
eq2endo[*it] = i;
equation_reordered[i] = i;
variable_reordered[*it] = i;
}
for (jacob_map::const_iterator it = static_jacobian_arg.begin(); it != static_jacobian_arg.end(); it ++)
IM[it->first.first * n + endo2eq[it->first.second]] = true;
bool something_has_been_done = true;
prologue = 0;
int k = 0;
// Find the prologue equations and place first the AR(1) shock equations first
while (something_has_been_done)
{
int tmp_prologue = prologue;
something_has_been_done = false;
for(int i = prologue;i < n; i++)
{
int nze = 0;
for(int j = tmp_prologue; j < n; j++)
if(IM[i * n + j])
{
nze ++;
k = j;
}
if(nze == 1)
{
for(int j = 0; j < n; j++)
{
bool tmp_bool = IM[tmp_prologue * n + j];
IM[tmp_prologue * n + j] = IM[i * n + j];
IM[i * n + j] = tmp_bool;
}
int tmp = equation_reordered[tmp_prologue];
equation_reordered[tmp_prologue] = equation_reordered[i];
equation_reordered[i] = tmp;
for(int j = 0; j < n; j++)
{
bool tmp_bool = IM[j * n + tmp_prologue];
IM[j * n + tmp_prologue] = IM[j * n + k];
IM[j * n + k] = tmp_bool;
}
tmp = variable_reordered[tmp_prologue];
variable_reordered[tmp_prologue] = variable_reordered[k];
variable_reordered[k] = tmp;
tmp_prologue++;
something_has_been_done = true;
}
}
prologue = tmp_prologue;
}
something_has_been_done = true;
epilogue = 0;
// Find the epilogue equations
while (something_has_been_done)
{
int tmp_epilogue = epilogue;
something_has_been_done = false;
for(int i = prologue;i < n - (int) epilogue; i++)
{
int nze = 0;
for(int j = prologue; j < n - tmp_epilogue; j++)
if(IM[j * n + i])
{
nze ++;
k = j;
}
if(nze == 1)
{
for(int j = 0; j < n; j++)
{
bool tmp_bool = IM[(n - 1 - tmp_epilogue) * n + j];
IM[(n - 1 - tmp_epilogue) * n + j] = IM[k * n + j];
IM[k * n + j] = tmp_bool;
}
int tmp = equation_reordered[n - 1 - tmp_epilogue];
equation_reordered[n - 1 - tmp_epilogue] = equation_reordered[k];
equation_reordered[k] = tmp;
for(int j = 0; j < n; j++)
{
bool tmp_bool = IM[j * n + n - 1 - tmp_epilogue];
IM[j * n + n - 1 - tmp_epilogue] = IM[j * n + i];
IM[j * n + i] = tmp_bool;
}
tmp = variable_reordered[n - 1 - tmp_epilogue];
variable_reordered[n - 1 - tmp_epilogue] = variable_reordered[i];
variable_reordered[i] = tmp;
tmp_epilogue++;
something_has_been_done = true;
}
}
epilogue = tmp_epilogue;
}
free(IM);
}
t_equation_type_and_normalized_equation
ModelTree::equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int mfs)
{
NodeID lhs, rhs;
ostringstream tmp_output;
BinaryOpNode *eq_node;
ostringstream tmp_s;
temporary_terms_type temporary_terms;
EquationType Equation_Simulation_Type;
t_equation_type_and_normalized_equation V_Equation_Simulation_Type(equations.size());
for (unsigned int i = 0; i < equations.size(); i++)
{
temporary_terms.clear();
int eq = Index_Equ_IM[i];
int var = Index_Var_IM[i];
eq_node = equations[eq];
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
Equation_Simulation_Type = E_SOLVE;
tmp_s.str("");
tmp_output.str("");
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
tmp_s << "y(it_, " << Index_Var_IM[i]+1 << ")";
map<pair<int, pair<int, int> >, NodeID>::iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
pair<bool, NodeID> res;
if(derivative != first_order_endo_derivatives.end())
{
set<pair<int, int> > result;
derivative->second->collectEndogenous(result);
set<pair<int, int> >::const_iterator d_endo_variable = result.find(make_pair(var, 0));
//Determine whether the equation could be evaluated rather than to be solved
ostringstream tt("");
derivative->second->writeOutput(tt, oMatlabDynamicModelSparse, temporary_terms);
if (tmp_output.str() == tmp_s.str() and tt.str()=="1")
{
Equation_Simulation_Type = E_EVALUATE;
}
else
{
vector<pair<int, pair<NodeID, NodeID> > > List_of_Op_RHS;
res = equations[eq]->normalizeEquation(var, List_of_Op_RHS);
if(mfs==2)
{
if(d_endo_variable == result.end() && res.second)
Equation_Simulation_Type = E_EVALUATE_S;
}
else if(mfs==3)
{
if(res.second) // The equation could be solved analytically
Equation_Simulation_Type = E_EVALUATE_S;
}
}
}
V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second));
}
return (V_Equation_Simulation_Type);
}
void
ModelTree::getVariableLeadLagByBlock(dynamic_jacob_map &dynamic_jacobian, vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_lag_lead_vector &equation_lead_lag, t_lag_lead_vector &variable_lead_lag, vector<int> equation_reordered, vector<int> variable_reordered) const
{
int nb_endo = symbol_table.endo_nbr();
variable_lead_lag = t_lag_lead_vector(nb_endo , make_pair(0,0));
equation_lead_lag = t_lag_lead_vector(nb_endo , make_pair(0,0));
vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
for (int i = 0; i < nb_endo; i++)
{
if (i < prologue)
{
variable_blck[variable_reordered[i]] = i;
equation_blck[equation_reordered[i]] = i;
}
else if (i < (int)components_set.size() + prologue)
{
variable_blck[variable_reordered[i]] = components_set[i-prologue] + prologue;
equation_blck[equation_reordered[i]] = components_set[i-prologue] + prologue;
}
else
{
variable_blck[variable_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
equation_blck[equation_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
}
}
for (dynamic_jacob_map::const_iterator it = dynamic_jacobian.begin(); it != dynamic_jacobian.end(); it++)
{
int lag = it->first.first;
int j_1 = it->first.second.second;
int i_1 = it->first.second.second;
if (variable_blck[i_1] == equation_blck[j_1])
{
if (lag > variable_lead_lag[i_1].second)
variable_lead_lag[i_1] = make_pair(variable_lead_lag[i_1].first, lag);
if (lag < -variable_lead_lag[i_1].first)
variable_lead_lag[i_1] = make_pair(-lag, variable_lead_lag[i_1].second);
if (lag > equation_lead_lag[j_1].second)
equation_lead_lag[j_1] = make_pair(equation_lead_lag[j_1].first, lag);
if (lag < -equation_lead_lag[j_1].first)
equation_lead_lag[j_1] = make_pair(-lag, equation_lead_lag[j_1].second);
}
}
}
void
ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian, int prologue, int epilogue, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, t_equation_type_and_normalized_equation &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const
{
int nb_var = variable_reordered.size();
int n = nb_var - prologue - epilogue;
typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
GraphvizDigraph G2(n);
vector<int> reverse_equation_reordered(nb_var), reverse_variable_reordered(nb_var);
for(int i=0; i<nb_var; i++)
{
reverse_equation_reordered[equation_reordered[i]] = i;
reverse_variable_reordered[variable_reordered[i]] = i;
}
for(jacob_map::const_iterator it = static_jacobian.begin(); it != static_jacobian.end(); it++)
if( reverse_equation_reordered[it->first.first]>=prologue && reverse_equation_reordered[it->first.first]<nb_var - epilogue
&& reverse_variable_reordered[it->first.second]>=prologue && reverse_variable_reordered[it->first.second]<nb_var - epilogue
&& it->first.first != endo2eq[it->first.second])
add_edge(reverse_equation_reordered[it->first.first]-prologue, reverse_equation_reordered[endo2eq[it->first.second]]-prologue, G2);
vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
// Compute strongly connected components
int num = strong_components(G2, &endo2block[0]);
blocks = vector<pair<int, int> >(num, make_pair(0, 0));
// Create directed acyclic graph associated to the strongly connected components
typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
DirectedGraph dag(num);
for (unsigned int i = 0;i < num_vertices(G2);i++)
{
GraphvizDigraph::out_edge_iterator it_out, out_end;
GraphvizDigraph::vertex_descriptor vi = vertex(i, G2);
for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
{
int t_b = endo2block[target(*it_out, G2)];
int s_b = endo2block[source(*it_out, G2)];
if (s_b != t_b)
add_edge(s_b, t_b, dag);
}
}
// Compute topological sort of DAG (ordered list of unordered SCC)
deque<int> ordered2unordered;
topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order
// Construct mapping from unordered SCC to ordered SCC
vector<int> unordered2ordered(num);
for(int i = 0; i < num; i++)
unordered2ordered[ordered2unordered[i]] = i;
//This vector contains for each block:
// - first set = equations belonging to the block,
// - second set = the feeback variables,
// - third vector = the reordered non-feedback variables.
vector<pair<set<int>, pair<set<int>, vector<int> > > > components_set(num);
for (unsigned int i = 0; i < endo2block.size(); i++)
{
endo2block[i] = unordered2ordered[endo2block[i]];
blocks[endo2block[i]].first++;
components_set[endo2block[i]].first.insert(i);
}
t_lag_lead_vector equation_lag_lead, variable_lag_lead;
getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, prologue, epilogue, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
vector<int> tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
int order = prologue;
//Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
if(select_feedback_variable)
{
for (int i = 0; i < n; i++)
if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
or variable_lag_lead[variable_reordered[i+prologue]].second > 0 or variable_lag_lead[variable_reordered[i+prologue]].first > 0
or equation_lag_lead[equation_reordered[i+prologue]].second > 0 or equation_lag_lead[equation_reordered[i+prologue]].first > 0
or mfs == 0)
add_edge(i, i, G2);
}
else
{
for (int i = 0; i < n; i++)
if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
add_edge(i, i, G2);
}
//For each block, the minimum set of feedback variable is computed
// and the non-feedback variables are reordered to get
// a sub-recursive block without feedback variables
for (int i = 0; i < num; i++)
{
AdjacencyList_type G = GraphvizDigraph_2_AdjacencyList(G2, components_set[i].first);
set<int> feed_back_vertices;
//Print(G);
AdjacencyList_type G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, G);
components_set[i].second.first = feed_back_vertices;
blocks[i].second = feed_back_vertices.size();
vector<int> Reordered_Vertice;
Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
//First we have the recursive equations conditional on feedback variables
for (vector<int>::iterator its = Reordered_Vertice.begin(); its != Reordered_Vertice.end(); its++)
{
equation_reordered[order] = tmp_equation_reordered[*its+prologue];
variable_reordered[order] = tmp_variable_reordered[*its+prologue];
order++;
}
components_set[i].second.second = Reordered_Vertice;
//Second we have the equations related to the feedback variables
for (set<int>::iterator its = feed_back_vertices.begin(); its != feed_back_vertices.end(); its++)
{
equation_reordered[order] = tmp_equation_reordered[v_index[vertex(*its, G)]+prologue];
variable_reordered[order] = tmp_variable_reordered[v_index[vertex(*its, G)]+prologue];
order++;
}
}
inv_equation_reordered = vector<int>(nb_var);
inv_variable_reordered = vector<int>(nb_var);
for(int i = 0; i < nb_var ; i++)
{
inv_variable_reordered[variable_reordered[i]] = i;
inv_equation_reordered[equation_reordered[i]] = i;
}
}
void ModelTree::printBlockDecomposition(vector<pair<int, int> > blocks)
{
int largest_block = 0;
int Nb_SimulBlocks = 0;
int Nb_feedback_variable = 0;
unsigned int Nb_TotalBlocks = getNbBlocks();
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)
{
Nb_SimulBlocks++;
int size = getBlockSize(block);
if (size > largest_block)
{
largest_block = size;
Nb_feedback_variable = blocks[Nb_SimulBlocks-1].second;
}
}
}
int Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks;
cout << Nb_TotalBlocks << " block(s) found:" << endl
<< " " << Nb_RecursBlocks << " recursive block(s) and " << Nb_SimulBlocks << " simultaneous block(s)." << endl
<< " the largest simultaneous block has " << largest_block << " equation(s)" << endl
<< " and " << Nb_feedback_variable << " feedback variable(s)." << endl;
}
t_block_type_firstequation_size_mfs
ModelTree::reduceBlocksAndTypeDetermination(dynamic_jacob_map &dynamic_jacobian, int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_equation_type_and_normalized_equation &Equation_Type, vector<int> &variable_reordered, vector<int> &equation_reordered)
{
int i = 0;
int count_equ = 0, blck_count_simult = 0;
int Blck_Size, MFS_Size;
int Lead, Lag;
t_block_type_firstequation_size_mfs block_type_size_mfs;
BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
int eq = 0;
for (i = 0; i < prologue+(int) blocks.size()+epilogue; i++)
{
int first_count_equ = count_equ;
if (i < prologue)
{
Blck_Size = 1;
MFS_Size = 1;
}
else if (i < prologue+(int) blocks.size())
{
Blck_Size = blocks[blck_count_simult].first;
MFS_Size = blocks[blck_count_simult].second;
blck_count_simult++;
}
else if (i < prologue+(int) blocks.size()+epilogue)
{
Blck_Size = 1;
MFS_Size = 1;
}
Lag = Lead = 0;
set<pair<int, int> > endo;
for(count_equ = first_count_equ; count_equ < Blck_Size+first_count_equ; count_equ++)
{
endo.clear();
equations[equation_reordered[count_equ]]->collectEndogenous(endo);
for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
{
int curr_variable = it->first;
int curr_lag = it->second;
vector<int>::const_iterator it = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
if(it != variable_reordered.begin()+(first_count_equ+Blck_Size))
if (dynamic_jacobian.find(make_pair(curr_lag, make_pair(equation_reordered[count_equ], curr_variable))) != dynamic_jacobian.end())
{
if (curr_lag > Lead)
Lead = curr_lag;
else if (-curr_lag > Lag)
Lag = -curr_lag;
}
}
}
if ((Lag > 0) && (Lead > 0))
{
if (Blck_Size == 1)
Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
else
Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE;
}
else if (Blck_Size > 1)
{
if (Lead > 0)
Simulation_Type = SOLVE_BACKWARD_COMPLETE;
else
Simulation_Type = SOLVE_FORWARD_COMPLETE;
}
else
{
if (Lead > 0)
Simulation_Type = SOLVE_BACKWARD_SIMPLE;
else
Simulation_Type = SOLVE_FORWARD_SIMPLE;
}
if (Blck_Size == 1)
{
if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE or Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
{
if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
Simulation_Type = EVALUATE_BACKWARD;
else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
Simulation_Type = EVALUATE_FORWARD;
}
if (i > 0)
{
if ((prev_Type == EVALUATE_FORWARD and Simulation_Type == EVALUATE_FORWARD)
or (prev_Type == EVALUATE_BACKWARD and Simulation_Type == EVALUATE_BACKWARD))
{
//merge the current block with the previous one
BlockSimulationType c_Type = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.first;
int c_Size = (block_type_size_mfs[block_type_size_mfs.size()-1]).second.first;
int first_equation = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.second;
block_type_size_mfs[block_type_size_mfs.size()-1] = make_pair(make_pair(c_Type, first_equation), make_pair(++c_Size, block_type_size_mfs[block_type_size_mfs.size()-1].second.second));
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] = make_pair(Lag, Lead);
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
}
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
}
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
}
prev_Type = Simulation_Type;
eq += Blck_Size;
}
return (block_type_size_mfs);
}
vector<bool>
ModelTree::BlockLinear(t_blocks_derivatives &blocks_derivatives, vector<int> &variable_reordered)
{
unsigned int nb_blocks = getNbBlocks();
vector<bool> blocks_linear(nb_blocks, true);
for (unsigned int block = 0;block < nb_blocks; block++)
{
BlockSimulationType simulation_type = getBlockSimulationType(block);
int block_size = getBlockSize(block);
t_block_derivatives_equation_variable_laglead_nodeid derivatives_block = blocks_derivatives[block];
int first_variable_position = getBlockFirstEquation(block);
if (simulation_type==SOLVE_BACKWARD_COMPLETE || simulation_type==SOLVE_FORWARD_COMPLETE)
{
for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
{
int lag = it->second.first;
if (lag == 0)
{
NodeID Id = it->second.second;
set<pair<int, int> > endogenous;
Id->collectEndogenous(endogenous);
if (endogenous.size() > 0)
{
for (int l=0;l<block_size;l++)
{
if (endogenous.find(make_pair(variable_reordered[first_variable_position+l], 0)) != endogenous.end())
{
blocks_linear[block] = false;
goto the_end;
}
}
}
}
}
}
else if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE)
{
for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
{
int lag = it->second.first;
NodeID Id = it->second.second;//
set<pair<int, int> > endogenous;
Id->collectEndogenous(endogenous);
if (endogenous.size() > 0)
{
for (int l=0;l<block_size;l++)
{
if (endogenous.find(make_pair(variable_reordered[first_variable_position+l], lag)) != endogenous.end())
{
blocks_linear[block] = false;
goto the_end;
}
}
}
}
}
the_end:;
}
return(blocks_linear);
}
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg) :

177
ModelTree.hh
View File

@ -30,9 +30,28 @@ using namespace std;
#include "DataTree.hh"
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
typedef vector<pair<EquationType, NodeID > > t_equation_type_and_normalized_equation;
//! Vector describing variables: max_lag in the block, max_lead in the block
typedef vector<pair< int, int> > t_lag_lead_vector;
//! for each block contains pair< pair<Simulation_Type, first_equation>, pair < Block_Size, Recursive_part_Size > >
typedef vector<pair< pair< BlockSimulationType, int> , pair<int, int> > > t_block_type_firstequation_size_mfs;
//! for a block contains derivatives pair< pair<block_equation_number, block_variable_number> , pair<lead_lag, NodeID> >
typedef vector< pair<pair<int, int> , pair< int, NodeID > > > t_block_derivatives_equation_variable_laglead_nodeid;
//! for all blocks derivatives description
typedef vector<t_block_derivatives_equation_variable_laglead_nodeid> t_blocks_derivatives;
//! Shared code for static and dynamic models
class ModelTree : public DataTree
{
friend class DynamicModel;
friend class StaticModel;
protected:
//! Stores declared and generated auxiliary equations
vector<BinaryOpNode *> equations;
@ -100,6 +119,96 @@ protected:
//! Writes LaTeX model file
void writeLatexModelFile(const string &filename, ExprNodeOutputType output_type) const;
//! Sparse matrix of double to store the values of the Jacobian
/*! First index is equation number, second index is endogenous type specific ID */
typedef map<pair<int, int>, double> jacob_map;
//! Sparse matrix of double to store the values of the Jacobian
/*! First index is lag, second index is equation number, third index is endogenous type specific ID */
typedef map<pair<int, pair<int, int> >, NodeID> dynamic_jacob_map;
//! Normalization of equations
/*! Maps endogenous type specific IDs to equation numbers */
vector<int> endo2eq;
//! number of equation in the prologue and in the epilogue
unsigned int epilogue, prologue;
//! for each block contains pair< max_lag, max_lead>
t_lag_lead_vector block_lag_lead;
//! Exception thrown when normalization fails
class NormalizationException
{
public:
//! A variable missing from the maximum cardinal matching
int symb_id;
NormalizationException(int symb_id_arg) : symb_id(symb_id_arg) { }
};
//! Compute the matching between endogenous and variable using the jacobian contemporaneous_jacobian
/*! \param endo_eqs_incidence A set indicating which endogenous appear in which equation. First element of pairs is equation number, second is type specific endo ID */
void computeNormalization(const set<pair<int, int> > &endo_eqs_incidence) throw (NormalizationException);
//! Try to compute the matching between endogenous and variable using a decreasing cutoff
/*! applied to the jacobian contemporaneous_jacobian and stop when a matching is found.*/
/*! if no matching is found with a cutoff close to zero an error message is printout */
void computePossiblySingularNormalization(const jacob_map &contemporaneous_jacobian, bool try_symbolic);
//! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS)
void computeNormalizedEquations(multimap<int, int> &endo2eqs) const;
//! Evaluate the jacobian and suppress all the elements below the cutoff
void evaluateAndReduceJacobian(const eval_context_type &eval_context, jacob_map &contemporaneous_jacobian, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian, double cutoff, bool verbose);
//! Search the equations and variables belonging to the prologue and the epilogue of the model
void computePrologueAndEpilogue(jacob_map &static_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, unsigned int &prologue, unsigned int &epilogue);
//! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
t_equation_type_and_normalized_equation equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int mfs);
//! Compute the block decomposition and for a non-recusive block find the minimum feedback set
void computeBlockDecompositionAndFeedbackVariablesForEachBlock(jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, t_equation_type_and_normalized_equation &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const;
//! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
t_block_type_firstequation_size_mfs reduceBlocksAndTypeDetermination(dynamic_jacob_map &dynamic_jacobian, int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_equation_type_and_normalized_equation &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! Determine the maximum number of lead and lag for the endogenous variable in a bloc
void getVariableLeadLagByBlock(dynamic_jacob_map &dynamic_jacobian, vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_lag_lead_vector &equation_lead_lag, t_lag_lead_vector &variable_lead_lag, vector<int> equation_reordered, vector<int> variable_reordered) const;
//! Print an abstract of the block structure of the model
void printBlockDecomposition(vector<pair<int, int> > blocks);
//! Determine for each block if it is linear or not
vector<bool> BlockLinear(t_blocks_derivatives &blocks_derivatives, vector<int> &variable_reordered);
virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException) = 0;
virtual int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException) = 0;
virtual int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException) = 0;
//! Determine the simulation type of each block
virtual BlockSimulationType getBlockSimulationType(int block_number) const = 0;
//! Return the number of blocks
virtual unsigned int getNbBlocks() const = 0;
//! Return the first equation number of a block
virtual unsigned int getBlockFirstEquation(int block_number) const = 0;
//! Return the size of the block block_number
virtual unsigned int getBlockSize(int block_number) const = 0;
//! Return the number of feedback variable of the block block_number
virtual unsigned int getBlockMfs(int block_number) const = 0;
//! Return the maximum lag in a block
virtual unsigned int getBlockMaxLag(int block_number) const = 0;
//! Return the maximum lead in a block
virtual unsigned int getBlockMaxLead(int block_number) const = 0;
//! Return the type of equation (equation_number) belonging to the block block_number
virtual EquationType getBlockEquationType(int block_number, int equation_number) const = 0;
//! Return true if the equation has been normalized
virtual bool isBlockEquationRenormalized(int block_number, int equation_number) const = 0;
//! Return the NodeID of the equation equation_number belonging to the block block_number
virtual NodeID getBlockEquationNodeID(int block_number, int equation_number) const = 0;
//! Return the NodeID of the renormalized equation equation_number belonging to the block block_number
virtual NodeID getBlockEquationRenormalizedNodeID(int block_number, int equation_number) const = 0;
//! Return the original number of equation equation_number belonging to the block block_number
virtual int getBlockEquationID(int block_number, int equation_number) const = 0;
//! Return the original number of variable variable_number belonging to the block block_number
virtual int getBlockVariableID(int block_number, int variable_number) const = 0;
//! Return the position of equation_number in the block number belonging to the block block_number
virtual int getBlockInitialEquationID(int block_number, int equation_number) const = 0;
//! Return the position of variable_number in the block number belonging to the block block_number
virtual int getBlockInitialVariableID(int block_number, int variable_number) const = 0;
public:
ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
//! Declare a node as an equation of the model
@ -110,6 +219,74 @@ public:
void addAuxEquation(NodeID eq);
//! Returns the number of equations in the model
int equation_number() const;
inline static std::string c_Equation_Type(int type)
{
char c_Equation_Type[4][13]=
{
"E_UNKNOWN ",
"E_EVALUATE ",
"E_EVALUATE_S",
"E_SOLVE "
};
return(c_Equation_Type[type]);
};
inline static std::string BlockType0(BlockType type)
{
switch (type)
{
case SIMULTANS:
return ("SIMULTANEOUS TIME SEPARABLE ");
break;
case PROLOGUE:
return ("PROLOGUE ");
break;
case EPILOGUE:
return ("EPILOGUE ");
break;
case SIMULTAN:
return ("SIMULTANEOUS TIME UNSEPARABLE");
break;
default:
return ("UNKNOWN ");
break;
}
};
inline static std::string BlockSim(int type)
{
switch (type)
{
case EVALUATE_FORWARD:
return ("EVALUATE FORWARD ");
break;
case EVALUATE_BACKWARD:
return ("EVALUATE BACKWARD ");
break;
case SOLVE_FORWARD_SIMPLE:
return ("SOLVE FORWARD SIMPLE ");
break;
case SOLVE_BACKWARD_SIMPLE:
return ("SOLVE BACKWARD SIMPLE ");
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
return ("SOLVE TWO BOUNDARIES SIMPLE ");
break;
case SOLVE_FORWARD_COMPLETE:
return ("SOLVE FORWARD COMPLETE ");
break;
case SOLVE_BACKWARD_COMPLETE:
return ("SOLVE BACKWARD COMPLETE ");
break;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
return ("SOLVE TWO BOUNDARIES COMPLETE");
break;
default:
return ("UNKNOWN ");
break;
}
};
};
#endif

7
ParsingDriver.cc
View File

@ -377,7 +377,7 @@ ParsingDriver::cutoff(string *value)
{
double val = atof(value->c_str());
mod_file->dynamic_model.cutoff = val;
mod_file->static_dll_model.cutoff = val;
mod_file->static_model.cutoff = val;
delete value;
}
@ -386,7 +386,7 @@ ParsingDriver::mfs(string *value)
{
int val = atoi(value->c_str());
mod_file->dynamic_model.mfs = val;
mod_file->static_dll_model.mfs = val;
mod_file->static_model.mfs = val;
delete value;
}
@ -731,9 +731,6 @@ ParsingDriver::option_num(const string &name_option, const string &opt)
if (options_list.num_options.find(name_option) != options_list.num_options.end())
error("option " + name_option + " declared twice");
if ((name_option == "periods") && mod_file->block)
mod_file->dynamic_model.block_triangular.periods = atoi(opt.c_str());
options_list.num_options[name_option] = opt;
}

1050
StaticDllModel.cc
View File

File diff suppressed because it is too large Load Diff

137
StaticDllModel.hh
View File

@ -1,137 +0,0 @@
/*
* Copyright (C) 2003-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _STATICDLLMODEL_HH
#define _STATICDLLMODEL_HH
using namespace std;
#include <fstream>
#include "StaticModel.hh"
#include "BlockTriangular.hh"
//! Stores a static model
class StaticDllModel : public ModelTree
{
private:
typedef map<pair<int, int>, int> deriv_id_table_t;
//! Maps a pair (symbol_id, lag) to a deriv ID
deriv_id_table_t deriv_id_table;
//! Maps a deriv ID to a pair (symbol_id, lag)
vector<pair<int, int> > inv_deriv_id_table;
//! Temporary terms for the file containing parameters dervicatives
temporary_terms_type params_derivs_temporary_terms;
typedef map< pair< int, pair< int, int> >, NodeID> first_chain_rule_derivatives_type;
first_chain_rule_derivatives_type first_chain_rule_derivatives;
//! Writes static model file (Matlab version)
void writeStaticMFile(const string &static_basename) const;
//! Writes static model file (C version)
/*! \todo add third derivatives handling */
void writeStaticCFile(const string &static_basename) const;
//! Writes the Block reordred structure of the model in M output
void writeModelEquationsOrdered_M(Model_Block *ModelBlock, const string &static_basename) const;
//! Writes the code of the Block reordred structure of the model in virtual machine bytecode
void writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, map_idx_type map_idx) const;
//! Computes jacobian and prepares for equation normalization
/*! Using values from initval/endval blocks and parameter initializations:
- computes the jacobian for the model w.r. to contemporaneous variables
- removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff)
*/
void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m, bool dynamic);
void BlockLinear(Model_Block *ModelBlock);
map_idx_type map_idx;
//! Build The incidence matrix form the modeltree
void BuildIncidenceMatrix();
void computeTemporaryTermsOrdered(Model_Block *ModelBlock);
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, map_idx_type &map_idx) const;
//! Get the type corresponding to a derivation ID
SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the lag corresponding to a derivation ID
int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the symbol ID corresponding to a derivation ID
int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Compute the column indices of the static Jacobian
void computeStatJacobianCols();
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(Model_Block *ModelBlock);
//! Collect only the first derivatives
map<pair<int, pair<int, int> >, NodeID> collect_first_order_derivatives_endogenous();
//! Helper for writing the Jacobian elements in MATLAB and C
/*! Writes either (i+1,j+1) or [i+j*no_eq] */
void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Helper for writing the sparse Hessian elements in MATLAB and C
/*! Writes either (i+1,j+1) or [i+j*NNZDerivatives[1]] */
void hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Write chain rule derivative of a recursive equation w.r. to a variable
void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
public:
StaticDllModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
//! Absolute value under which a number is considered to be zero
double cutoff;
//! Compute the minimum feedback set in the static model:
/*! 0 : all endogenous variables are considered as feedback variables
1 : the variables belonging to a non linear equation are considered as feedback variables
2 : the variables belonging to a non normalizable non linear equation are considered as feedback variables
default value = 0 */
int mfs;
//! the file containing the model and the derivatives code
ofstream code_file;
//! Execute computations (variable sorting + derivation)
/*!
\param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
\param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true)
\param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true)
\param paramsDerivatives whether 2nd derivatives w.r. to a pair (endo/exo/exo_det, parameter) should be computed (implies jacobianExo = true)
\param eval_context evaluation context for normalization
\param no_tmp_terms if true, no temporary terms will be computed in the static files
*/
void computingPass(const eval_context_type &eval_context, bool no_tmp_terms, bool block);
//! Complete set to block decompose the model
BlockTriangular block_triangular;
//! Adds informations for simulation in a binary file
void Write_Inf_To_Bin_File(const string &static_basename, const string &bin_basename,
const int &num, int &u_count_int, bool &file_open) const;
//! Writes static model file
void writeStaticFile(const string &basename, bool block) const;
//! Writes LaTeX file with the equations of the static model
void writeLatexFile(const string &basename) const;
//! Writes initializations in oo_.steady_state for the auxiliary variables
void writeAuxVarInitval(ostream &output) const;
virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
};
#endif

1502
StaticModel.cc
View File

File diff suppressed because it is too large Load Diff

217
StaticModel.hh
View File

@ -17,86 +17,166 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _STATICMODEL_HH
#define _STATICMODEL_HH
#ifndef _STATIC_MODEL_HH
#define _STATIC_MODEL_HH
using namespace std;
#include <fstream>
#include "ModelTree.hh"
//! Stores a static model
/*! Derivation IDs are allocated only for endogenous, and are equal to symbol ID in that case */
class StaticModel : public ModelTree
{
private:
//! Normalization of equations
/*! Maps endogenous type specific IDs to equation numbers */
vector<int> endo2eq;
typedef map<pair<int, int>, int> deriv_id_table_t;
//! Maps a pair (symbol_id, lag) to a deriv ID
deriv_id_table_t deriv_id_table;
//! Maps a deriv ID to a pair (symbol_id, lag)
vector<pair<int, int> > inv_deriv_id_table;
//! Block decomposition of the model
/*! List of blocks in topological order. Lists the set of endogenous type specific IDs for each block. */
vector<set<int> > blocks;
//! Temporary terms for the file containing parameters dervicatives
temporary_terms_type params_derivs_temporary_terms;
//! Minimum feedback set for each block
/*! Elements of blocksMFS are subset of elements of blocks */
vector<set<int> > blocksMFS;
//! Temporary terms for block decomposed models
vector<vector<temporary_terms_type> > v_temporary_terms;
//! Variables not in minimum feedback set for each block, sorted in topological order
/*! This is the set difference blocks - blocksMFS. The variables are sorted in topological order. */
vector<vector<int> > blocksRecursive;
vector<temporary_terms_inuse_type> v_temporary_terms_inuse;
//! Jacobian for matrix restricted to MFS
/*! Maps a pair (equation ID, endogenous type specific ID) to the derivative expression. Stores only non-null derivatives. */
map<pair<int, int>, NodeID> blocksMFSJacobian;
typedef map< pair< int, pair< int, int> >, NodeID> first_chain_rule_derivatives_type;
first_chain_rule_derivatives_type first_chain_rule_derivatives;
//! Writes static model file (standard Matlab version)
void writeStaticMFile(ostream &output, const string &func_name) const;
void writeStaticMFile(const string &static_basename) const;
//! Writes static model file (block+MFS version)
void writeStaticBlockMFSFile(ostream &output, const string &func_name) const;
//! Writes the static function calling the block to solve (Matlab version)
void writeStaticBlockMFSFile(const string &basename) const;
//! Computes normalization of the static model
void computeNormalization();
//! Computes blocks of the static model, sorted in topological order
/*! Must be called after computeNormalization() */
void computeSortedBlockDecomposition();
//! Writes static model file (C version)
/*! \todo add third derivatives handling */
void writeStaticCFile(const string &static_basename) const;
//! For each block of the static model, computes minimum feedback set (MFS)
/*! Must be called after computeSortedBlockDecomposition() */
void computeMFS();
//! Writes the Block reordred structure of the model in M output
void writeModelEquationsOrdered_M(const string &dynamic_basename) const;
//! For each block of the static model, computes resursive variables (those not in minimum feedback set), and sort them in topological order
/*! Must be called after computeMFS() */
void computeSortedRecursive();
//! Writes the code of the Block reordred structure of the model in virtual machine bytecode
void writeModelEquationsCodeOrdered(const string file_name, const string bin_basename, map_idx_type map_idx) const;
//! Computes derivatives of each MFS
/*! Must be called after computeSortedRecursive() */
void computeBlockMFSJacobian();
//! Computes the list of equations which are already in normalized form
/*! Returns a multimap mapping endogenous which are normalized (represented by their type specific ID) to the equation(s) which define it */
void computeNormalizedEquations(multimap<int, int> &endo2eqs) const;
//! Helper for writing model local variables in block+MFS mode
/*!
Write the definition of model local variables which are used in expr, except those in local_var_written.
Add these variables to local_var_written at the end.
//! Computes jacobian and prepares for equation normalization
/*! Using values from initval/endval blocks and parameter initializations:
- computes the jacobian for the model w.r. to contemporaneous variables
- removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff)
*/
void writeLocalVars(ostream &output, NodeID expr, set<int> &local_var_written) const;
void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m, bool dynamic);
map_idx_type map_idx;
void computeTemporaryTermsOrdered();
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, map_idx_type &map_idx) const;
//! Get the type corresponding to a derivation ID
virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the lag corresponding to a derivation ID
virtual int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the symbol ID corresponding to a derivation ID
virtual int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Compute the column indices of the static Jacobian
void computeStatJacobianCols();
//! return a map on the block jacobian
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> get_Derivatives(int block);
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(t_blocks_derivatives &blocks_derivatives);
//! Collect only the first derivatives
map<pair<int, pair<int, int> >, NodeID> collect_first_order_derivatives_endogenous();
//! Helper for writing the Jacobian elements in MATLAB and C
/*! Writes either (i+1,j+1) or [i+j*no_eq] */
void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Helper for writing the sparse Hessian elements in MATLAB and C
/*! Writes either (i+1,j+1) or [i+j*NNZDerivatives[1]] */
void hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Write chain rule derivative of a recursive equation w.r. to a variable
void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
//! Collecte the derivatives w.r. to endogenous of the block, to endogenous of previouys blocks and to exogenous
void collect_block_first_order_derivatives();
protected:
//! Indicate if the temporary terms are computed for the overall model (true) or not (false). Default value true
bool global_temporary_terms;
//! vector of block reordered variables and equations
vector<int> equation_reordered, variable_reordered, inv_equation_reordered, inv_variable_reordered;
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
t_equation_type_and_normalized_equation equation_type_and_normalized_equation;
//! for each block contains pair< Simulation_Type, pair < Block_Size, Recursive_part_Size > >
t_block_type_firstequation_size_mfs block_type_firstequation_size_mfs;
//! for all blocks derivatives description
t_blocks_derivatives blocks_derivatives;
//! The jacobian without the elements below the cutoff
dynamic_jacob_map dynamic_jacobian;
//! Vector indicating if the block is linear in endogenous variable (true) or not (false)
vector<bool> blocks_linear;
//! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), NodeID>
typedef map<pair< int, pair<int, int> >, NodeID> t_derivative;
//! Vector of derivative for each blocks
vector<t_derivative> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
//!List for each block and for each lag-leag all the other endogenous variables and exogenous variables
typedef set<int> t_var;
typedef map<int, t_var> t_lag_var;
vector<t_lag_var> other_endo_block, exo_block, exo_det_block;
//!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
vector<pair<int, int> > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
public:
StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
//! Execute computations (derivation)
/*!
\param block whether block decomposition and minimum feedback set should be computed
\param hessian whether Hessian (w.r. to endogenous only) should be computed
\param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
void computingPass(bool block, bool hessian, bool no_tmp_terms);
//! Writes information on block decomposition when relevant
void writeOutput(ostream &output, bool block) const;
//! Absolute value under which a number is considered to be zero
double cutoff;
//! Compute the minimum feedback set in the static model:
/*! 0 : all endogenous variables are considered as feedback variables
1 : the variables belonging to a non linear equation are considered as feedback variables
2 : the variables belonging to a non normalizable non linear equation are considered as feedback variables
default value = 0 */
int mfs;
//! the file containing the model and the derivatives code
ofstream code_file;
//! Execute computations (variable sorting + derivation)
/*!
\param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
\param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true)
\param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true)
\param paramsDerivatives whether 2nd derivatives w.r. to a pair (endo/exo/exo_det, parameter) should be computed (implies jacobianExo = true)
\param eval_context evaluation context for normalization
\param no_tmp_terms if true, no temporary terms will be computed in the static files
*/
void computingPass(const eval_context_type &eval_context, bool no_tmp_terms, bool hessian, bool block);
//! Adds informations for simulation in a binary file
void Write_Inf_To_Bin_File(const string &static_basename, const string &bin_basename, const int &num,
int &u_count_int, bool &file_open) const;
//! Writes static model file
void writeStaticFile(const string &basename, bool block) const;
void writeStaticFile(const string &basename, bool block, bool bytecode) const;
//! Writes LaTeX file with the equations of the static model
void writeLatexFile(const string &basename) const;
@ -105,6 +185,39 @@ public:
void writeAuxVarInitval(ostream &output) const;
virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
//! Return the number of blocks
virtual unsigned int getNbBlocks() const {return(block_type_firstequation_size_mfs.size());};
//! Determine the simulation type of each block
virtual BlockSimulationType getBlockSimulationType(int block_number) const {return(block_type_firstequation_size_mfs[block_number].first.first);};
//! Return the first equation number of a block
virtual unsigned int getBlockFirstEquation(int block_number) const {return(block_type_firstequation_size_mfs[block_number].first.second);};
//! Return the size of the block block_number
virtual unsigned int getBlockSize(int block_number) const {return(block_type_firstequation_size_mfs[block_number].second.first);};
//! Return the number of feedback variable of the block block_number
virtual unsigned int getBlockMfs(int block_number) const {return(block_type_firstequation_size_mfs[block_number].second.second);};
//! Return the maximum lag in a block
virtual unsigned int getBlockMaxLag(int block_number) const {return(block_lag_lead[block_number].first);};
//! Return the maximum lead in a block
virtual unsigned int getBlockMaxLead(int block_number) const {return(block_lag_lead[block_number].second);};
//! Return the type of equation (equation_number) belonging to the block block_number
virtual EquationType getBlockEquationType(int block_number, int equation_number) const {return( equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first);};
//! Return true if the equation has been normalized
virtual bool isBlockEquationRenormalized(int block_number, int equation_number) const {return( equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);};
//! Return the NodeID of the equation equation_number belonging to the block block_number
virtual NodeID getBlockEquationNodeID(int block_number, int equation_number) const {return( equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);};
//! Return the NodeID of the renormalized equation equation_number belonging to the block block_number
virtual NodeID getBlockEquationRenormalizedNodeID(int block_number, int equation_number) const {return( equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);};
//! Return the original number of equation equation_number belonging to the block block_number
virtual int getBlockEquationID(int block_number, int equation_number) const {return( equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]);};
//! Return the original number of variable variable_number belonging to the block block_number
virtual int getBlockVariableID(int block_number, int variable_number) const {return( variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);};
//! Return the position of equation_number in the block number belonging to the block block_number
virtual int getBlockInitialEquationID(int block_number, int equation_number) const {return((int)inv_equation_reordered[equation_number] - (int)block_type_firstequation_size_mfs[block_number].first.second);};
//! Return the position of variable_number in the block number belonging to the block block_number
virtual int getBlockInitialVariableID(int block_number, int variable_number) const {return((int)inv_variable_reordered[variable_number] - (int)block_type_firstequation_size_mfs[block_number].first.second);};
};
#endif