Merge branch 'remots/mytarget/non_linear_decomposition' into 'master'

Adds a new model option "linear_decomposition" that creates a block composed of…

See merge request Dynare/preprocessor!1

src/ComputingTasks.cc

@@ -683,6 +683,58 @@ ForecastStatement::writeJsonOutput(ostream &output) const
   output << "}";
 }
 
+DetCondForecast::DetCondForecast(const SymbolList &symbol_list_arg,
+                                 const OptionsList &options_list_arg,
+                                 const bool linear_decomposition_arg) :
+     options_list(options_list_arg),
+     symbol_list(symbol_list_arg),
+     linear_decomposition(linear_decomposition_arg)
+{
+
+}
+
+void
+DetCondForecast::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
+{
+  options_list.writeOutput(output);
+  if (linear_decomposition)
+    {
+      output << "first_order_solution_to_compute = 1;" << endl;
+      output << "if exist('oo_')" << endl;
+      output << "  if isfield(oo_, 'dr')" << endl;
+      output << "    if isfield(oo_.dr, 'ghx') && isfield(oo_.dr, 'ghu') && isfield(oo_.dr, 'state_var') && isfield(oo_.dr, 'order_var')" << endl;
+      output << "      first_order_solution_to_compute = 0;" << endl;
+      output << "    end;" << endl;
+      output << "  end;" << endl;
+      output << "end;" << endl;
+      output << "if first_order_solution_to_compute" << endl;
+      output << " fprintf('%s','Computing the first order solution ...');" << endl;
+      output << " options_.nograph = 1;" << endl;
+      output << " options_.order = 1;" << endl;
+      output << " options_.noprint = 1;" << endl;
+      output << " options_.nocorr = 1;" << endl;
+      output << " options_.nomoments = 1;" << endl;
+      output << " options_.nodecomposition = 1;" << endl;
+      output << " options_.nofunctions = 1;" << endl;
+      output << " options_.irf = 0;" << endl;
+      output << " tmp_periods = options_.periods;" << endl;
+      output << " options_.periods = 0;" << endl;
+      output << " var_list_ = char();" << endl;
+      output << " info = stoch_simul(var_list_);" << endl;
+      output << " fprintf('%s\\n','done');" << endl;
+      output << " options_.periods = tmp_periods;" << endl;
+      output << "end;" << endl;
+    }
+  vector<string> symbols = symbol_list.get_symbols();
+  if (symbols.size() > 0)
+    output << symbols[1] << " = det_cond_forecast(" ;
+  for (unsigned int i = 0; i < symbols.size() - 1; i++)
+    output << symbols[i] << ", ";
+  if (symbols.size() > 0)
+    output << symbols[symbols.size() - 1];
+  output << ");" << endl;
+}
+
 RamseyModelStatement::RamseyModelStatement(OptionsList options_list_arg) :
   options_list(move(options_list_arg))
 {

src/ComputingTasks.hh

@@ -95,6 +95,21 @@ public:
   void writeJsonOutput(ostream &output) const override;
 };
 
+
+class DetCondForecast : public Statement
+{
+private:
+  const OptionsList options_list;
+  const SymbolList symbol_list;
+  const bool linear_decomposition;
+public:
+  DetCondForecast(const SymbolList &symbol_list_arg,
+                  const OptionsList &options_list_arg,
+                  const bool linear_decompositiontion_arg);
+  //virtual void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings);
+  virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const;
+};
+
 class ModelInfoStatement : public Statement
 {
 private:

src/DynamicModel.cc

@@ -1043,7 +1043,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
 }
 
 void
-DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx) const
+DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx, const bool linear_decomposition) const
 {
   struct Uff_l
   {
@@ -1068,10 +1068,12 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
   deriv_node_temp_terms_t tef_terms;
   vector<int> feedback_variables;
   bool file_open = false;
-
+  string main_name;
   boost::filesystem::create_directories(basename + "/model/bytecode");
-
-  string main_name = basename + "/model/bytecode/dynamic.cod";
+  if (linear_decomposition)
+    main_name = basename + "/model/bytecode/non_linear.cod";
+  else
+    main_name = basename + "/model/bytecode/dynamic.cod";
   code_file.open(main_name, ios::out | ios::binary | ios::ate);
   if (!code_file.is_open())
     {
@@ -1104,7 +1106,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
           || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
         {
           Write_Inf_To_Bin_File_Block(basename, block, u_count_int, file_open,
-                                      simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE);
+                                      simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE, linear_decomposition);
           file_open = true;
         }
       map<pair<int, pair<int, int>>, expr_t> tmp_block_endo_derivative;
@@ -1187,13 +1189,16 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
                                );
       fbeginblock.write(code_file, instruction_number);
 
+      if (linear_decomposition)
+        compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
+
       // The equations
       for (i = 0; i < (int) block_size; i++)
         {
           //The Temporary terms
           temporary_terms_t tt2;
           tt2.clear();
-          if (v_temporary_terms[block][i].size())
+          if (v_temporary_terms[block][i].size() && !linear_decomposition)
             {
               for (auto it : v_temporary_terms[block][i])
                 {
@@ -1715,11 +1720,17 @@ DynamicModel::setNonZeroHessianEquations(map<int, string> &eqs)
 
 void
 DynamicModel::Write_Inf_To_Bin_File_Block(const string &basename, const int &num,
-                                          int &u_count_int, bool &file_open, bool is_two_boundaries) const
+                                          int &u_count_int, bool &file_open, bool is_two_boundaries, const bool linear_decomposition) const
 {
   int j;
   std::ofstream SaveCode;
-  string filename = basename + "/model/bytecode/dynamic.bin";
+  string filename;
+
+  if (!linear_decomposition)
+    filename = basename + "/model/bytecode/dynamic.bin";
+  else
+    filename = basename + "/model/bytecode/non_linear.bin";
+
   if (file_open)
     SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate);
   else
@@ -2814,7 +2825,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
 }
 
 void
-DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const
+DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const
 {
   /* Writing initialisation for M_.lead_lag_incidence matrix
      M_.lead_lag_incidence is a matrix with as many columns as there are
@@ -2956,7 +2967,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
         }
 
   //In case of sparse model, writes the block_decomposition structure of the model
-  if (block_decomposition)
+  if (block_decomposition || linear_decomposition)
     {
       vector<int> state_equ;
       int count_lead_lag_incidence = 0;
@@ -4179,7 +4190,7 @@ DynamicModel::substitutePacExpectation()
 void
 DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
                             const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll,
-                            bool bytecode, const bool nopreprocessoroutput)
+                            bool bytecode, const bool nopreprocessoroutput, bool linear_decomposition)
 {
   assert(jacobianExo || !(hessian || thirdDerivatives || paramsDerivsOrder));
 
@@ -4203,8 +4214,15 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
   // Launch computations
   if (!nopreprocessoroutput)
-    cout << "Computing dynamic model derivatives:" << endl
-         << " - order 1" << endl;
+    {
+      if (linear_decomposition)
+        cout << "Computing nonlinear dynamic model derivatives:" << endl
+             << " - order 1" << endl;
+      else
+        cout << "Computing dynamic model derivatives:" << endl
+             << " - order 1" << endl;
+    }
+
   computeJacobian(vars);
 
   if (hessian)
@@ -4231,13 +4249,51 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
       computeThirdDerivatives(vars);
     }
 
-  if (block)
+  jacob_map_t contemporaneous_jacobian, static_jacobian;
+  map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives;
+  // for each block contains pair<Size, Feddback_variable>
+  vector<pair<int, int> > blocks;
+  vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
+
+  if (linear_decomposition)
     {
-      vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
-      jacob_map_t contemporaneous_jacobian, static_jacobian;
+      map<pair<int, pair<int, int>>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
+      is_equation_linear = equationLinear(first_order_endo_derivatives);
+
+      evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
+
+      if (!computeNaturalNormalization())
+        computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
 
-      // for each block contains pair<Size, Feddback_variable>
-      vector<pair<int, int>> blocks;
+      lag_lead_vector_t equation_lag_lead, variable_lag_lead;
+
+      blocks = select_non_linear_equations_and_variables(is_equation_linear, dynamic_jacobian, equation_reordered, variable_reordered,
+                                                         inv_equation_reordered, inv_variable_reordered,
+                                                         equation_lag_lead, variable_lag_lead,
+                                                         n_static, n_forward, n_backward, n_mixed);
+
+      equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, 0);
+      prologue = 0;
+      epilogue = 0;
+
+      block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, linear_decomposition);
+
+      computeChainRuleJacobian(blocks_derivatives);
+
+      blocks_linear = BlockLinear(blocks_derivatives, variable_reordered);
+
+      collect_block_first_order_derivatives();
+
+      collectBlockVariables();
+
+      global_temporary_terms = true;
+      if (!no_tmp_terms)
+        computeTemporaryTermsOrdered();
+
+    }
+
+  if (block)
+    {
 
       evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
 
@@ -4245,7 +4301,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
       computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
 
-      map<pair<int, pair<int, int>>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
+      first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
 
       equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
@@ -4256,7 +4312,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
       computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed);
 
-      block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type);
+      block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, linear_decomposition);
 
       printBlockDecomposition(blocks);
 
@@ -4685,12 +4741,16 @@ DynamicModel::collectBlockVariables()
 }
 
 void
-DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order, bool julia) const
+DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, int order, bool julia) const
 {
   if (block && bytecode)
-    writeModelEquationsCode_Block(basename, map_idx);
+    writeModelEquationsCode_Block(basename, map_idx, linear_decomposition);
   else if (!block && bytecode)
-    writeModelEquationsCode(basename, map_idx);
+    {
+      if (linear_decomposition)
+        writeModelEquationsCode_Block(basename, map_idx, linear_decomposition);
+      writeModelEquationsCode(basename, map_idx);
+    }
   else if (block && !bytecode)
     writeSparseDynamicMFile(basename);
   else if (use_dll)
@@ -4891,6 +4951,15 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model, const boo
     addEquation(neweq, -1);
 }
 
+void
+DynamicModel::toNonlinearPart(DynamicModel &non_linear_equations_dynamic_model) const
+{
+  // Convert model local variables (need to be done first)
+  for (map<int, expr_t>::const_iterator it = local_variables_table.begin();
+       it != local_variables_table.end(); it++)
+    non_linear_equations_dynamic_model.AddLocalVariable(it->first, it->second);
+}
+
 void
 DynamicModel::toStatic(StaticModel &static_model) const
 {

src/DynamicModel.hh

@@ -118,7 +118,7 @@ private:
   //! Writes the Block reordred structure of the model in M output
   void writeModelEquationsOrdered_M(const string &basename) const;
   //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
-  void writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx) const;
+  void writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx, const bool linear_decomposition) const;
   //! Writes the code of the model in virtual machine bytecode
   void writeModelEquationsCode(const string &basename, const map_idx_t &map_idx) const;
 
@@ -280,9 +280,9 @@ public:
     \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
   */
   void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
-                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, const bool nopreprocessoroutput);
+                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, const bool nopreprocessoroutput, bool linear_decomposition);
   //! Writes model initialization and lead/lag incidence matrix to output
-  void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
+  void writeOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
 
   //! Write JSON AST
   void writeJsonAST(ostream &output) const;
@@ -351,9 +351,9 @@ public:
 
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File_Block(const string &basename,
-                                   const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
+                                   const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries, const bool linear_decomposition) const;
   //! Writes dynamic model file
-  void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order, bool julia) const;
+  void writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, int order, bool julia) const;
   //! Writes file containing parameters derivatives
   void writeParamsDerivativesFile(const string &basename, bool julia) const;
 
@@ -361,6 +361,11 @@ public:
   /*! It assumes that the static model given in argument has just been allocated */
   void toStatic(StaticModel &static_model) const;
 
+  
+  //! Converts to nonlinear model (only the equations)
+  /*! It assumes that the nonlinear model given in argument has just been allocated */
+  void toNonlinearPart(DynamicModel &non_linear_equations_dynamic_model) const;
+
   //! Find endogenous variables not used in model
   set<int> findUnusedEndogenous();
   //! Find exogenous variables not used in model

src/DynareBison.yy

@@ -81,7 +81,7 @@ class ParsingDriver;
 %token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED PROPOSAL_DISTRIBUTION REALTIME VINTAGE
 %token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
 %token COMMA CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED INITIAL_CONDITION_DECOMPOSITION
-%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
+%token DATAFILE FILE SERIES DET_COND_FORECAST DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
 %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR
 %token FILENAME DIRNAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME OSR_PARAMS_BOUNDS KEEP_KALMAN_ALGO_IF_SINGULARITY_IS_DETECTED
 %token <string> FLOAT_NUMBER DATES
@@ -93,7 +93,7 @@ class ParsingDriver;
 %token <string> INT_NUMBER
 %token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS IRF_PLOT_THRESHOLD IRF_CALIBRATION
 %token FAST_KALMAN_FILTER KALMAN_ALGO KALMAN_TOL DIFFUSE_KALMAN_TOL SUBSAMPLES OPTIONS TOLF TOLX PLOT_INIT_DATE PLOT_END_DATE
-%token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_RESULTS_AFTER_LOAD_MH LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LOGDATA LYAPUNOV LINEAR_APPROXIMATION
+%token LAPLACE LIK_ALGO LIK_INIT LINEAR LINEAR_DECOMPOSITION LOAD_IDENT_FILES LOAD_MH_FILE LOAD_RESULTS_AFTER_LOAD_MH LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LOGDATA LYAPUNOV LINEAR_APPROXIMATION
 %token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT
 %token MFS MH_CONF_SIG MH_DROP MH_INIT_SCALE MH_JSCALE MH_TUNE_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER POSTERIOR_MAX_SUBSAMPLE_DRAWS MIN MINIMAL_SOLVING_PERIODS
 %token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
@@ -301,6 +301,7 @@ statement : parameters
           | gmm_estimation
           | smm_estimation
           | shock_groups
+          | det_cond_forecast
           | var_expectation_model
           ;
 
@@ -970,6 +971,7 @@ model_options : BLOCK { driver.block(); }
               | DIFFERENTIATE_FORWARD_VARS EQUAL '(' symbol_list ')' { driver.differentiate_forward_vars_some(); }
               | o_linear
               | PARALLEL_LOCAL_FILES EQUAL '(' parallel_local_filename_list ')'
+              | LINEAR_DECOMPOSITION { driver.linear_decomposition(); }
               ;
 
 model_options_list : model_options_list COMMA model_options
@@ -1477,6 +1479,16 @@ prior_posterior_function_options : o_function
                                  | o_sampling_draws
                                  ;
 
+det_cond_forecast : DET_COND_FORECAST '(' symbol ')' ';'
+                    { driver.det_cond_forecast_linear_decomposition($3); }
+                  | DET_COND_FORECAST '(' symbol COMMA symbol ')' ';'
+                    { driver.det_cond_forecast_linear_decomposition($3,$5); }
+                  | DET_COND_FORECAST '(' symbol COMMA simul_options_list ')' ';'
+                    { driver.det_cond_forecast_linear_decomposition($3); }
+                  | DET_COND_FORECAST '(' symbol COMMA symbol COMMA simul_options_list ')' ';'
+                    { driver.det_cond_forecast_linear_decomposition($3,$5); }
+                  ;
+
 simul : SIMUL ';'
         { driver.simul(); }
       | SIMUL '(' simul_options_list ')' ';'

src/DynareFlex.ll

@@ -186,6 +186,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
 <INITIAL>smoother2histval {BEGIN DYNARE_STATEMENT; return token::SMOOTHER2HISTVAL;}
 <INITIAL>perfect_foresight_setup {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_SETUP;}
 <INITIAL>perfect_foresight_solver {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_SOLVER;}
