NumericalInitialization.cc 9.43 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
sebastien's avatar
sebastien committed
2
 * Copyright (C) 2003-2009 Dynare Team
sebastien's avatar
sebastien committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * This file is part of Dynare.
 *
 * Dynare is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Dynare is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 */

20
21
#include <iostream>
#include <fstream>
22
#include <cstdlib>
23

24
#include "NumericalInitialization.hh"
sebastien's avatar
sebastien committed
25
26

InitParamStatement::InitParamStatement(int symb_id_arg,
sebastien's avatar
sebastien committed
27
                                       const NodeID param_value_arg,
sebastien's avatar
sebastien committed
28
                                       const SymbolTable &symbol_table_arg) :
sebastien's avatar
sebastien committed
29
  symb_id(symb_id_arg),
sebastien's avatar
sebastien committed
30
31
  param_value(param_value_arg),
  symbol_table(symbol_table_arg)
32
33
34
{
}

sebastien's avatar
sebastien committed
35
void
sebastien's avatar
sebastien committed
36
InitParamStatement::writeOutput(ostream &output, const string &basename) const
37
{
sebastien's avatar
sebastien committed
38
  int id = symbol_table.getTypeSpecificID(symb_id) + 1;
sebastien's avatar
sebastien committed
39
40
41
  output << "M_.params( " << id << " ) = ";
  param_value->writeOutput(output);
  output << ";" << endl;
sebastien's avatar
sebastien committed
42
  output << symbol_table.getName(symb_id) << " = M_.params( " << id << " );\n";
43
44
}

sebastien's avatar
sebastien committed
45
46
void
InitParamStatement::fillEvalContext(eval_context_type &eval_context) const
47
{
sebastien's avatar
sebastien committed
48
49
50
51
52
53
54
55
  try
    {
      eval_context[symb_id] = param_value->eval(eval_context);
    }
  catch(ExprNode::EvalException &e)
    {
      // Do nothing
    }
56
57
}

sebastien's avatar
sebastien committed
58
59
60
61
InitOrEndValStatement::InitOrEndValStatement(const init_values_type &init_values_arg,
                                             const SymbolTable &symbol_table_arg) :
  init_values(init_values_arg),
  symbol_table(symbol_table_arg)
62
63
64
{
}

sebastien's avatar
sebastien committed
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
void
InitOrEndValStatement::fillEvalContext(eval_context_type &eval_context) const
{
  for(init_values_type::const_iterator it = init_values.begin();
      it != init_values.end(); it++)
    {
      try
        {
          eval_context[it->first] = (it->second)->eval(eval_context);
        }
      catch(ExprNode::EvalException &e)
        {
          // Do nothing
        }
    }
}

sebastien's avatar
sebastien committed
82
83
void
InitOrEndValStatement::writeInitValues(ostream &output) const
84
{
sebastien's avatar
sebastien committed
85
86
  for(init_values_type::const_iterator it = init_values.begin();
      it != init_values.end(); it++)
87
    {
sebastien's avatar
sebastien committed
88
      const int symb_id = it->first;
sebastien's avatar
sebastien committed
89
      const NodeID expression = it->second;
sebastien's avatar
sebastien committed
90

sebastien's avatar
sebastien committed
91
92
      SymbolType type = symbol_table.getType(symb_id);
      int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
sebastien's avatar
sebastien committed
93
94

      if (type == eEndogenous)
sebastien's avatar
sebastien committed
95
        output << "oo_.steady_state";
sebastien's avatar
sebastien committed
96
      else if (type == eExogenous)
sebastien's avatar
sebastien committed
97
        output << "oo_.exo_steady_state";
sebastien's avatar
sebastien committed
98
      else if (type == eExogenousDet)
sebastien's avatar
sebastien committed
99
100
        output << "oo_.exo_det_steady_state";

sebastien's avatar
sebastien committed
101
      output << "( " << tsid << " ) = ";
sebastien's avatar
sebastien committed
102
103
      expression->writeOutput(output);
      output << ";" << endl;
104
105
106
    }
}

sebastien's avatar
sebastien committed
107
InitValStatement::InitValStatement(const init_values_type &init_values_arg,
108
109
                                   const SymbolTable &symbol_table_arg) :
  InitOrEndValStatement(init_values_arg, symbol_table_arg)
110
{
sebastien's avatar
sebastien committed
111
}
112

sebastien's avatar
sebastien committed
113
void
sebastien's avatar
sebastien committed
114
InitValStatement::writeOutput(ostream &output, const string &basename) const
sebastien's avatar
sebastien committed
115
{
sebastien's avatar
sebastien committed
116
117
118
  output << "%" << endl
         << "% INITVAL instructions" << endl
         << "%" << endl;
119
  // Writing initval block to set initial values for variables
sebastien's avatar
trunk:    
sebastien committed
120
  output << "options_.initval_file = 0;" << endl;
121

sebastien's avatar
sebastien committed
122
  writeInitValues(output);
123

124
  output << "oo_.endo_simul=[oo_.steady_state*ones(1,M_.maximum_lag)];\n";
sebastien's avatar
sebastien committed
125
126
127
128
129
130
  output << "if M_.exo_nbr > 0;\n";
  output << "\too_.exo_simul = [ones(M_.maximum_lag,1)*oo_.exo_steady_state'];\n";
  output <<"end;\n";
  output << "if M_.exo_det_nbr > 0;\n";
  output << "\too_.exo_det_simul = [ones(M_.maximum_lag,1)*oo_.exo_det_steady_state'];\n";
  output <<"end;\n";
131
132
}

sebastien's avatar
sebastien committed
133
134
135
136

EndValStatement::EndValStatement(const init_values_type &init_values_arg,
                                 const SymbolTable &symbol_table_arg) :
  InitOrEndValStatement(init_values_arg, symbol_table_arg)
137
138
139
{
}

sebastien's avatar
sebastien committed
140
141

void
sebastien's avatar
sebastien committed
142
EndValStatement::writeOutput(ostream &output, const string &basename) const
143
{
sebastien's avatar
sebastien committed
144
145
146
  output << "%" << endl
         << "% ENDVAL instructions" << endl
         << "%" << endl;
147
  // Writing endval block to set terminal values for variables
sebastien's avatar
sebastien committed
148
  output << "ys0_= oo_.steady_state;" << endl
sebastien's avatar
trunk:    
sebastien committed
149
         << "ex0_ = oo_.exo_steady_state;" << endl;
150

sebastien's avatar
sebastien committed
151
  writeInitValues(output);
152
153
}

sebastien's avatar
sebastien committed
154
155
156
157
HistValStatement::HistValStatement(const hist_values_type &hist_values_arg,
                                   const SymbolTable &symbol_table_arg) :
  hist_values(hist_values_arg),
  symbol_table(symbol_table_arg)
158
159
160
{
}

sebastien's avatar
sebastien committed
161
void
sebastien's avatar
sebastien committed
162
HistValStatement::writeOutput(ostream &output, const string &basename) const
163
{
sebastien's avatar
sebastien committed
164
165
166
  output << "%" << endl
         << "% HISTVAL instructions" << endl
         << "%" << endl;
167

sebastien's avatar
sebastien committed
168
169
  for(hist_values_type::const_iterator it = hist_values.begin();
      it != hist_values.end(); it++)
170
    {
sebastien's avatar
sebastien committed
171
      const int &symb_id = it->first.first;
sebastien's avatar
sebastien committed
172
      const int &lag = it->first.second;
sebastien's avatar
sebastien committed
173
      const NodeID expression = it->second;
sebastien's avatar
sebastien committed
174

sebastien's avatar
sebastien committed
175
176
      SymbolType type = symbol_table.getType(symb_id);
      int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
sebastien's avatar
sebastien committed
177
178

      if (type == eEndogenous)
sebastien's avatar
sebastien committed
179
        output << "oo_.endo_simul( " << tsid << ", M_.maximum_lag + " << lag << ") = ";
sebastien's avatar
sebastien committed
180
      else if (type == eExogenous)
sebastien's avatar
sebastien committed
181
        output << "oo_.exo_simul( M_.maximum_lag + " << lag << ", " << tsid << " ) = ";
sebastien's avatar
sebastien committed
182
      else if (type != eExogenousDet)
sebastien's avatar
sebastien committed
183
        output << "oo_.exo_det_simul( M_.maximum_lag + " << lag  << ", " << tsid << " ) = ";
sebastien's avatar
sebastien committed
184
185
186

      expression->writeOutput(output);
      output << ";" << endl;
187
188
    }
}
189

sebastien's avatar
sebastien committed
190
191
InitvalFileStatement::InitvalFileStatement(const string &filename_arg) :
  filename(filename_arg)
michel's avatar
michel committed
192
193
194
{
}

sebastien's avatar
sebastien committed
195
196
void
InitvalFileStatement::writeOutput(ostream &output, const string &basename) const
michel's avatar
michel committed
197
{
sebastien's avatar
sebastien committed
198
199
200
201
202
  output << "%" << endl
         << "% INITVAL_FILE statement" << endl
         << "%" << endl
         << "options_.initval_file = 1;" << endl
         << "initvalf('" << filename << "');" << endl;
michel's avatar
michel committed
203
204
}

