diff --git a/BlockTriangular.cc b/BlockTriangular.cc index 1810c3f6..54201d96 100644 --- a/BlockTriangular.cc +++ b/BlockTriangular.cc @@ -120,6 +120,9 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog typedef adjacency_list BipartiteGraph; + + if(verbose == 2) + cout << "trying to normlized even in singular case\n"; /* Vertices 0 to n-1 are for endogenous (using type specific ID) Vertices n to 2*n-1 are for equations (using equation no.) @@ -148,7 +151,7 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog { if (verbose == 1) { - cerr << "WARNING: Could not normalize dynamic model. Variable " + cout << "WARNING: Could not normalize dynamic model. Variable " << symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin())) << " is not in the maximum cardinality matching. Trying to find a singular normalization." << endl; //exit(EXIT_FAILURE); @@ -169,6 +172,7 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog memcpy(SIM, IM, equation_number*equation_number*sizeof(bool)); for (int i = 0; i < n; i++) { + //printf("match equation %4d with variable %4d \n", mate_map[i] - n, i); Index_Equ_IM[i + prologue] = Index_Equ_IM_tmp[mate_map[i] - n + prologue]; for (int k = 0; k < n; k++) IM[(i+prologue)*equation_number+k +prologue] = SIM[(mate_map[i]-n+prologue)*equation_number+k +prologue]; @@ -337,6 +341,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block bool *Cur_IM; bool *IM, OK; int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo, Lag_Other_Endo, Lead_Other_Endo; + //cout << "block " << count_Block << " size " << size << " SimType=" << BlockSim(SimType) << "\n"; ModelBlock->Periods = periods; ModelBlock->Block_List[count_Block].is_linear = true; @@ -725,35 +730,35 @@ BlockTriangular::Equation_Type_determination(vector &equations, lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); tmp_s << "y(it_, " << Index_Var_IM[i]+1 << ")"; map >, NodeID>::iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0))); - set > result; pair res; - derivative->second->collectEndogenous(result); - set >::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); - //cout << tt.str().c_str() << " tmp_output.str()=" << tmp_output.str() << " tmp_s.str()=" << tmp_s.str(); - if (tmp_output.str() == tmp_s.str() and tt.str()=="1") + if(derivative != first_order_endo_derivatives.end()) { - Equation_Simulation_Type = E_EVALUATE; - //cout << " E_EVALUATE "; - } - else - { - vector > > List_of_Op_RHS; - res = equations[eq]->normalizeEquation(var, List_of_Op_RHS); - if(mfs==2) + set > result; + derivative->second->collectEndogenous(result); + set >::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") { - if(d_endo_variable == result.end() && res.second) - Equation_Simulation_Type = E_EVALUATE_S; + Equation_Simulation_Type = E_EVALUATE; } - else if(mfs==3) + else { - if(res.second) // The equation could be solved analytically - Equation_Simulation_Type = E_EVALUATE_S; + vector > > 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; + } } } - //cout << " " << c_Equation_Type(Equation_Simulation_Type) << endl; V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, dynamic_cast(res.second)); } return (V_Equation_Simulation_Type); @@ -1000,10 +1005,10 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, Compute_Normalization(IM, n, prologue, epilogue, 2, IM_0, Index_Equ_IM); } } - 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); + free(SIM0); V_Equation_Type = Equation_Type_determination(equations, first_order_endo_derivatives, Index_Var_IM, Index_Equ_IM, mfs); diff --git a/DynamicModel.cc b/DynamicModel.cc index 1451f7a9..4727a1e5 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -2424,12 +2424,16 @@ DynamicModel::toStaticDll(StaticDllModel &static_model) const // Convert model local variables (need to be done first) for (map::const_iterator it = local_variables_table.begin(); it != local_variables_table.end(); it++) + { static_model.AddLocalVariable(symbol_table.getName(it->first), it->second->toStatic(static_model)); + } // Convert equations for (vector::const_iterator it = equations.begin(); it != equations.end(); it++) + { static_model.addEquation((*it)->toStatic(static_model)); + } } void diff --git a/ExprNode.cc b/ExprNode.cc index 87f4bc39..7226eac6 100644 --- a/ExprNode.cc +++ b/ExprNode.cc @@ -52,12 +52,10 @@ ExprNode::getDerivative(int deriv_id) { if (!preparedForDerivation) prepareForDerivation(); - // Return zero if derivative is necessarily null (using symbolic a priori) set::const_iterator it = non_null_derivatives.find(deriv_id); if (it == non_null_derivatives.end()) return datatree.Zero; - // If derivative is stored in cache, use the cached value, otherwise compute it (and cache it) map::const_iterator it2 = derivatives.find(deriv_id); if (it2 != derivatives.end()) @@ -308,7 +306,7 @@ VariableNode::prepareForDerivation() { if (preparedForDerivation) return; - + int der_id; preparedForDerivation = true; // Fill in non_null_derivatives @@ -319,12 +317,25 @@ VariableNode::prepareForDerivation() case eExogenousDet: case eParameter: // For a variable or a parameter, the only non-null derivative is with respect to itself - non_null_derivatives.insert(datatree.getDerivID(symb_id, lag)); + //non_null_derivatives.insert(datatree.getDerivID(symb_id, lag)); + der_id = datatree.getDerivID(symb_id, lag); + non_null_derivatives = set(&der_id, &der_id+1); break; case eModelLocalVariable: datatree.local_variables_table[symb_id]->prepareForDerivation(); // Non null derivatives are those of the value of the local parameter non_null_derivatives = datatree.local_variables_table[symb_id]->non_null_derivatives; + /*for(set::const_iterator it=datatree.local_variables_table[symb_id]->non_null_derivatives.begin(); it != datatree.local_variables_table[symb_id]->non_null_derivatives.end(); it++) + non_null_derivatives.insert(*it);*/ + //cout << "symb_id=" << symb_id << " " << datatree.symbol_table.getName(symb_id) << " " << datatree.symbol_table.maxID(); + //non_null_derivatives = datatree.local_variables_table[symb_id]->non_null_derivatives; + /*cout << " non_null_derivatives.size()= " << non_null_derivatives.size() << endl; + if(non_null_derivatives.size()>0) + for(set::const_iterator it=non_null_derivatives.begin(); it != non_null_derivatives.end(); it++) + { + cout << "symb_id=" << symb_id << " *it = " << *it << "\n" ; + //cout << datatree.symbol_table.getName(staticdllmodel.getSymbIDByDerivID(*it)) << "(" << staticdllmodel.getLagByDerivID(*it) << ")\n"; + }*/ break; case eModFileLocalVariable: // Such a variable is never derived @@ -349,6 +360,7 @@ VariableNode::computeDerivative(int deriv_id) else return datatree.Zero; case eModelLocalVariable: + //datatree.local_variables_table[symb_id]->writeOutput(cout); return datatree.local_variables_table[symb_id]->getDerivative(deriv_id); case eModFileLocalVariable: cerr << "ModFileLocalVariable is not derivable" << endl; @@ -628,7 +640,7 @@ VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_ break; case eModelLocalVariable: case eModFileLocalVariable: - //cout << "eModelLocalVariable=" << symb_id << "\n"; + //cout << "eModelLocalVariable=" << tsid << "lhs_rhs=" << lhs_rhs << "\n"; datatree.local_variables_table[symb_id]->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic); break; case eUnknownFunction: diff --git a/ExprNode.hh b/ExprNode.hh index 08bafae8..903e5a87 100644 --- a/ExprNode.hh +++ b/ExprNode.hh @@ -241,7 +241,7 @@ public: typedef map subst_table_t; //! Creates auxiliary 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 diff --git a/ModelTree.cc b/ModelTree.cc index ac19cc4c..eecfc3d7 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -53,13 +53,15 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag, void ModelTree::computeJacobian(const set &vars) { + int i = 0; for(set::const_iterator it = vars.begin(); - it != vars.end(); it++) + it != vars.end(); it++, i++) for (int eq = 0; eq < (int) equations.size(); eq++) { NodeID d1 = equations[eq]->getDerivative(*it); if (d1 == Zero) continue; + //printf("eq=%4d var=%4d i=%4d\n",eq,*it,i); first_derivatives[make_pair(eq, *it)] = d1; ++NNZDerivatives[0]; } diff --git a/StaticDllModel.cc b/StaticDllModel.cc index b135122c..a2395f45 100644 --- a/StaticDllModel.cc +++ b/StaticDllModel.cc @@ -47,7 +47,8 @@ void StaticDllModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const { //first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag))); - first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag))); + //first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag))); + first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, symbol_table.getID(eEndogenous, symb_id))); if (it != first_derivatives.end()) (it->second)->compile(code_file, false, temporary_terms, map_idx, false, false); else @@ -645,8 +646,8 @@ StaticDllModel::evaluateJacobian(const eval_context_type &eval_context, jacob_ma it != first_derivatives.end(); it++) { //cout << "it->first.second=" << it->first.second << " variable_table.getSymbolID(it->first.second)=" << variable_table.getSymbolID(it->first.second) << " Type=" << variable_table.getType(it->first.second) << " eEndogenous=" << eEndogenous << " eExogenous=" << eExogenous << " variable_table.getLag(it->first.second)=" << variable_table.getLag(it->first.second) << "\n"; - if (getTypeByDerivID(it->first.second) == eEndogenous) - { + /*if (getTypeByDerivID(it->first.second) == eEndogenous) + {*/ NodeID Id = it->second; double val = 0; try @@ -655,14 +656,17 @@ StaticDllModel::evaluateJacobian(const eval_context_type &eval_context, jacob_ma } catch (ExprNode::EvalException &e) { - cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ") [" << getSymbIDByDerivID(it->first.second) << "] !" << endl; + //cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ") [" << getSymbIDByDerivID(it->first.second) << "] !" << endl; + cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(it->first.second) << "(" << 0 << ") [" << it->first.second << "] !" << endl; Id->writeOutput(cout, oMatlabStaticModelSparse, temporary_terms); cout << "\n"; - cerr << "StaticDllModel::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ")!" << endl; + //cerr << "StaticDllModel::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ")!" << endl; + cerr << "StaticDllModel::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(it->first.second) << "(" << 0 << ")!" << endl; } int eq=it->first.first; - int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it->first.second));///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second); - int k1 = getLagByDerivID(it->first.second); + //int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it->first.second));///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second); + int var = symbol_table.getTypeSpecificID(it->first.second);///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second); + int k1 = 0;//getLagByDerivID(it->first.second); if (a_variable_lag!=k1) { IM=block_triangular.incidencematrix.Get_IM(k1, eEndogenous); @@ -673,22 +677,26 @@ StaticDllModel::evaluateJacobian(const eval_context_type &eval_context, jacob_ma j++; (*j_m)[make_pair(eq,var)]+=val; } - if (IM[eq*symbol_table.endo_nbr()+var] && (fabs(val) < cutoff)) + /*if (IM[eq*symbol_table.endo_nbr()+var] && (fabs(val) < cutoff)) { - if (block_triangular.bt_verbose) + //if (block_triangular.bt_verbose) cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")\n"; block_triangular.incidencematrix.unfill_IM(eq, var, k1, eEndogenous); i++; - } - } + }*/ + /*}*/ } //Get ride of the elements of the incidence matrix equal to Zero IM=block_triangular.incidencematrix.Get_IM(0, eEndogenous); - for (int i=0;i0) { cout << i << " elements among " << first_derivatives.size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded\n"; @@ -711,7 +719,8 @@ StaticDllModel::BlockLinear(Model_Block *ModelBlock) int eq=ModelBlock->Block_List[j].IM_lead_lag[ll].Equ_Index[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[ll].Var_Index[i]; //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,0))); - first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var),0))); + //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var),0))); + first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,symbol_table.getID(eEndogenous, var))); if (it!= first_derivatives.end()) { NodeID Id = it->second; @@ -741,7 +750,8 @@ StaticDllModel::BlockLinear(Model_Block *ModelBlock) 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]; //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,k1))); - first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var),k1))); + //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var),k1))); + first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,symbol_table.getID(eEndogenous, var))); NodeID Id = it->second; if (it!= first_derivatives.end()) { @@ -775,14 +785,17 @@ StaticDllModel::collect_first_order_derivatives_endogenous() for (first_derivatives_type::iterator it2 = first_derivatives.begin(); it2 != first_derivatives.end(); it2++) { - if (getTypeByDerivID(it2->first.second)==eEndogenous) - { + //if (getTypeByDerivID(it2->first.second)==eEndogenous) + /*if (it2->first.second)==eEndogenous) + {*/ int eq = it2->first.first; - int var=symbol_table.getTypeSpecificID(getSymbIDByDerivID(it2->first.second)); - int lag=getLagByDerivID(it2->first.second); + //int var=symbol_table.getTypeSpecificID(getSymbIDByDerivID(it2->first.second)); + int var=symbol_table.getTypeSpecificID(it2->first.second); + //int lag=getLagByDerivID(it2->first.second); + int lag = 0; //if (lag==0) endo_derivatives[make_pair(eq, make_pair(var, lag))] = it2->second; - } + //} } return endo_derivatives; } @@ -796,26 +809,27 @@ StaticDllModel::computingPass(const eval_context_type &eval_context, bool no_tmp // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic set vars; - for (deriv_id_table_t::const_iterator it = deriv_id_table.begin(); + /*for (deriv_id_table_t::const_iterator it = deriv_id_table.begin(); it != deriv_id_table.end(); it++) { SymbolType type = symbol_table.getType(it->first.first); if (type == eEndogenous) vars.insert(it->second); - } + }*/ + for(int i = 0; i < symbol_table.endo_nbr(); i++) + vars.insert(symbol_table.getID(eEndogenous, i)); // Launch computations cout << "Computing static model derivatives:" << endl << " - order 1" << endl; computeJacobian(vars); //cout << "mode=" << mode << " eSparseDLLMode=" << eSparseDLLMode << " eSparseMode=" << eSparseMode << "\n"; - BuildIncidenceMatrix(); + jacob_map j_m; evaluateJacobian(eval_context, &j_m, true); - if (block_triangular.bt_verbose) { cout << "The gross incidence matrix \n"; @@ -914,9 +928,11 @@ StaticDllModel::computeChainRuleJacobian(Model_Block *ModelBlock) 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]; + //recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = ModelBlock->Block_List[blck].Equation_Normalized[i]; + recursive_variables[symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i])] = 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]]; + //recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]]; + recursive_variables[symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i])] = equations[ModelBlock->Block_List[blck].Equation[i]]; } map >, pair >, int> Derivatives = block_triangular.get_Derivatives(ModelBlock, blck); @@ -934,18 +950,20 @@ StaticDllModel::computeChainRuleJacobian(Model_Block *ModelBlock) int varr = it_l.second.second; if(Deriv_type == 0) { - first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))]; + //first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))]; + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, symbol_table.getID(eEndogenous, varr))]; } else if (Deriv_type == 1) { - first_chain_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); + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(symbol_table.getID(eEndogenous, varr), recursive_variables); } else if (Deriv_type == 2) { if(ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S && eqBlock_List[blck].Nb_Recursives) - first_chain_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); + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(symbol_table.getID(eEndogenous, varr), recursive_variables); else - first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables); + //first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables); + first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(symbol_table.getID(eEndogenous, varr), recursive_variables); } ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr))); } @@ -957,9 +975,11 @@ StaticDllModel::computeChainRuleJacobian(Model_Block *ModelBlock) 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]; + //recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = ModelBlock->Block_List[blck].Equation_Normalized[i]; + recursive_variables[symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i])] = 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]]; + //recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]]; + recursive_variables[symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i])] = equations[ModelBlock->Block_List[blck].Equation[i]]; } for(int eq = ModelBlock->Block_List[blck].Nb_Recursives; eq < ModelBlock->Block_List[blck].Size; eq++) { @@ -967,7 +987,8 @@ StaticDllModel::computeChainRuleJacobian(Model_Block *ModelBlock) 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); + //NodeID d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables); + NodeID d1 = equations[eqr]->getChainRuleDerivative(symbol_table.getID(eEndogenous, varr), recursive_variables); if (d1 == Zero) continue; first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1; diff --git a/StaticModel.cc b/StaticModel.cc index d848b739..15981f62 100644 --- a/StaticModel.cc +++ b/StaticModel.cc @@ -250,6 +250,8 @@ StaticModel::computeNormalization() #if 1 bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]); + for (int i = 0; i < n; i++) + printf("match equation %4d with variable %4d \n", mate_map[i] - n, i); #else // Alternative way to compute normalization, by giving an initial matching using natural normalizations fill(mate_map.begin(), mate_map.end(), graph_traits::null_vertex()); @@ -404,11 +406,11 @@ StaticModel::computeSortedBlockDecomposition() for(int i = 0; i < n; i++) blocks[unordered2ordered[endo2block[i]]].insert(i); -#ifdef DEBUG +//#ifdef DEBUG cout << "Found " << m << " blocks" << endl; for(int i = 0; i < m; i++) cout << " Block " << i << " of size " << blocks[i].size() << endl; -#endif +//#endif } void