diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 92fcc25326e3db2a87c18901c4030d215ff31456..9f206d63860a66981fd0bf287b37129439f3224f 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3972,7 +3972,56 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode
else if (julia)
writeDynamicJuliaFile(basename);
else
- writeDynamicMFile(t_basename);
+ {
+ writeDynamicMFile(t_basename);
+ writeSetAuxiliaryVariables(t_basename, julia);
+ }
+}
+
+void
+DynamicModel::writeSetAuxiliaryVariables(const string &basename, const bool julia) const
+{
+ ostringstream output_func_body;
+ writeAuxVarRecursiveDefinitions(output_func_body, oMatlabDseries);
+
+ if (output_func_body.str().empty())
+ return;
+
+ string func_name = basename + "_set_auxiliary_variables_as_dseries";
+ string filename = julia ? func_name + ".jl" : func_name + ".m";
+ string comment = julia ? "#" : "%";
+
+ 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 ds = " << func_name + "(ds, params)" << endl
+ << comment << endl
+ << comment << " Status : Computes Auxiliary variables of the dynamic model and returns a dseries" << endl
+ << comment << endl
+ << comment << " Warning : this file is generated automatically by Dynare" << endl
+ << comment << " from model file (.mod)" << endl << endl
+ << output_func_body.str();
+}
+
+void
+DynamicModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const
+{
+ deriv_node_temp_terms_t tef_terms;
+ temporary_terms_t temporary_terms;
+ for (int i = 0; i < (int) aux_equations.size(); i++)
+ if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
+ dynamic_cast<ExprNode *>(aux_equations[i])->writeExternalFunctionOutput(output, output_type,
+ temporary_terms, tef_terms);
+ for (int i = 0; i < (int) aux_equations.size(); i++)
+ {
+ dynamic_cast<ExprNode *>(aux_equations[i])->writeOutput(output, output_type, temporary_terms, tef_terms);
+ output << ";" << endl;
+ }
}
void
@@ -4848,6 +4897,39 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
aux_equations[i] = substeq;
}
+ // Substitute in diff_aux_equations
+ // Without this loop, the auxiliary equations in equations
+ // will diverge from those in diff_aux_equations
+ for (int i = 0; i < (int) diff_aux_equations.size(); i++)
+ {
+ expr_t subst;
+ switch (type)
+ {
+ case avEndoLead:
+ subst = diff_aux_equations[i]->substituteEndoLeadGreaterThanTwo(subst_table,
+ neweqs, deterministic_model);
+ break;
+ case avEndoLag:
+ subst = diff_aux_equations[i]->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+ break;
+ case avExoLead:
+ subst = diff_aux_equations[i]->substituteExoLead(subst_table, neweqs, deterministic_model);
+ break;
+ case avExoLag:
+ subst = diff_aux_equations[i]->substituteExoLag(subst_table, neweqs);
+ break;
+ case avDiffForward:
+ subst = diff_aux_equations[i]->differentiateForwardVars(subset, subst_table, neweqs);
+ break;
+ default:
+ cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
+ exit(EXIT_FAILURE);
+ }
+ BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(subst);
+ assert(substeq != NULL);
+ diff_aux_equations[i] = substeq;
+ }
+
// Add new equations
for (int i = 0; i < (int) neweqs.size(); i++)
addEquation(neweqs[i], -1);
@@ -4928,12 +5010,18 @@ DynamicModel::substituteDiff(StaticModel &static_model)
for (int i = 0; i < (int) neweqs.size(); i++)
addEquation(neweqs[i], -1);
- copy(neweqs.begin(), neweqs.end(), back_inserter(aux_equations));
+ copy(neweqs.begin(), neweqs.end(), back_inserter(diff_aux_equations));
if (diff_subst_table.size() > 0)
cout << "Substitution of Diff operator: added " << neweqs.size() << " auxiliary variables and equations." << endl;
}
+void
+DynamicModel::combineDiffAuxEquations()
+{
+ copy(diff_aux_equations.begin(), diff_aux_equations.end(), back_inserter(aux_equations));
+}
+
void
DynamicModel::substituteExpectation(bool partial_information_model)
{
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index dda2a2aaba63f5767580884a90d6b6fe1b869b31..ad41e03f5f9c7b2dead60ab4a4bb761296cf3a91 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -116,6 +116,9 @@ private:
//! Writes the code of the model in virtual machine bytecode
void writeModelEquationsCode(string &file_name, const string &bin_basename, const map_idx_t &map_idx) const;
+ void writeSetAuxiliaryVariables(const string &basename, const bool julia) const;
+ void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const;
+
//! Computes jacobian and prepares for equation normalization
/*! Using values from initval/endval blocks and parameter initializations:
- computes the jacobian for the model w.r. to contemporaneous variables
@@ -394,6 +397,9 @@ public:
//! Substitutes diff operator
void substituteDiff(StaticModel &static_model);
+ //! Adds contents of diff_aux_equations to the back of aux_equations
+ void combineDiffAuxEquations();
+
//! Fill var_expectation_functions_to_write
void fillVarExpectationFunctionsToWrite();
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index d0c28be33e9a61508f215e81481cde6d83932134..fa0cdf4275616fce7c82c59b91fa8de01a423af5 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -850,6 +850,11 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case oCSteadyStateFile:
output << "ys_[" << tsid << "]";
break;
+ case oMatlabDseries:
+ output << "ds." << datatree.symbol_table.getName(symb_id);
+ if (lag != 0)
+ output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag << RIGHT_ARRAY_SUBSCRIPT(output_type);
+ break;
default:
cerr << "VariableNode::writeOutput: should not reach this point" << endl;
exit(EXIT_FAILURE);
@@ -902,6 +907,11 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case oCSteadyStateFile:
output << "exo_[" << i - 1 << "]";
break;
+ case oMatlabDseries:
+ output << "ds." << datatree.symbol_table.getName(symb_id);
+ if (lag != 0)
+ output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag << RIGHT_ARRAY_SUBSCRIPT(output_type);
+ break;
default:
cerr << "VariableNode::writeOutput: should not reach this point" << endl;
exit(EXIT_FAILURE);
@@ -954,6 +964,11 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case oCSteadyStateFile:
output << "exo_[" << i - 1 << "]";
break;
+ case oMatlabDseries:
+ output << "ds." << datatree.symbol_table.getName(symb_id);
+ if (lag != 0)
+ output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag << RIGHT_ARRAY_SUBSCRIPT(output_type);
+ break;
default:
cerr << "VariableNode::writeOutput: should not reach this point" << endl;
exit(EXIT_FAILURE);
@@ -2907,40 +2922,55 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
diff_table_t::iterator it = diff_table.find(sthis);
assert(it != diff_table.end() && it->second[-arg->maxLag()] == this);
- int origin_arg_max_lag;
- VariableNode *origin_aux_var;
+ int last_arg_max_lag;
+ VariableNode *last_aux_var;
for (map<int, expr_t>::reverse_iterator rit = it->second.rbegin();
rit != it->second.rend(); rit++)
{
expr_t argsubst = dynamic_cast<UnaryOpNode *>(rit->second)->
get_arg()->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
-
int symb_id;
VariableNode *vn = dynamic_cast<VariableNode *>(argsubst);
- if (vn != NULL)
- symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst, vn->get_symb_id(), vn->get_lag());
- else
- symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
-
if (rit == it->second.rbegin())
{
+ if (vn != NULL)
+ symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst, vn->get_symb_id(), vn->get_lag());
+ else
+ symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
+
// make originating aux var & equation
- origin_arg_max_lag = rit->first;
- origin_aux_var = datatree.AddVariable(symb_id, 0);
+ last_arg_max_lag = rit->first;
+ last_aux_var = datatree.AddVariable(symb_id, 0);
//ORIG_AUX_DIFF = argsubst - argsubst(-1)
- neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(origin_aux_var,
+ neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(last_aux_var,
datatree.AddMinus(argsubst,
argsubst->decreaseLeadsLags(1)))));
- subst_table[rit->second] = dynamic_cast<VariableNode *>(origin_aux_var);
+ subst_table[rit->second] = dynamic_cast<VariableNode *>(last_aux_var);
}
else
{
- // just add equation of form: AUX_DIFF = ORIG_AUX_DIFF(origin_arg_max_lag - rit->first)
- VariableNode *new_aux_var = datatree.AddVariable(symb_id, 0);
- neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(new_aux_var,
- origin_aux_var->decreaseLeadsLags(origin_arg_max_lag
- - rit->first))));
+ // just add equation of form: AUX_DIFF = ORIG_AUX_DIFF(last_arg_max_lag - rit->first)
+ VariableNode *new_aux_var;
+ for (int i = last_arg_max_lag; i > rit->first; i--)
+ {
+ if (vn == NULL)
+ symb_id = datatree.symbol_table.addDiffAuxiliaryVar(new_aux_var->idx, new_aux_var);
+ else
+ if (i == last_arg_max_lag)
+ symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst,
+ vn->get_symb_id(), i - 1);
+ else
+ symb_id = datatree.symbol_table.addDiffAuxiliaryVar(new_aux_var->idx, new_aux_var,
+ vn->get_symb_id(), i - 1);
+
+
+ new_aux_var = datatree.AddVariable(symb_id, 0);
+ neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(new_aux_var,
+ last_aux_var->decreaseLeadsLags(1))));
+ last_aux_var = new_aux_var;
+ }
subst_table[rit->second] = dynamic_cast<VariableNode *>(new_aux_var);
+ last_arg_max_lag = rit->first;
}
}
return const_cast<VariableNode *>(subst_table[this]);
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 052c148bc816515da77e7cb6f2e93454df0fcb42..1ea6dbb5a48020dd1078e70b9776a724862ec312 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -83,7 +83,8 @@ enum ExprNodeOutputType
oJuliaDynamicSteadyStateOperator, //!< Julia code, dynamic model, inside a steady state operator
oSteadyStateFile, //!< Matlab code, in the generated steady state file
oCSteadyStateFile, //!< C code, in the generated steady state file
- oJuliaSteadyStateFile //!< Julia code, in the generated steady state file
+ oJuliaSteadyStateFile, //!< Julia code, in the generated steady state file
+ oMatlabDseries //!< Matlab code for dseries
};
#define IS_MATLAB(output_type) ((output_type) == oMatlabStaticModel \
@@ -93,7 +94,8 @@ enum ExprNodeOutputType
|| (output_type) == oMatlabDynamicModelSparse \
|| (output_type) == oMatlabDynamicSteadyStateOperator \
|| (output_type) == oMatlabDynamicSparseSteadyStateOperator \
- || (output_type) == oSteadyStateFile)
+ || (output_type) == oSteadyStateFile \
+ || (output_type) == oMatlabDseries)
#define IS_JULIA(output_type) ((output_type) == oJuliaStaticModel \
|| (output_type) == oJuliaDynamicModel \
diff --git a/src/ModFile.cc b/src/ModFile.cc
index df36512cc309b254e981147daf5ee8b93757c7d9..b5fd81d87ba1cd0f32bf056b72a0c45d89c5abef 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -479,6 +479,8 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
dynamic_model.substituteEndoLagGreaterThanTwo(true);
}
+ dynamic_model.combineDiffAuxEquations();
+
if (differentiate_forward_vars)
dynamic_model.differentiateForwardVars(differentiate_forward_vars_subset);
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index be6b66cbbf7646be6c35e1fa1e2414eadf836bd8..312169a1f90bfc7d3af1d8d3f2f14f76cfe6f1ed 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -60,6 +60,12 @@ protected:
//! Only stores generated auxiliary equations, in an order meaningful for evaluation
deque<BinaryOpNode *> aux_equations;
+ //! Temporarily stores aux equations for diff operator
+ //! This is a hack to make diff aux vars work: they must be substituted before we consider VARs
+ //! But leads lags can't be substituted until after we consider VARs
+ //! The diff_aux_vars may depend on aux_vars created for leads lags, as in
+ //! diff(x(-1)) => x(-1) - x(-2) => x(-1) - aux_var(-1); aux_var = x(-1);
+ deque<BinaryOpNode *> diff_aux_equations;
//! Stores equation tags
vector<pair<int, pair<string, string> > > equation_tags;