205
HomotopyStatement::HomotopyStatement(const homotopy_values_type &homotopy_values_arg,
206
                                     const SymbolTable &symbol_table_arg) :
207
208
209
210
211
212
213
214
  homotopy_values(homotopy_values_arg),
  symbol_table(symbol_table_arg)
{
}

void
HomotopyStatement::writeOutput(ostream &output, const string &basename) const
{
sebastien's avatar
sebastien committed
215
216
217
  output << "%" << endl
         << "% HOMOTOPY_SETUP instructions" << endl
         << "%" << endl
218
         << "options_.homotopy_values = [];" << endl;
219
220
221
222

  for(homotopy_values_type::const_iterator it = homotopy_values.begin();
      it != homotopy_values.end(); it++)
    {
sebastien's avatar
sebastien committed
223
      const int &symb_id = it->first;
224
225
226
      const NodeID expression1 = it->second.first;
      const NodeID expression2 = it->second.second;

sebastien's avatar
sebastien committed
227
228
      const SymbolType type = symbol_table.getType(symb_id);
      const int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
229

sebastien's avatar
sebastien committed
230
      output << "options_.homotopy_values = vertcat(options_.homotopy_values, [ " << type << ", " << tsid << ", ";
231
232
233
234
      if (expression1 != NULL)
        expression1->writeOutput(output);
      else
        output << "NaN";
235
236
      output << ", ";
      expression2->writeOutput(output);
237
      output << "]);" << endl;
238
239
    }
}
sebastien's avatar
sebastien committed
240
241
242
243
244
245
246
247
248
249
250
251

SaveParamsAndSteadyStateStatement::SaveParamsAndSteadyStateStatement(const string &filename_arg) :
  filename(filename_arg)
{
}

void
SaveParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename) const
{
  output << "save_params_and_steady_state('" << filename << "');" << endl;
}

252
LoadParamsAndSteadyStateStatement::LoadParamsAndSteadyStateStatement(const string &filename,
sebastien's avatar
sebastien committed
253
                                                                     const SymbolTable &symbol_table_arg) :
254
  symbol_table(symbol_table_arg)
sebastien's avatar
sebastien committed
255
{
256
  cout << "Reading " << filename << "." << endl;
sebastien's avatar
sebastien committed
257

258
259
  ifstream f;
  f.open(filename.c_str(), ios::in);
260
  if (f.fail())
261
262
263
264
    {
      cerr << "ERROR: Can't open " << filename << endl;
      exit(EXIT_FAILURE);
    }
sebastien's avatar
sebastien committed
265

266
  while(true)
sebastien's avatar
sebastien committed
267
    {
268
269
270
271
      string symb_name, value;
      f >> symb_name >> value;
      if (f.eof())
        break;
sebastien's avatar
sebastien committed
272
273
274

      try
        {
275
276
          int symb_id = symbol_table.getID(symb_name);
          content[symb_id] = value;
sebastien's avatar
sebastien committed
277
278
279
        }
      catch(SymbolTable::UnknownSymbolNameException &e)
        {
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
          cerr << "WARNING: Unknown symbol " << symb_name << " in " << filename << endl;
        }
    }
}

void
LoadParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename) const
{
  for(map<int, string>::const_iterator it = content.begin();
      it != content.end(); it++)
    {
      switch(symbol_table.getType(it->first))
        {
        case eParameter:
          output << "M_.params";
          break;
        case eEndogenous:
          output << "oo_.steady_state";
          break;
        case eExogenous:
          output << "oo_.exo_steady_state";
          break;
        case eExogenousDet:
          output << "oo_.exo_det_steady_state";
          break;
        default:
          cerr << "ERROR: Unsupported variable type for " << symbol_table.getName(it->first) << " in load_params_and_steady_state" << endl;
          exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
308
        }
309
310
311

      int tsid = symbol_table.getTypeSpecificID(it->first) + 1;
      output << "(" << tsid << ") = " << it->second << ";" << endl;
sebastien's avatar
sebastien committed
312
313
    }
}
314
315
316
317
318
319
320
321

void
LoadParamsAndSteadyStateStatement::fillEvalContext(eval_context_type &eval_context) const
{
  for(map<int, string>::const_iterator it = content.begin();
      it != content.end(); it++)
    eval_context[it->first] = atof(it->second.c_str());
}