diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 2954c023810dc6d62cd36675371054ef81ebe97f..ce800ac8da04b8dc6cf7cf5b78e53f95d7ef0318 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -907,8 +907,7 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct)
void
PlannerObjectiveStatement::computingPass()
{
- model_tree->computeStaticHessian = true;
- model_tree->computingPass(false);
+ model_tree->computingPass(true, false);
}
void
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index a6d334dcfd0785262aeb2805dfcf4e8382f6aae8..a04e0ef106ffe58efd2ae31f80d30e236ccf3820 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -34,17 +34,13 @@
DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg) :
ModelTree(symbol_table_arg, num_constants_arg),
- var_endo_nbr(0),
max_lag(0), max_lead(0),
max_endo_lag(0), max_endo_lead(0),
max_exo_lag(0), max_exo_lead(0),
max_exo_det_lag(0), max_exo_det_lead(0),
+ dynJacobianColsNbr(0),
cutoff(1e-15),
markowitz(0.7),
- computeJacobian(false),
- computeJacobianExo(false),
- computeHessian(false),
- computeThirdDerivatives(false),
block_triangular(symbol_table_arg)
{
}
@@ -94,7 +90,7 @@ DynamicModel::BuildIncidenceMatrix()
}
void
-DynamicModel::computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock)
+DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
{
map<NodeID, pair<int, int> > first_occurence;
map<NodeID, int> reference_count;
@@ -1157,7 +1153,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
<< " {" << endl
<< " /* Set the output pointer to the output matrix g1. */" << endl
- << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << getDynJacobianColsNbr() << ", mxREAL);" << endl
+ << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr << ", mxREAL);" << endl
<< " /* Create a C pointer to a copy of the output matrix g1. */" << endl
<< " g1 = mxGetPr(plhs[1]);" << endl
<< " }" << endl
@@ -1166,7 +1162,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
<< " if (nlhs >= 3)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix g2. */" << endl
- << " plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << getDynJacobianColsNbr()*getDynJacobianColsNbr()
+ << " plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr*dynJacobianColsNbr
<< ", mxREAL);" << endl
<< " /* Create a C pointer to a copy of the output matrix g1. */" << endl
<< " g2 = mxGetPr(plhs[2]);" << endl
@@ -1282,62 +1278,77 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
int i, k, Nb_SGE=0;
bool skip_head, open_par=false;
- if (computeJacobian || computeJacobianExo || computeHessian)
+
+ mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n";
+ mDynamicModelFile << " global oo_ options_ M_ ;\n";
+ mDynamicModelFile << " g2=[];g3=[];\n";
+ //Temporary variables declaration
+ OK=true;
+ ostringstream tmp_output;
+ for (temporary_terms_type::const_iterator it = temporary_terms.begin();
+ it != temporary_terms.end(); it++)
+ {
+ if (OK)
+ OK=false;
+ else
+ tmp_output << " ";
+ (*it)->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms);
+ }
+ if (tmp_output.str().length()>0)
+ mDynamicModelFile << " global " << tmp_output.str() << " M_ ;\n";
+
+ mDynamicModelFile << " T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n";
+ tmp_output.str("");
+ for (temporary_terms_type::const_iterator it = temporary_terms.begin();
+ it != temporary_terms.end(); it++)
+ {
+ tmp_output << " ";
+ (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
+ tmp_output << "=T_init;\n";
+ }
+ if (tmp_output.str().length()>0)
+ mDynamicModelFile << tmp_output.str();
+
+ mDynamicModelFile << " y_kmin=M_.maximum_lag;\n";
+ mDynamicModelFile << " y_kmax=M_.maximum_lead;\n";
+ mDynamicModelFile << " y_size=M_.endo_nbr;\n";
+ mDynamicModelFile << " if(length(varargin)>0)\n";
+ mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n";
+ mDynamicModelFile << " params=varargin{3};\n";
+ mDynamicModelFile << " it_=varargin{4};\n";
+ /*i = symbol_table.endo_nbr*(variable_table.max_endo_lag+variable_table.max_endo_lead+1)+
+ symbol_table.exo_nbr*(variable_table.max_exo_lag+variable_table.max_exo_lead+1);
+ mDynamicModelFile << " g1=spalloc(" << symbol_table.endo_nbr << ", " << i << ", " << i*symbol_table.endo_nbr << ");\n";*/
+ mDynamicModelFile << " Per_u_=0;\n";
+ mDynamicModelFile << " Per_y_=it_*y_size;\n";
+ mDynamicModelFile << " y=varargin{1};\n";
+ mDynamicModelFile << " ys=y(it_,:);\n";
+ mDynamicModelFile << " x=varargin{2};\n";
+ prev_Simulation_Type=-1;
+ tmp.str("");
+ tmp_eq.str("");
+ for (int count_call=1, i = 0;i < block_triangular.ModelBlock->Size;i++, count_call++)
{
- mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n";
- mDynamicModelFile << " global oo_ options_ M_ ;\n";
- mDynamicModelFile << " g2=[];g3=[];\n";
- //Temporary variables declaration
- OK=true;
- ostringstream tmp_output;
- for (temporary_terms_type::const_iterator it = temporary_terms.begin();
- it != temporary_terms.end(); it++)
+ k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
+ if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) &&
+ ((prev_Simulation_Type==EVALUATE_FORWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FORWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R)
+ || (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)))
{
- if (OK)
- OK=false;
- else
- tmp_output << " ";
- (*it)->writeOutput(tmp_output, oMatlabStaticModelSparse, temporary_terms);
+ mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
+ tmp_eq.str("");
+ mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
+ tmp.str("");
+ mDynamicModelFile << tmp1.str();
+ tmp1.str("");
}
- if (tmp_output.str().length()>0)
- mDynamicModelFile << " global " << tmp_output.str() << " M_ ;\n";
-
- mDynamicModelFile << " T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n";
- tmp_output.str("");
- for (temporary_terms_type::const_iterator it = temporary_terms.begin();
- it != temporary_terms.end(); it++)
+ for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
- tmp_output << " ";
- (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms);
- tmp_output << "=T_init;\n";
+ tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
}
- if (tmp_output.str().length()>0)
- mDynamicModelFile << tmp_output.str();
-
- mDynamicModelFile << " y_kmin=M_.maximum_lag;\n";
- mDynamicModelFile << " y_kmax=M_.maximum_lead;\n";
- mDynamicModelFile << " y_size=M_.endo_nbr;\n";
- mDynamicModelFile << " if(length(varargin)>0)\n";
- mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n";
- mDynamicModelFile << " params=varargin{3};\n";
- mDynamicModelFile << " it_=varargin{4};\n";
- /*i = symbol_table.endo_nbr*(variable_table.max_endo_lag+variable_table.max_endo_lead+1)+
- symbol_table.exo_nbr*(variable_table.max_exo_lag+variable_table.max_exo_lead+1);
- mDynamicModelFile << " g1=spalloc(" << symbol_table.endo_nbr << ", " << i << ", " << i*symbol_table.endo_nbr << ");\n";*/
- mDynamicModelFile << " Per_u_=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " y=varargin{1};\n";
- mDynamicModelFile << " ys=y(it_,:);\n";
- mDynamicModelFile << " x=varargin{2};\n";
- prev_Simulation_Type=-1;
- tmp.str("");
- tmp_eq.str("");
- for (int count_call=1, i = 0;i < block_triangular.ModelBlock->Size;i++, count_call++)
+ if (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)
{
- k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) &&
- ((prev_Simulation_Type==EVALUATE_FORWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FORWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R)
- || (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)))
+ if (i==block_triangular.ModelBlock->Size-1)
{
mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
tmp_eq.str("");
@@ -1346,116 +1357,100 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << tmp1.str();
tmp1.str("");
}
- for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
- {
- tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
- tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
- }
- if (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)
- {
- if (i==block_triangular.ModelBlock->Size-1)
- {
- mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
- tmp_eq.str("");
- mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
- tmp.str("");
- mDynamicModelFile << tmp1.str();
- tmp1.str("");
- }
- }
- if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
- (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))
- skip_head=true;
- else
- skip_head=false;
- switch (k)
- {
- case EVALUATE_FORWARD:
- case EVALUATE_BACKWARD:
- case EVALUATE_FORWARD_R:
- case EVALUATE_BACKWARD_R:
- if (!skip_head)
- {
- tmp1 << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 1, it_-1, 1);\n";
- tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n";
- }
- break;
- case SOLVE_FORWARD_SIMPLE:
- case SOLVE_BACKWARD_SIMPLE:
- mDynamicModelFile << " y_index_eq = " << block_triangular.ModelBlock->Block_List[i].Equation[0]+1 << ";\n";
- mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n";
- mDynamicModelFile << " residual(y_index_eq)=r;\n";
- tmp_eq.str("");
- tmp.str("");
- break;
- case SOLVE_FORWARD_COMPLETE:
- case SOLVE_BACKWARD_COMPLETE:
- mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
- mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n";
- mDynamicModelFile << " residual(y_index_eq)=r;\n";
- break;
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- case SOLVE_TWO_BOUNDARIES_SIMPLE:
- int j;
- mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
- tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
- mDynamicModelFile << " y_index = [";
- for (j=0;j<tmp_i;j++)
- for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
- {
- mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1+j*symbol_table.endo_nbr();
- }
- int tmp_ix=block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1;
- for (j=0;j<tmp_ix;j++)
- for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].nb_exo;ik++)
- mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Exogenous[ik]+1+j*symbol_table.exo_nbr()+symbol_table.endo_nbr()*tmp_i;
- mDynamicModelFile << " ];\n";
- tmp.str("");
- tmp_eq.str("");
- //mDynamicModelFile << " ga = [];\n";
- j = block_triangular.ModelBlock->Block_List[i].Size*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1)
- + block_triangular.ModelBlock->Block_List[i].nb_exo*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1);
- /*mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " <<
- block_triangular.ModelBlock->Block_List[i].Size*j << ");\n";*/
- tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
- mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_-" << max_lag << ", 1, " << max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n";
- /*if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead)
- mDynamicModelFile << " g1(y_index_eq,y_index) = ga;\n";
- else
- mDynamicModelFile << " g1(y_index_eq,y_index) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";*/
- mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
- break;
- }
- prev_Simulation_Type=k;
}
- if (tmp1.str().length())
+ if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
+ (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))
+ skip_head=true;
+ else
+ skip_head=false;
+ switch (k)
{
- mDynamicModelFile << tmp1.str();
- tmp1.str("");
+ case EVALUATE_FORWARD:
+ case EVALUATE_BACKWARD:
+ case EVALUATE_FORWARD_R:
+ case EVALUATE_BACKWARD_R:
+ if (!skip_head)
+ {
+ tmp1 << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 1, it_-1, 1);\n";
+ tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n";
+ }
+ break;
+ case SOLVE_FORWARD_SIMPLE:
+ case SOLVE_BACKWARD_SIMPLE:
+ mDynamicModelFile << " y_index_eq = " << block_triangular.ModelBlock->Block_List[i].Equation[0]+1 << ";\n";
+ mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n";
+ mDynamicModelFile << " residual(y_index_eq)=r;\n";
+ tmp_eq.str("");
+ tmp.str("");
+ break;
+ case SOLVE_FORWARD_COMPLETE:
+ case SOLVE_BACKWARD_COMPLETE:
+ mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
+ mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n";
+ mDynamicModelFile << " residual(y_index_eq)=r;\n";
+ break;
+ case SOLVE_TWO_BOUNDARIES_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_SIMPLE:
+ int j;
+ mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
+ tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
+ mDynamicModelFile << " y_index = [";
+ for (j=0;j<tmp_i;j++)
+ for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
+ {
+ mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1+j*symbol_table.endo_nbr();
+ }
+ int tmp_ix=block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1;
+ for (j=0;j<tmp_ix;j++)
+ for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].nb_exo;ik++)
+ mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Exogenous[ik]+1+j*symbol_table.exo_nbr()+symbol_table.endo_nbr()*tmp_i;
+ mDynamicModelFile << " ];\n";
+ tmp.str("");
+ tmp_eq.str("");
+ //mDynamicModelFile << " ga = [];\n";
+ j = block_triangular.ModelBlock->Block_List[i].Size*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1)
+ + block_triangular.ModelBlock->Block_List[i].nb_exo*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1);
+ /*mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " <<
+ block_triangular.ModelBlock->Block_List[i].Size*j << ");\n";*/
+ tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
+ mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_-" << max_lag << ", 1, " << max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n";
+ /*if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead)
+ mDynamicModelFile << " g1(y_index_eq,y_index) = ga;\n";
+ else
+ mDynamicModelFile << " g1(y_index_eq,y_index) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";*/
+ mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
+ break;
}
- mDynamicModelFile << " varargout{1}=residual;\n";
- mDynamicModelFile << " varargout{2}=dr;\n";
- mDynamicModelFile << " return;\n";
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n";
- mDynamicModelFile << " if(options_.simulation_method==0)\n";
- mDynamicModelFile << " mthd='Sparse LU';\n";
- mDynamicModelFile << " elseif(options_.simulation_method==2)\n";
- mDynamicModelFile << " mthd='GMRES';\n";
- mDynamicModelFile << " elseif(options_.simulation_method==3)\n";
- mDynamicModelFile << " mthd='BICGSTAB';\n";
- mDynamicModelFile << " else\n";
- mDynamicModelFile << " mthd='UNKNOWN';\n";
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " disp (['-----------------------------------------------------']) ;\n";
- mDynamicModelFile << " disp (['MODEL SIMULATION: (method=' mthd ')']) ;\n";
- mDynamicModelFile << " fprintf('\\n') ;\n";
- mDynamicModelFile << " periods=options_.periods;\n";
- mDynamicModelFile << " maxit_=options_.maxit_;\n";
- mDynamicModelFile << " solve_tolf=options_.solve_tolf;\n";
- mDynamicModelFile << " y=oo_.endo_simul';\n";
- mDynamicModelFile << " x=oo_.exo_simul;\n";
+ prev_Simulation_Type=k;
+ }
+ if (tmp1.str().length())
+ {
+ mDynamicModelFile << tmp1.str();
+ tmp1.str("");
}
+ mDynamicModelFile << " varargout{1}=residual;\n";
+ mDynamicModelFile << " varargout{2}=dr;\n";
+ mDynamicModelFile << " return;\n";
+ mDynamicModelFile << " end;\n";
+ mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n";
+ mDynamicModelFile << " if(options_.simulation_method==0)\n";
+ mDynamicModelFile << " mthd='Sparse LU';\n";
+ mDynamicModelFile << " elseif(options_.simulation_method==2)\n";
+ mDynamicModelFile << " mthd='GMRES';\n";
+ mDynamicModelFile << " elseif(options_.simulation_method==3)\n";
+ mDynamicModelFile << " mthd='BICGSTAB';\n";
+ mDynamicModelFile << " else\n";
+ mDynamicModelFile << " mthd='UNKNOWN';\n";
+ mDynamicModelFile << " end;\n";
+ mDynamicModelFile << " disp (['-----------------------------------------------------']) ;\n";
+ mDynamicModelFile << " disp (['MODEL SIMULATION: (method=' mthd ')']) ;\n";
+ mDynamicModelFile << " fprintf('\\n') ;\n";
+ mDynamicModelFile << " periods=options_.periods;\n";
+ mDynamicModelFile << " maxit_=options_.maxit_;\n";
+ mDynamicModelFile << " solve_tolf=options_.solve_tolf;\n";
+ mDynamicModelFile << " y=oo_.endo_simul';\n";
+ mDynamicModelFile << " x=oo_.exo_simul;\n";
+
prev_Simulation_Type=-1;
mDynamicModelFile << " params=M_.params;\n";
mDynamicModelFile << " oo_.deterministic_simulation.status = 0;\n";
@@ -1629,105 +1624,98 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
writeModelEquations(model_output, output_type);
int nrows = equations.size();
- int nvars = getDynJacobianColsNbr();
- int nvars_sq = nvars * nvars;
+ int hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr;
// Writing Jacobian
- if (computeJacobian || computeJacobianExo)
- for (first_derivatives_type::const_iterator it = first_derivatives.begin();
- it != first_derivatives.end(); it++)
- {
- int eq = it->first.first;
- int var = it->first.second;
- NodeID d1 = it->second;
-
- if (computeJacobianExo || getTypeByDerivID(var) == eEndogenous)
- {
- ostringstream g1;
- g1 << " g1";
- matrixHelper(g1, eq, getDynJacobianCol(var), output_type);
+ for (first_derivatives_type::const_iterator it = first_derivatives.begin();
+ it != first_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int var = it->first.second;
+ NodeID d1 = it->second;
- jacobian_output << g1.str() << "=" << g1.str() << "+";
- d1->writeOutput(jacobian_output, output_type, temporary_terms);
- jacobian_output << ";" << endl;
- }
- }
+ ostringstream g1;
+ g1 << " g1";
+ matrixHelper(g1, eq, getDynJacobianCol(var), output_type);
+
+ jacobian_output << g1.str() << "=" << g1.str() << "+";
+ d1->writeOutput(jacobian_output, output_type, temporary_terms);
+ jacobian_output << ";" << endl;
+ }
// Writing Hessian
- if (computeHessian)
- for (second_derivatives_type::const_iterator it = second_derivatives.begin();
- it != second_derivatives.end(); it++)
- {
- int eq = it->first.first;
- int var1 = it->first.second.first;
- int var2 = it->first.second.second;
- NodeID d2 = it->second;
-
- int id1 = getDynJacobianCol(var1);
- int id2 = getDynJacobianCol(var2);
-
- int col_nb = id1*nvars+id2;
- int col_nb_sym = id2*nvars+id1;
-
- hessian_output << " g2";
- matrixHelper(hessian_output, eq, col_nb, output_type);
- hessian_output << " = ";
- d2->writeOutput(hessian_output, output_type, temporary_terms);
- hessian_output << ";" << endl;
-
- // Treating symetric elements
- if (id1 != id2)
- {
- lsymetric << " g2";
- matrixHelper(lsymetric, eq, col_nb_sym, output_type);
- lsymetric << " = " << "g2";
- matrixHelper(lsymetric, eq, col_nb, output_type);
- lsymetric << ";" << endl;
- }
- }
+ for (second_derivatives_type::const_iterator it = second_derivatives.begin();
+ it != second_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int var1 = it->first.second.first;
+ int var2 = it->first.second.second;
+ NodeID d2 = it->second;
+
+ int id1 = getDynJacobianCol(var1);
+ int id2 = getDynJacobianCol(var2);
+
+ int col_nb = id1 * dynJacobianColsNbr + id2;
+ int col_nb_sym = id2 * dynJacobianColsNbr + id1;
+
+ hessian_output << " g2";
+ matrixHelper(hessian_output, eq, col_nb, output_type);
+ hessian_output << " = ";
+ d2->writeOutput(hessian_output, output_type, temporary_terms);
+ hessian_output << ";" << endl;
+
+ // Treating symetric elements
+ if (id1 != id2)
+ {
+ lsymetric << " g2";
+ matrixHelper(lsymetric, eq, col_nb_sym, output_type);
+ lsymetric << " = " << "g2";
+ matrixHelper(lsymetric, eq, col_nb, output_type);
+ lsymetric << ";" << endl;
+ }
+ }
// Writing third derivatives
- if (computeThirdDerivatives)
- for (third_derivatives_type::const_iterator it = third_derivatives.begin();
- it != third_derivatives.end(); it++)
- {
- int eq = it->first.first;
- int var1 = it->first.second.first;
- int var2 = it->first.second.second.first;
- int var3 = it->first.second.second.second;
- NodeID d3 = it->second;
-
- int id1 = getDynJacobianCol(var1);
- int id2 = getDynJacobianCol(var2);
- int id3 = getDynJacobianCol(var3);
-
- // Reference column number for the g3 matrix
- int ref_col = id1 * nvars_sq + id2 * nvars + id3;
-
- third_derivatives_output << " g3";
- matrixHelper(third_derivatives_output, eq, ref_col, output_type);
- third_derivatives_output << " = ";
- d3->writeOutput(third_derivatives_output, output_type, temporary_terms);
- third_derivatives_output << ";" << endl;
-
- // Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal)
- set<int> cols;
- cols.insert(id1 * nvars_sq + id3 * nvars + id2);
- cols.insert(id2 * nvars_sq + id1 * nvars + id3);
- cols.insert(id2 * nvars_sq + id3 * nvars + id1);
- cols.insert(id3 * nvars_sq + id1 * nvars + id2);
- cols.insert(id3 * nvars_sq + id2 * nvars + id1);
-
- for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
- if (*it2 != ref_col)
- {
- third_derivatives_output << " g3";
- matrixHelper(third_derivatives_output, eq, *it2, output_type);
- third_derivatives_output << " = " << "g3";
- matrixHelper(third_derivatives_output, eq, ref_col, output_type);
- third_derivatives_output << ";" << endl;
- }
- }
+ for (third_derivatives_type::const_iterator it = third_derivatives.begin();
+ it != third_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int var1 = it->first.second.first;
+ int var2 = it->first.second.second.first;
+ int var3 = it->first.second.second.second;
+ NodeID d3 = it->second;
+
+ int id1 = getDynJacobianCol(var1);
+ int id2 = getDynJacobianCol(var2);
+ int id3 = getDynJacobianCol(var3);
+
+ // Reference column number for the g3 matrix
+ int ref_col = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
+
+ third_derivatives_output << " g3";
+ matrixHelper(third_derivatives_output, eq, ref_col, output_type);
+ third_derivatives_output << " = ";
+ d3->writeOutput(third_derivatives_output, output_type, temporary_terms);
+ third_derivatives_output << ";" << endl;
+
+ // Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal)
+ set<int> cols;
+ cols.insert(id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2);
+ cols.insert(id2 * hessianColsNbr + id1 * dynJacobianColsNbr + id3);
+ cols.insert(id2 * hessianColsNbr + id3 * dynJacobianColsNbr + id1);
+ cols.insert(id3 * hessianColsNbr + id1 * dynJacobianColsNbr + id2);
+ cols.insert(id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1);
+
+ for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
+ if (*it2 != ref_col)
+ {
+ third_derivatives_output << " g3";
+ matrixHelper(third_derivatives_output, eq, *it2, output_type);
+ third_derivatives_output << " = " << "g3";
+ matrixHelper(third_derivatives_output, eq, ref_col, output_type);
+ third_derivatives_output << ";" << endl;
+ }
+ }
if (mode == eStandardMode)
{
@@ -1736,27 +1724,23 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
<< "%" << endl
<< endl
<< "residual = zeros(" << nrows << ", 1);" << endl
- << model_output.str();
-
- if (computeJacobian || computeJacobianExo)
- {
- // Writing initialization instruction for matrix g1
- DynamicOutput << "if nargout >= 2," << endl
- << " g1 = zeros(" << nrows << ", " << nvars << ");" << endl
- << endl
- << "%" << endl
- << "% Jacobian matrix" << endl
- << "%" << endl
- << endl
- << jacobian_output.str()
- << "end" << endl;
- }
- if (computeHessian)
+ << model_output.str()
+ // Writing initialization instruction for matrix g1
+ << "if nargout >= 2," << endl
+ << " g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");" << endl
+ << endl
+ << "%" << endl
+ << "% Jacobian matrix" << endl
+ << "%" << endl
+ << endl
+ << jacobian_output.str()
+ << "end" << endl;
+
+ if (second_derivatives.size())
{
// Writing initialization instruction for matrix g2
- int ncols = nvars_sq;
DynamicOutput << "if nargout >= 3," << endl
- << " g2 = sparse([],[],[], " << nrows << ", " << ncols << ", " << 5*ncols << ");" << endl
+ << " g2 = sparse([],[],[], " << nrows << ", " << hessianColsNbr << ", " << 5*hessianColsNbr << ");" << endl
<< endl
<< "%" << endl
<< "% Hessian matrix" << endl
@@ -1766,9 +1750,9 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
<< lsymetric.str()
<< "end;" << endl;
}
- if (computeThirdDerivatives)
+ if (third_derivatives.size())
{
- int ncols = nvars_sq * nvars;
+ int ncols = hessianColsNbr * dynJacobianColsNbr;
DynamicOutput << "if nargout >= 4," << endl
<< " g3 = sparse([],[],[], " << nrows << ", " << ncols << ", " << 5*ncols << ");" << endl
<< endl
@@ -1787,19 +1771,16 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
<< " double lhs, rhs;" << endl
<< endl
<< " /* Residual equations */" << endl
- << model_output.str();
-
- if (computeJacobian || computeJacobianExo)
- {
- DynamicOutput << " /* Jacobian */" << endl
- << " if (g1 == NULL)" << endl
- << " return;" << endl
- << " else" << endl
- << " {" << endl
- << jacobian_output.str()
- << " }" << endl;
- }
- if (computeHessian)
+ << model_output.str()
+ << " /* Jacobian */" << endl
+ << " if (g1 == NULL)" << endl
+ << " return;" << endl
+ << " else" << endl
+ << " {" << endl
+ << jacobian_output.str()
+ << " }" << endl;
+
+ if (second_derivatives.size())
{
DynamicOutput << " /* Hessian for endogenous and exogenous variables */" << endl
<< " if (g2 == NULL)" << endl
@@ -2169,22 +2150,43 @@ DynamicModel::BlockLinear(Model_Block *ModelBlock)
}
void
-DynamicModel::computingPass(const eval_context_type &eval_context, bool no_tmp_terms)
+DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives,
+ const eval_context_type &eval_context, bool no_tmp_terms)
{
+ if (!jacobianExo && (hessian || thirdDerivatives))
+ {
+ cerr << "DynamicModel::computingPass: computing 2nd or 3rd order derivatives imply computing 1st derivatives w.r. to exogenous" << endl;
+ exit(EXIT_FAILURE);
+ }
+
// Computes dynamic jacobian columns
- computeDynJacobianCols();
+ computeDynJacobianCols(jacobianExo);
- // Determine derivation order
- int order = 1;
- if (computeThirdDerivatives)
- order = 3;
- else if (computeHessian)
- order = 2;
+ // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic
+ set<int> vars;
+ for(int i = 0; i < getDerivIDNbr(); i++)
+ {
+ SymbolType type = getTypeByDerivID(i);
+ if (type == eEndogenous || (jacobianExo && (type == eExogenous || type == eExogenousDet)))
+ vars.insert(i);
+ }
// Launch computations
- cout << "Computing dynamic model derivatives at order " << order << "...";
- derive(order);
- cout << " done" << endl;
+ cout << "Computing dynamic model derivatives:" << endl
+ << " - order 1" << endl;
+ computeJacobian(vars);
+
+ if (hessian)
+ {
+ cout << " - order 2" << endl;
+ computeHessian(vars);
+ }
+
+ if (thirdDerivatives)
+ {
+ cout << " - order 3" << endl;
+ computeThirdDerivatives(vars);
+ }
if (mode == eSparseDLLMode || mode == eSparseMode)
{
@@ -2203,11 +2205,11 @@ DynamicModel::computingPass(const eval_context_type &eval_context, bool no_tmp_t
block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations);
BlockLinear(block_triangular.ModelBlock);
if (!no_tmp_terms)
- computeTemporaryTermsOrdered(order, block_triangular.ModelBlock);
+ computeTemporaryTermsOrdered(block_triangular.ModelBlock);
}
else
if (!no_tmp_terms)
- computeTemporaryTerms(order);
+ computeTemporaryTerms();
}
void
@@ -2261,10 +2263,6 @@ DynamicModel::toStatic(StaticModel &static_model) const
int
DynamicModel::computeDerivID(int symb_id, int lag)
{
- /* NOTE: we can't use computeJacobianExo here to decide
- which variables to include, since this boolean is not
- yet set at this point. */
-
// Setting maximum and minimum lags
if (max_lead < lag)
max_lead = lag;
@@ -2310,12 +2308,12 @@ DynamicModel::computeDerivID(int symb_id, int lag)
inv_deriv_id_table.push_back(key);
if (type == eEndogenous)
- var_endo_nbr++;
+ dynJacobianColsNbr++;
return deriv_id;
}
-int
+SymbolType
DynamicModel::getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException)
{
return symbol_table.getType(getSymbIDByDerivID(deriv_id));
@@ -2356,7 +2354,7 @@ DynamicModel::getDerivIDNbr() const
}
void
-DynamicModel::computeDynJacobianCols()
+DynamicModel::computeDynJacobianCols(bool jacobianExo)
{
/* Sort the dynamic endogenous variables by lexicographic order over (lag, type_specific_symbol_id)
and fill the dynamic columns for exogenous and exogenous deterministic */
@@ -2377,10 +2375,14 @@ DynamicModel::computeDynJacobianCols()
ordered_dyn_endo[make_pair(lag, tsid)] = deriv_id;
break;
case eExogenous:
- dyn_jacobian_cols_table[deriv_id] = var_endo_nbr + tsid;
+ // At this point, dynJacobianColsNbr contains the number of dynamic endogenous
+ if (jacobianExo)
+ dyn_jacobian_cols_table[deriv_id] = dynJacobianColsNbr + tsid;
break;
case eExogenousDet:
- dyn_jacobian_cols_table[deriv_id] = var_endo_nbr + symbol_table.exo_nbr() + tsid;
+ // At this point, dynJacobianColsNbr contains the number of dynamic endogenous
+ if (jacobianExo)
+ dyn_jacobian_cols_table[deriv_id] = dynJacobianColsNbr + symbol_table.exo_nbr() + tsid;
break;
default:
// Shut up GCC
@@ -2393,6 +2395,10 @@ DynamicModel::computeDynJacobianCols()
for(map<pair<int, int>, int>::const_iterator it = ordered_dyn_endo.begin();
it != ordered_dyn_endo.end(); it++)
dyn_jacobian_cols_table[it->second] = sorted_id++;
+
+ // Set final value for dynJacobianColsNbr
+ if (jacobianExo)
+ dynJacobianColsNbr += symbol_table.exo_nbr() + symbol_table.exo_det_nbr();
}
int
@@ -2405,11 +2411,3 @@ DynamicModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDExcepti
return it->second;
}
-int
-DynamicModel::getDynJacobianColsNbr() const
-{
- if (computeJacobianExo)
- return var_endo_nbr + symbol_table.exo_nbr() + symbol_table.exo_det_nbr();
- else
- return var_endo_nbr;
-}
diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh
index c246c0e9bbdcb0f26d5450835cd2ab63120222a8..871909a8b8e61b85e2e5b39fcc094a3b4bc25efe 100644
--- a/preprocessor/DynamicModel.hh
+++ b/preprocessor/DynamicModel.hh
@@ -39,10 +39,6 @@ private:
/*! Contains only endogenous, exogenous and exogenous deterministic */
map<int, int> dyn_jacobian_cols_table;
- //! Number of dynamic endogenous variables inside the model block
- /*! Set by computeDerivID() */
- int var_endo_nbr;
-
//! Maximum lag and lead over all types of variables (positive values)
/*! Set by computeDerivID() */
int max_lag, max_lead;
@@ -56,6 +52,10 @@ private:
/*! Set by computeDerivID() */
int max_exo_det_lag, max_exo_det_lead;
+ //! Number of columns of dynamic jacobian
+ /*! Set by computeDerivID() and computeDynJacobianCols() */
+ int dynJacobianColsNbr;
+
//! Writes dynamic model file (Matlab version)
void writeDynamicMFile(const string &dynamic_basename) const;
//! Writes dynamic model file (C version)
@@ -82,21 +82,19 @@ private:
//! Build The incidence matrix form the modeltree
void BuildIncidenceMatrix();
- void computeTemporaryTermsOrdered(int order, Model_Block *ModelBlock);
+ void computeTemporaryTermsOrdered(Model_Block *ModelBlock);
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, ExprNodeOutputType output_type, map_idx_type &map_idx) const;
virtual int computeDerivID(int symb_id, int lag);
- //! Get the number of columns of dynamic jacobian
- int getDynJacobianColsNbr() const;
//! Get the type corresponding to a derivation ID
- int getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
+ SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the lag corresponding to a derivation ID
int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Get the symbol ID corresponding to a derivation ID
int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Compute the column indices of the dynamic Jacobian
- void computeDynJacobianCols();
+ void computeDynJacobianCols(bool jacobianExo);
public:
DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
@@ -109,18 +107,16 @@ public:
double markowitz;
//! the file containing the model and the derivatives code
ofstream code_file;
- //! Whether dynamic Jacobian (w.r. to endogenous) should be written
- bool computeJacobian;
- //! Whether dynamic Jacobian (w.r. to endogenous and exogenous) should be written
- bool computeJacobianExo;
- //! Whether dynamic Hessian (w.r. to endogenous and exogenous) should be written
- bool computeHessian;
- //! Whether dynamic third order derivatives (w.r. to endogenous and exogenous) should be written
- bool computeThirdDerivatives;
//! Execute computations (variable sorting + derivation)
- /*! You must set computeJacobian, computeJacobianExo, computeHessian and computeThirdDerivatives to correct values before calling this function
- \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files */
- void computingPass(const eval_context_type &eval_context, bool no_tmp_terms);
+ /*!
+ \param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
+ \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true)
+ \param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true)
+ \param eval_context evaluation context for normalization
+ \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
+ */
+ void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives,
+ const eval_context_type &eval_context, bool no_tmp_terms);
//! Writes model initialization and lead/lag incidence matrix to output
void writeOutput(ostream &output) const;
//! Complete set to block decompose the model
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index b3ae75732c16c7b591cc744cced6351e9af843ae..7d1d75fc3644d0067ad6293730d4753cf7e2ca74 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -144,11 +144,12 @@ ModFile::computingPass(bool no_tmp_terms)
{
// Compute static model and its derivatives
dynamic_model.toStatic(static_model);
- static_model.computingPass(no_tmp_terms);
+ static_model.computingPass(false, no_tmp_terms);
// Set things to compute for dynamic model
+
if (mod_file_struct.simul_present)
- dynamic_model.computeJacobian = true;
+ dynamic_model.computingPass(false, false, false, global_eval_context, no_tmp_terms);
else
{
if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3)
@@ -156,14 +157,12 @@ ModFile::computingPass(bool no_tmp_terms)
cerr << "Incorrect order option..." << endl;
exit(EXIT_FAILURE);
}
- dynamic_model.computeJacobianExo = true;
- if (mod_file_struct.order_option >= 2)
- dynamic_model.computeHessian = true;
- if (mod_file_struct.order_option == 3)
- dynamic_model.computeThirdDerivatives = true;
+ bool hessian = mod_file_struct.order_option >= 2;
+ bool thirdDerivatives = mod_file_struct.order_option == 3;
+ dynamic_model.computingPass(true, hessian, thirdDerivatives, global_eval_context, no_tmp_terms);
}
- dynamic_model.computingPass(global_eval_context, no_tmp_terms);
}
+
for(vector<Statement *>::iterator it = statements.begin();
it != statements.end(); it++)
(*it)->computingPass();
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 1e2f392e14029952cc9a3b5256c8090d174f1970..75919b1cbabe81eafe98db23915ffbb6cdf5a102 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -48,65 +48,77 @@ ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
}
void
-ModelTree::derive(int order)
+ModelTree::computeJacobian(const set<int> &vars)
{
-
- for (int var = 0; var < getDerivIDNbr(); var++)
+ for(set<int>::const_iterator it = vars.begin();
+ it != vars.end(); it++)
for (int eq = 0; eq < (int) equations.size(); eq++)
{
- NodeID d1 = equations[eq]->getDerivative(var);
+ NodeID d1 = equations[eq]->getDerivative(*it);
if (d1 == Zero)
continue;
- first_derivatives[make_pair(eq, var)] = d1;
+ first_derivatives[make_pair(eq, *it)] = d1;
}
+}
- if (order >= 2)
+void
+ModelTree::computeHessian(const set<int> &vars)
+{
+ for (first_derivatives_type::const_iterator it = first_derivatives.begin();
+ it != first_derivatives.end(); it++)
{
- for (first_derivatives_type::const_iterator it = first_derivatives.begin();
- it != first_derivatives.end(); it++)
+ int eq = it->first.first;
+ int var1 = it->first.second;
+ NodeID d1 = it->second;
+
+ // Store only second derivatives with var2 <= var1
+ for(set<int>::const_iterator it2 = vars.begin();
+ it2 != vars.end(); it2++)
{
- int eq = it->first.first;
- int var1 = it->first.second;
- NodeID d1 = it->second;
-
- // Store only second derivatives with var2 <= var1
- for (int var2 = 0; var2 <= var1; var2++)
- {
- NodeID d2 = d1->getDerivative(var2);
- if (d2 == Zero)
- continue;
- second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
- }
+ int var2 = *it2;
+ if (var2 > var1)
+ continue;
+
+ NodeID d2 = d1->getDerivative(var2);
+ if (d2 == Zero)
+ continue;
+ second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
}
}
+}
- if (order >= 3)
+void
+ModelTree::computeThirdDerivatives(const set<int> &vars)
+{
+ for (second_derivatives_type::const_iterator it = second_derivatives.begin();
+ it != second_derivatives.end(); it++)
{
- for (second_derivatives_type::const_iterator it = second_derivatives.begin();
- it != second_derivatives.end(); it++)
+ int eq = it->first.first;
+
+ int var1 = it->first.second.first;
+ int var2 = it->first.second.second;
+ // By construction, var2 <= var1
+
+ NodeID d2 = it->second;
+
+ // Store only third derivatives such that var3 <= var2 <= var1
+ for(set<int>::const_iterator it2 = vars.begin();
+ it2 != vars.end(); it2++)
{
- int eq = it->first.first;
-
- int var1 = it->first.second.first;
- int var2 = it->first.second.second;
- // By construction, var2 <= var1
-
- NodeID d2 = it->second;
-
- // Store only third derivatives such that var3 <= var2 <= var1
- for (int var3 = 0; var3 <= var2; var3++)
- {
- NodeID d3 = d2->getDerivative(var3);
- if (d3 == Zero)
- continue;
- third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
- }
+ int var3 = *it2;
+ if (var3 > var2)
+ continue;
+
+ NodeID d3 = d2->getDerivative(var3);
+ if (d3 == Zero)
+ continue;
+ third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
}
}
}
void
-ModelTree::computeTemporaryTerms(int order)
+ModelTree::computeTemporaryTerms()
{
map<NodeID, int> reference_count;
temporary_terms.clear();
@@ -121,15 +133,13 @@ ModelTree::computeTemporaryTerms(int order)
it != first_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
- if (order >= 2)
- for (second_derivatives_type::iterator it = second_derivatives.begin();
- it != second_derivatives.end(); it++)
- it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
+ for (second_derivatives_type::iterator it = second_derivatives.begin();
+ it != second_derivatives.end(); it++)
+ it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
- if (order >= 3)
- for (third_derivatives_type::iterator it = third_derivatives.begin();
- it != third_derivatives.end(); it++)
- it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
+ for (third_derivatives_type::iterator it = third_derivatives.begin();
+ it != third_derivatives.end(); it++)
+ it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
}
void
diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh
index 6e6cb44b58147562598ecd6ca467002c0fe2e5dd..278623ec2ca01e5d21674a3e30371999df412b11 100644
--- a/preprocessor/ModelTree.hh
+++ b/preprocessor/ModelTree.hh
@@ -74,12 +74,20 @@ protected:
//! Temporary terms (those which will be noted Txxxx)
temporary_terms_type temporary_terms;
- //! Computes derivatives of ModelTree
- void derive(int order);
+ //! Computes 1st derivatives
+ /*! \param vars the derivation IDs w.r. to which compute the derivatives */
+ void computeJacobian(const set<int> &vars);
+ //! Computes 2nd derivatives
+ /*! \param vars the derivation IDs w.r. to which derive the 1st derivatives */
+ void computeHessian(const set<int> &vars);
+ //! Computes 3rd derivatives
+ /*! \param vars the derivation IDs w.r. to which derive the 2nd derivatives */
+ void computeThirdDerivatives(const set<int> &vars);
+
//! 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;
- //! Computes temporary terms
- void computeTemporaryTerms(int order);
+ //! Computes temporary terms (for all equations and derivatives)
+ void computeTemporaryTerms();
//! Writes temporary terms
void writeTemporaryTerms(ostream &output, ExprNodeOutputType output_type) const;
//! Writes model local variables
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 8cc372575726881891eb5be5c1e47db03984b6c5..90aa477532e68ae65a9c9c88a0a80b8125ef0685 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -23,8 +23,7 @@
StaticModel::StaticModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg) :
- ModelTree(symbol_table_arg, num_constants_arg),
- computeStaticHessian(false)
+ ModelTree(symbol_table_arg, num_constants_arg)
{
}
@@ -151,38 +150,37 @@ StaticModel::writeStaticModel(ostream &StaticOutput) const
jacobian_output << ";" << endl;
}
- // Write Hessian w.r. to endogenous only
- if (computeStaticHessian)
- for (second_derivatives_type::const_iterator it = second_derivatives.begin();
- it != second_derivatives.end(); it++)
- {
- int eq = it->first.first;
- int symb_id1 = inv_deriv_id_table[it->first.second.first];
- int symb_id2 = inv_deriv_id_table[it->first.second.second];
- NodeID d2 = it->second;
-
- int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
- int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
-
- int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
- int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
-
- hessian_output << " g2";
- matrixHelper(hessian_output, eq, col_nb, output_type);
- hessian_output << " = ";
- d2->writeOutput(hessian_output, output_type, temporary_terms);
- hessian_output << ";" << endl;
-
- // Treating symetric elements
- if (symb_id1 != symb_id2)
- {
- lsymetric << " g2";
- matrixHelper(lsymetric, eq, col_nb_sym, output_type);
- lsymetric << " = " << "g2";
- matrixHelper(lsymetric, eq, col_nb, output_type);
- lsymetric << ";" << endl;
- }
- }
+ // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
+ for (second_derivatives_type::const_iterator it = second_derivatives.begin();
+ it != second_derivatives.end(); it++)
+ {
+ int eq = it->first.first;
+ int symb_id1 = inv_deriv_id_table[it->first.second.first];
+ int symb_id2 = inv_deriv_id_table[it->first.second.second];
+ NodeID d2 = it->second;
+
+ int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
+ int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
+
+ int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
+ int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
+
+ hessian_output << " g2";
+ matrixHelper(hessian_output, eq, col_nb, output_type);
+ hessian_output << " = ";
+ d2->writeOutput(hessian_output, output_type, temporary_terms);
+ hessian_output << ";" << endl;
+
+ // Treating symetric elements
+ if (symb_id1 != symb_id2)
+ {
+ lsymetric << " g2";
+ matrixHelper(lsymetric, eq, col_nb_sym, output_type);
+ lsymetric << " = " << "g2";
+ matrixHelper(lsymetric, eq, col_nb, output_type);
+ lsymetric << ";" << endl;
+ }
+ }
// Writing ouputs
if (mode != eDLLMode)
@@ -208,9 +206,11 @@ StaticModel::writeStaticModel(ostream &StaticOutput) const
<< " g1 = real(g1)+2*imag(g1);" << endl
<< " end" << endl
<< "end" << endl;
- if (computeStaticHessian)
+
+ // If 2nd order derivatives have been computed
+ if (second_derivatives.size())
{
- StaticOutput << "if nargout >= 3,\n";
+ StaticOutput << "if nargout >= 3," << endl;
// Writing initialization instruction for matrix g2
int ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
StaticOutput << " g2 = sparse([],[],[], " << equations.size() << ", " << ncols << ", " << 5*ncols << ");" << endl
@@ -266,20 +266,26 @@ StaticModel::writeStaticFile(const string &basename) const
}
void
-StaticModel::computingPass(bool no_tmp_terms)
+StaticModel::computingPass(bool hessian, bool no_tmp_terms)
{
- // Determine derivation order
- int order = 1;
- if (computeStaticHessian)
- order = 2;
+ // Compute derivatives w.r. to all derivation IDs (i.e. all endogenous)
+ set<int> vars;
+ for(int i = 0; i < getDerivIDNbr(); i++)
+ vars.insert(i);
// Launch computations
- cout << "Computing static model derivatives at order " << order << "...";
- derive(order);
- cout << " done" << endl;
+ cout << "Computing static model derivatives:" << endl
+ << " - order 1" << endl;
+ computeJacobian(vars);
+
+ if (hessian)
+ {
+ cout << " - order 2" << endl;
+ computeHessian(vars);
+ }
if (!no_tmp_terms)
- computeTemporaryTerms(order);
+ computeTemporaryTerms();
}
int
diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh
index 407e3b53243ceb1f7e7b25c2bfab979e776c480f..50882b941cc9cf0105cca8d90ec6f0f2be53f089 100644
--- a/preprocessor/StaticModel.hh
+++ b/preprocessor/StaticModel.hh
@@ -46,12 +46,11 @@ private:
public:
StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
- //! Whether static Hessian (w.r. to endogenous only) should be written
- bool computeStaticHessian;
//! Execute computations (derivation)
/*! You must set computeStaticHessian before calling this function
+ \param hessian whether Hessian (w.r. to endogenous only) should be computed
\param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
- void computingPass(bool no_tmp_terms);
+ void computingPass(bool hessian, bool no_tmp_terms);
//! Writes static model file
void writeStaticFile(const string &basename) const;