1 /*
2 * Copyright © 2006 Ondra Kamenik
3 * Copyright © 2019 Dynare Team
4 *
5 * This file is part of Dynare.
6 *
7 * Dynare is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * Dynare is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 #include "planner_builder.hh"
22 #include "dynare_exception.hh"
23 #include "dynare_model.hh"
24
25 #include <cmath>
26 #include <utility>
27
28 using namespace ogdyn;
29
30 const IntegerMatrix &
operator =(const IntegerMatrix & im)31 IntegerMatrix::operator=(const IntegerMatrix &im)
32 {
33 if (nr != im.nr || nc != im.nc)
34 throw DynareException(__FILE__, __LINE__,
35 "Matrices have different dimensions in IntegerMatrix::operator=");
36 std::copy_n(im.data.get(), nr*nc, data.get());
37 return *this;
38 }
39
40 const IntegerArray3 &
operator =(const IntegerArray3 & ia3)41 IntegerArray3::operator=(const IntegerArray3 &ia3)
42 {
43 if (n1 != ia3.n1 || n2 != ia3.n2 || n3 != ia3.n3)
44 throw DynareException(__FILE__, __LINE__,
45 "Arrays have different dimensions in IntegerArray3::operator=");
46 std::copy_n(ia3.data.get(), n1*n2*n3, data.get());
47 return *this;
48 }
49
PlannerBuilder(DynareModel & m,const Tvarset & yyset,Teqset ffset)50 PlannerBuilder::PlannerBuilder(DynareModel &m, const Tvarset &yyset,
51 Teqset ffset)
52 : yset(), fset(std::move(ffset)), model(m),
53 tb(model.t_plobjective), tbeta(model.t_pldiscount),
54 maxlead(model.atoms.get_maxlead()),
55 minlag(model.atoms.get_minlag()),
56 diff_b(yyset.size(), 1-minlag),
57 diff_f(yyset.size(), fset.size(), 1+maxlead-minlag),
58 static_atoms(),
59 static_tree(),
60 diff_b_static(yyset.size(), 1-minlag),
61 diff_f_static(yyset.size(), fset.size(), 1+maxlead-minlag)
62 {
63 info.num_new_terms -= model.getParser().getTree().get_num_op();
64
65 fill_yset(m.atoms.get_name_storage(), yyset);
66
67 add_derivatives_of_b();
68 add_derivatives_of_f();
69 shift_derivatives_of_b();
70 shift_derivatives_of_f();
71 beta_multiply_b();
72 beta_multiply_f();
73 make_static_version();
74 lagrange_mult_f();
75 form_equations();
76
77 info.num_new_terms += model.getParser().getTree().get_num_op();
78 }
79
PlannerBuilder(const PlannerBuilder & pb,ogdyn::DynareModel & m)80 PlannerBuilder::PlannerBuilder(const PlannerBuilder &pb, ogdyn::DynareModel &m)
81 : yset(), fset(pb.fset), model(m),
82 tb(pb.tb), tbeta(pb.tbeta),
83 maxlead(pb.maxlead), minlag(pb.minlag),
84 diff_b(pb.diff_b), diff_f(pb.diff_f),
85 static_atoms(pb.static_atoms),
86 static_tree(pb.static_tree),
87 diff_b_static(pb.diff_b_static),
88 diff_f_static(pb.diff_f_static),
89 aux_map(), static_aux_map()
90 {
91 fill_yset(m.atoms.get_name_storage(), pb.yset);
92 fill_aux_map(m.atoms.get_name_storage(), pb.aux_map, pb.static_aux_map);
93 }
94
95 void
add_derivatives_of_b()96 PlannerBuilder::add_derivatives_of_b()
97 {
98 int yi = 0;
99 for (auto yname = yset.begin(); yname != yset.end(); ++yname, yi++)
100 for (int ll = minlag; ll <= 0; ll++)
101 {
102 int yt = model.atoms.index(*yname, ll);
103 if (yt != -1)
104 diff_b(yi, ll-minlag) = model.eqs.add_derivative(tb, yt);
105 else
106 diff_b(yi, ll-minlag) = ogp::OperationTree::zero;
107 }
108 }
109
110 void
add_derivatives_of_f()111 PlannerBuilder::add_derivatives_of_f()
112 {
113 int yi = 0;
114 for (auto yname = yset.begin(); yname != yset.end(); ++yname, yi++)
115 for (unsigned int fi = 0; fi < fset.size(); fi++)
116 for (int ll = minlag; ll <= maxlead; ll++)
117 {
118 int yt = model.atoms.index(*yname, ll);
119 if (yt != -1)
120 diff_f(yi, fi, ll-minlag)
121 = model.eqs.add_derivative(model.eqs.formula(fset[fi]), yt);
122 else
123 diff_f(yi, fi, ll-minlag) = ogp::OperationTree::zero;
124 }
125 }
126
127 void
shift_derivatives_of_b()128 PlannerBuilder::shift_derivatives_of_b()
129 {
130 map<int, int> subst;
131 for (int yi = 0; yi < diff_b.nrows(); yi++)
132 for (int ll = minlag; ll < 0; ll++)
133 if (diff_b(yi, ll-minlag) != ogp::OperationTree::zero)
134 {
135 model.variable_shift_map(model.eqs.nulary_of_term(diff_b(yi, ll-minlag)),
136 -ll, subst);
137 diff_b(yi, ll-minlag) = model.eqs.add_substitution(diff_b(yi, ll-minlag), subst);
138 }
139 }
140
141 void
shift_derivatives_of_f()142 PlannerBuilder::shift_derivatives_of_f()
143 {
144 map<int, int> subst;
145 for (int yi = 0; yi < diff_f.dim1(); yi++)
146 for (int fi = 0; fi < diff_f.dim2(); fi++)
147 {
148 // first do it leads which are put under expectation before t: no problem
149 for (int ll = 0; ll <= maxlead; ll++)
150 if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
151 {
152 model.variable_shift_map(model.eqs.nulary_of_term(diff_f(yi, fi, ll-minlag)),
153 -ll, subst);
154 diff_f(yi, fi, ll-minlag)
155 = model.eqs.add_substitution(diff_f(yi, fi, ll-minlag), subst);
156 }
157 /* now do it for lags, these are put as leads under expectations after
158 time t, so we have to introduce auxiliary variables at time t, and
159 make leads of them here */
160 for (int ll = minlag; ll < 0; ll++)
161 {
162 int ft = diff_f(yi, fi, ll-minlag);
163 if (ft != ogp::OperationTree::zero)
164 {
165 /* if the ft term has a lead, than we need to introduce an
166 auxiliary variable zₜ, define it as ₜ[ft] and put z_{t-ll}
167 to the equation. Otherwise, we just put leaded ft to the
168 equation directly. */
169 int ft_maxlead, ft_minlag;
170 model.termspan(ft, ft_maxlead, ft_minlag);
171 if (ft_maxlead > 0)
172 {
173 // make an auxiliary variable
174 std::string name;
175 name = "AUX_" + std::to_string(yi) + '_' + std::to_string(fset[fi]) + '_' + std::to_string(-ll);
176 model.atoms.register_uniq_endo(name);
177 info.num_aux_variables++;
178 int taux = model.eqs.add_nulary(name);
179 name = "AUX_" + std::to_string(yi) + '_' + std::to_string(fset[fi]) + '_' + std::to_string(-ll) + '(' + std::to_string(-ll) + ')';
180 int taux_leaded = model.eqs.add_nulary(name);
181 // put aux_leaded to the equation
182 diff_f(yi, fi, ll-minlag) = taux_leaded;
183 // save auxiliary variable and the term
184 aux_map.emplace(model.atoms.name(taux), ft);
185 }
186 else
187 {
188 /* no auxiliary variable is needed and the term ft can be
189 leaded in place */
190 model.variable_shift_map(model.eqs.nulary_of_term(ft), -ll, subst);
191 diff_f(yi, fi, ll-minlag) = model.eqs.add_substitution(ft, subst);
192 }
193 }
194 }
195 }
196 }
197
198 void
beta_multiply_b()199 PlannerBuilder::beta_multiply_b()
200 {
201 int beta_pow = ogp::OperationTree::one;
202 for (int ll = 0; ll >= minlag; ll--,
203 beta_pow = model.eqs.add_binary(ogp::code_t::TIMES, beta_pow, tbeta))
204 for (int yi = 0; yi < diff_b.nrows(); yi++)
205 if (diff_b(yi, ll-minlag) != ogp::OperationTree::zero)
206 diff_b(yi, ll-minlag)
207 = model.eqs.add_binary(ogp::code_t::TIMES, beta_pow, diff_b(yi, ll-minlag));
208 }
209
210 void
beta_multiply_f()211 PlannerBuilder::beta_multiply_f()
212 {
213 int beta_pow = ogp::OperationTree::one;
214 for (int ll = 0; ll <= maxlead; ll++,
215 beta_pow = model.eqs.add_binary(ogp::code_t::DIVIDE, beta_pow, tbeta))
216 for (int yi = 0; yi < diff_f.dim1(); yi++)
217 for (int fi = 0; fi < diff_f.dim2(); fi++)
218 if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
219 diff_f(yi, fi, ll-minlag)
220 = model.eqs.add_binary(ogp::code_t::TIMES, beta_pow, diff_f(yi, fi, ll-minlag));
221
222 beta_pow = ogp::OperationTree::one;
223 for (int ll = 0; ll >= minlag; ll--,
224 beta_pow = model.eqs.add_binary(ogp::code_t::TIMES, beta_pow, tbeta))
225 for (int yi = 0; yi < diff_f.dim1(); yi++)
226 for (int fi = 0; fi < diff_f.dim2(); fi++)
227 if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
228 diff_f(yi, fi, ll-minlag)
229 = model.eqs.add_binary(ogp::code_t::TIMES, beta_pow, diff_f(yi, fi, ll-minlag));
230 }
231
232 void
make_static_version()233 PlannerBuilder::make_static_version()
234 {
235 // map holding substitutions from dynamic to static
236 ogp::StaticFineAtoms::Tintintmap tmap;
237
238 // fill static atoms with outer ordering
239 static_atoms.import_atoms(model.atoms, static_tree, tmap);
240
241 // go through diff_b and fill diff_b_static
242 for (int ll = minlag; ll <= 0; ll++)
243 for (int yi = 0; yi < diff_b.nrows(); yi++)
244 diff_b_static(yi, ll-minlag)
245 = static_tree.add_substitution(diff_b(yi, ll-minlag),
246 tmap, model.eqs.getTree());
247
248 // go through diff_f and fill diff_f_static
249 for (int ll = minlag; ll <= maxlead; ll++)
250 for (int yi = 0; yi < diff_f.dim1(); yi++)
251 for (int fi = 0; fi < diff_f.dim2(); fi++)
252 diff_f_static(yi, fi, ll-minlag)
253 = static_tree.add_substitution(diff_f(yi, fi, ll-minlag),
254 tmap, model.eqs.getTree());
255
256 // go through aux_map and fill static_aux_map
257 for (const auto &it : aux_map)
258 {
259 int tstatic = static_tree.add_substitution(it.second, tmap, model.eqs.getTree());
260 static_aux_map.emplace(it.first, tstatic);
261 }
262 }
263
264 void
lagrange_mult_f()265 PlannerBuilder::lagrange_mult_f()
266 {
267 // register multipliers
268 std::string mult_name;
269 for (int fi = 0; fi < diff_f.dim2(); fi++)
270 {
271 mult_name = "MULT" + std::to_string(fset[fi]);
272 model.atoms.register_uniq_endo(mult_name);
273 info.num_lagrange_mults++;
274 }
275 // multiply with the multipliers
276 for (int yi = 0; yi < diff_f.dim1(); yi++)
277 for (int fi = 0; fi < diff_f.dim2(); fi++)
278 for (int ll = minlag; ll <= maxlead; ll++)
279 if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
280 {
281 mult_name = "MULT" + std::to_string(fset[fi]) + '(' + std::to_string(-ll) + ')';
282 int tm = model.eqs.add_nulary(mult_name);
283 diff_f(yi, fi, ll-minlag)
284 = model.eqs.add_binary(ogp::code_t::TIMES, tm, diff_f(yi, fi, ll-minlag));
285 }
286 }
287
288 void
form_equations()289 PlannerBuilder::form_equations()
290 {
291 // add planner’s FOCs
292 for (int yi = 0; yi < diff_f.dim1(); yi++)
293 {
294 int eq = ogp::OperationTree::zero;
295 for (int ll = minlag; ll <= 0; ll++)
296 eq = model.eqs.add_binary(ogp::code_t::PLUS, eq, diff_b(yi, ll-minlag));
297 for (int fi = 0; fi < diff_f.dim2(); fi++)
298 for (int ll = minlag; ll <= maxlead; ll++)
299 eq = model.eqs.add_binary(ogp::code_t::PLUS, eq, diff_f(yi, fi, ll-minlag));
300 model.eqs.add_formula(eq);
301 }
302
303 // add equations for auxiliary variables
304 for (const auto &it : aux_map)
305 {
306 int t = model.atoms.index(it.first, 0);
307 model.eqs.add_formula(model.eqs.add_binary(ogp::code_t::MINUS, t, it.second));
308 }
309 }
310
311 void
fill_yset(const ogp::NameStorage & ns,const PlannerBuilder::Tvarset & yyset)312 PlannerBuilder::fill_yset(const ogp::NameStorage &ns,
313 const PlannerBuilder::Tvarset &yyset)
314 {
315 for (auto it : yyset)
316 yset.insert(it);
317 }
318
319 void
fill_aux_map(const ogp::NameStorage & ns,const Tsubstmap & aaux_map,const Tsubstmap & astatic_aux_map)320 PlannerBuilder::fill_aux_map(const ogp::NameStorage &ns, const Tsubstmap &aaux_map,
321 const Tsubstmap &astatic_aux_map)
322 {
323 // fill aux_map
324 for (auto it : aaux_map)
325 aux_map.insert(it);
326
327 // fill static_aux_map
328 for (auto it : astatic_aux_map)
329 static_aux_map.insert(it);
330 }
331
MultInitSS(const PlannerBuilder & pb,const Vector & pvals,Vector & yy)332 MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy)
333 : builder(pb), b(builder.diff_b_static.nrows()),
334 F(builder.diff_f_static.dim1(), builder.diff_f_static.dim2())
335 {
336 b.zeros();
337 F.zeros();
338
339 // first evaluate substitutions (auxiliary variables) from the builder
340 ogdyn::DynareStaticSteadySubstitutions dss(builder.model.atoms, builder.static_atoms,
341 builder.static_tree,
342 builder.static_aux_map, pvals, yy);
343
344 /* gather all the terms from builder.diff_b_static and builder.diff_f_static
345 to the vector, the ordering is important, since the index of this vector
346 will have to be decoded to the position in b and F. */
347 vector<int> terms;
348 for (int yi = 0; yi < builder.diff_b_static.nrows(); yi++)
349 for (int l = 0; l < builder.diff_b_static.ncols(); l++)
350 terms.push_back(builder.diff_b_static(yi, l));
351 for (int yi = 0; yi < builder.diff_f_static.dim1(); yi++)
352 for (int fi = 0; fi < builder.diff_f_static.dim2(); fi++)
353 for (int l = 0; l < builder.diff_f_static.dim3(); l++)
354 terms.push_back(builder.diff_f_static(yi, fi, l));
355
356 /* evaluate the terms, it will call a series of load(i,res), which sum the
357 results through lags/leads to b and F */
358 DynareStaticSteadyAtomValues dssav(builder.model.atoms, builder.static_atoms, pvals, yy);
359 ogp::FormulaCustomEvaluator fe(builder.static_tree, terms);
360 fe.eval(dssav, *this);
361
362 // solve overdetermined system b+F*lambda=0 using SVD decomposition
363 SVDDecomp decomp(F);
364 Vector lambda(builder.diff_f_static.dim2());
365 decomp.solve(b, lambda);
366 lambda.mult(-1);
367
368 // take values of lambda and put it to yy
369 for (int fi = 0; fi < builder.diff_f_static.dim2(); fi++)
370 {
371 std::string mult_name = "MULT" + std::to_string(builder.fset[fi]);
372 int iouter = builder.model.atoms.name2outer_endo(mult_name);
373 int iy = builder.model.atoms.outer2y_endo()[iouter];
374 if (!std::isfinite(yy[iy]))
375 yy[iy] = lambda[fi];
376
377 /* go through all substitutions of the multiplier and set them as well */
378 if (builder.model.atom_substs)
379 {
380 const ogp::AtomSubstitutions::Toldnamemap &old2new
381 = builder.model.atom_substs->get_old2new();
382 auto it = old2new.find(mult_name);
383 if (it != old2new.end())
384 {
385 const ogp::AtomSubstitutions::Tshiftnameset &sset = it->second;
386 for (const auto &itt : sset)
387 {
388 const std::string &newname = itt.first;
389 int iouter = builder.model.atoms.name2outer_endo(newname);
390 int iy = builder.model.atoms.outer2y_endo()[iouter];
391 if (!std::isfinite(yy[iy]))
392 yy[iy] = lambda[fi];
393 }
394 }
395 }
396 }
397 }
398
399 void
load(int i,double res)400 MultInitSS::load(int i, double res)
401 {
402 /* we can afford it, since the evaluator sets res to exact zero if the term
403 is zero */
404 if (res == 0)
405 return;
406 // decode i and add to either b or F
407 if (i < builder.diff_b_static.nrows()*builder.diff_b_static.ncols())
408 // add to b
409 b[i / builder.diff_b_static.ncols()] += res;
410 else
411 {
412 // add to F
413 i -= builder.diff_b_static.nrows()*builder.diff_b_static.ncols();
414 int yifi = i / builder.diff_f_static.dim3();
415 int yi = yifi / builder.diff_f_static.dim2();
416 int fi = yifi % builder.diff_f_static.dim2();
417 F.get(yi, fi) += res;
418 }
419 }
420