Commit ec215e2a authored by Houtan Bastani's avatar Houtan Bastani
Browse files

separate temporary terms: WIP

parent 7cbd7b6e
......@@ -449,6 +449,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
else
sps = "";
// The equations
temporary_terms_idxs_t temporary_terms_idxs;
for (unsigned int i = 0; i < block_size; i++)
{
temporary_terms_t tt2;
......@@ -460,12 +461,12 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
it != v_temporary_terms[block][i].end(); it++)
{
if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
(*it)->writeExternalFunctionOutput(output, local_output_type, tt2, tef_terms);
(*it)->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
output << " " << sps;
(*it)->writeOutput(output, local_output_type, local_temporary_terms, tef_terms);
(*it)->writeOutput(output, local_output_type, local_temporary_terms, temporary_terms_idxs, tef_terms);
output << " = ";
(*it)->writeOutput(output, local_output_type, tt2, tef_terms);
(*it)->writeOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
// Insert current node into tt2
tt2.insert(*it);
output << ";" << endl;
......@@ -1522,49 +1523,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
void
DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
{
string filename = dynamic_basename + ".m";
ofstream mDynamicModelFile;
mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
if (!mDynamicModelFile.is_open())
{
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x, params, steady_state, it_)" << endl
<< "%" << endl
<< "% Status : Computes dynamic model for Dynare" << endl
<< "%" << endl
<< "% Inputs :" << endl
<< "% y [#dynamic variables by 1] double vector of endogenous variables in the order stored" << endl
<< "% in M_.lead_lag_incidence; see the Manual" << endl
<< "% x [nperiods by M_.exo_nbr] double matrix of exogenous variables (in declaration order)" << endl
<< "% for all simulation periods" << endl
<< "% steady_state [M_.endo_nbr by 1] double vector of steady state values" << endl
<< "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
<< "% it_ scalar double time period for exogenous variables for which to evaluate the model" << endl
<< "%" << endl
<< "% Outputs:" << endl
<< "% residual [M_.endo_nbr by 1] double vector of residuals of the dynamic model equations in order of " << endl
<< "% declaration of the equations." << endl
<< "% Dynare may prepend auxiliary equations, see M_.aux_vars" << endl
<< "% g1 [M_.endo_nbr by #dynamic variables] double Jacobian matrix of the dynamic model equations;" << endl
<< "% rows: equations in order of declaration" << endl
<< "% columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names" << endl
<< "% g2 [M_.endo_nbr by (#dynamic variables)^2] double Hessian matrix of the dynamic model equations;" << endl
<< "% rows: equations in order of declaration" << endl
<< "% columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names" << endl
<< "% g3 [M_.endo_nbr by (#dynamic variables)^3] double Third order derivative matrix of the dynamic model equations;" << endl
<< "% rows: equations in order of declaration" << endl
<< "% columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names" << endl
<< "%" << endl
<< "%" << endl
<< "% Warning : this file is generated automatically by Dynare" << endl
<< "% from model file (.mod)" << endl << endl;
writeDynamicModel(mDynamicModelFile, false, false);
mDynamicModelFile << "end" << endl; // Close *_dynamic function
mDynamicModelFile.close();
writeDynamicModel(dynamic_basename, false, false);
}
void
......@@ -1600,26 +1559,7 @@ DynamicModel::writeVarExpectationCalls(ostream &output) const
void
DynamicModel::writeDynamicJuliaFile(const string &basename) const
{
string filename = basename + "Dynamic.jl";
ofstream output;
output.open(filename.c_str(), ios::out | ios::binary);
if (!output.is_open())
{
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "module " << basename << "Dynamic" << endl
<< "#" << endl
<< "# NB: this file was automatically generated by Dynare" << endl
<< "# from " << basename << ".mod" << endl
<< "#" << endl
<< "using Utils" << endl << endl
<< "export dynamic!" << endl << endl;
writeDynamicModel(output, false, true);
output << "end" << endl;
output.close();
writeDynamicModel(basename, false, true);
}
void
......@@ -2201,200 +2141,374 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
chdir("..");
}
void
DynamicModel::writeWrapperFunctions(const string &basename, const string &ending) const
{
string name;
if (ending == "g1")
name = basename + "_resid_g1";
else if (ending == "g2")
name= basename + "_resid_g1_g2";
else if (ending == "g3")
name = basename + "_resid_g1_g2_g3";
string filename = name + ".m";
ofstream output;
output.open(filename.c_str(), ios::out | ios::binary);
if (!output.is_open())
{
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
if (ending == "g1")
output << "function [residual, g1] = " << name << "(T, y, x, params, steady_state, it_)" << endl
<< "% function [residual, g1] = " << name << "(T, y, x, params, steady_state, it_)" << endl;
else if (ending == "g2")
output << "function [residual, g1, g2] = " << name << "(T, y, x, params, steady_state, it_)" << endl
<< "% function [residual, g1, g2] = " << name << "(T, y, x, params, steady_state, it_)" << endl;
else if (ending == "g3")
output << "function [residual, g1, g2, g3] = " << name << "(T, y, x, params, steady_state, it_)" << endl
<< "% function [residual, g1, g2, g3] = " << name << "(T, y, x, params, steady_state, it_)" << endl;
output << "%" << endl
<< "% Wrapper function automatically created by Dynare" << endl
<< "%" << endl << endl;
if (ending == "g1")
output << " T = " << basename + "_" + ending + "_tt(T, y, x, params, steady_state, it_);" << endl
<< " residual = " << basename + "_resid(T, y, x, params, steady_state, it_, false);" << endl
<< " g1 = " << basename + "_g1(T, y, x, params, steady_state, it_, false);" << endl;
else if (ending == "g2")
output << " T = " << basename + "_" + ending + "_tt(T, y, x, params, steady_state, it_);" << endl
<< " [residual, g1] = " << basename + "_resid_g1(T, y, x, params, steady_state, it_, false);" << endl
<< " g2 = " << basename + "_g2(T, y, x, params, steady_state, it_, false);" << endl;
else if (ending == "g3")
output << " T = " << basename + "_" + ending + "_tt(T, y, x, params, steady_state, it_);" << endl
<< " [residual, g1, g2] = " << basename + "_resid_g1(T, y, x, params, steady_state, it_, false);" << endl
<< " g3 = " << basename + "_g3(T, y, x, params, steady_state, it_, false);" << endl;
output << endl << "end" << endl;
output.close();
}
void
DynamicModel::writeDynamicModelHelper(const string &name, const string &retvalname,
const string &name_tt, size_t ttlen,
const string &previous_tt_name,
const ostringstream &init_s,
const ostringstream &end_s,
const ostringstream &s, const ostringstream &s_tt) const
{
string filename = name_tt + ".m";
ofstream output;
output.open(filename.c_str(), ios::out | ios::binary);
if (!output.is_open())
{
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "function T = " << name_tt << "(T, y, x, params, steady_state, it_)" << endl
<< "% function T = " << name_tt << "(T, y, x, params, steady_state, it_)" << endl
<< "%" << endl
<< "% File created by Dynare Preprocessor from .mod file" << endl
<< "%" << endl
<< "% Inputs:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms to be filled by function" << endl
<< "% y [#dynamic variables by 1] double vector of endogenous variables in the order stored" << endl
<< "% in M_.lead_lag_incidence; see the Manual" << endl
<< "% x [nperiods by M_.exo_nbr] double matrix of exogenous variables (in declaration order)" << endl
<< "% for all simulation periods" << endl
<< "% steady_state [M_.endo_nbr by 1] double vector of steady state values" << endl
<< "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
<< "% it_ scalar double time period for exogenous variables for which" << endl
<< "% to evaluate the model" << endl
<< "%" << endl
<< "% Output:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms" << endl
<< "%" << endl << endl
<< "assert(length(T) >= " << ttlen << ");" << endl
<< endl;
if (!previous_tt_name.empty())
output << "T = " << previous_tt_name << "(T, y, x, params, steady_state, it_);" << endl << endl;
output << s_tt.str() << endl
<< "end" << endl;
output.close();
filename = name + ".m";
output.open(filename.c_str(), ios::out | ios::binary);
if (!output.is_open())
{
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "function " << retvalname << " = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "% function " << retvalname << " = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "%" << endl
<< "% File created by Dynare Preprocessor from .mod file" << endl
<< "%" << endl
<< "% Inputs:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms to be filled by function" << endl
<< "% y [#dynamic variables by 1] double vector of endogenous variables in the order stored" << endl
<< "% in M_.lead_lag_incidence; see the Manual" << endl
<< "% x [nperiods by M_.exo_nbr] double matrix of exogenous variables (in declaration order)" << endl
<< "% for all simulation periods" << endl
<< "% steady_state [M_.endo_nbr by 1] double vector of steady state values" << endl
<< "% params [M_.param_nbr by 1] double vector of parameter values in declaration order" << endl
<< "% it_ scalar double time period for exogenous variables for which" << endl
<< "% to evaluate the model" << endl
<< "% T_flag boolean boolean flag saying whether or not to calculate temporary terms" << endl
<< "%" << endl
<< "% Output:" << endl
<< "% " << retvalname << endl
<< "%" << endl << endl;
if (!name_tt.empty())
output << "if T_flag" << endl
<< " T = " << name_tt << "(T, y, x, params, steady_state, it_);" << endl
<< "end" << endl;
output << init_s.str() << endl
<< s.str()
<< end_s.str() << endl
<< "end" << endl;
output.close();
}
void
DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const
{
ostringstream model_local_vars_output; // Used for storing model local vars
ostringstream model_output; // Used for storing model temp vars and equations
ostringstream jacobian_output; // Used for storing jacobian equations
ostringstream hessian_output; // Used for storing Hessian equations
ostringstream third_derivatives_output; // Used for storing third order derivatives equations
writeDynamicModel("", DynamicOutput, use_dll, julia);
}
void
DynamicModel::writeDynamicModel(const string &dynamic_basename, bool use_dll, bool julia) const
{
ofstream DynamicOutput;
writeDynamicModel(dynamic_basename, DynamicOutput, use_dll, julia);
}
void
DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &DynamicOutput, bool use_dll, bool julia) const
{
ostringstream model_tt_output; // Used for storing model temp vars
ostringstream model_output; // Used for storing model equations
ostringstream jacobian_tt_output; // Used for storing jacobian temp vars
ostringstream jacobian_output; // Used for storing jacobian equations
ostringstream hessian_tt_output; // Used for storing Hessian temp vars
ostringstream hessian_output; // Used for storing Hessian equations
ostringstream third_derivatives_tt_output; // Used for storing third order derivatives temp terms
ostringstream third_derivatives_output; // Used for storing third order derivatives equations
ExprNodeOutputType output_type = (use_dll ? oCDynamicModel :
julia ? oJuliaDynamicModel : oMatlabDynamicModel);
deriv_node_temp_terms_t tef_terms;
temporary_terms_t temp_term_empty;
temporary_terms_t temp_term_union = temporary_terms_res;
temporary_terms_t temp_term_union_m_1;
writeModelLocalVariables(model_local_vars_output, output_type, tef_terms);
writeTemporaryTerms(temporary_terms_res, temp_term_union_m_1, model_output, output_type, tef_terms);
writeModelEquations(model_output, output_type);
temporary_terms_t temp_term_union;
temporary_terms_idxs_t temp_term_union_idxs;
for (map<expr_t, expr_t>::const_iterator it = temporary_terms_mlv.begin();
it != temporary_terms_mlv.end(); it++)
temp_term_union.insert(it->first);
writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_mlv,
temporary_terms_mlv_idxs,
model_tt_output, output_type, tef_terms);
temp_term_union_idxs = temporary_terms_mlv_idxs;
writeTemporaryTerms(temporary_terms_res, temporary_terms_res_idxs,
temp_term_union, temp_term_union_idxs,
model_tt_output, output_type, tef_terms);
temp_term_union.insert(temporary_terms_res.begin(), temporary_terms_res.end());
temp_term_union_idxs.insert(temp_term_union_idxs.end(),
temporary_terms_res_idxs.begin(), temporary_terms_res_idxs.end());
writeModelEquations(model_output, output_type, temp_term_union, temp_term_union_idxs);
int nrows = equations.size();
int hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr;
// Writing Jacobian
temp_term_union_m_1 = temp_term_union;
temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
if (!first_derivatives.empty())
if (julia)
writeTemporaryTerms(temp_term_union, temp_term_empty, jacobian_output, output_type, tef_terms);
else
writeTemporaryTerms(temp_term_union, temp_term_union_m_1, jacobian_output, output_type, tef_terms);
for (first_derivatives_t::const_iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
int eq = it->first.first;
int var = it->first.second;
expr_t d1 = it->second;
writeTemporaryTerms(temporary_terms_g1, temporary_terms_g1_idxs,
temp_term_union, temp_term_union_idxs,
jacobian_tt_output, output_type, tef_terms);
temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
temp_term_union_idxs.insert(temp_term_union_idxs.end(),
temporary_terms_g1_idxs.begin(), temporary_terms_g1_idxs.end());
jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
jacobian_output << "=";
d1->writeOutput(jacobian_output, output_type, temp_term_union, tef_terms);
jacobian_output << ";" << endl;
for (first_derivatives_t::const_iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
int eq = it->first.first;
int var = it->first.second;
expr_t d1 = it->second;
jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
jacobian_output << "=";
d1->writeOutput(jacobian_output, output_type,
temp_term_union, temp_term_union_idxs, tef_terms);
jacobian_output << ";" << endl;
}
}
// Writing Hessian
temp_term_union_m_1 = temp_term_union;
temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
if (!second_derivatives.empty())
if (julia)
writeTemporaryTerms(temp_term_union, temp_term_empty, hessian_output, output_type, tef_terms);
else
writeTemporaryTerms(temp_term_union, temp_term_union_m_1, hessian_output, output_type, tef_terms);
int k = 0; // Keep the line of a 2nd derivative in v2
for (second_derivatives_t::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;
expr_t d2 = it->second;
writeTemporaryTerms(temporary_terms_g2, temporary_terms_g2_idxs,
temp_term_union, temp_term_union_idxs,
hessian_tt_output, output_type, tef_terms);
temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
temp_term_union_idxs.insert(temp_term_union_idxs.end(),
temporary_terms_g2_idxs.begin(), temporary_terms_g2_idxs.end());
int id1 = getDynJacobianCol(var1);
int id2 = getDynJacobianCol(var2);
int k = 0; // Keep the line of a 2nd derivative in v2
for (second_derivatives_t::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;
expr_t d2 = it->second;
int col_nb = id1 * dynJacobianColsNbr + id2;
int col_nb_sym = id2 * dynJacobianColsNbr + id1;
int id1 = getDynJacobianCol(var1);
int id2 = getDynJacobianCol(var2);
ostringstream for_sym;
if (output_type == oJuliaDynamicModel)
{
for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
hessian_output << " @inbounds " << for_sym.str() << " = ";
d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
hessian_output << endl;
}
else
{
sparseHelper(2, hessian_output, k, 0, output_type);
hessian_output << "=" << eq + 1 << ";" << endl;
int col_nb = id1 * dynJacobianColsNbr + id2;
int col_nb_sym = id2 * dynJacobianColsNbr + id1;
sparseHelper(2, hessian_output, k, 1, output_type);
hessian_output << "=" << col_nb + 1 << ";" << endl;
ostringstream for_sym;
if (output_type == oJuliaDynamicModel)
{
for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
hessian_output << " @inbounds " << for_sym.str() << " = ";
d2->writeOutput(hessian_output, output_type, temp_term_union, temp_term_union_idxs, tef_terms);
hessian_output << endl;
}
else
{
sparseHelper(2, hessian_output, k, 0, output_type);
hessian_output << "=" << eq + 1 << ";" << endl;
sparseHelper(2, hessian_output, k, 2, output_type);
hessian_output << "=";
d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
hessian_output << ";" << endl;
sparseHelper(2, hessian_output, k, 1, output_type);
hessian_output << "=" << col_nb + 1 << ";" << endl;
k++;
}
sparseHelper(2, hessian_output, k, 2, output_type);
hessian_output << "=";
d2->writeOutput(hessian_output, output_type, temp_term_union, temp_term_union_idxs, tef_terms);
hessian_output << ";" << endl;
// Treating symetric elements
if (id1 != id2)
if (output_type == oJuliaDynamicModel)
hessian_output << " @inbounds g2[" << eq + 1 << "," << col_nb_sym + 1 << "] = "
<< for_sym.str() << endl;
else
{
sparseHelper(2, hessian_output, k, 0, output_type);
hessian_output << "=" << eq + 1 << ";" << endl;
k++;
}
// Treating symetric elements
if (id1 != id2)
if (output_type == oJuliaDynamicModel)
hessian_output << " @inbounds g2[" << eq + 1 << "," << col_nb_sym + 1 << "] = "
<< for_sym.str() << endl;
else
{
sparseHelper(2, hessian_output, k, 0, output_type);
hessian_output << "=" << eq + 1 << ";" << endl;
sparseHelper(2, hessian_output, k, 1, output_type);
hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
sparseHelper(2, hessian_output, k, 1, output_type);
hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
sparseHelper(2, hessian_output, k, 2, output_type);
hessian_output << "=";
sparseHelper(2, hessian_output, k-1, 2, output_type);
hessian_output << ";" << endl;
sparseHelper(2, hessian_output, k, 2, output_type);
hessian_output << "=";
sparseHelper(2, hessian_output, k-1, 2, output_type);
hessian_output << ";" << endl;
k++;
}
k++;
}
}
}
// Writing third derivatives
temp_term_union_m_1 = temp_term_union;
temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
if (!third_derivatives.empty())
if (julia)
writeTemporaryTerms(temp_term_union, temp_term_empty, third_derivatives_output, output_type, tef_terms);
else
writeTemporaryTerms(temp_term_union, temp_term_union_m_1, third_derivatives_output, output_type, tef_terms);
k = 0; // Keep the line of a 3rd derivative in v3
for (third_derivatives_t::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;
expr_t d3 = it->second;
int id1 = getDynJacobianCol(var1);
int id2 = getDynJacobianCol(var2);
int id3 = getDynJacobianCol(var3);
writeTemporaryTerms(temporary_terms_g3, temporary_terms_g3_idxs,
temp_term_union, temp_term_union_idxs,
third_derivatives_tt_output, output_type, tef_terms);
temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
temp_term_union_idxs.insert(temp_term_union_idxs.end(),
temporary_terms_g3_idxs.begin(), temporary_terms_g3_idxs.end());
// Reference column number for the g3 matrix
int ref_col = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
ostringstream for_sym;
if (output_type == oJuliaDynamicModel)
{
for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
third_derivatives_output << " @inbounds " << for_sym.str() << " = ";
d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
third_derivatives_output << endl;
}
else
int k = 0; // Keep the line of a 3rd derivative in v3
for (third_derivatives_t::const_iterator it = third_derivatives.begin();
it != third_derivatives.end(); it++)
{
sparseHelper(3, third_derivatives_output, k, 0, output_type);
third_derivatives_output << "=" << eq + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k, 1, output_type);
third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
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;
expr_t d3 = it->second;
sparseHelper(3, third_derivatives_output, k, 2, output_type);
third_derivatives_output << "=";
d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
third_derivatives_output << ";" << endl;
}
int id1 = getDynJacobianCol(var1);
int id2 = getDynJacobianCol(var2);
int id3 = getDynJacobianCol(var3);
// 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);
// Reference column number for the g3 matrix
int ref_col = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
int k2 = 1; // Keeps the offset of the permutation relative to k
for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
if (*it2 != ref_col)
ostringstream for_sym;
if (output_type == oJuliaDynamicModel)
third_derivatives_output << " @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
<< for_sym.str() << endl;
{
for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
third_derivatives_output << " @inbounds " << for_sym.str() << " = ";
d3->writeOutput(third_derivatives_output, output_type, temp_term_union, temp_term_union_idxs, tef_terms);
third_derivatives_output << endl;
}
else
{
sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
sparseHelper(3, third_derivatives_output, k, 0, output_type);
third_derivatives_output << "=" << eq + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k, 1, output_type);
third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
third_derivatives_output << "=";
sparseHelper(3, third_derivatives_output, k, 2, output_type);
third_derivatives_output << "=";
d3->writeOutput(third_derivatives_output, output_type, temp_term_union, temp_term_union_idxs, tef_terms);
third_derivatives_output << ";" << endl;
k2++;
}
k += k2;
// 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);
int k2 = 1; // Keeps the offset of the permutation relative to k
for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
if (*it2 != ref_col)
if (output_type == oJuliaDynamicModel)
third_derivatives_output << " @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
<< for_sym.str() << endl;
else
{
sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
third_derivatives_output << "=" << eq + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
third_derivatives_output << "=";
sparseHelper(3, third_derivatives_output, k, 2, output_type);
third_derivatives_output << ";" << endl;
k2++;
}
k += k2;
}
}
if (output_type == oMatlabDynamicModel)
......@@ -2403,192 +2517,306 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
map<string, string> tmp_paren_vars;
bool message_printed = false;
fixNestedParenthesis(model_output, tmp_paren_vars, message_printed);
fixNestedParenthesis(model_local_vars_output, tmp_paren_vars, message_printed);
fixNestedParenthesis(model_tt_output, tmp_paren_vars, message_printed);
fixNestedParenthesis(jacobian_output, tmp_paren_vars, message_printed);
fixNestedParenthesis(jacobian_tt_output, tmp_paren_vars, message_printed);
fixNestedParenthesis(hessian_output, tmp_paren_vars, message_printed);
fixNestedParenthesis(hessian_tt_output, tmp_paren_vars, message_printed);
fixNestedParenthesis(third_derivatives_output, tmp_paren_vars, message_printed);
DynamicOutput << "%" << endl
<< "% Model equations" << endl
<< "%" << endl
<< endl;
writeVarExpectationCalls(DynamicOutput);
DynamicOutput << "residual = zeros(" << nrows << ", 1);" << endl
<< model_local_vars_output.str()
<< model_output.str()
// Writing initialization instruction for matrix g1
<< "if nargout >= 2," << endl
<< " g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");" << endl
<< endl
<< " %" << endl
<< " % Jacobian matrix" << endl