1 /*
2  * Copyright © 2006-2011 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 #ifndef PLANNER_BUILDER_H
22 #define PLANNER_BUILDER_H
23 
24 #include <unordered_set>
25 #include <map>
26 #include <vector>
27 #include <memory>
28 #include <algorithm>
29 
30 #include "parser/cc/static_fine_atoms.hh"
31 #include "dynare_atoms.hh"
32 #include "GeneralMatrix.hh"
33 
34 namespace ogdyn
35 {
36   using std::unordered_set;
37   using std::map;
38   using std::vector;
39 
40   /** This is a two dimensional array of integers. Nothing
41    * difficult. */
42   class IntegerMatrix
43   {
44   protected:
45     /** Number of rows. */
46     int nr;
47     /** Number of columns. */
48     int nc;
49     /** The pointer to the data. */
50     std::unique_ptr<int[]> data;
51   public:
52     /** Construct uninitialized array. */
IntegerMatrix(int nrr,int ncc)53     IntegerMatrix(int nrr, int ncc)
54       : nr(nrr), nc(ncc), data(std::make_unique<int[]>(nr*nc))
55     {
56     }
57     /** Copy constructor. */
IntegerMatrix(const IntegerMatrix & im)58     IntegerMatrix(const IntegerMatrix &im)
59       : nr(im.nr), nc(im.nc), data(std::make_unique<int[]>(nr*nc))
60     {
61       std::copy_n(im.data.get(), nr*nc, data.get());
62     }
63     /** Assignment operator. It can only assing array with the
64      * same dimensions. */
65     const IntegerMatrix &operator=(const IntegerMatrix &im);
66     int &
operator ()(int i,int j)67     operator()(int i, int j)
68     {
69       return data[i+j*nr];
70     }
71     const int &
operator ()(int i,int j) const72     operator()(int i, int j) const
73     {
74       return data[i+j*nr];
75     }
76     int
nrows() const77     nrows() const
78     {
79       return nr;
80     }
81     int
ncols() const82     ncols() const
83     {
84       return nc;
85     }
86   };
87 
88   /** The three dimensional array of integers. Nothing difficult. */
89   class IntegerArray3
90   {
91   protected:
92     /** First dimension. */
93     int n1;
94     /** Second dimension. */
95     int n2;
96     /** Third dimension. */
97     int n3;
98     /** The data. */
99     std::unique_ptr<int[]> data;
100   public:
101     /** Constrcut unitialized array. */
IntegerArray3(int nn1,int nn2,int nn3)102     IntegerArray3(int nn1, int nn2, int nn3)
103       : n1(nn1), n2(nn2), n3(nn3), data(std::make_unique<int[]>(n1*n2*n3))
104     {
105     }
106     /** Copy constructor. */
IntegerArray3(const IntegerArray3 & ia3)107     IntegerArray3(const IntegerArray3 &ia3)
108       : n1(ia3.n1), n2(ia3.n2), n3(ia3.n3), data(std::make_unique<int[]>(n1*n2*n3))
109     {
110       std::copy_n(ia3.data.get(), n1*n2*n3, data.get());
111     }
112     /** Assignment operator assigning the arrays with the same dimensions. */
113     const IntegerArray3 &operator=(const IntegerArray3 &ia3);
114     int &
operator ()(int i,int j,int k)115     operator()(int i, int j, int k)
116     {
117       return data[i+j*n1+k*n1*n2];
118     }
119     const int &
operator ()(int i,int j,int k) const120     operator()(int i, int j, int k) const
121     {
122       return data[i+j*n1+k*n1*n2];
123     }
124     int
dim1() const125     dim1() const
126     {
127       return n1;
128     }
129     int
dim2() const130     dim2() const
131     {
132       return n2;
133     }
134     int
dim3() const135     dim3() const
136     {
137       return n3;
138     }
139   };
140 
141   /** This struct encapsulates information about the building of a
142    * planner's problem. */
143   struct PlannerInfo
144   {
145     int num_lagrange_mults{0};
146     int num_aux_variables{0};
147     int num_new_terms{0};
148   };
149 
150   class MultInitSS;
151   class DynareModel;
152 
153   /** This class builds the first order conditions of the social
154    * planner problem with constraints being the equations in the
155    * model. The model is non-const parameter to the constructor
156    * which adds appropriate FOCs to the system. It also allows for
157    * an estimation of the lagrange multipliers given all other
158    * endogenous variables of the static system. For this purpose we
159    * need to create static atoms and static versions of all the tree
160    * index matrices. The algorithm and algebra are documented in
161    * dynare++-ramsey.pdf. */
162   class PlannerBuilder
163   {
164     friend class MultInitSS;
165   public:
166     /** Type for a set of variable names. */
167     using Tvarset = unordered_set<string>;
168     /** Type for a set of equations. An equation is identified by
169      * an index to an equation in the equation vector given by
170      * DynareModel::eqs. The tree index of the i-th formula is
171      * retrieved as DynareModel::egs.formula(i). */
172     using Teqset = vector<int>;
173   protected:
174     /** This is a set of variables wrt which the planner
175      * optimizes. These could be all endogenous variables, but it
176      * is beneficial to exclude all variables which are
177      * deterministic transformations of past exogenous variables,
178      * since the planner cannot influence them. This could save a
179      * few equations. This is not changed after it is constructed,
180      * but it is constructed manually, so it cannot be declared as
181      * const. */
182     Tvarset yset;
183     /** These are the equation indices constituing the constraints
184      * for the planner. Again, it is beneficial to exclude all
185      * equations defining exogenous variables excluded from
186      * yset. */
187     const Teqset fset;
188     /** Reference to the model. */
189     ogdyn::DynareModel &model;
190     /** Tree index of the planner objective. */
191     int tb;
192     /** Tree index of the planner discount parameter. */
193     int tbeta;
194     /** The maximum lead in the model including the planner's
195      * objective before building the planner's FOCs. */
196     const int maxlead;
197     /** The minimum lag in the model including the planner's objective
198      * before building the planner's FOCs. */
199     const int minlag;
200     /** Tree indices of formulas in the planner FOCs involving
201      * derivatives of the planner's objective. Rows correspond to the
202      * endogenous variables, columns correspond to lags in the
203      * objective function. The contents of the matrix will evolve as
204      * the algorithm proceeds. */
205     IntegerMatrix diff_b;
206     /** Tree indices of formulas in the planner FOCs involving
207      * derivatives of the model equations (constraints). The first
208      * dimension corresponds to endogenous variables, the second to
209      * the constraints, the third to lags or leads of endogenous
210      * variables in the constraints. The contents of the array will
211      * evolve as the algorithm proceeds.*/
212     IntegerArray3 diff_f;
213     /** Static version of the model atoms. It is needed to build
214      * static version of diff_b and diff_f. */
215     ogp::StaticFineAtoms static_atoms;
216     /** Static version of all the trees of diff_b and diff_f build
217      * over static_atoms. */
218     ogp::OperationTree static_tree;
219     /** Tree indices of static version of diff_b over static_atoms and static_tree. */
220     IntegerMatrix diff_b_static;
221     /** Tree indices of static version of diff_f over static_atoms
222      * and static_tree. This member is created before calling
223      * lagrange_mult_f(), so it does not contain the
224      * multiplication with the lagrange multipliers. */
225     IntegerArray3 diff_f_static;
226     /** Auxiliary variables mapping. During the algorithm, some
227      * auxiliary variables for the terms might be created, so we
228      * remember their names and tree indices of the terms. This
229      * maps a name to the tree index of an expression equal to the
230      * auxiliary variable at time zero. The auxiliary variables
231      * names point to the dynamic atoms storage, tree inidices to
232      * the dynamic model tree. */
233     Tsubstmap aux_map;
234     /** Static version of aux_map. The names point to static_atoms
235      * storage, the tree indices to the static_tree. */
236     Tsubstmap static_aux_map;
237     /** Information about the number of various things. */
238     PlannerInfo info;
239   public:
240     /** Build the planner problem for the given model optimizing
241      * through the given endogenous variables with the given
242      * constraints. We allow for a selection of a subset of
243      * equations and variables in order to eliminate exogenous
244      * predetermined process which cannot be influenced by the
245      * social planner. */
246     PlannerBuilder(ogdyn::DynareModel &m, const Tvarset &yyset,
247                    Teqset ffset);
248     /** Construct a copy of the builder with provided model, which
249      * is supposed to be the copy of the model in the builder. */
250     PlannerBuilder(const PlannerBuilder &pb, ogdyn::DynareModel &m);
251     /** Avoid copying from only PlannerBuilder. */
252     PlannerBuilder(const PlannerBuilder &pb) = delete;
253     /** Return the information. */
254     const PlannerInfo &
get_info() const255     get_info() const
256     {
257       return info;
258     }
259   protected:
260     /** Differentiate the planner objective wrt endogenous
261      * variables with different lags. */
262     void add_derivatives_of_b();
263     /** Differentiate the constraints wrt endogenous variables
264      * with different lags and leads. */
265     void add_derivatives_of_f();
266     /** Shift derivatives of diff_b. */
267     void shift_derivatives_of_b();
268     /** Shift derivatives of diff_ff. */
269     void shift_derivatives_of_f();
270     /** Multiply with the discount factor terms in diff_b. */
271     void beta_multiply_b();
272     /** Multiply with the discount factor terms in diff_f. */
273     void beta_multiply_f();
274     /** Fill static_atoms and static_tree and build diff_b_static,
275      * diff_f_static and aux_map_static with static versions of diff_b,
276      * diff_f and aux_map. */
277     void make_static_version();
278     /** Multiply diff_f with Langrange multipliers. */
279     void lagrange_mult_f();
280     /** Add the equations to the mode, including equation for auxiliary variables. */
281     void form_equations();
282   private:
283     /** Fill yset for a given yyset and given name storage. */
284     void fill_yset(const ogp::NameStorage &ns, const Tvarset &yyset);
285     /** Fill aux_map and aux_map_static for a given aaux_map and
286      * aaux_map_static for a given storage of dynamic atoms (used
287      * for aux_map) and static atoms storage from this object for
288      * aux_map_static. */
289     void fill_aux_map(const ogp::NameStorage &ns, const Tsubstmap &aaux_map,
290                       const Tsubstmap &astatic_aux_map);
291   };
292 
293   /** This class only calculates for the given initial guess of
294    * endogenous variables, initial guess of the Langrange
295    * multipliers of the social planner problem yielding the least
296    * square error. It is used by just calling its constructor. The
297    * constructor takes non-const reference to the vector of
298    * endogenous variables, calculates lambdas and put the values of
299    * lambdas to the vector. The algbera is found in
300    * dynare++-ramsey.pdf.
301    *
302    * The code can be run only after the parsing has been finished in
303    * atoms. */
304   class MultInitSS : public ogp::FormulaEvalLoader
305   {
306   protected:
307     /** The constant reference to the builder. */
308     const PlannerBuilder &builder;
309     /** The constant term of the problem. Its length is the number
310      * of endogenous variable wrt the planner optimizes. */
311     Vector b;
312     /** The matrix of the overdetermined problem. The number of
313      * rows is equal to the number of endogenous variables wrt
314      * which the planner optimizes, the number of columns is equal
315      * to the number of Langrange multipliers which is equal to
316      * the number of constraints which is smaller than the number
317      * of endogenous variables. Hence the system b+F*lambda=0 is
318      * overdetermined. */
319     GeneralMatrix F;
320   public:
321     /** The constructor of the object which does everything. Its
322      * main goal is to update yy. Note that if an item of yy
323      * corresponding to a lagrange multiplier is already set, it
324      * is not reset. */
325     MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy);
326     /** This loads evaluated parts of b or F and decodes i and
327      * advances b or F depending on the decoded i. The decoding is
328      * dependent on the way how the terms of builder.diff_b and
329      * builder.diff_f_save have been put the the
330      * ogp::FormulaCustomEvaluator. This is documented in the code
331      * of the constructor. */
332     void load(int i, double res) override;
333   };
334 };
335 
336 #endif
337