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

- Complete implementation of feedback variables in dynamic model with sparse option 

- Normalization of linear in endogenous variable equations 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2834 ac1d8469-bf42-47a9-8791-bf33cf982152
issue#70
ferhat 2009-07-10 15:05:09 +00:00
parent 81b6823bb3
commit 8421712f80
7 changed files with 453 additions and 312 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

113
BlockTriangular.cc
View File

@ -167,8 +167,11 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog
return check;
}
t_vtype
BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue) const
BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) const
{
int nb_endo = symbol_table.endo_nbr();
vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
@ -190,9 +193,9 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int
variable_blck[Index_Var_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
equation_blck[Index_Equ_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
}
//cout << "equation_blck[" << Index_Equ_IM[i] << "]=" << equation_blck[Index_Equ_IM[i]] << " variable_blck[" << Index_Var_IM[i] << "] = " << variable_blck[Index_Var_IM[i]] << "\n";
Variable_Type[i] = make_pair(0, 0);
}
equation_lead_lag = Variable_Type;
for (int k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
{
bool *Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
@ -203,19 +206,22 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int
int i_1 = Index_Var_IM[i];
for (int j = 0; j < nb_endo; j++)
{
if (Cur_IM[i_1 + Index_Equ_IM[ j] * nb_endo] and variable_blck[i_1] == equation_blck[Index_Equ_IM[ j]])
int j_l = Index_Equ_IM[ j];
if (Cur_IM[i_1 + Index_Equ_IM[ j] * nb_endo] and variable_blck[i_1] == equation_blck[j_l])
{
if (k > Variable_Type[i_1].second)
Variable_Type[i_1] = make_pair(Variable_Type[i_1].first, k);
if (k < -Variable_Type[i_1].first)
Variable_Type[i_1] = make_pair(-k, Variable_Type[i_1].second);
if (k > equation_lead_lag[j_l].second)
equation_lead_lag[j_l] = make_pair(equation_lead_lag[j_l].first, k);
if (k < -equation_lead_lag[j_l].first)
equation_lead_lag[j_l] = make_pair(-k, equation_lead_lag[j_l].second);
}
}
}
}
}
/*for(int i = 0; i < nb_endo; i++)
cout << "Variable_Type[" << Index_Equ_IM[i] << "].first = " << Variable_Type[Index_Equ_IM[i]].first << " Variable_Type[" << Index_Equ_IM[i] << "].second=" << Variable_Type[Index_Equ_IM[i]].second << "\n";*/
return (Variable_Type);
}
@ -255,7 +261,9 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
components_set[component[i]].first.insert(i);
}
V_Variable_Type = Get_Variable_LeadLag_By_Block(component, num, prologue, epilogue);
t_vtype equation_lead_lag;
V_Variable_Type = Get_Variable_LeadLag_By_Block(component, num, prologue, epilogue, equation_lead_lag);
vector<int> tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM);
int order = prologue;
@ -265,7 +273,8 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
//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
for (int i = 0; i < n; i++)
if (Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or V_Variable_Type[Index_Var_IM[i+prologue]].second >= 0 or V_Variable_Type[Index_Var_IM[i+prologue]].first > 0)
if (Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or V_Variable_Type[Index_Var_IM[i+prologue]].second > 0 or V_Variable_Type[Index_Var_IM[i+prologue]].first > 0
or equation_lead_lag[Index_Equ_IM[i+prologue]].second > 0 or equation_lead_lag[Index_Equ_IM[i+prologue]].first > 0)
add_edge(i, i, G2);
//For each block, the minimum set of feedback variable is computed
@ -297,13 +306,14 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
Index_Var_IM[order] = tmp_Index_Var_IM[v_index[vertex(*its, G)]+prologue];
order++;
}
}
free(AMp);
free(SIM);
}
void
BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block *ModelBlock, t_etype &Equation_Type, int recurs_Size)
BlockTriangular::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)
{
int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1, li;
int *tmp_size, *tmp_size_other_endo, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_other_endo, *tmp_exo, tmp_nb_other_endo, tmp_nb_exo, nb_lead_lag_endo;
@ -318,11 +328,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].Type = type;
ModelBlock->Block_List[count_Block].Nb_Recursives = recurs_Size;
ModelBlock->Block_List[count_Block].Temporary_InUse = new temporary_terms_inuse_type();
ModelBlock->Block_List[count_Block].Chaine_Rule_Derivatives = new chaine_rule_derivatives_type();
ModelBlock->Block_List[count_Block].Temporary_InUse->clear();
ModelBlock->Block_List[count_Block].Simulation_Type = SimType;
ModelBlock->Block_List[count_Block].Equation = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
ModelBlock->Block_List[count_Block].Equation_Type = (EquationType *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(EquationType));
ModelBlock->Block_List[count_Block].Equation_Normalized = (NodeID *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(NodeID));
ModelBlock->Block_List[count_Block].Equation_Normalized = (NodeID*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(NodeID));
ModelBlock->Block_List[count_Block].Variable = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation = (temporary_terms_type **) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(temporary_terms_type));
ModelBlock->Block_List[count_Block].Own_Derivative = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
@ -668,6 +679,7 @@ BlockTriangular::Free_Block(Model_Block *ModelBlock) const
delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i];
free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation);
delete (ModelBlock->Block_List[blk].Temporary_InUse);
delete ModelBlock->Block_List[blk].Chaine_Rule_Derivatives;
}
free(ModelBlock->Block_List);
free(ModelBlock);
@ -735,7 +747,7 @@ BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations,
BlockTriangular::Recompute_Deriavtives_Respect_to_Feedback_Variables(
*/
t_type
BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type)
BlockTriangular::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)
{
int i = 0;
int count_equ = 0, blck_count_simult = 0;
@ -842,6 +854,66 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
return (Type);
}
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int>
BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
{
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int> Derivatives;
Derivatives.clear();
int nb_endo = symbol_table.endo_nbr();
/*ModelBlock.Block_List[Blck].first_order_determinstic_simulation_derivatives = new*/
for(int lag = -ModelBlock->Block_List[blck].Max_Lag; lag <= ModelBlock->Block_List[blck].Max_Lead; lag++)
{
bool *IM=incidencematrix.Get_IM(lag, eEndogenous);
if(IM)
{
for(int eq = 0; eq < ModelBlock->Block_List[blck].Size; eq++)
{
int eqr = ModelBlock->Block_List[blck].Equation[eq];
for(int var = 0; var < ModelBlock->Block_List[blck].Size; var++)
{
int varr = ModelBlock->Block_List[blck].Variable[var];
/*cout << "IM=" << IM << "\n";
cout << "varr=" << varr << " eqr=" << eqr << " lag=" << lag << "\n";*/
if(IM[varr+eqr*nb_endo])
{
/*if(eq<ModelBlock->Block_List[blck].Nb_Recursives and var<ModelBlock->Block_List[blck].Nb_Recursives)
{*/
bool OK = true;
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int>::const_iterator its = Derivatives.find(make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag)));
if(its!=Derivatives.end())
{
if(its->second == 2)
OK=false;
}
if(OK)
{
if (ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_List[blck].Nb_Recursives)
Derivatives[make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag))] = 1;
else
Derivatives[make_pair(make_pair(eqr, eq),make_pair(make_pair(varr, var), lag))] = 0;
}
/*}
else if(eq<ModelBlock->Block_List[blck].Nb_Recursives and var<ModelBlock->Block_List[blck].Nb_Recursives)*/
if(var<ModelBlock->Block_List[blck].Nb_Recursives)
{
int eqs = ModelBlock->Block_List[blck].Equation[var];
for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; vars<ModelBlock->Block_List[blck].Size; vars++)
{
int varrs = ModelBlock->Block_List[blck].Variable[vars];
if(Derivatives.find(make_pair(make_pair(eqs, var), make_pair(make_pair(varrs, vars), lag)))!=Derivatives.end())
Derivatives[make_pair(make_pair(eqr, eq),make_pair(make_pair(varrs, vars), lag))] = 2;
}
}
}
}
}
}
}
return(Derivatives);
}
void
BlockTriangular::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 &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives)
{
@ -850,6 +922,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
int count_Block, count_Equ;
bool *SIM0, *SIM00;
SIM0 = (bool *) malloc(n * n * sizeof(bool));
memcpy(SIM0, IM_0, n*n*sizeof(bool));
Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0);
@ -873,7 +946,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
//double bi=1e-13;
int suppressed = 0;
vector<int> Index_Equ_IM_save(Index_Equ_IM);
while (!OK && bi > 1e-14)
while (!OK && bi > 1e-19)
{
int suppress = 0;
Index_Equ_IM = Index_Equ_IM_save;
@ -907,7 +980,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
//bi/=1.07;
bi /= 2;
counted++;
if (bi > 1e-14)
if (bi > 1e-19)
free(SIM00);
}
if (!OK)
@ -921,12 +994,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
if (prologue+epilogue < n)
Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(IM, n, prologue, epilogue, Index_Equ_IM, Index_Var_IM, blocks, V_Equation_Type, false);
t_type Type = Reduce_Blocks_and_type_determination(prologue, epilogue, blocks, equations, V_Equation_Type);
t_type Type = Reduce_Blocks_and_type_determination(prologue, epilogue, blocks, equations, V_Equation_Type, Index_Var_IM, Index_Equ_IM);
i = 0;
j = 0;
Nb_SimulBlocks = 0;
int Nb_fv = 0;
int Nb_feedback_variable = 0;
for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
{
if (it->first == SOLVE_FORWARD_COMPLETE || it->first == SOLVE_BACKWARD_COMPLETE || it->first == SOLVE_TWO_BOUNDARIES_COMPLETE)
@ -935,7 +1008,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
if (it->second.first > j)
{
j = it->second.first;
Nb_fv = blocks[Nb_SimulBlocks-1].second;
Nb_feedback_variable = blocks[Nb_SimulBlocks-1].second;
}
}
}
@ -945,7 +1018,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
cout << Nb_TotalBlocks << " block(s) found:\n";
cout << " " << Nb_RecursBlocks << " recursive block(s) and " << blocks.size() << " simultaneous block(s). \n";
cout << " the largest simultaneous block has " << j << " equation(s)\n"
<<" and " << Nb_fv << " feedback variable(s).\n";
<<" and " << Nb_feedback_variable << " feedback variable(s).\n";
ModelBlock->Size = Nb_TotalBlocks;
ModelBlock->Periods = periods;
@ -953,6 +1026,8 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
count_Equ = count_Block = 0;
for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
{
if (count_Equ < prologue)
@ -961,10 +1036,12 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
if (it->second.first == 1)
Btype = PROLOGUE;
else
Btype = SIMULTANS;
{
Btype = SIMULTANS;
}
else
Btype = EPILOGUE;
Allocate_Block(it->second.first, &count_Equ, count_Block++, Btype, it->first, ModelBlock, V_Equation_Type, it->second.second);
Allocate_Block(it->second.first, &count_Equ, count_Block++, Btype, it->first, ModelBlock, V_Equation_Type, it->second.second, Index_Var_IM, Index_Equ_IM);
}
}

