Commit 1a058fea authored by ferhat's avatar ferhat
Browse files

- Bugs correction in the new block decomposition (incorporating the feedback variables)

- First draft of DynamicModel.cc files with feedback variables. 
TODO : 
- reduction of the Jacobian matrix 
- symbolic normalization of equations
- application to the binary code evaluation (simulate.dll).

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2678 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 2d846f0d
This diff is collapsed.
......@@ -27,8 +27,7 @@
#include "ModelNormalization.hh"
#include "ModelBlocks.hh"
#include "IncidenceMatrix.hh"
#include "Modeltree.hh"
......@@ -37,7 +36,9 @@
//! Sparse matrix of double to store the values of the Jacobian
typedef map<pair<int ,int >,double> jacob_map;
typedef vector<pair<BlockSimulationType, int> > t_type;
typedef vector<pair<BlockSimulationType, pair<int, int> > > t_type;
typedef vector<pair<EquationType, int> > t_etype;
//! Creates the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition
class BlockTriangular
......@@ -46,13 +47,15 @@ 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);
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);
//! Finds a matching between equations and endogenous variables
bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool 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, bool verbose_) const;
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_) const;
//! determine the type of each equation of the model (couble evaluated or need to be solved)
t_etype Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, int >, NodeID> &first_cur_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! 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_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type);
public:
const SymbolTable &symbol_table;
BlockTriangular(const SymbolTable &symbol_table_arg);
......@@ -63,8 +66,8 @@ public:
Blocks blocks;
Normalization normalization;
IncidenceMatrix incidencematrix;
void Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m, vector<BinaryOpNode *> equations);
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);
void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, int >, NodeID> &first_cur_endo_derivatives);
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, int >, NodeID> &first_cur_endo_derivatives);
vector<int> Index_Equ_IM;
vector<int> Index_Var_IM;
int prologue, epilogue;
......@@ -128,5 +131,17 @@ public:
break;
}
};
inline static std::string c_Equation_Type(int type)
{
char c_Equation_Type[5][13]=
{
"E_UNKNOWN ",
"E_EVALUATE ",
"E_EVALUATE_R",
"E_EVALUATE_S",
"E_SOLVE "
};
return(c_Equation_Type[type]);
};
};
#endif
......@@ -49,6 +49,17 @@ enum BlockType
SIMULTAN = 3 //<! Simultaneous time unseparable block
};
enum EquationType
{
E_UNKNOWN, //!< Unknown equation type
E_EVALUATE, //!< Simple evaluation, normalized variable on left-hand side
E_EVALUATE_R, //!< Simple evaluation, normalized variable on right-hand side
E_EVALUATE_S, //!< Simple evaluation, normalize using the first order derivative which does not involve the normalized variable
E_SOLVE //!< No simple evaluation of the equation, it has to be solved
};
enum BlockSimulationType
{
UNKNOWN, //!< Unknown simulation type
......
......@@ -280,7 +280,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " if(jacobian_eval)\n";
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size*(1+ModelBlock->Block_List[j].Max_Lag_Endo+ModelBlock->Block_List[j].Max_Lead_Endo) << ", " << nze << ");\n";
//output << " g1_x=spalloc(" << ModelBlock->Block_List[j].Size << ", " << (ModelBlock->Block_List[j].nb_exo + ModelBlock->Block_List[j].nb_exo_det)*(1+ModelBlock->Block_List[j].Max_Lag_Exo+ModelBlock->Block_List[j].Max_Lead_Exo) << ", " << nze_exo << ");\n";
output << " g1_o=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].nb_other_endo*(1+ModelBlock->Block_List[j].Max_Lag_Other_Endo+ModelBlock->Block_List[j].Max_Lead_Other_Endo) << ", " << nze_other_endo << ");\n";
//output << " g1_o=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].nb_other_endo*(1+ModelBlock->Block_List[j].Max_Lag_Other_Endo+ModelBlock->Block_List[j].Max_Lead_Other_Endo) << ", " << nze_other_endo << ");\n";
output << " end;\n";
}
else
......@@ -357,15 +357,30 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
evaluation:
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
output << " ";
output << tmp_output.str();
output << " = ";
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
{
output << tmp_output.str();
output << " = ";
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
}
else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_R)
{
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << " = ";
output << tmp_output.str();
output << "; %reversed " << ModelBlock->Block_List[j].Equation_Type[i] << " \n";
}
else
{
cerr << "Type missmatch for equation " << ModelBlock->Block_List[j].Equation[i]+1 << "\n";
exit(-1);
}
output << ";\n";
break;
case EVALUATE_BACKWARD_R:
break; /*case EVALUATE_BACKWARD_R:
case EVALUATE_FORWARD_R:
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
......@@ -374,21 +389,25 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " = ";
lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << ";\n";
break;
break;*/
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
if(i<ModelBlock->Block_List[j].Nb_Recursives)
goto evaluation;
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
output << " " << "residual(" << i+1 << ") = (";
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
output << " " << "residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ") = (";
goto end;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
if(i<ModelBlock->Block_List[j].Nb_Recursives)
goto evaluation;
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << "+Per_J_) = -residual(" << i+1 << ", it_)";
output << " residual(" << i+1 << ", it_) = (";
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_) = -residual(" << i+1 << ", it_)";
output << " residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ", it_) = (";
goto end;
default:
end:
......@@ -2164,6 +2183,31 @@ DynamicModel::BlockLinear(Model_Block *ModelBlock)
}
map<pair<int, int >, NodeID>
DynamicModel::collect_first_order_derivatives_current_endogenous()
{
map<pair<int, int >, NodeID> curr_endo_derivatives;
for (first_derivatives_type::iterator it2 = first_derivatives.begin();
it2 != first_derivatives.end(); it2++)
{
if(getTypeByDerivID(it2->first.second)==eEndogenous)
{
int eq = it2->first.first;
int var=symbol_table.getTypeSpecificID(getSymbIDByDerivID(it2->first.second));
int lag=getLagByDerivID(it2->first.second);
/*cout << "eq=" << eq << " var=" << symbol_table.getName(var) << " (" << var << ") lag=" << lag;
ExprNodeOutputType output_type=oMatlabDynamicModelSparse;
const temporary_terms_type temporary_terms;
cout << " derivative = ";
(it2->second)->writeOutput(cout, output_type, temporary_terms);
cout << "\n";*/
if(lag==0)
curr_endo_derivatives[make_pair(eq, var)] = it2->second;
}
}
return curr_endo_derivatives;
}
void
......@@ -2225,9 +2269,10 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
cout << "The gross incidence matrix \n";
block_triangular.incidencematrix.Print_IM(eEndogenous);
}
block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations);
t_etype equation_simulation_type;
map<pair<int, int >, NodeID> first_cur_endo_derivatives = collect_first_order_derivatives_current_endogenous();
block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations, equation_simulation_type, first_cur_endo_derivatives);
BlockLinear(block_triangular.ModelBlock);
if (!no_tmp_terms)
computeTemporaryTermsOrdered(block_triangular.ModelBlock);
......
......@@ -109,7 +109,8 @@ private:
void computeParamsDerivatives();
//! Computes temporary terms for the file containing parameters derivatives
void computeParamsDerivativesTemporaryTerms();
//! Collect only the first derivatives
map<pair<int, int >, NodeID> collect_first_order_derivatives_current_endogenous();
public:
DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
......
......@@ -25,7 +25,9 @@ using namespace std;
#include <cstdlib>
#include <cstring>
#ifndef PACKAGE_VERSION
#define PACKAGE_VERSION 4.
#endif
#include "macro/MacroDriver.hh"
/* Prototype for second part of main function
......
......@@ -383,7 +383,7 @@ struct IM_compact
//! One block of the model
struct Block
{
int Size, Sized, nb_exo, nb_exo_det, nb_other_endo;
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;
......@@ -392,7 +392,8 @@ struct Block
int Max_Lag_Exo, Max_Lead_Exo;
bool is_linear;
int *Equation, *Own_Derivative;
int *Variable, *Other_Endogenous, *Exogenous;
EquationType *Equation_Type;
int *Variable, *Other_Endogenous, *Exogenous, *Equation_Type_Var;
temporary_terms_type **Temporary_Terms_in_Equation;
//temporary_terms_type *Temporary_terms;
temporary_terms_inuse_type *Temporary_InUse;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment