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