10
BlockTriangular.hh
View File

@ -44,6 +44,8 @@ typedef vector<pair< int, int> > t_vtype;
typedef set<int> temporary_terms_inuse_type;
typedef vector<pair< pair<int, int>, pair<int, pair<int, int> > > > chaine_rule_derivatives_type;
//! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
struct IM_compact
{
@ -71,6 +73,7 @@ struct Block
//temporary_terms_type *Temporary_terms;
temporary_terms_inuse_type *Temporary_InUse;
IM_compact *IM_lead_lag;
chaine_rule_derivatives_type *Chaine_Rule_Derivatives;
int Code_Start, Code_Length;
};
@ -92,7 +95,7 @@ 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);
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, 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
@ -100,9 +103,9 @@ private:
//! 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);
//! 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);
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) const;
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;
@ -115,6 +118,7 @@ public:
//! Frees the Model structure describing the content of each block
void Free_Block(Model_Block* ModelBlock) const;
map<pair<pair<int, int>, pair<pair<int, 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);

603
DynamicModel.cc
View File

@ -105,6 +105,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
ostringstream tmp_output;
BinaryOpNode *eq_node;
first_derivatives_type::const_iterator it;
first_chaine_rule_derivatives_type::const_iterator it_chr;
ostringstream tmp_s;
temporary_terms.clear();
@ -132,6 +133,16 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
//if(it!=first_derivatives.end())
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
{
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
eq=it.first.second;
var=it.second.second.first;
lag=it.second.second.second;
it_chr=first_chaine_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
//it_chr->second->writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
it_chr->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
@ -187,6 +198,16 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
//if(it!=first_derivatives.end())
it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
}
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
{
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
eq=it.first.second;
var=it.second.second.first;
lag=it.second.second.second;
it_chr=first_chaine_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
//it_chr->second->writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
it_chr->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
@ -240,13 +261,13 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
ofstream output;
//temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
int nze, nze_exo, nze_other_endo;
map<int, NodeID> recursive_variables;
//map<int, NodeID> recursive_variables;
vector<int> feedback_variables;
//----------------------------------------------------------------------
//For each block
for (j = 0;j < ModelBlock->Size;j++)
{
recursive_variables.clear();
//recursive_variables.clear();
feedback_variables.clear();
//For a block composed of a single equation determines wether we have to evaluate or to solve the equation
nze = nze_exo = nze_other_endo = 0;
@ -388,14 +409,17 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
tmp_output.str("");
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
if((ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD) and (i<ModelBlock->Block_List[j].Nb_Recursives))
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
else
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
switch (ModelBlock->Block_List[j].Simulation_Type)
{
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 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
evaluation: if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
output << " ";
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
{
@ -418,9 +442,19 @@ evaluation:
{
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << "\n ";
temporary_terms_type tt2;
tt2.clear();
ModelBlock->Block_List[j].Equation_Normalized[i]->writeOutput(output , oMatlabDynamicModelSparse, temporary_terms/*tt2*/);
/*temporary_terms_type tt2;
tt2.clear();*/
tmp_output.str("");
eq_node = (BinaryOpNode *)ModelBlock->Block_List[j].Equation_Normalized[i];
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
if(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD)
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
else
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
output << tmp_output.str();
output << " = ";
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
}
}
else
@ -436,10 +470,10 @@ evaluation:
case SOLVE_FORWARD_COMPLETE:
if (i<ModelBlock->Block_List[j].Nb_Recursives)
{
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
/*if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = ModelBlock->Block_List[j].Equation_Normalized[i];
else
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];*/
goto evaluation;
}
feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
@ -451,10 +485,10 @@ evaluation:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
if (i<ModelBlock->Block_List[j].Nb_Recursives)
{
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
/*if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = ModelBlock->Block_List[j].Equation_Normalized[i];
else
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];*/
goto evaluation;
}
feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
@ -618,79 +652,44 @@ end:
m=ModelBlock->Block_List[j].Max_Lag;
//cout << "\nDerivatives in Block " << j << "\n";
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
{
//Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
int eqr=it.first.first;
int eq=it.first.second;
int varr=it.second.first;
int var=it.second.second.first;
k=it.second.second.second;
/*for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
bool derivative_exist;
ostringstream tmp_output;
if (eqr<ModelBlock->Block_List[j].Nb_Recursives)
{
if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
{
/*tmp_output << " g1_tmp_r(" << eqr+1 << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
NodeID tmp_n;
if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE)
tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
else
tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];
int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),0);
NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables);
ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << " %1 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k
<< ") " << var+1
<< ", equation=" << eq+1 << endl;*/
}
}
else
{
if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
{
{*/
output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
/*writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);*/
/*if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, 0, 0, recursive_variables, feedback_variables);
else
derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, 0, 0, recursive_variables, feedback_variables);
//if (derivative_exist)
output << tmp_output.str() << ";";*/
NodeID tmp_n;
//if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
/*else
tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];*/
//cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
//cout << "derivaive eq=" << eq << " var=" << var << " k0=" << k << "\n";
int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),0);
NodeID ChaineRule_Derivative = tmp_n->getChainRuleDerivative(deriv_id, recursive_variables);
ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
/*}
}
/*if (eqr<ModelBlock->Block_List[j].Nb_Recursives or varr<ModelBlock->Block_List[j].Nb_Recursives)
output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), 0, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
<< "(0) " << var+1
<< ", equation=" << eq+1 << endl;*/
}
}*/
output << " end;\n";
//output << " ya = y;\n";
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
output << " if ~jacobian_eval" << endl;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
k=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
@ -698,163 +697,86 @@ end:
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
bool derivative_exist;
ostringstream tmp_output;
//cout << "ModelBlock->Block_List[" << j << "].Nb_Recursives=" << ModelBlock->Block_List[j].Nb_Recursives << "\n";
if (eqr<ModelBlock->Block_List[j].Nb_Recursives)
{
/*if (varr<ModelBlock->Block_List[j].Nb_Recursives)
{
if (k==0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
<< ", " << varr+1
<< "+" << ModelBlock->Block_List[j].Size*ModelBlock->Block_List[j].Max_Lag << ")*y(it_, " << var+1 << ")";
else if (k>0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
<< ", " << varr+1
<< "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ")*y(it_+" << k << ", " << var+1 << ")";
else if (k<0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
<< ", " << varr+1
<< "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ")*y(it_" << k << ", " << var+1 << ")";
if (k==0)
tmp_output << " g1_tmp_r(" << eqr+1 << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = ";
else if (k>0)
tmp_output << " g1_tmp_r(" << eqr+1 << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
else if (k<0)
tmp_output << " g1_tmp_r(" << eqr+1 << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE)
derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->get_arg2()->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
else
{
BinaryOpNode* tt = (BinaryOpNode*)ModelBlock->Block_List[j].Equation_Normalized[eqr];
derivative_exist = tt->get_arg2()->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
}
if (derivative_exist)
output << tmp_output.str() << ";";
else
{
//output << "1" << ";";
if (ModelBlock->Block_List[j].Equation_Type[eqr] != E_EVALUATE)
{
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];*/
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
{
//Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
int eqr=it.first.first;
int eq=it.first.second;
int varr=it.second.first;
int var=it.second.second.first;
k=it.second.second.second;
}
}
//writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << " %1 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << " derivative_exist=" << derivative_exist << " varr+1=" << varr+1 << endl;
}*/
}
else
//bool derivative_exist;
ostringstream tmp_output;
/*if (eqr>=ModelBlock->Block_List[j].Nb_Recursives)
{
if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
{*/
/*for(int equation = ModelBlock->Block_List[j].Nb_Recursives; equation<ModelBlock->Block_List[j].Size; equation++)
{
int eq = ModelBlock->Block_List[j].Equation[equation];
int eqr = equation - ModelBlock->Block_List[j].Nb_Recursives;
for(int variable = ModelBlock->Block_List[j].Nb_Recursives; variable<ModelBlock->Block_List[j].Size; variable++)
{
if (k==0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_K_)*y(it_, " << var+1 << ")";
else if (k==1)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_y_)*y(it_+1, " << var+1 << ")";
else if (k>0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")";
else if (k<0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")";
if (k==0)
tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_K_) = ";
else if (k==1)
tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_y_) = ";
else if (k>0)
tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_+" << k-1 << ")) = ";
else if (k<0)
tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_" << k-1 << ")) = ";
/*NodeID tmp_n;
if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
else
tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];*/
/*int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),k);
//cout << "-------------------------------------------------------------------------------------\n";
//cout << "derivaive eq=" << eq << " var=" << var << " k=" << k << "\n";
//output << " " << tmp_output.str();
map<int, NodeID> recursive_variables_save(recursive_variables);
NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables, var, k);
recursive_variables = recursive_variables_save;
//ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
*/
/*if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
else
derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
if (derivative_exist)
output << tmp_output.str() << ";";*/
int var = ModelBlock->Block_List[j].Variable[variable];
int varr = variable - ModelBlock->Block_List[j].Nb_Recursives;*/
//cout << "eqr=" << eqr << " varr=" << varr;
//cout << "k=" << k << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << " ModelBlock->Block_List[j].Equation[eq]=" << ModelBlock->Block_List[j].Equation[eq] << "\n";
if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives)
{
if (k==0)
Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_K_)*y(it_, " << varr+1 << ")";
else if (k==1)
Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_y_)*y(it_+1, " << varr+1 << ")";
else if (k>0)
Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << varr+1 << ")";
else if (k<0)
Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << varr+1 << ")";
if (k==0)
tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< var+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_K_) = ";
else if (k==1)
tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< var+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_y_) = ";
else if (k>0)
tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< var+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_+" << k-1 << ")) = ";
else if (k<0)
tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< var+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_" << k-1 << ")) = ";
output << " " << tmp_output.str();
writeChaineRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
/*output << tmp_output.str();*/
//output << ";";
//output << "\n%";
output << tmp_output.str();
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
/*else
{
if (k==0)
tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = ";
else if (k>0)
tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
else if (k<0)
tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
else
derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
if (derivative_exist)
output << tmp_output.str() << ";";
if (k==0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag
<< ")*y(it_, " << var+1 << ")";
else if (k>0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag
<< ")*y(it_+" << k << ", " << var+1 << ")";
else if (k<0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag
<< ")*y(it_" << k << ", " << var+1 << ")";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}*/
}
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
<< "(" << k << ") " << varr+1
<< ", equation=" << eqr+1 << endl;
}
//cout << " done\n";
/* }
}*/
#ifdef CONDITION
output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
output << " condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
#endif
}
//}
}
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
@ -975,12 +897,10 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
bool lhs_rhs_done;
Uff Uf[symbol_table.endo_nbr()];
map<NodeID, int> reference_count;
map<int,int> ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number;
//map<int,int> ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number;
vector<int> feedback_variables;
int prev_Simulation_Type=-1;
//SymbolicGaussElimination SGE;
bool file_open=false;
//temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
//----------------------------------------------------------------------
string main_name=file_name;
main_name+=".cod";
code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate );
@ -993,52 +913,26 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
code_file.write(&FDIMT, sizeof(FDIMT));
k=temporary_terms.size();
code_file.write(reinterpret_cast<char *>(&k),sizeof(k));
//search for successive and identical blocks
i=k=k0=0;
ModelBlock_Aggregated_Count=-1;
for (j = 0;j < ModelBlock->Size;j++)
{
if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type)
&& (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
/*||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R */))
{
}
else
{
k=k0=0;
ModelBlock_Aggregated_Count++;
}
k0+=ModelBlock->Block_List[j].Size;
ModelBlock_Aggregated_Number[ModelBlock_Aggregated_Count]=k0;
ModelBlock_Aggregated_Size[ModelBlock_Aggregated_Count]=++k;
prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
}
ModelBlock_Aggregated_Count++;
ModelBlock_Aggregated_Count = ModelBlock->Size;
//For each block
j=0;
for (k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++)
for (j = 0; j < ModelBlock->Size ;j++)
{
k1=j;
if (k0>0)
feedback_variables.clear();
if (j>0)
code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK));
v=ModelBlock_Aggregated_Number[k0];
v=ModelBlock->Block_List[j].Size;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
v=ModelBlock->Block_List[j].Simulation_Type;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
for (k=0; k<ModelBlock_Aggregated_Size[k0]; k++)
for (i=0; i < ModelBlock->Block_List[j].Size;i++)
{
for (i=0; i < ModelBlock->Block_List[j].Size;i++)
{
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i]));
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i]));
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i]));
}
j++;
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i]));
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i]));
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i]));
}
j=k1;
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
{
@ -1051,8 +945,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
v=block_triangular.ModelBlock->Block_List[j].Max_Lead;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
//if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
//{
int u_count_int=0;
Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open,
ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE);
@ -1061,8 +953,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
file_open=true;
//}
}
for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++)
{
//For a block composed of a single equation determines whether we have to evaluate or to solve the equation
if (ModelBlock->Block_List[j].Size==1)
{
@ -1076,7 +966,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
// The equations
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
//ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0);
//The Temporary terms
temporary_terms_type tt2;
#ifdef DEBUGC
@ -1116,6 +1005,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
}
switch (ModelBlock->Block_List[j].Simulation_Type)
{
evaluation:
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
@ -1133,19 +1023,13 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
lhs->compile(code_file, true, temporary_terms, map_idx);
}
break;
/*case EVALUATE_BACKWARD_R:
case EVALUATE_FORWARD_R:
lhs->compile(code_file, false, temporary_terms, map_idx);
rhs->compile(code_file, true, temporary_terms, map_idx);
break;*/
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
v=ModelBlock->Block_List[j].Equation[i];
Uf[v].eqr=i;
Uf[v].Ufl=NULL;
goto end;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
if (i<ModelBlock->Block_List[j].Nb_Recursives)
goto evaluation;
feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
v=ModelBlock->Block_List[j].Equation[i];
Uf[v].eqr=i;
Uf[v].Ufl=NULL;
@ -1167,6 +1051,7 @@ end:
}
code_file.write(&FENDEQU, sizeof(FENDEQU));
// The Jacobian if we have to solve the block
bool feedback_variable = (feedback_variables.size()>0);
if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD
&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD
/*&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R
@ -1183,30 +1068,73 @@ end:
break;
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
m=ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
if(feedback_variable)
{
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int v=ModelBlock->Block_List[j].Equation[eqr];
if (!Uf[v].Ufl)
int u = feedback_variables.size();
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
{
Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl_First=Uf[v].Ufl;
//Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
int eqr=it.first.first;
int eq=it.first.second;
int varr=it.second.first;
int var=it.second.second.first;
int v=ModelBlock->Block_List[j].Equation[eqr];
k=it.second.second.second;
if (!Uf[v].Ufl)
{
Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl_First=Uf[v].Ufl;
}
else
{
Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl=Uf[v].Ufl->pNext;
}
Uf[v].Ufl->pNext=NULL;
Uf[v].Ufl->u=u;
Uf[v].Ufl->var=var;
compileDerivative(code_file, eq, var, 0, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
u++;
/*output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k
<< ") " << var+1
<< ", equation=" << eq+1 << endl;*/
}
else
}
else
{
m=ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl=Uf[v].Ufl->pNext;
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int v=ModelBlock->Block_List[j].Equation[eqr];
if (!Uf[v].Ufl)
{
Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl_First=Uf[v].Ufl;
}
else
{
Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl=Uf[v].Ufl->pNext;
}
Uf[v].Ufl->pNext=NULL;
Uf[v].Ufl->u=u;
Uf[v].Ufl->var=var;
compileDerivative(code_file, eq, var, 0, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
}
Uf[v].Ufl->pNext=NULL;
Uf[v].Ufl->u=u;
Uf[v].Ufl->var=var;
compileDerivative(code_file, eq, var, 0, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
}
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
@ -1340,8 +1268,6 @@ end:
prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
}
j++;
}
}
code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
code_file.write(&FEND, sizeof(FEND));
@ -1993,7 +1919,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
hessian_output << "v2";
hessianHelper(hessian_output, k, 2, output_type);
hessian_output << "=v2";
hessian_output << "=v2";
hessianHelper(hessian_output, k-1, 2, output_type);
hessian_output << ";" << endl;
@ -2041,7 +1967,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
k += k2;
}
if (mode == eStandardMode)
{
DynamicOutput << "%" << endl
@ -2597,6 +2523,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
BlockLinear(block_triangular.ModelBlock);
if (!no_tmp_terms)
computeTemporaryTermsOrdered(block_triangular.ModelBlock);
computeChaineRuleJacobian(block_triangular.ModelBlock);
}
else
if (!no_tmp_terms)
@ -2817,6 +2745,107 @@ DynamicModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDExcepti
return it->second;
}
void
DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock)
{
//cout << "computeChaineRuleJacobian\n";
//clock_t t1 = clock();
map<int, NodeID> recursive_variables;
first_chaine_rule_derivatives.clear();
for(int blck = 0; blck<ModelBlock->Size; blck++)
{
//cout << "blck=" << blck << "\n";
recursive_variables.clear();
if (ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
{
//cout << "SOLVE_TWO_BOUNDARIES_COMPLETE \n";
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->clear();
for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++)
{
if (ModelBlock->Block_List[blck].Equation_Type[i] == E_EVALUATE_S)
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = ModelBlock->Block_List[blck].Equation_Normalized[i];
else
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]];
}
//cout << "After recursive_alloc\n";
map<pair<pair<int, int >, pair<pair<int, int>,int> > , int > Derivatives = block_triangular.get_Derivatives(ModelBlock, blck);
for(map<pair<pair<int, int >, pair<pair<int, int>,int> > , int >::const_iterator it = Derivatives.begin(); it != Derivatives.end(); it++)
{
int eqr = it->first.first.first;
int eq = it->first.first.second;
int varr = it->first.second.first.first;
int var = it->first.second.first.second;
int lag = it->first.second.second;
int Deriv_type = it->second;
if(Deriv_type == 0)
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))];
else if (Deriv_type == 1)
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
else if (Deriv_type == 2)
{
if(ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_List[blck].Nb_Recursives)
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
else
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
}
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag))));
}
/*for(int lag = -ModelBlock->Block_List[blck].Max_Lag; lag <= ModelBlock->Block_List[blck].Max_Lead; lag++)
{
for(int eq = 0; eq < ModelBlock->Block_List[blck].Size; eq++)
{
int eqr = ModelBlock->Block_List[blck].Equation[eq];
for(int var = ModelBlock->Block_List[blck].Nb_Recursives; var < ModelBlock->Block_List[blck].Size; var++)
{
int varr = ModelBlock->Block_List[blck].Variable[var];
NodeID d1 = equations[eqr]->getChaineRuleDerivative(recursive_variables, varr, lag);
if (d1 == Zero)
continue;
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = d1;
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag))));
}
}
}*/
}
else if( ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_SIMPLE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_SIMPLE
or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_COMPLETE)
{
//cout << "SOLVE_FORWARD_SIMPLE \n";
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->clear();
for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++)
{
if (ModelBlock->Block_List[blck].Equation_Type[i] == E_EVALUATE_S)
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = ModelBlock->Block_List[blck].Equation_Normalized[i];
else
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]];
}
for(int eq = ModelBlock->Block_List[blck].Nb_Recursives; eq < ModelBlock->Block_List[blck].Size; eq++)
{
int eqr = ModelBlock->Block_List[blck].Equation[eq];
for(int var = ModelBlock->Block_List[blck].Nb_Recursives; var < ModelBlock->Block_List[blck].Size; var++)
{
int varr = ModelBlock->Block_List[blck].Variable[var];
NodeID d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
if (d1 == Zero)
continue;
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1;
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, 0))));
}
}
}
}
//cout << "elapsed time in milliseconds = " << 1000.0*(double(clock()) - double(t1))/double(CLOCKS_PER_SEC) << "\n";
}
void
DynamicModel::computeParamsDerivatives()
{

2
DynamicModel.hh
View File

@ -116,6 +116,8 @@ private:
int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Compute the column indices of the dynamic Jacobian
void computeDynJacobianCols(bool jacobianExo);
//! Computes chaine rule derivatives of the Jacobian w.r. to endogenous variables
void computeChaineRuleJacobian(Model_Block *ModelBlock);
//! Computes derivatives of the Jacobian w.r. to parameters
void computeParamsDerivatives();
//! Computes temporary terms for the file containing parameters derivatives

18
MinimumFeedbackSet.cc
View File

@ -26,6 +26,8 @@ namespace MFS
void
Suppress(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G)
{
/*clear all in and out edges of vertex_to_eliminate
and remove vertex_to_eliminate from the graph*/
clear_vertex(vertex_to_eliminate, G);
remove_vertex(vertex_to_eliminate, G);
}
@ -41,6 +43,7 @@ namespace MFS
void
Eliminate(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G)
{
/*before the vertex i suppression replace all edges e_k_i and e_i_j by e_k_j*/
if (in_degree (vertex_to_eliminate, G) > 0 && out_degree (vertex_to_eliminate, G) > 0)
{
AdjacencyList_type::in_edge_iterator it_in, in_end;
@ -66,11 +69,14 @@ namespace MFS
color[u] = gray_color;
graph_traits<AdjacencyList_type>::out_edge_iterator vi, vi_end;
for (tie(vi, vi_end) = out_edges(u, g); vi != vi_end; ++vi)
if (color[target(*vi, g)] == white_color && has_cycle_dfs(g, target(*vi, g), color, circuit_stack))
if (color[target(*vi, g)] == white_color)
{
// cycle detected, return immediately
circuit_stack.push_back(v_index[target(*vi, g)]);
return true;
if (has_cycle_dfs(g, target(*vi, g), color, circuit_stack))
{
// cycle detected, return immediately
circuit_stack.push_back(v_index[target(*vi, g)]);
return true;
}
}
else if (color[target(*vi, g)] == gray_color)
{
@ -209,6 +215,8 @@ namespace MFS
vector_vertex_descriptor
Collect_Doublet(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G)
{
/*collect all doublet (for each edge e_i_k there is an edge e_k_i with k!=i) in the graph
and return the vector of doublet*/
AdjacencyList_type::in_edge_iterator it_in, in_end;
AdjacencyList_type::out_edge_iterator it_out, out_end;
vector<AdjacencyList_type::vertex_descriptor> Doublet;
@ -223,6 +231,7 @@ namespace MFS
bool
Vertex_Belong_to_a_Clique(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G)
{
/*Detect all the clique (all vertex in a clique are related to each other) in the graph*/
vector<AdjacencyList_type::vertex_descriptor> liste;
bool agree = true;
AdjacencyList_type::in_edge_iterator it_in, in_end;
@ -262,6 +271,7 @@ namespace MFS
bool
Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_type& G)
{
/*Graph reduction: eliminating purely intermediate variables or variables outside of any circuit*/
bool something_has_been_done = false;
bool not_a_loop;
int i;

14
ModelTree.cc
View File

@ -65,6 +65,20 @@ ModelTree::computeJacobian(const set<int> &vars)
}
}
void
ModelTree::writeChaineRuleDerivative(ostream &output, int eq, int var, int lag,
ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chaine_rule_derivatives.find(make_pair(eq, make_pair(var, lag)));
if (it != first_chaine_rule_derivatives.end())
(it->second)->writeOutput(output, output_type, temporary_terms);
else
output << 0;
}
void
ModelTree::computeHessian(const set<int> &vars)
{

5
ModelTree.hh
View File

@ -47,6 +47,9 @@ protected:
*/
first_derivatives_type first_derivatives;
typedef map< pair< int, pair< int, int> >, NodeID> first_chaine_rule_derivatives_type;
first_chaine_rule_derivatives_type first_chaine_rule_derivatives;
typedef map<pair<int, pair<int, int> >, NodeID> second_derivatives_type;
//! Second order derivatives
/*! First index is equation number, second and third are variables w.r. to which is computed the derivative.
@ -80,6 +83,8 @@ protected:
//! Write derivative of an equation w.r. to a variable
void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
//! Write chaine rule derivative of a recursive equation w.r. to a variable
void writeChaineRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
//! Computes temporary terms (for all equations and derivatives)
void computeTemporaryTerms(bool is_matlab);
//! Writes temporary terms