+<INITIAL>det_cond_forecast {BEGIN DYNARE_STATEMENT; return token::DET_COND_FORECAST;}
 
 <DYNARE_STATEMENT>; {
   if (!sigma_e)
@@ -776,6 +777,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
 <DYNARE_BLOCK>use_dll {return token::USE_DLL;}
 <DYNARE_BLOCK>block {return token::BLOCK;}
 <DYNARE_BLOCK>bytecode {return token::BYTECODE;}
+<DYNARE_BLOCK>linear_decomposition {return token::LINEAR_DECOMPOSITION;}
 <DYNARE_BLOCK>all_values_required {return token::ALL_VALUES_REQUIRED;}
 <DYNARE_BLOCK>no_static {return token::NO_STATIC;}
 <DYNARE_BLOCK>differentiate_forward_vars {return token::DIFFERENTIATE_FORWARD_VARS;}

src/ModFile.cc

@@ -43,6 +43,8 @@ ModFile::ModFile(WarningConsolidation &warnings_arg)
                                        trend_component_model_table, var_model_table),
     orig_ramsey_dynamic_model(symbol_table, num_constants, external_functions_table,
                               trend_component_model_table, var_model_table),
+    non_linear_equations_dynamic_model(symbol_table, num_constants, external_functions_table,
+                  trend_component_model_table, var_model_table),
     epilogue(symbol_table, num_constants, external_functions_table,
              trend_component_model_table, var_model_table),
     static_model(symbol_table, num_constants, external_functions_table),
@@ -191,9 +193,20 @@ ModFile::checkPass(bool nostrict, bool stochastic)
       exit(EXIT_FAILURE);
     }
 
-  if (use_dll && (block || byte_code))
+  if (use_dll && (block || byte_code || linear_decomposition))
+    {
+      cerr << "ERROR: In 'model' block, 'use_dll' option is not compatible with 'block', 'bytecode' or 'linear_decomposition'" << endl;
+      exit(EXIT_FAILURE);
+    }
+  if (block && linear_decomposition)
+    {
+      cerr << "ERROR: In 'model' block, 'block' option is not compatible with 'linear_decomposition'" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  if (!byte_code && linear_decomposition)
     {
-      cerr << "ERROR: In 'model' block, 'use_dll' option is not compatible with 'block' or 'bytecode'" << endl;
+      cerr << "ERROR: For the moment in 'model' block, 'linear_decomposition' option is compatible only with 'bytecode' option" << endl;
       exit(EXIT_FAILURE);
     }
 
@@ -683,6 +696,12 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
 
       // Compute static model and its derivatives
       dynamic_model.toStatic(static_model);
+	  if (linear_decomposition)
+        {
+          dynamic_model.cloneDynamic(non_linear_equations_dynamic_model);
+          non_linear_equations_dynamic_model.set_cutoff_to_zero();
+          non_linear_equations_dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
+        }
       if (!no_static)
         {
           if (mod_file_struct.stoch_simul_present
@@ -707,7 +726,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
           || mod_file_struct.calib_smoother_present)
         {
           if (mod_file_struct.perfect_foresight_solver_present)
-            dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput);
+            dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
           else
             {
               if (mod_file_struct.stoch_simul_present
@@ -732,13 +751,13 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
               int paramsDerivsOrder = 0;
               if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
                 paramsDerivsOrder = params_derivs_order;
-              dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput);
+              dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
               if (linear && mod_file_struct.ramsey_model_present)
-                orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput);
+                orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
             }
         }
       else // No computing task requested, compute derivatives up to 2nd order by default
-        dynamic_model.computingPass(true, true, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput);
+        dynamic_model.computingPass(true, true, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
 
       map<int, string> eqs;
       if (mod_file_struct.ramsey_model_present)
@@ -764,7 +783,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
   for (auto & statement : statements)
     statement->computingPass();
 
-  epilogue.computingPass(true, true, false, 0, global_eval_context, true, false, false, false, true);
+  epilogue.computingPass(true, true, false, 0, global_eval_context, true, false, false, false, true, false);
 }
 
 void
@@ -903,7 +922,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
 
   mOutputFile << "options_.block=" << block << ";" << endl
               << "options_.bytecode=" << byte_code << ";" << endl
-              << "options_.use_dll=" << use_dll << ";" << endl;
+              << "options_.use_dll=" << use_dll << ";" << endl
+              << "options_.linear_decomposition=" << linear_decomposition << ";" << endl;
 
   if (parallel_local_files.size() > 0)
     {
@@ -992,7 +1012,9 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
 
   if (dynamic_model.equation_number() > 0)
     {
-      dynamic_model.writeOutput(mOutputFile, basename, block, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false);
+      if (linear_decomposition)
+        non_linear_equations_dynamic_model.writeOutput(mOutputFile, basename, block, true, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false);
+      dynamic_model.writeOutput(mOutputFile, basename, block, false, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false);
       if (!no_static)
         static_model.writeOutput(mOutputFile, block);
     }
@@ -1066,7 +1088,14 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
               static_model.writeParamsDerivativesFile(basename, false);
             }
 
-          dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option, false);
+          if (linear_decomposition)
+            {
+              non_linear_equations_dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll, mod_file_struct.order_option, false);
+              non_linear_equations_dynamic_model.writeParamsDerivativesFile(basename, false);
+            }
+
+          dynamic_model.writeDynamicFile(basename, block, false, byte_code, use_dll, mod_file_struct.order_option, false);
+
           dynamic_model.writeParamsDerivativesFile(basename, false);
         }
 
@@ -1179,7 +1208,7 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output,
 
   if (dynamic_model.equation_number() > 0)
     {
-      dynamic_model.writeOutput(jlOutputFile, basename, false, false, false,
+      dynamic_model.writeOutput(jlOutputFile, basename, false, false, false, false,
                                 mod_file_struct.order_option,
                                 mod_file_struct.estimation_present, false, true);
       if (!no_static)
@@ -1187,7 +1216,7 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output,
           static_model.writeStaticFile(basename, false, false, false, true);
           static_model.writeParamsDerivativesFile(basename, true);
         }
-      dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll,
+      dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll,
                                      mod_file_struct.order_option, true);
       dynamic_model.writeParamsDerivativesFile(basename, true);
     }

src/ModFile.hh

@@ -67,6 +67,8 @@ public:
   DynamicModel ramsey_FOC_equations_dynamic_model;
   //! A copy of the original model, used to test model linearity under ramsey problem
   DynamicModel orig_ramsey_dynamic_model;
+  //! A copy of the Dynamic model, containing only the non-linear equations w.r. to endogenous variables
+  DynamicModel non_linear_equations_dynamic_model;
   //! Epilogue model, as declared in the "epilogue" block
   Epilogue epilogue;
   //! Static model, as derived from the "model" block when leads and lags have been removed
@@ -100,6 +102,9 @@ public:
     with a lead */
   vector<string> differentiate_forward_vars_subset;
 
+  //! Is the model is block-decomposed according the linear and the non-linear equations
+  bool linear_decomposition;
+
   //! Are nonstationary variables present ?
   bool nonstationary_variables;
 

src/ModelTree.cc

@@ -316,6 +316,105 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
     }
 }
 
+vector<pair<int, int> >
+ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_linear, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered,
+                                                     vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered,
+                                                     lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead,
+                                                     vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed)
+{
+  vector<int> eq2endo(equations.size(), 0);
+  /*equation_reordered.resize(equations.size());
+  variable_reordered.resize(equations.size());*/
+  unsigned int num = 0;
+  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++)
+    if (!is_equation_linear[*it])
+      num++;
+  vector<int> endo2block = vector<int>(endo2eq.size(), 1);
+  vector<pair<set<int>, pair<set<int>, vector<int> > > > components_set(num);
+  int i = 0, j = 0;
+  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, j++)
+    {
+      if (!is_equation_linear[*it])
+        {
+          equation_reordered[i] = *it;
+          variable_reordered[i] = j;
+          endo2block[j] = 0;
+          components_set[endo2block[j]].first.insert(i);
+          /*cout << " -----------------------------------------" << endl;
+          cout.flush();
+          cout << "equation_reordered[" << *it << "] = " << i << endl;
+          cout.flush();
+          cout << "variable_reordered[" << j << "] = " << i << endl;
+          cout.flush();
+          cout << "components_set[" << endo2block[j] << "].first[" << i << "] = " << i << endl;