2009-04-14 16:39:53 +02:00
/*
2017-01-05 15:19:13 +01:00
* Copyright ( C ) 2003 - 2017 Dynare Team
2009-04-14 16:39:53 +02:00
*
* 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/>.
*/
2009-12-16 14:21:31 +01:00
# include <iostream>
# include <cmath>
2009-04-14 16:47:57 +02:00
# include <cstdlib>
2009-04-27 19:15:14 +02:00
# include <cassert>
2009-12-16 14:21:31 +01:00
# include <cstdio>
# include <cerrno>
2009-04-28 19:11:48 +02:00
# include <algorithm>
2009-12-16 14:21:31 +01:00
# include "StaticModel.hh"
// For mkdir() and chdir()
# ifdef _WIN32
# include <direct.h>
# else
# include <unistd.h>
# include <sys / stat.h>
# include <sys / types.h>
2009-06-30 17:07:09 +02:00
# endif
2009-04-28 19:11:48 +02:00
2009-12-16 14:21:31 +01:00
StaticModel : : StaticModel ( SymbolTable & symbol_table_arg ,
2010-02-22 17:33:38 +01:00
NumericalConstants & num_constants_arg ,
ExternalFunctionsTable & external_functions_table_arg ) :
ModelTree ( symbol_table_arg , num_constants_arg , external_functions_table_arg ) ,
2011-06-22 11:56:07 +02:00
global_temporary_terms ( true )
2009-12-16 14:21:31 +01:00
{
}
2009-04-14 16:47:57 +02:00
2009-12-16 14:21:31 +01:00
void
2010-10-11 19:21:32 +02:00
StaticModel : : compileDerivative ( ofstream & code_file , unsigned int & instruction_number , int eq , int symb_id , map_idx_t & map_idx , temporary_terms_t temporary_terms ) const
2009-12-16 14:21:31 +01:00
{
2012-12-01 23:31:37 +01:00
first_derivatives_t : : const_iterator it = first_derivatives . find ( make_pair ( eq , getDerivID ( symbol_table . getID ( eEndogenous , symb_id ) , 0 ) ) ) ;
2009-12-16 14:21:31 +01:00
if ( it ! = first_derivatives . end ( ) )
2010-07-23 11:20:24 +02:00
( it - > second ) - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , false , false ) ;
2009-12-16 14:21:31 +01:00
else
{
FLDZ_ fldz ;
2010-07-23 11:20:24 +02:00
fldz . write ( code_file , instruction_number ) ;
2009-12-16 14:21:31 +01:00
}
}
2009-04-14 16:39:53 +02:00
2009-12-16 14:21:31 +01:00
void
2010-10-11 19:21:32 +02:00
StaticModel : : compileChainRuleDerivative ( ofstream & code_file , unsigned int & instruction_number , int eqr , int varr , int lag , map_idx_t & map_idx , temporary_terms_t temporary_terms ) const
2009-12-16 14:21:31 +01:00
{
2010-09-16 19:18:45 +02:00
map < pair < int , pair < int , int > > , expr_t > : : const_iterator it = first_chain_rule_derivatives . find ( make_pair ( eqr , make_pair ( varr , lag ) ) ) ;
2009-12-16 14:21:31 +01:00
if ( it ! = first_chain_rule_derivatives . end ( ) )
2010-07-23 11:20:24 +02:00
( it - > second ) - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , false , false ) ;
2009-12-16 14:21:31 +01:00
else
{
FLDZ_ fldz ;
2010-07-23 11:20:24 +02:00
fldz . write ( code_file , instruction_number ) ;
2009-12-16 14:21:31 +01:00
}
}
void
StaticModel : : computeTemporaryTermsOrdered ( )
{
2010-09-16 19:18:45 +02:00
map < expr_t , pair < int , int > > first_occurence ;
map < expr_t , int > reference_count ;
2009-12-16 14:21:31 +01:00
BinaryOpNode * eq_node ;
2010-09-16 19:00:48 +02:00
first_derivatives_t : : const_iterator it ;
first_chain_rule_derivatives_t : : const_iterator it_chr ;
2009-12-16 14:21:31 +01:00
ostringstream tmp_s ;
v_temporary_terms . clear ( ) ;
map_idx . clear ( ) ;
unsigned int nb_blocks = getNbBlocks ( ) ;
2010-09-16 19:00:48 +02:00
v_temporary_terms = vector < vector < temporary_terms_t > > ( nb_blocks ) ;
2010-10-11 19:21:32 +02:00
v_temporary_terms_local = vector < vector < temporary_terms_t > > ( nb_blocks ) ;
2009-12-16 14:21:31 +01:00
2010-09-16 19:00:48 +02:00
v_temporary_terms_inuse = vector < temporary_terms_inuse_t > ( nb_blocks ) ;
2009-12-16 14:21:31 +01:00
2010-10-11 19:21:32 +02:00
map_idx2 = vector < map_idx_t > ( nb_blocks ) ;
2009-12-16 14:21:31 +01:00
temporary_terms . clear ( ) ;
2010-10-11 19:21:32 +02:00
//local temporay terms
for ( unsigned int block = 0 ; block < nb_blocks ; block + + )
2009-12-16 14:21:31 +01:00
{
2010-10-11 19:21:32 +02:00
map < expr_t , int > reference_count_local ;
reference_count_local . clear ( ) ;
map < expr_t , pair < int , int > > first_occurence_local ;
first_occurence_local . clear ( ) ;
temporary_terms_t temporary_terms_l ;
temporary_terms_l . clear ( ) ;
2009-12-16 14:21:31 +01:00
2010-10-11 19:21:32 +02:00
unsigned int block_size = getBlockSize ( block ) ;
unsigned int block_nb_mfs = getBlockMfs ( block ) ;
unsigned int block_nb_recursives = block_size - block_nb_mfs ;
v_temporary_terms_local [ block ] = vector < temporary_terms_t > ( block_size ) ;
for ( unsigned int i = 0 ; i < block_size ; i + + )
{
if ( i < block_nb_recursives & & isBlockEquationRenormalized ( block , i ) )
getBlockEquationRenormalizedExpr ( block , i ) - > computeTemporaryTerms ( reference_count_local , temporary_terms_l , first_occurence_local , block , v_temporary_terms_local , i ) ;
else
2009-12-16 14:21:31 +01:00
{
2010-10-11 19:21:32 +02:00
eq_node = ( BinaryOpNode * ) getBlockEquationExpr ( block , i ) ;
eq_node - > computeTemporaryTerms ( reference_count_local , temporary_terms_l , first_occurence_local , block , v_temporary_terms_local , i ) ;
2009-12-16 14:21:31 +01:00
}
}
2010-10-11 19:21:32 +02:00
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = blocks_derivatives [ block ] . begin ( ) ; it ! = ( blocks_derivatives [ block ] ) . end ( ) ; it + + )
{
expr_t id = it - > second . second ;
id - > computeTemporaryTerms ( reference_count_local , temporary_terms_l , first_occurence_local , block , v_temporary_terms_local , block_size - 1 ) ;
}
set < int > temporary_terms_in_use ;
temporary_terms_in_use . clear ( ) ;
v_temporary_terms_inuse [ block ] = temporary_terms_in_use ;
computeTemporaryTermsMapping ( temporary_terms_l , map_idx2 [ block ] ) ;
2009-12-16 14:21:31 +01:00
}
2010-10-11 19:21:32 +02:00
// global temporay terms
for ( unsigned int block = 0 ; block < nb_blocks ; block + + )
2009-12-16 14:21:31 +01:00
{
2010-10-11 19:21:32 +02:00
// Compute the temporary terms reordered
unsigned int block_size = getBlockSize ( block ) ;
unsigned int block_nb_mfs = getBlockMfs ( block ) ;
unsigned int block_nb_recursives = block_size - block_nb_mfs ;
v_temporary_terms [ block ] = vector < temporary_terms_t > ( block_size ) ;
for ( unsigned int i = 0 ; i < block_size ; i + + )
2009-12-16 14:21:31 +01:00
{
2010-10-11 19:21:32 +02:00
if ( i < block_nb_recursives & & isBlockEquationRenormalized ( block , i ) )
getBlockEquationRenormalizedExpr ( block , i ) - > computeTemporaryTerms ( reference_count , temporary_terms , first_occurence , block , v_temporary_terms , i ) ;
else
2009-12-16 14:21:31 +01:00
{
2010-10-11 19:21:32 +02:00
eq_node = ( BinaryOpNode * ) getBlockEquationExpr ( block , i ) ;
eq_node - > computeTemporaryTerms ( reference_count , temporary_terms , first_occurence , block , v_temporary_terms , i ) ;
2009-12-16 14:21:31 +01:00
}
}
2010-10-11 19:21:32 +02:00
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = blocks_derivatives [ block ] . begin ( ) ; it ! = ( blocks_derivatives [ block ] ) . end ( ) ; it + + )
2009-12-16 14:21:31 +01:00
{
2010-10-11 19:21:32 +02:00
expr_t id = it - > second . second ;
id - > computeTemporaryTerms ( reference_count , temporary_terms , first_occurence , block , v_temporary_terms , block_size - 1 ) ;
}
}
for ( unsigned int block = 0 ; block < nb_blocks ; block + + )
{
// Collecte the temporary terms reordered
unsigned int block_size = getBlockSize ( block ) ;
unsigned int block_nb_mfs = getBlockMfs ( block ) ;
unsigned int block_nb_recursives = block_size - block_nb_mfs ;
set < int > temporary_terms_in_use ;
for ( unsigned int i = 0 ; i < block_size ; i + + )
{
if ( i < block_nb_recursives & & isBlockEquationRenormalized ( block , i ) )
getBlockEquationRenormalizedExpr ( block , i ) - > collectTemporary_terms ( temporary_terms , temporary_terms_in_use , block ) ;
else
2009-12-16 14:21:31 +01:00
{
2010-10-11 19:21:32 +02:00
eq_node = ( BinaryOpNode * ) getBlockEquationExpr ( block , i ) ;
eq_node - > collectTemporary_terms ( temporary_terms , temporary_terms_in_use , block ) ;
2009-12-16 14:21:31 +01:00
}
}
2010-10-11 19:21:32 +02:00
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = blocks_derivatives [ block ] . begin ( ) ; it ! = ( blocks_derivatives [ block ] ) . end ( ) ; it + + )
{
expr_t id = it - > second . second ;
id - > collectTemporary_terms ( temporary_terms , temporary_terms_in_use , block ) ;
}
for ( int i = 0 ; i < ( int ) getBlockSize ( block ) ; i + + )
for ( temporary_terms_t : : const_iterator it = v_temporary_terms [ block ] [ i ] . begin ( ) ;
it ! = v_temporary_terms [ block ] [ i ] . end ( ) ; it + + )
( * it ) - > collectTemporary_terms ( temporary_terms , temporary_terms_in_use , block ) ;
v_temporary_terms_inuse [ block ] = temporary_terms_in_use ;
2009-12-16 14:21:31 +01:00
}
2010-10-11 19:21:32 +02:00
computeTemporaryTermsMapping ( temporary_terms , map_idx ) ;
2010-01-22 11:03:29 +01:00
}
void
2010-10-11 19:21:32 +02:00
StaticModel : : computeTemporaryTermsMapping ( temporary_terms_t & temporary_terms , map_idx_t & map_idx )
2010-01-22 11:03:29 +01:00
{
2009-12-16 14:21:31 +01:00
// Add a mapping form node ID to temporary terms order
2009-12-16 18:13:23 +01:00
int j = 0 ;
2010-09-16 19:00:48 +02:00
for ( temporary_terms_t : : const_iterator it = temporary_terms . begin ( ) ;
2011-02-04 16:25:38 +01:00
it ! = temporary_terms . end ( ) ; it + + )
2009-12-16 18:13:23 +01:00
map_idx [ ( * it ) - > idx ] = j + + ;
2009-12-16 14:21:31 +01:00
}
void
StaticModel : : writeModelEquationsOrdered_M ( const string & static_basename ) const
2009-12-16 18:13:23 +01:00
{
string tmp_s , sps ;
ostringstream tmp_output , tmp1_output , global_output ;
2010-09-16 19:18:45 +02:00
expr_t lhs = NULL , rhs = NULL ;
2009-12-16 18:13:23 +01:00
BinaryOpNode * eq_node ;
2010-09-16 19:18:45 +02:00
map < expr_t , int > reference_count ;
2010-09-16 19:00:48 +02:00
temporary_terms_t local_temporary_terms ;
2009-12-16 18:13:23 +01:00
ofstream output ;
vector < int > feedback_variables ;
2010-12-10 11:50:27 +01:00
deriv_node_temp_terms_t tef_terms ;
2009-12-16 18:13:23 +01:00
ExprNodeOutputType local_output_type ;
2009-12-16 14:21:31 +01:00
2010-09-16 16:57:35 +02:00
local_output_type = oMatlabStaticModelSparse ;
2009-12-16 18:13:23 +01:00
if ( global_temporary_terms )
2010-09-16 16:57:35 +02:00
local_temporary_terms = temporary_terms ;
2009-12-16 18:13:23 +01:00
//----------------------------------------------------------------------
//For each block
for ( unsigned int block = 0 ; block < getNbBlocks ( ) ; block + + )
{
//recursive_variables.clear();
feedback_variables . clear ( ) ;
//For a block composed of a single equation determines wether we have to evaluate or to solve the equation
BlockSimulationType simulation_type = getBlockSimulationType ( block ) ;
unsigned int block_size = getBlockSize ( block ) ;
unsigned int block_mfs = getBlockMfs ( block ) ;
unsigned int block_recursive = block_size - block_mfs ;
tmp1_output . str ( " " ) ;
tmp1_output < < static_basename < < " _ " < < block + 1 < < " .m " ;
output . open ( tmp1_output . str ( ) . c_str ( ) , ios : : out | ios : : binary ) ;
output < < " % \n " ;
output < < " % " < < tmp1_output . str ( ) < < " : Computes static model for Dynare \n " ;
output < < " % \n " ;
output < < " % Warning : this file is generated automatically by Dynare \n " ;
output < < " % from model file (.mod) \n \n " ;
output < < " %/ \n " ;
if ( simulation_type = = EVALUATE_BACKWARD | | simulation_type = = EVALUATE_FORWARD )
output < < " function y = " < < static_basename < < " _ " < < block + 1 < < " (y, x, params) \n " ;
else
output < < " function [residual, y, g1] = " < < static_basename < < " _ " < < block + 1 < < " (y, x, params) \n " ;
BlockType block_type ;
if ( simulation_type = = SOLVE_FORWARD_COMPLETE | | simulation_type = = SOLVE_BACKWARD_COMPLETE )
block_type = SIMULTANS ;
else if ( ( simulation_type = = SOLVE_FORWARD_SIMPLE | | simulation_type = = SOLVE_BACKWARD_SIMPLE
| | simulation_type = = EVALUATE_BACKWARD | | simulation_type = = EVALUATE_FORWARD )
& & getBlockFirstEquation ( block ) < prologue )
block_type = PROLOGUE ;
else if ( ( simulation_type = = SOLVE_FORWARD_SIMPLE | | simulation_type = = SOLVE_BACKWARD_SIMPLE
| | simulation_type = = EVALUATE_BACKWARD | | simulation_type = = EVALUATE_FORWARD )
& & getBlockFirstEquation ( block ) > = equations . size ( ) - epilogue )
block_type = EPILOGUE ;
else
block_type = SIMULTANS ;
output < < " % //////////////////////////////////////////////////////////////////////// " < < endl
< < " % // " < < string ( " Block " ) . substr ( int ( log10 ( block + 1 ) ) ) < < block + 1 < < " " < < BlockType0 ( block_type )
< < " // " < < endl
< < " % // Simulation type "
< < BlockSim ( simulation_type ) < < " // " < < endl
< < " % //////////////////////////////////////////////////////////////////////// " < < endl ;
output < < " global options_; " < < endl ;
//The Temporary terms
if ( simulation_type ! = EVALUATE_BACKWARD & & simulation_type ! = EVALUATE_FORWARD )
2010-10-22 16:20:28 +02:00
output < < " g1 = spalloc( " < < block_mfs < < " , " < < block_mfs < < " , " < < derivative_endo [ block ] . size ( ) < < " ); " < < endl ;
2009-12-16 18:13:23 +01:00
if ( v_temporary_terms_inuse [ block ] . size ( ) )
{
tmp_output . str ( " " ) ;
2010-09-16 19:00:48 +02:00
for ( temporary_terms_inuse_t : : const_iterator it = v_temporary_terms_inuse [ block ] . begin ( ) ;
2009-12-16 18:13:23 +01:00
it ! = v_temporary_terms_inuse [ block ] . end ( ) ; it + + )
tmp_output < < " T " < < * it ;
output < < " global " < < tmp_output . str ( ) < < " ; \n " ;
}
if ( simulation_type ! = EVALUATE_BACKWARD & & simulation_type ! = EVALUATE_FORWARD )
output < < " residual=zeros( " < < block_mfs < < " ,1); \n " ;
// The equations
for ( unsigned int i = 0 ; i < block_size ; i + + )
{
if ( ! global_temporary_terms )
local_temporary_terms = v_temporary_terms [ block ] [ i ] ;
2010-09-16 19:00:48 +02:00
temporary_terms_t tt2 ;
2009-12-16 18:13:23 +01:00
tt2 . clear ( ) ;
if ( v_temporary_terms [ block ] . size ( ) )
{
output < < " " < < " % //Temporary variables " < < endl ;
2010-09-16 19:00:48 +02:00
for ( temporary_terms_t : : const_iterator it = v_temporary_terms [ block ] [ i ] . begin ( ) ;
2009-12-16 18:13:23 +01:00
it ! = v_temporary_terms [ block ] [ i ] . end ( ) ; it + + )
{
2014-03-13 11:22:00 +01:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( * it ) ! = NULL )
2010-12-10 11:50:27 +01:00
( * it ) - > writeExternalFunctionOutput ( output , local_output_type , tt2 , tef_terms ) ;
2009-12-16 18:13:23 +01:00
output < < " " < < sps ;
2010-12-10 11:50:27 +01:00
( * it ) - > writeOutput ( output , local_output_type , local_temporary_terms , tef_terms ) ;
2009-12-16 18:13:23 +01:00
output < < " = " ;
2010-12-10 11:50:27 +01:00
( * it ) - > writeOutput ( output , local_output_type , tt2 , tef_terms ) ;
2009-12-16 18:13:23 +01:00
// Insert current node into tt2
tt2 . insert ( * it ) ;
output < < " ; " < < endl ;
}
}
int variable_ID = getBlockVariableID ( block , i ) ;
int equation_ID = getBlockEquationID ( block , i ) ;
EquationType equ_type = getBlockEquationType ( block , i ) ;
string sModel = symbol_table . getName ( symbol_table . getID ( eEndogenous , variable_ID ) ) ;
2010-09-16 19:18:45 +02:00
eq_node = ( BinaryOpNode * ) getBlockEquationExpr ( block , i ) ;
2009-12-16 18:13:23 +01:00
lhs = eq_node - > get_arg1 ( ) ;
rhs = eq_node - > get_arg2 ( ) ;
tmp_output . str ( " " ) ;
lhs - > writeOutput ( tmp_output , local_output_type , local_temporary_terms ) ;
switch ( simulation_type )
{
case EVALUATE_BACKWARD :
case EVALUATE_FORWARD :
evaluation :
output < < " % equation " < < getBlockEquationID ( block , i ) + 1 < < " variable : " < < sModel
< < " ( " < < variable_ID + 1 < < " ) " < < c_Equation_Type ( equ_type ) < < endl ;
output < < " " ;
if ( equ_type = = E_EVALUATE )
{
output < < tmp_output . str ( ) ;
output < < " = " ;
rhs - > writeOutput ( output , local_output_type , local_temporary_terms ) ;
}
else if ( equ_type = = E_EVALUATE_S )
{
output < < " % " < < tmp_output . str ( ) ;
output < < " = " ;
if ( isBlockEquationRenormalized ( block , i ) )
{
rhs - > writeOutput ( output , local_output_type , local_temporary_terms ) ;
output < < " \n " ;
tmp_output . str ( " " ) ;
2010-09-16 19:18:45 +02:00
eq_node = ( BinaryOpNode * ) getBlockEquationRenormalizedExpr ( block , i ) ;
2009-12-16 18:13:23 +01:00
lhs = eq_node - > get_arg1 ( ) ;
rhs = eq_node - > get_arg2 ( ) ;
lhs - > writeOutput ( output , local_output_type , local_temporary_terms ) ;
output < < " = " ;
rhs - > writeOutput ( output , local_output_type , local_temporary_terms ) ;
}
}
else
{
cerr < < " Type missmatch for equation " < < equation_ID + 1 < < " \n " ;
exit ( EXIT_FAILURE ) ;
}
output < < " ; \n " ;
break ;
case SOLVE_BACKWARD_SIMPLE :
case SOLVE_FORWARD_SIMPLE :
case SOLVE_BACKWARD_COMPLETE :
case SOLVE_FORWARD_COMPLETE :
if ( i < block_recursive )
goto evaluation ;
feedback_variables . push_back ( variable_ID ) ;
output < < " % equation " < < equation_ID + 1 < < " variable : " < < sModel
< < " ( " < < variable_ID + 1 < < " ) " < < c_Equation_Type ( equ_type ) < < endl ;
output < < " " < < " residual( " < < i + 1 - block_recursive < < " ) = ( " ;
goto end ;
default :
end :
output < < tmp_output . str ( ) ;
output < < " ) - ( " ;
rhs - > writeOutput ( output , local_output_type , local_temporary_terms ) ;
output < < " ); \n " ;
}
}
// The Jacobian if we have to solve the block
if ( simulation_type = = SOLVE_BACKWARD_SIMPLE | | simulation_type = = SOLVE_FORWARD_SIMPLE
| | simulation_type = = SOLVE_BACKWARD_COMPLETE | | simulation_type = = SOLVE_FORWARD_COMPLETE )
output < < " " < < sps < < " % Jacobian " < < endl ;
switch ( simulation_type )
{
case SOLVE_BACKWARD_SIMPLE :
case SOLVE_FORWARD_SIMPLE :
case SOLVE_BACKWARD_COMPLETE :
case SOLVE_FORWARD_COMPLETE :
2010-09-16 19:00:48 +02:00
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = blocks_derivatives [ block ] . begin ( ) ; it ! = ( blocks_derivatives [ block ] ) . end ( ) ; it + + )
2009-12-16 18:13:23 +01:00
{
unsigned int eq = it - > first . first ;
unsigned int var = it - > first . second ;
unsigned int eqr = getBlockEquationID ( block , eq ) ;
unsigned int varr = getBlockVariableID ( block , var ) ;
2010-09-16 19:18:45 +02:00
expr_t id = it - > second . second ;
2009-12-16 18:13:23 +01:00
output < < " g1( " < < eq + 1 - block_recursive < < " , " < < var + 1 - block_recursive < < " ) = " ;
id - > writeOutput ( output , local_output_type , local_temporary_terms ) ;
output < < " ; % variable= " < < symbol_table . getName ( symbol_table . getID ( eEndogenous , varr ) )
< < " ( " < < 0
< < " ) " < < varr + 1
< < " , equation= " < < eqr + 1 < < endl ;
}
break ;
default :
break ;
}
2010-12-17 14:21:35 +01:00
output < < " end " < < endl ;
2009-12-16 18:13:23 +01:00
output . close ( ) ;
}
}
2009-12-16 14:21:31 +01:00
void
2010-09-16 19:00:48 +02:00
StaticModel : : writeModelEquationsCode ( const string file_name , const string bin_basename , map_idx_t map_idx ) const
2010-01-22 11:03:29 +01:00
{
ostringstream tmp_output ;
ofstream code_file ;
2010-07-23 11:20:24 +02:00
unsigned int instruction_number = 0 ;
2010-01-22 11:03:29 +01:00
bool file_open = false ;
string main_name = file_name ;
main_name + = " .cod " ;
code_file . open ( main_name . c_str ( ) , ios : : out | ios : : binary | ios : : ate ) ;
if ( ! code_file . is_open ( ) )
{
2017-03-01 12:44:31 +01:00
cerr < < " Error : Can't open file \" " < < main_name < < " \" for writing " < < endl ;
2010-01-22 11:03:29 +01:00
exit ( EXIT_FAILURE ) ;
}
int count_u ;
int u_count_int = 0 ;
Write_Inf_To_Bin_File ( file_name , u_count_int , file_open , false , symbol_table . endo_nbr ( ) ) ;
file_open = true ;
//Temporary variables declaration
2010-10-11 19:21:32 +02:00
FDIMST_ fdimst ( temporary_terms . size ( ) ) ;
fdimst . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
FBEGINBLOCK_ fbeginblock ( symbol_table . endo_nbr ( ) ,
SOLVE_FORWARD_COMPLETE ,
0 ,
symbol_table . endo_nbr ( ) ,
variable_reordered ,
equation_reordered ,
false ,
symbol_table . endo_nbr ( ) ,
0 ,
0 ,
2010-07-23 11:20:24 +02:00
u_count_int ,
symbol_table . endo_nbr ( )
2010-01-22 11:03:29 +01:00
) ;
2010-07-23 11:20:24 +02:00
fbeginblock . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
// Add a mapping form node ID to temporary terms order
int j = 0 ;
2010-09-16 19:00:48 +02:00
for ( temporary_terms_t : : const_iterator it = temporary_terms . begin ( ) ;
2010-01-22 11:03:29 +01:00
it ! = temporary_terms . end ( ) ; it + + )
map_idx [ ( * it ) - > idx ] = j + + ;
2010-07-23 11:20:24 +02:00
compileTemporaryTerms ( code_file , instruction_number , temporary_terms , map_idx , false , false ) ;
2010-01-22 11:03:29 +01:00
2010-07-23 11:20:24 +02:00
compileModelEquations ( code_file , instruction_number , temporary_terms , map_idx , false , false ) ;
2010-01-22 11:03:29 +01:00
FENDEQU_ fendequ ;
2010-07-23 11:20:24 +02:00
fendequ . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
2010-10-18 17:28:21 +02:00
// Get the current code_file position and jump if eval = true
streampos pos1 = code_file . tellp ( ) ;
FJMPIFEVAL_ fjmp_if_eval ( 0 ) ;
fjmp_if_eval . write ( code_file , instruction_number ) ;
int prev_instruction_number = instruction_number ;
2010-01-22 11:03:29 +01:00
vector < vector < pair < int , int > > > derivatives ;
derivatives . resize ( symbol_table . endo_nbr ( ) ) ;
count_u = symbol_table . endo_nbr ( ) ;
2010-09-16 19:00:48 +02:00
for ( first_derivatives_t : : const_iterator it = first_derivatives . begin ( ) ;
2010-01-22 11:03:29 +01:00
it ! = first_derivatives . end ( ) ; it + + )
{
int deriv_id = it - > first . second ;
if ( getTypeByDerivID ( deriv_id ) = = eEndogenous )
{
2010-09-16 19:18:45 +02:00
expr_t d1 = it - > second ;
2010-01-22 11:03:29 +01:00
unsigned int eq = it - > first . first ;
int symb = getSymbIDByDerivID ( deriv_id ) ;
unsigned int var = symbol_table . getTypeSpecificID ( symb ) ;
2010-01-22 17:42:08 +01:00
FNUMEXPR_ fnumexpr ( FirstEndoDerivative , eq , var ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
if ( ! derivatives [ eq ] . size ( ) )
derivatives [ eq ] . clear ( ) ;
derivatives [ eq ] . push_back ( make_pair ( var , count_u ) ) ;
2010-07-23 11:20:24 +02:00
d1 - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , false , false ) ;
2010-01-22 11:03:29 +01:00
FSTPSU_ fstpsu ( count_u ) ;
2010-07-23 11:20:24 +02:00
fstpsu . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
count_u + + ;
}
}
for ( int i = 0 ; i < symbol_table . endo_nbr ( ) ; i + + )
{
FLDR_ fldr ( i ) ;
2010-07-23 11:20:24 +02:00
fldr . write ( code_file , instruction_number ) ;
2010-11-20 16:06:07 +01:00
if ( derivatives [ i ] . size ( ) )
2010-01-22 11:03:29 +01:00
{
2011-02-04 16:25:38 +01:00
for ( vector < pair < int , int > > : : const_iterator it = derivatives [ i ] . begin ( ) ;
it ! = derivatives [ i ] . end ( ) ; it + + )
2010-01-22 11:03:29 +01:00
{
2010-11-20 16:06:07 +01:00
FLDSU_ fldsu ( it - > second ) ;
fldsu . write ( code_file , instruction_number ) ;
FLDSV_ fldsv ( eEndogenous , it - > first ) ;
fldsv . write ( code_file , instruction_number ) ;
FBINARY_ fbinary ( oTimes ) ;
2010-07-23 11:20:24 +02:00
fbinary . write ( code_file , instruction_number ) ;
2010-11-20 16:06:07 +01:00
if ( it ! = derivatives [ i ] . begin ( ) )
{
FBINARY_ fbinary ( oPlus ) ;
fbinary . write ( code_file , instruction_number ) ;
}
2010-01-22 11:03:29 +01:00
}
2010-11-20 16:06:07 +01:00
FBINARY_ fbinary ( oMinus ) ;
fbinary . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
}
FSTPSU_ fstpsu ( i ) ;
2010-07-23 11:20:24 +02:00
fstpsu . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
}
2010-10-18 17:28:21 +02:00
// Get the current code_file position and jump = true
streampos pos2 = code_file . tellp ( ) ;
FJMP_ fjmp ( 0 ) ;
fjmp . write ( code_file , instruction_number ) ;
// Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
streampos pos3 = code_file . tellp ( ) ;
code_file . seekp ( pos1 ) ;
FJMPIFEVAL_ fjmp_if_eval1 ( instruction_number - prev_instruction_number ) ;
fjmp_if_eval1 . write ( code_file , instruction_number ) ;
code_file . seekp ( pos3 ) ;
2011-02-04 16:25:38 +01:00
prev_instruction_number = instruction_number ;
2010-10-18 17:28:21 +02:00
temporary_terms_t tt2 ;
tt2 . clear ( ) ;
temporary_terms_t tt3 ;
tt3 . clear ( ) ;
// The Jacobian if we have to solve the block determinsitic bloc
for ( first_derivatives_t : : const_iterator it = first_derivatives . begin ( ) ;
it ! = first_derivatives . end ( ) ; it + + )
{
int deriv_id = it - > first . second ;
if ( getTypeByDerivID ( deriv_id ) = = eEndogenous )
{
expr_t d1 = it - > second ;
unsigned int eq = it - > first . first ;
int symb = getSymbIDByDerivID ( deriv_id ) ;
unsigned int var = symbol_table . getTypeSpecificID ( symb ) ;
FNUMEXPR_ fnumexpr ( FirstEndoDerivative , eq , var ) ;
fnumexpr . write ( code_file , instruction_number ) ;
if ( ! derivatives [ eq ] . size ( ) )
derivatives [ eq ] . clear ( ) ;
derivatives [ eq ] . push_back ( make_pair ( var , count_u ) ) ;
d1 - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , false , false ) ;
2011-02-04 16:25:38 +01:00
FSTPG2_ fstpg2 ( eq , var ) ;
2010-10-18 17:28:21 +02:00
fstpg2 . write ( code_file , instruction_number ) ;
}
}
// Set codefile position to previous JMP_ and set the number of instructions to jump
pos1 = code_file . tellp ( ) ;
code_file . seekp ( pos2 ) ;
FJMP_ fjmp1 ( instruction_number - prev_instruction_number ) ;
fjmp1 . write ( code_file , instruction_number ) ;
code_file . seekp ( pos1 ) ;
2010-01-22 11:03:29 +01:00
FENDBLOCK_ fendblock ;
2010-07-23 11:20:24 +02:00
fendblock . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
FEND_ fend ;
2010-07-23 11:20:24 +02:00
fend . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
code_file . close ( ) ;
}
void
2010-10-11 19:21:32 +02:00
StaticModel : : writeModelEquationsCode_Block ( const string file_name , const string bin_basename , map_idx_t map_idx , vector < map_idx_t > map_idx2 ) const
2009-12-16 18:13:23 +01:00
{
struct Uff_l
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
int u , var , lag ;
Uff_l * pNext ;
} ;
2009-12-16 14:21:31 +01:00
2009-12-16 18:13:23 +01:00
struct Uff
{
Uff_l * Ufl , * Ufl_First ;
} ;
int i , v ;
string tmp_s ;
ostringstream tmp_output ;
ofstream code_file ;
2010-07-23 11:20:24 +02:00
unsigned int instruction_number = 0 ;
2010-09-16 19:18:45 +02:00
expr_t lhs = NULL , rhs = NULL ;
2009-12-16 18:13:23 +01:00
BinaryOpNode * eq_node ;
Uff Uf [ symbol_table . endo_nbr ( ) ] ;
2010-09-16 19:18:45 +02:00
map < expr_t , int > reference_count ;
2009-12-16 18:13:23 +01:00
vector < int > feedback_variables ;
2010-12-10 11:50:27 +01:00
deriv_node_temp_terms_t tef_terms ;
2009-12-16 18:13:23 +01:00
bool file_open = false ;
string main_name = file_name ;
main_name + = " .cod " ;
code_file . open ( main_name . c_str ( ) , ios : : out | ios : : binary | ios : : ate ) ;
if ( ! code_file . is_open ( ) )
{
2017-03-01 12:44:31 +01:00
cerr < < " Error : Can't open file \" " < < main_name < < " \" for writing " < < endl ;
2009-12-16 18:13:23 +01:00
exit ( EXIT_FAILURE ) ;
}
//Temporary variables declaration
2010-10-11 19:21:32 +02:00
FDIMST_ fdimst ( temporary_terms . size ( ) ) ;
fdimst . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
for ( unsigned int block = 0 ; block < getNbBlocks ( ) ; block + + )
{
feedback_variables . clear ( ) ;
if ( block > 0 )
{
FENDBLOCK_ fendblock ;
2010-07-23 11:20:24 +02:00
fendblock . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
}
int count_u ;
int u_count_int = 0 ;
BlockSimulationType simulation_type = getBlockSimulationType ( block ) ;
unsigned int block_size = getBlockSize ( block ) ;
unsigned int block_mfs = getBlockMfs ( block ) ;
unsigned int block_recursive = block_size - block_mfs ;
if ( simulation_type = = SOLVE_TWO_BOUNDARIES_SIMPLE | | simulation_type = = SOLVE_TWO_BOUNDARIES_COMPLETE
| | simulation_type = = SOLVE_BACKWARD_COMPLETE | | simulation_type = = SOLVE_FORWARD_COMPLETE )
{
2010-01-22 11:03:29 +01:00
Write_Inf_To_Bin_File_Block ( file_name , bin_basename , block , u_count_int , file_open ) ;
2009-12-16 18:13:23 +01:00
file_open = true ;
}
FBEGINBLOCK_ fbeginblock ( block_mfs ,
simulation_type ,
getBlockFirstEquation ( block ) ,
block_size ,
variable_reordered ,
equation_reordered ,
blocks_linear [ block ] ,
symbol_table . endo_nbr ( ) ,
0 ,
0 ,
2010-07-23 11:20:24 +02:00
u_count_int ,
2011-02-04 16:25:38 +01:00
/*symbol_table.endo_nbr()*/ block_size
2009-12-16 18:13:23 +01:00
) ;
2010-10-11 19:21:32 +02:00
2010-07-23 11:20:24 +02:00
fbeginblock . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
2010-12-31 16:46:09 +01:00
// Get the current code_file position and jump if eval = true
streampos pos1 = code_file . tellp ( ) ;
FJMPIFEVAL_ fjmp_if_eval ( 0 ) ;
fjmp_if_eval . write ( code_file , instruction_number ) ;
int prev_instruction_number = instruction_number ;
2009-12-16 18:13:23 +01:00
for ( i = 0 ; i < ( int ) block_size ; i + + )
{
//The Temporary terms
2010-09-16 19:00:48 +02:00
temporary_terms_t tt2 ;
2009-12-16 18:13:23 +01:00
tt2 . clear ( ) ;
if ( v_temporary_terms [ block ] . size ( ) )
{
2010-09-16 19:00:48 +02:00
for ( temporary_terms_t : : const_iterator it = v_temporary_terms [ block ] [ i ] . begin ( ) ;
2009-12-16 18:13:23 +01:00
it ! = v_temporary_terms [ block ] [ i ] . end ( ) ; it + + )
{
2014-03-13 11:22:00 +01:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( * it ) ! = NULL )
2010-12-10 11:50:27 +01:00
( * it ) - > compileExternalFunctionOutput ( code_file , instruction_number , false , tt2 , map_idx , false , false , tef_terms ) ;
2011-02-04 16:25:38 +01:00
FNUMEXPR_ fnumexpr ( TemporaryTerm , ( int ) ( map_idx . find ( ( * it ) - > idx ) - > second ) ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-12-10 11:50:27 +01:00
( * it ) - > compile ( code_file , instruction_number , false , tt2 , map_idx , false , false , tef_terms ) ;
2011-02-04 16:25:38 +01:00
FSTPST_ fstpst ( ( int ) ( map_idx . find ( ( * it ) - > idx ) - > second ) ) ;
2010-07-23 11:20:24 +02:00
fstpst . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
// Insert current node into tt2
tt2 . insert ( * it ) ;
}
}
2010-10-11 19:21:32 +02:00
// The equations
2009-12-16 18:13:23 +01:00
int variable_ID , equation_ID ;
EquationType equ_type ;
switch ( simulation_type )
{
evaluation :
case EVALUATE_BACKWARD :
case EVALUATE_FORWARD :
equ_type = getBlockEquationType ( block , i ) ;
2010-01-22 17:42:08 +01:00
{
FNUMEXPR_ fnumexpr ( ModelEquation , getBlockEquationID ( block , i ) ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-01-22 17:42:08 +01:00
}
2009-12-16 18:13:23 +01:00
if ( equ_type = = E_EVALUATE )
{
2010-09-16 19:18:45 +02:00
eq_node = ( BinaryOpNode * ) getBlockEquationExpr ( block , i ) ;
2009-12-16 18:13:23 +01:00
lhs = eq_node - > get_arg1 ( ) ;
rhs = eq_node - > get_arg2 ( ) ;
2010-07-23 11:20:24 +02:00
rhs - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , false , false ) ;
lhs - > compile ( code_file , instruction_number , true , temporary_terms , map_idx , false , false ) ;
2009-12-16 18:13:23 +01:00
}
else if ( equ_type = = E_EVALUATE_S )
{
2010-09-16 19:18:45 +02:00
eq_node = ( BinaryOpNode * ) getBlockEquationRenormalizedExpr ( block , i ) ;
2009-12-16 18:13:23 +01:00
lhs = eq_node - > get_arg1 ( ) ;
rhs = eq_node - > get_arg2 ( ) ;
2010-07-23 11:20:24 +02:00
rhs - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , false , false ) ;
lhs - > compile ( code_file , instruction_number , true , temporary_terms , map_idx , false , false ) ;
2009-12-16 18:13:23 +01:00
}
break ;
case SOLVE_BACKWARD_COMPLETE :
case SOLVE_FORWARD_COMPLETE :
if ( i < ( int ) block_recursive )
goto evaluation ;
variable_ID = getBlockVariableID ( block , i ) ;
equation_ID = getBlockEquationID ( block , i ) ;
feedback_variables . push_back ( variable_ID ) ;
Uf [ equation_ID ] . Ufl = NULL ;
goto end ;
default :
end :
2010-01-22 17:42:08 +01:00
FNUMEXPR_ fnumexpr ( ModelEquation , getBlockEquationID ( block , i ) ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-09-16 19:18:45 +02:00
eq_node = ( BinaryOpNode * ) getBlockEquationExpr ( block , i ) ;
2009-12-16 18:13:23 +01:00
lhs = eq_node - > get_arg1 ( ) ;
rhs = eq_node - > get_arg2 ( ) ;
2010-07-23 11:20:24 +02:00
lhs - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , false , false ) ;
rhs - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , false , false ) ;
2009-12-16 18:13:23 +01:00
FBINARY_ fbinary ( oMinus ) ;
2010-07-23 11:20:24 +02:00
fbinary . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
FSTPR_ fstpr ( i - block_recursive ) ;
2010-07-23 11:20:24 +02:00
fstpr . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
}
}
FENDEQU_ fendequ ;
2010-07-23 11:20:24 +02:00
fendequ . write ( code_file , instruction_number ) ;
2010-12-10 11:50:27 +01:00
2009-12-16 18:13:23 +01:00
// The Jacobian if we have to solve the block
if ( simulation_type ! = EVALUATE_BACKWARD
& & simulation_type ! = EVALUATE_FORWARD )
{
switch ( simulation_type )
{
case SOLVE_BACKWARD_SIMPLE :
case SOLVE_FORWARD_SIMPLE :
2010-01-22 17:42:08 +01:00
{
FNUMEXPR_ fnumexpr ( FirstEndoDerivative , 0 , 0 ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-01-22 17:42:08 +01:00
}
2010-10-11 19:21:32 +02:00
compileDerivative ( code_file , instruction_number , getBlockEquationID ( block , 0 ) , getBlockVariableID ( block , 0 ) , map_idx , temporary_terms ) ;
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
FSTPG_ fstpg ( 0 ) ;
2010-07-23 11:20:24 +02:00
fstpg . write ( code_file , instruction_number ) ;
2009-12-16 14:21:31 +01:00
}
2009-12-16 18:13:23 +01:00
break ;
2009-12-16 14:21:31 +01:00
2009-12-16 18:13:23 +01:00
case SOLVE_BACKWARD_COMPLETE :
case SOLVE_FORWARD_COMPLETE :
count_u = feedback_variables . size ( ) ;
2010-09-16 19:00:48 +02:00
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = blocks_derivatives [ block ] . begin ( ) ; it ! = ( blocks_derivatives [ block ] ) . end ( ) ; it + + )
2009-12-16 18:13:23 +01:00
{
unsigned int eq = it - > first . first ;
unsigned int var = it - > first . second ;
unsigned int eqr = getBlockEquationID ( block , eq ) ;
unsigned int varr = getBlockVariableID ( block , var ) ;
2010-04-28 16:03:32 +02:00
if ( eq > = block_recursive & & var > = block_recursive )
2009-12-16 18:13:23 +01:00
{
if ( ! Uf [ eqr ] . Ufl )
{
Uf [ eqr ] . Ufl = ( Uff_l * ) malloc ( sizeof ( Uff_l ) ) ;
Uf [ eqr ] . Ufl_First = Uf [ eqr ] . Ufl ;
}
else
{
Uf [ eqr ] . Ufl - > pNext = ( Uff_l * ) malloc ( sizeof ( Uff_l ) ) ;
Uf [ eqr ] . Ufl = Uf [ eqr ] . Ufl - > pNext ;
}
Uf [ eqr ] . Ufl - > pNext = NULL ;
Uf [ eqr ] . Ufl - > u = count_u ;
Uf [ eqr ] . Ufl - > var = varr ;
2010-01-22 17:42:08 +01:00
FNUMEXPR_ fnumexpr ( FirstEndoDerivative , eqr , varr ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-10-11 19:21:32 +02:00
compileChainRuleDerivative ( code_file , instruction_number , eqr , varr , 0 , map_idx , temporary_terms ) ;
2009-12-16 18:13:23 +01:00
FSTPSU_ fstpsu ( count_u ) ;
2010-07-23 11:20:24 +02:00
fstpsu . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
count_u + + ;
}
}
for ( i = 0 ; i < ( int ) block_size ; i + + )
{
if ( i > = ( int ) block_recursive )
{
FLDR_ fldr ( i - block_recursive ) ;
2010-07-23 11:20:24 +02:00
fldr . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
FLDZ_ fldz ;
2010-07-23 11:20:24 +02:00
fldz . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
v = getBlockEquationID ( block , i ) ;
for ( Uf [ v ] . Ufl = Uf [ v ] . Ufl_First ; Uf [ v ] . Ufl ; Uf [ v ] . Ufl = Uf [ v ] . Ufl - > pNext )
{
FLDSU_ fldsu ( Uf [ v ] . Ufl - > u ) ;
2010-07-23 11:20:24 +02:00
fldsu . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
FLDSV_ fldsv ( eEndogenous , Uf [ v ] . Ufl - > var ) ;
2010-07-23 11:20:24 +02:00
fldsv . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
FBINARY_ fbinary ( oTimes ) ;
2010-07-23 11:20:24 +02:00
fbinary . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
FCUML_ fcuml ;
2010-07-23 11:20:24 +02:00
fcuml . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
}
Uf [ v ] . Ufl = Uf [ v ] . Ufl_First ;
while ( Uf [ v ] . Ufl )
{
Uf [ v ] . Ufl_First = Uf [ v ] . Ufl - > pNext ;
free ( Uf [ v ] . Ufl ) ;
Uf [ v ] . Ufl = Uf [ v ] . Ufl_First ;
}
FBINARY_ fbinary ( oMinus ) ;
2010-07-23 11:20:24 +02:00
fbinary . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
FSTPSU_ fstpsu ( i - block_recursive ) ;
2010-07-23 11:20:24 +02:00
fstpsu . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
}
}
break ;
default :
break ;
}
}
2010-10-11 19:21:32 +02:00
// Get the current code_file position and jump = true
streampos pos2 = code_file . tellp ( ) ;
FJMP_ fjmp ( 0 ) ;
fjmp . write ( code_file , instruction_number ) ;
// Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
streampos pos3 = code_file . tellp ( ) ;
code_file . seekp ( pos1 ) ;
FJMPIFEVAL_ fjmp_if_eval1 ( instruction_number - prev_instruction_number ) ;
fjmp_if_eval1 . write ( code_file , instruction_number ) ;
code_file . seekp ( pos3 ) ;
2011-02-04 16:25:38 +01:00
prev_instruction_number = instruction_number ;
2010-10-11 19:21:32 +02:00
temporary_terms_t tt2 ;
tt2 . clear ( ) ;
temporary_terms_t tt3 ;
tt3 . clear ( ) ;
2010-12-31 16:46:09 +01:00
deriv_node_temp_terms_t tef_terms2 ;
2010-10-11 19:21:32 +02:00
for ( i = 0 ; i < ( int ) block_size ; i + + )
{
if ( v_temporary_terms_local [ block ] . size ( ) )
{
for ( temporary_terms_t : : const_iterator it = v_temporary_terms_local [ block ] [ i ] . begin ( ) ;
it ! = v_temporary_terms_local [ block ] [ i ] . end ( ) ; it + + )
{
2014-03-13 11:22:00 +01:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( * it ) ! = NULL )
2010-12-31 16:46:09 +01:00
( * it ) - > compileExternalFunctionOutput ( code_file , instruction_number , false , tt3 , map_idx2 [ block ] , false , false , tef_terms2 ) ;
2010-12-10 11:50:27 +01:00
2011-02-04 16:25:38 +01:00
FNUMEXPR_ fnumexpr ( TemporaryTerm , ( int ) ( map_idx2 [ block ] . find ( ( * it ) - > idx ) - > second ) ) ;
2010-10-11 19:21:32 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-12-10 11:50:27 +01:00
( * it ) - > compile ( code_file , instruction_number , false , tt3 , map_idx2 [ block ] , false , false , tef_terms ) ;
2011-02-04 16:25:38 +01:00
FSTPST_ fstpst ( ( int ) ( map_idx2 [ block ] . find ( ( * it ) - > idx ) - > second ) ) ;
2010-10-11 19:21:32 +02:00
fstpst . write ( code_file , instruction_number ) ;
// Insert current node into tt2
tt3 . insert ( * it ) ;
tt2 . insert ( * it ) ;
}
}
// The equations
int variable_ID , equation_ID ;
EquationType equ_type ;
switch ( simulation_type )
{
evaluation_l :
case EVALUATE_BACKWARD :
case EVALUATE_FORWARD :
equ_type = getBlockEquationType ( block , i ) ;
{
FNUMEXPR_ fnumexpr ( ModelEquation , getBlockEquationID ( block , i ) ) ;
fnumexpr . write ( code_file , instruction_number ) ;
}
if ( equ_type = = E_EVALUATE )
{
eq_node = ( BinaryOpNode * ) getBlockEquationExpr ( block , i ) ;
lhs = eq_node - > get_arg1 ( ) ;
rhs = eq_node - > get_arg2 ( ) ;
2010-12-31 16:46:09 +01:00
rhs - > compile ( code_file , instruction_number , false , tt2 , map_idx2 [ block ] , false , false ) ;
lhs - > compile ( code_file , instruction_number , true , tt2 , map_idx2 [ block ] , false , false ) ;
2010-10-11 19:21:32 +02:00
}
else if ( equ_type = = E_EVALUATE_S )
{
eq_node = ( BinaryOpNode * ) getBlockEquationRenormalizedExpr ( block , i ) ;
lhs = eq_node - > get_arg1 ( ) ;
rhs = eq_node - > get_arg2 ( ) ;
2010-12-31 16:46:09 +01:00
rhs - > compile ( code_file , instruction_number , false , tt2 , map_idx2 [ block ] , false , false ) ;
lhs - > compile ( code_file , instruction_number , true , tt2 , map_idx2 [ block ] , false , false ) ;
2010-10-11 19:21:32 +02:00
}
break ;
case SOLVE_BACKWARD_COMPLETE :
case SOLVE_FORWARD_COMPLETE :
if ( i < ( int ) block_recursive )
goto evaluation_l ;
variable_ID = getBlockVariableID ( block , i ) ;
equation_ID = getBlockEquationID ( block , i ) ;
feedback_variables . push_back ( variable_ID ) ;
Uf [ equation_ID ] . Ufl = NULL ;
goto end_l ;
default :
end_l :
FNUMEXPR_ fnumexpr ( ModelEquation , getBlockEquationID ( block , i ) ) ;
fnumexpr . write ( code_file , instruction_number ) ;
eq_node = ( BinaryOpNode * ) getBlockEquationExpr ( block , i ) ;
lhs = eq_node - > get_arg1 ( ) ;
rhs = eq_node - > get_arg2 ( ) ;
2010-12-31 16:46:09 +01:00
lhs - > compile ( code_file , instruction_number , false , tt2 , map_idx2 [ block ] , false , false ) ;
rhs - > compile ( code_file , instruction_number , false , tt2 , map_idx2 [ block ] , false , false ) ;
2010-10-11 19:21:32 +02:00
FBINARY_ fbinary ( oMinus ) ;
fbinary . write ( code_file , instruction_number ) ;
FSTPR_ fstpr ( i - block_recursive ) ;
fstpr . write ( code_file , instruction_number ) ;
}
}
FENDEQU_ fendequ_l ;
fendequ_l . write ( code_file , instruction_number ) ;
// The Jacobian if we have to solve the block determinsitic bloc
switch ( simulation_type )
{
case SOLVE_BACKWARD_SIMPLE :
case SOLVE_FORWARD_SIMPLE :
{
FNUMEXPR_ fnumexpr ( FirstEndoDerivative , 0 , 0 ) ;
fnumexpr . write ( code_file , instruction_number ) ;
}
2011-02-04 16:25:38 +01:00
compileDerivative ( code_file , instruction_number , getBlockEquationID ( block , 0 ) , getBlockVariableID ( block , 0 ) , map_idx2 [ block ] , tt2 /*temporary_terms*/ ) ;
2010-10-11 19:21:32 +02:00
{
2011-02-04 16:25:38 +01:00
FSTPG2_ fstpg2 ( 0 , 0 ) ;
2010-10-11 19:21:32 +02:00
fstpg2 . write ( code_file , instruction_number ) ;
}
break ;
case EVALUATE_BACKWARD :
case EVALUATE_FORWARD :
case SOLVE_BACKWARD_COMPLETE :
case SOLVE_FORWARD_COMPLETE :
count_u = feedback_variables . size ( ) ;
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = blocks_derivatives [ block ] . begin ( ) ; it ! = ( blocks_derivatives [ block ] ) . end ( ) ; it + + )
{
unsigned int eq = it - > first . first ;
unsigned int var = it - > first . second ;
unsigned int eqr = getBlockEquationID ( block , eq ) ;
unsigned int varr = getBlockVariableID ( block , var ) ;
FNUMEXPR_ fnumexpr ( FirstEndoDerivative , eqr , varr , 0 ) ;
fnumexpr . write ( code_file , instruction_number ) ;
2011-02-04 16:25:38 +01:00
compileChainRuleDerivative ( code_file , instruction_number , eqr , varr , 0 , map_idx2 [ block ] , tt2 /*temporary_terms*/ ) ;
2010-10-11 19:21:32 +02:00
2011-02-04 16:25:38 +01:00
FSTPG2_ fstpg2 ( eq , var ) ;
2010-10-11 19:21:32 +02:00
fstpg2 . write ( code_file , instruction_number ) ;
}
break ;
default :
break ;
}
// Set codefile position to previous JMP_ and set the number of instructions to jump
pos1 = code_file . tellp ( ) ;
code_file . seekp ( pos2 ) ;
FJMP_ fjmp1 ( instruction_number - prev_instruction_number ) ;
fjmp1 . write ( code_file , instruction_number ) ;
code_file . seekp ( pos1 ) ;
2009-12-16 18:13:23 +01:00
}
FENDBLOCK_ fendblock ;
2010-07-23 11:20:24 +02:00
fendblock . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
FEND_ fend ;
2010-07-23 11:20:24 +02:00
fend . write ( code_file , instruction_number ) ;
2009-12-16 18:13:23 +01:00
code_file . close ( ) ;
}
2009-12-16 14:21:31 +01:00
void
2010-01-22 11:03:29 +01:00
StaticModel : : Write_Inf_To_Bin_File_Block ( const string & static_basename , const string & bin_basename , const int & num ,
2011-02-04 16:25:38 +01:00
int & u_count_int , bool & file_open ) const
2009-12-16 18:13:23 +01:00
{
int j ;
std : : ofstream SaveCode ;
if ( file_open )
SaveCode . open ( ( bin_basename + " _static.bin " ) . c_str ( ) , ios : : out | ios : : in | ios : : binary | ios : : ate ) ;
else
SaveCode . open ( ( bin_basename + " _static.bin " ) . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! SaveCode . is_open ( ) )
{
2017-03-01 12:44:31 +01:00
cerr < < " Error : Can't open file \" " < < bin_basename < < " _static.bin \" for writing " < < endl ;
2009-12-16 18:13:23 +01:00
exit ( EXIT_FAILURE ) ;
}
u_count_int = 0 ;
unsigned int block_size = getBlockSize ( num ) ;
unsigned int block_mfs = getBlockMfs ( num ) ;
unsigned int block_recursive = block_size - block_mfs ;
2010-09-16 19:00:48 +02:00
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = blocks_derivatives [ num ] . begin ( ) ; it ! = ( blocks_derivatives [ num ] ) . end ( ) ; it + + )
2009-12-16 18:13:23 +01:00
{
unsigned int eq = it - > first . first ;
unsigned int var = it - > first . second ;
int lag = 0 ;
2010-04-28 16:03:32 +02:00
if ( eq > = block_recursive & & var > = block_recursive )
2009-12-16 18:13:23 +01:00
{
int v = eq - block_recursive ;
SaveCode . write ( reinterpret_cast < char * > ( & v ) , sizeof ( v ) ) ;
int varr = var - block_recursive ;
SaveCode . write ( reinterpret_cast < char * > ( & varr ) , sizeof ( varr ) ) ;
SaveCode . write ( reinterpret_cast < char * > ( & lag ) , sizeof ( lag ) ) ;
int u = u_count_int + block_mfs ;
SaveCode . write ( reinterpret_cast < char * > ( & u ) , sizeof ( u ) ) ;
u_count_int + + ;
}
}
for ( j = block_recursive ; j < ( int ) block_size ; j + + )
{
unsigned int varr = getBlockVariableID ( num , j ) ;
SaveCode . write ( reinterpret_cast < char * > ( & varr ) , sizeof ( varr ) ) ;
}
for ( j = block_recursive ; j < ( int ) block_size ; j + + )
{
unsigned int eqr = getBlockEquationID ( num , j ) ;
SaveCode . write ( reinterpret_cast < char * > ( & eqr ) , sizeof ( eqr ) ) ;
}
SaveCode . close ( ) ;
}
2009-12-16 14:21:31 +01:00
2010-09-16 19:18:45 +02:00
map < pair < int , pair < int , int > > , expr_t >
2009-12-16 14:21:31 +01:00
StaticModel : : collect_first_order_derivatives_endogenous ( )
2009-04-14 16:39:53 +02:00
{
2010-09-16 19:18:45 +02:00
map < pair < int , pair < int , int > > , expr_t > endo_derivatives ;
2010-09-16 19:00:48 +02:00
for ( first_derivatives_t : : iterator it2 = first_derivatives . begin ( ) ;
2009-12-16 14:21:31 +01:00
it2 ! = first_derivatives . end ( ) ; it2 + + )
{
2009-12-16 18:13:23 +01:00
if ( getTypeByDerivID ( it2 - > first . second ) = = eEndogenous )
2009-12-16 14:21:31 +01:00
{
int eq = it2 - > first . first ;
2012-12-01 23:31:37 +01:00
int var = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( it2 - > first . second ) ) ;
2009-12-16 14:21:31 +01:00
int lag = 0 ;
endo_derivatives [ make_pair ( eq , make_pair ( var , lag ) ) ] = it2 - > second ;
}
}
2009-12-16 18:13:23 +01:00
return endo_derivatives ;
2009-12-16 14:21:31 +01:00
}
void
2016-05-18 12:26:19 +02:00
StaticModel : : computingPass ( const eval_context_t & eval_context , bool no_tmp_terms , bool hessian , bool thirdDerivatives , int paramsDerivsOrder , bool block , bool bytecode )
2009-12-16 14:21:31 +01:00
{
2011-06-22 11:34:25 +02:00
initializeVariablesAndEquations ( ) ;
2016-03-18 14:55:14 +01:00
vector < BinaryOpNode * > neweqs ;
2016-03-21 11:51:48 +01:00
for ( unsigned int eq = 0 ; eq < equations . size ( ) - aux_equations . size ( ) ; eq + + )
{
expr_t eq_tmp = equations [ eq ] - > substituteStaticAuxiliaryVariable ( ) ;
neweqs . push_back ( dynamic_cast < BinaryOpNode * > ( eq_tmp - > toStatic ( * this ) ) ) ;
}
2016-03-18 14:55:14 +01:00
2016-03-21 20:42:35 +01:00
for ( unsigned int eq = 0 ; eq < aux_equations . size ( ) ; eq + + )
{
expr_t eq_tmp = aux_equations [ eq ] - > substituteStaticAuxiliaryDefinition ( ) ;
neweqs . push_back ( dynamic_cast < BinaryOpNode * > ( eq_tmp - > toStatic ( * this ) ) ) ;
}
2016-03-18 14:55:14 +01:00
equations . clear ( ) ;
copy ( neweqs . begin ( ) , neweqs . end ( ) , back_inserter ( equations ) ) ;
2009-12-16 18:13:23 +01:00
// Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic
2009-12-16 14:21:31 +01:00
set < int > vars ;
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < symbol_table . endo_nbr ( ) ; i + + )
2016-03-21 11:51:48 +01:00
{
int id = symbol_table . getID ( eEndogenous , i ) ;
2016-03-21 20:42:35 +01:00
// if (!symbol_table.isAuxiliaryVariableButNotMultiplier(id))
vars . insert ( getDerivID ( id , 0 ) ) ;
2016-03-21 11:51:48 +01:00
}
2009-12-16 14:21:31 +01:00
// Launch computations
cout < < " Computing static model derivatives: " < < endl
2009-12-16 18:13:23 +01:00
< < " - order 1 " < < endl ;
2009-12-16 14:21:31 +01:00
first_derivatives . clear ( ) ;
computeJacobian ( vars ) ;
if ( hessian )
{
cout < < " - order 2 " < < endl ;
computeHessian ( vars ) ;
}
2013-11-22 21:09:04 +01:00
if ( thirdDerivatives )
{
cout < < " - order 3 " < < endl ;
computeThirdDerivatives ( vars ) ;
}
2016-05-18 12:26:19 +02:00
if ( paramsDerivsOrder > 0 )
2012-11-29 18:07:48 +01:00
{
cout < < " - derivatives of Jacobian/Hessian w.r. to parameters " < < endl ;
2016-05-18 12:26:19 +02:00
computeParamsDerivatives ( paramsDerivsOrder ) ;
2012-11-29 18:07:48 +01:00
if ( ! no_tmp_terms )
computeParamsDerivativesTemporaryTerms ( ) ;
}
2009-12-16 14:21:31 +01:00
if ( block )
{
2010-09-16 19:00:48 +02:00
jacob_map_t contemporaneous_jacobian , static_jacobian ;
2010-07-23 11:20:24 +02:00
vector < unsigned int > n_static , n_forward , n_backward , n_mixed ;
2009-12-16 14:21:31 +01:00
// for each block contains pair<Size, Feddback_variable>
vector < pair < int , int > > blocks ;
evaluateAndReduceJacobian ( eval_context , contemporaneous_jacobian , static_jacobian , dynamic_jacobian , cutoff , false ) ;
2009-12-21 11:29:21 +01:00
computeNonSingularNormalization ( contemporaneous_jacobian , cutoff , static_jacobian , dynamic_jacobian ) ;
2009-12-16 14:21:31 +01:00
2010-09-16 17:51:50 +02:00
computePrologueAndEpilogue ( static_jacobian , equation_reordered , variable_reordered ) ;
2009-12-16 14:21:31 +01:00
2010-09-16 19:18:45 +02:00
map < pair < int , pair < int , int > > , expr_t > first_order_endo_derivatives = collect_first_order_derivatives_endogenous ( ) ;
2009-12-16 14:21:31 +01:00
2010-09-16 17:51:50 +02:00
equation_type_and_normalized_equation = equationTypeDetermination ( first_order_endo_derivatives , variable_reordered , equation_reordered , mfs ) ;
2009-12-16 14:21:31 +01:00
cout < < " Finding the optimal block decomposition of the model ... \n " ;
2010-07-23 11:20:24 +02:00
lag_lead_vector_t equation_lag_lead , variable_lag_lead ;
computeBlockDecompositionAndFeedbackVariablesForEachBlock ( static_jacobian , dynamic_jacobian , equation_reordered , variable_reordered , blocks , equation_type_and_normalized_equation , false , false , mfs , inv_equation_reordered , inv_variable_reordered , equation_lag_lead , variable_lag_lead , n_static , n_forward , n_backward , n_mixed ) ;
2009-12-16 14:21:31 +01:00
2010-07-23 11:20:24 +02:00
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 ) ;
2009-12-16 14:21:31 +01:00
printBlockDecomposition ( blocks ) ;
computeChainRuleJacobian ( blocks_derivatives ) ;
blocks_linear = BlockLinear ( blocks_derivatives , variable_reordered ) ;
collect_block_first_order_derivatives ( ) ;
global_temporary_terms = true ;
if ( ! no_tmp_terms )
computeTemporaryTermsOrdered ( ) ;
}
else
2010-01-22 11:03:29 +01:00
{
if ( ! no_tmp_terms )
{
computeTemporaryTerms ( true ) ;
if ( bytecode )
2010-10-11 19:21:32 +02:00
computeTemporaryTermsMapping ( temporary_terms , map_idx ) ;
2010-01-22 11:03:29 +01:00
}
}
2009-04-14 16:39:53 +02:00
}
void
2009-12-16 14:21:31 +01:00
StaticModel : : writeStaticMFile ( const string & func_name ) const
2009-04-14 16:39:53 +02:00
{
// Writing comments and function definition command
2009-12-16 14:21:31 +01:00
string filename = func_name + " _static.m " ;
ofstream output ;
output . open ( filename . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2013-11-22 21:09:04 +01:00
output < < " function [residual, g1, g2, g3] = " < < func_name + " _static(y, x, params) " < < endl
2009-07-07 16:20:48 +02:00
< < " % " < < endl
< < " % Status : Computes static model for Dynare " < < endl
< < " % " < < endl
2013-07-28 11:32:14 +02:00
< < " % Inputs : " < < endl
< < " % y [M_.endo_nbr by 1] double vector of endogenous variables in declaration order " < < endl
< < " % x [M_.exo_nbr by 1] double vector of exogenous variables in declaration order " < < endl
< < " % params [M_.param_nbr by 1] double vector of parameter values in declaration order " < < endl
< < " % " < < endl
< < " % Outputs: " < < endl
< < " % residual [M_.endo_nbr by 1] double vector of residuals of the static model equations " < < endl
2014-07-17 20:28:16 +02:00
< < " % in order of declaration of the equations. " < < endl
2016-03-29 15:25:51 +02:00
< < " % Dynare may prepend or append auxiliary equations, see M_.aux_vars " < < endl
2013-07-28 11:32:14 +02:00
< < " % g1 [M_.endo_nbr by M_.endo_nbr] double Jacobian matrix of the static model equations; " < < endl
2014-03-14 10:21:42 +01:00
< < " % columns: variables in declaration order " < < endl
< < " % rows: equations in order of declaration " < < endl
2013-07-28 11:32:14 +02:00
< < " % g2 [M_.endo_nbr by (M_.endo_nbr)^2] double Hessian matrix of the static model equations; " < < endl
2014-03-14 10:21:42 +01:00
< < " % columns: variables in declaration order " < < endl
< < " % rows: equations in order of declaration " < < endl
2013-11-22 21:09:04 +01:00
< < " % g3 [M_.endo_nbr by (M_.endo_nbr)^3] double Third derivatives matrix of the static model equations; " < < endl
< < " % columns: variables in declaration order " < < endl
< < " % rows: equations in order of declaration " < < endl
2013-07-28 11:32:14 +02:00
< < " % " < < endl
< < " % " < < endl
2009-07-07 16:20:48 +02:00
< < " % Warning : this file is generated automatically by Dynare " < < endl
2011-08-18 12:44:11 +02:00
< < " % from model file (.mod) " < < endl < < endl ;
2009-07-07 16:20:48 +02:00
2015-07-27 15:33:38 +02:00
writeStaticModel ( output , false , false ) ;
2011-08-18 12:44:11 +02:00
output < < " end " < < endl ;
output . close ( ) ;
}
void
2015-07-27 15:33:38 +02:00
StaticModel : : writeStaticModel ( ostream & StaticOutput , bool use_dll , bool julia ) const
2011-08-18 12:44:11 +02:00
{
2015-09-03 13:50:02 +02:00
ostringstream model_local_vars_output ; // Used for storing model local vars
ostringstream model_output ; // Used for storing model
ostringstream jacobian_output ; // Used for storing jacobian equations
ostringstream hessian_output ; // Used for storing Hessian equations
2013-11-22 21:09:04 +01:00
ostringstream third_derivatives_output ; // Used for storing third order derivatives equations
2015-08-20 14:41:15 +02:00
ostringstream for_sym ;
2015-07-27 15:33:38 +02:00
ExprNodeOutputType output_type = ( use_dll ? oCStaticModel :
julia ? oJuliaStaticModel : oMatlabStaticModel ) ;
2009-07-07 16:20:48 +02:00
2011-08-18 12:44:11 +02:00
deriv_node_temp_terms_t tef_terms ;
2016-06-01 15:28:11 +02:00
temporary_terms_t temp_term_empty ;
2015-09-03 13:50:02 +02:00
temporary_terms_t temp_term_union = temporary_terms_res ;
2016-06-01 15:28:11 +02:00
temporary_terms_t temp_term_union_m_1 ;
2009-07-07 16:20:48 +02:00
2015-09-03 13:50:02 +02:00
writeModelLocalVariables ( model_local_vars_output , output_type , tef_terms ) ;
2009-07-07 16:20:48 +02:00
2016-06-01 15:28:11 +02:00
writeTemporaryTerms ( temporary_terms_res , temp_term_union_m_1 , model_output , output_type , tef_terms ) ;
2015-09-03 13:50:02 +02:00
writeModelEquations ( model_output , output_type ) ;
2009-04-14 16:39:53 +02:00
2013-11-22 21:09:04 +01:00
int nrows = equations . size ( ) ;
int JacobianColsNbr = symbol_table . endo_nbr ( ) ;
int hessianColsNbr = JacobianColsNbr * JacobianColsNbr ;
2009-04-14 16:39:53 +02:00
// Write Jacobian w.r. to endogenous only
2016-06-01 15:28:11 +02:00
temp_term_union_m_1 = temp_term_union ;
2015-09-03 15:47:53 +02:00
temp_term_union . insert ( temporary_terms_g1 . begin ( ) , temporary_terms_g1 . end ( ) ) ;
2015-09-03 13:50:02 +02:00
if ( ! first_derivatives . empty ( ) )
if ( julia )
2016-06-01 15:28:11 +02:00
writeTemporaryTerms ( temp_term_union , temp_term_empty , jacobian_output , output_type , tef_terms ) ;
2015-09-03 13:50:02 +02:00
else
2016-06-01 15:28:11 +02:00
writeTemporaryTerms ( temp_term_union , temp_term_union_m_1 , jacobian_output , output_type , tef_terms ) ;
2010-09-16 19:00:48 +02:00
for ( first_derivatives_t : : const_iterator it = first_derivatives . begin ( ) ;
2009-04-14 16:39:53 +02:00
it ! = first_derivatives . end ( ) ; it + + )
{
int eq = it - > first . first ;
2012-12-01 23:31:37 +01:00
int symb_id = getSymbIDByDerivID ( it - > first . second ) ;
2010-09-16 19:18:45 +02:00
expr_t d1 = it - > second ;
2009-04-14 16:39:53 +02:00
2011-08-18 12:44:11 +02:00
jacobianHelper ( jacobian_output , eq , symbol_table . getTypeSpecificID ( symb_id ) , output_type ) ;
jacobian_output < < " = " ;
2015-09-03 13:50:02 +02:00
d1 - > writeOutput ( jacobian_output , output_type , temp_term_union , tef_terms ) ;
2011-08-18 12:44:11 +02:00
jacobian_output < < " ; " < < endl ;
2009-04-14 16:39:53 +02:00
}
2009-07-07 16:20:48 +02:00
int g2ncols = symbol_table . endo_nbr ( ) * symbol_table . endo_nbr ( ) ;
2011-08-18 12:44:11 +02:00
// Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
2016-06-01 15:28:11 +02:00
temp_term_union_m_1 = temp_term_union ;
2015-09-03 15:47:53 +02:00
temp_term_union . insert ( temporary_terms_g2 . begin ( ) , temporary_terms_g2 . end ( ) ) ;
2015-09-03 13:50:02 +02:00
if ( ! second_derivatives . empty ( ) )
if ( julia )
2016-06-01 15:28:11 +02:00
writeTemporaryTerms ( temp_term_union , temp_term_empty , hessian_output , output_type , tef_terms ) ;
2015-09-03 13:50:02 +02:00
else
2016-06-01 15:28:11 +02:00
writeTemporaryTerms ( temp_term_union , temp_term_union_m_1 , hessian_output , output_type , tef_terms ) ;
2011-08-18 12:44:11 +02:00
int k = 0 ; // Keep the line of a 2nd derivative in v2
for ( second_derivatives_t : : const_iterator it = second_derivatives . begin ( ) ;
it ! = second_derivatives . end ( ) ; it + + )
2009-04-20 12:48:54 +02:00
{
2011-08-18 12:44:11 +02:00
int eq = it - > first . first ;
2012-12-01 23:31:37 +01:00
int symb_id1 = getSymbIDByDerivID ( it - > first . second . first ) ;
int symb_id2 = getSymbIDByDerivID ( it - > first . second . second ) ;
2011-08-18 12:44:11 +02:00
expr_t d2 = it - > second ;
2009-04-20 12:48:54 +02:00
2011-08-18 12:44:11 +02:00
int tsid1 = symbol_table . getTypeSpecificID ( symb_id1 ) ;
int tsid2 = symbol_table . getTypeSpecificID ( symb_id2 ) ;
2009-07-07 16:20:48 +02:00
2011-08-18 12:44:11 +02:00
int col_nb = tsid1 * symbol_table . endo_nbr ( ) + tsid2 ;
int col_nb_sym = tsid2 * symbol_table . endo_nbr ( ) + tsid1 ;
2009-06-25 11:56:41 +02:00
2015-08-20 12:18:19 +02:00
if ( output_type = = oJuliaDynamicModel )
{
2015-08-20 14:02:22 +02:00
for_sym < < " g2[ " < < eq + 1 < < " , " < < col_nb + 1 < < " ] " ;
2015-08-20 12:18:19 +02:00
hessian_output < < " @inbounds " < < for_sym . str ( ) < < " = " ;
2015-09-03 13:50:02 +02:00
d2 - > writeOutput ( hessian_output , output_type , temp_term_union , tef_terms ) ;
2015-08-20 12:18:19 +02:00
hessian_output < < endl ;
}
else
2011-08-18 12:44:11 +02:00
{
2011-08-19 15:05:57 +02:00
sparseHelper ( 2 , hessian_output , k , 0 , output_type ) ;
2011-08-18 12:44:11 +02:00
hessian_output < < " = " < < eq + 1 < < " ; " < < endl ;
2011-08-19 15:05:57 +02:00
sparseHelper ( 2 , hessian_output , k , 1 , output_type ) ;
2015-08-20 12:18:19 +02:00
hessian_output < < " = " < < col_nb + 1 < < " ; " < < endl ;
2011-08-18 12:44:11 +02:00
2011-08-19 15:05:57 +02:00
sparseHelper ( 2 , hessian_output , k , 2 , output_type ) ;
2011-08-18 12:44:11 +02:00
hessian_output < < " = " ;
2015-09-03 13:50:02 +02:00
d2 - > writeOutput ( hessian_output , output_type , temp_term_union , tef_terms ) ;
2011-08-18 12:44:11 +02:00
hessian_output < < " ; " < < endl ;
k + + ;
2009-04-20 12:48:54 +02:00
}
2015-08-20 12:18:19 +02:00
// Treating symetric elements
if ( symb_id1 ! = symb_id2 )
if ( output_type = = oJuliaDynamicModel )
2015-08-20 14:02:22 +02:00
hessian_output < < " @inbounds g2[ " < < eq + 1 < < " , " < < col_nb_sym + 1 < < " ] = "
2015-08-20 12:18:19 +02:00
< < for_sym . str ( ) < < endl ;
else
{
sparseHelper ( 2 , hessian_output , k , 0 , output_type ) ;
hessian_output < < " = " < < eq + 1 < < " ; " < < endl ;
sparseHelper ( 2 , hessian_output , k , 1 , output_type ) ;
hessian_output < < " = " < < col_nb_sym + 1 < < " ; " < < endl ;
sparseHelper ( 2 , hessian_output , k , 2 , output_type ) ;
hessian_output < < " = " ;
sparseHelper ( 2 , hessian_output , k - 1 , 2 , output_type ) ;
hessian_output < < " ; " < < endl ;
k + + ;
}
2011-08-18 12:44:11 +02:00
}
2009-04-14 16:39:53 +02:00
2013-11-22 21:09:04 +01:00
// Writing third derivatives
2016-06-01 15:28:11 +02:00
temp_term_union_m_1 = temp_term_union ;
2015-09-03 15:47:53 +02:00
temp_term_union . insert ( temporary_terms_g3 . begin ( ) , temporary_terms_g3 . end ( ) ) ;
2015-09-03 13:50:02 +02:00
if ( ! third_derivatives . empty ( ) )
if ( julia )
2016-06-01 15:28:11 +02:00
writeTemporaryTerms ( temp_term_union , temp_term_empty , third_derivatives_output , output_type , tef_terms ) ;
2015-09-03 13:50:02 +02:00
else
2016-06-01 15:28:11 +02:00
writeTemporaryTerms ( temp_term_union , temp_term_union_m_1 , third_derivatives_output , output_type , tef_terms ) ;
2013-11-22 21:09:04 +01:00
k = 0 ; // Keep the line of a 3rd derivative in v3
for ( third_derivatives_t : : const_iterator it = third_derivatives . begin ( ) ;
it ! = third_derivatives . end ( ) ; it + + )
{
int eq = it - > first . first ;
int var1 = it - > first . second . first ;
int var2 = it - > first . second . second . first ;
int var3 = it - > first . second . second . second ;
expr_t d3 = it - > second ;
int id1 = getSymbIDByDerivID ( var1 ) ;
int id2 = getSymbIDByDerivID ( var2 ) ;
int id3 = getSymbIDByDerivID ( var3 ) ;
// Reference column number for the g3 matrix
int ref_col = id1 * hessianColsNbr + id2 * JacobianColsNbr + id3 ;
2015-08-20 14:41:15 +02:00
if ( output_type = = oJuliaDynamicModel )
{
for_sym < < " g3[ " < < eq + 1 < < " , " < < ref_col + 1 < < " ] " ;
third_derivatives_output < < " @inbounds " < < for_sym . str ( ) < < " = " ;
2015-09-03 13:50:02 +02:00
d3 - > writeOutput ( third_derivatives_output , output_type , temp_term_union , tef_terms ) ;
2015-08-20 14:41:15 +02:00
third_derivatives_output < < endl ;
}
else
{
sparseHelper ( 3 , third_derivatives_output , k , 0 , output_type ) ;
third_derivatives_output < < " = " < < eq + 1 < < " ; " < < endl ;
2013-11-22 21:09:04 +01:00
2015-08-20 14:41:15 +02:00
sparseHelper ( 3 , third_derivatives_output , k , 1 , output_type ) ;
third_derivatives_output < < " = " < < ref_col + 1 < < " ; " < < endl ;
2013-11-22 21:09:04 +01:00
2015-08-20 14:41:15 +02:00
sparseHelper ( 3 , third_derivatives_output , k , 2 , output_type ) ;
third_derivatives_output < < " = " ;
2015-09-03 13:50:02 +02:00
d3 - > writeOutput ( third_derivatives_output , output_type , temp_term_union , tef_terms ) ;
2015-08-20 14:41:15 +02:00
third_derivatives_output < < " ; " < < endl ;
}
2013-11-22 21:09:04 +01:00
2015-08-20 14:41:15 +02:00
// Compute the column numbers for the 5 other permutations of (id1,id2,id3)
// and store them in a set (to avoid duplicates if two indexes are equal)
2013-11-22 21:09:04 +01:00
set < int > cols ;
cols . insert ( id1 * hessianColsNbr + id3 * JacobianColsNbr + id2 ) ;
cols . insert ( id2 * hessianColsNbr + id1 * JacobianColsNbr + id3 ) ;
cols . insert ( id2 * hessianColsNbr + id3 * JacobianColsNbr + id1 ) ;
cols . insert ( id3 * hessianColsNbr + id1 * JacobianColsNbr + id2 ) ;
cols . insert ( id3 * hessianColsNbr + id2 * JacobianColsNbr + id1 ) ;
int k2 = 1 ; // Keeps the offset of the permutation relative to k
for ( set < int > : : iterator it2 = cols . begin ( ) ; it2 ! = cols . end ( ) ; it2 + + )
if ( * it2 ! = ref_col )
2015-08-20 14:41:15 +02:00
if ( output_type = = oJuliaDynamicModel )
third_derivatives_output < < " @inbounds g3[ " < < eq + 1 < < " , " < < * it2 + 1 < < " ] = "
< < for_sym . str ( ) < < endl ;
else
{
sparseHelper ( 3 , third_derivatives_output , k + k2 , 0 , output_type ) ;
third_derivatives_output < < " = " < < eq + 1 < < " ; " < < endl ;
2013-11-22 21:09:04 +01:00
2015-08-20 14:41:15 +02:00
sparseHelper ( 3 , third_derivatives_output , k + k2 , 1 , output_type ) ;
third_derivatives_output < < " = " < < * it2 + 1 < < " ; " < < endl ;
2013-11-22 21:09:04 +01:00
2015-08-20 14:41:15 +02:00
sparseHelper ( 3 , third_derivatives_output , k + k2 , 2 , output_type ) ;
third_derivatives_output < < " = " ;
sparseHelper ( 3 , third_derivatives_output , k , 2 , output_type ) ;
third_derivatives_output < < " ; " < < endl ;
2013-11-22 21:09:04 +01:00
2015-08-20 14:41:15 +02:00
k2 + + ;
}
2013-11-22 21:09:04 +01:00
k + = k2 ;
}
2015-08-18 15:41:27 +02:00
if ( output_type = = oMatlabStaticModel )
2011-08-18 12:44:11 +02:00
{
2016-12-28 14:02:50 +01:00
// Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201
2016-12-30 11:26:56 +01:00
map < string , string > tmp_paren_vars ;
2017-01-05 15:19:13 +01:00
bool message_printed = false ;
fixNestedParenthesis ( model_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( model_local_vars_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( jacobian_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( hessian_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( third_derivatives_output , tmp_paren_vars , message_printed ) ;
2016-12-28 14:02:50 +01:00
2011-08-18 12:44:11 +02:00
StaticOutput < < " residual = zeros( " < < equations . size ( ) < < " , 1); " < < endl < < endl
< < " % " < < endl
< < " % Model equations " < < endl
< < " % " < < endl < < endl
2015-09-03 13:50:02 +02:00
< < model_local_vars_output . str ( )
2011-08-18 12:44:11 +02:00
< < model_output . str ( )
< < " if ~isreal(residual) " < < endl
< < " residual = real(residual)+imag(residual).^2; " < < endl
< < " end " < < endl
< < " if nargout >= 2, " < < endl
< < " g1 = zeros( " < < equations . size ( ) < < " , " < < symbol_table . endo_nbr ( ) < < " ); " < < endl < < endl
< < " % " < < endl
< < " % Jacobian matrix " < < endl
< < " % " < < endl < < endl
< < jacobian_output . str ( )
< < " if ~isreal(g1) " < < endl
2012-05-18 16:19:07 +02:00
< < " g1 = real(g1)+2*imag(g1); " < < endl
2011-08-18 12:44:11 +02:00
< < " end " < < endl
< < " if nargout >= 3, " < < endl
< < " % " < < endl
< < " % Hessian matrix " < < endl
< < " % " < < endl
< < endl ;
if ( second_derivatives . size ( ) )
StaticOutput < < " v2 = zeros( " < < NNZDerivatives [ 1 ] < < " ,3); " < < endl
< < hessian_output . str ( )
< < " g2 = sparse(v2(:,1),v2(:,2),v2(:,3), " < < equations . size ( ) < < " , " < < g2ncols < < " ); " < < endl ;
else
StaticOutput < < " g2 = sparse([],[],[], " < < equations . size ( ) < < " , " < < g2ncols < < " ); " < < endl ;
2015-09-04 14:58:14 +02:00
2013-11-22 21:09:04 +01:00
// Initialize g3 matrix
StaticOutput < < " if nargout >= 4, " < < endl
< < " % " < < endl
< < " % Third order derivatives " < < endl
< < " % " < < endl
< < endl ;
int ncols = hessianColsNbr * JacobianColsNbr ;
if ( third_derivatives . size ( ) )
StaticOutput < < " v3 = zeros( " < < NNZDerivatives [ 2 ] < < " ,3); " < < endl
< < third_derivatives_output . str ( )
< < " g3 = sparse(v3(:,1),v3(:,2),v3(:,3), " < < nrows < < " , " < < ncols < < " ); " < < endl ;
else // Either 3rd derivatives is all zero, or we didn't compute it
StaticOutput < < " g3 = sparse([],[],[], " < < nrows < < " , " < < ncols < < " ); " < < endl ;
2015-09-04 14:58:14 +02:00
StaticOutput < < " end " < < endl
< < " end " < < endl
< < " end " < < endl ;
2015-07-27 15:33:38 +02:00
}
2015-08-18 15:41:27 +02:00
else if ( output_type = = oCStaticModel )
2011-08-18 12:44:11 +02:00
{
StaticOutput < < " void Static(double *y, double *x, int nb_row_x, double *params, double *residual, double *g1, double *v2) " < < endl
< < " { " < < endl
< < " double lhs, rhs; " < < endl
< < endl
< < " /* Residual equations */ " < < endl
2015-09-03 13:50:02 +02:00
< < model_local_vars_output . str ( )
2011-08-18 12:44:11 +02:00
< < model_output . str ( )
< < " /* Jacobian */ " < < endl
< < " if (g1 == NULL) " < < endl
< < " return; " < < endl
2015-09-04 14:58:14 +02:00
< < endl
2011-08-18 12:44:11 +02:00
< < jacobian_output . str ( )
2015-09-04 14:58:14 +02:00
< < endl ;
2011-08-18 12:44:11 +02:00
if ( second_derivatives . size ( ) )
StaticOutput < < " /* Hessian for endogenous and exogenous variables */ " < < endl
< < " if (v2 == NULL) " < < endl
< < " return; " < < endl
2015-09-04 14:58:14 +02:00
< < endl
2011-08-18 12:44:11 +02:00
< < hessian_output . str ( )
2015-09-04 14:58:14 +02:00
< < endl ;
2013-11-22 21:09:04 +01:00
if ( third_derivatives . size ( ) )
StaticOutput < < " /* Third derivatives for endogenous and exogenous variables */ " < < endl
2015-09-04 14:58:14 +02:00
< < " if (v3 == NULL) " < < endl
< < " return; " < < endl
< < endl
< < third_derivatives_output . str ( )
< < endl ;
2011-08-18 12:44:11 +02:00
}
2015-08-18 15:41:27 +02:00
else
{
2015-08-21 12:00:47 +02:00
ostringstream comments ;
comments < < " ## Function Arguments " < < endl
< < endl
< < " ## Input " < < endl
< < " 1 y: Array{Float64, length(model.endo), 1} Vector of endogenous variables in declaration order " < < endl
< < " 2 x: Array{Float64, length(model.exo), 1} Vector of exogenous variables in declaration order " < < endl
< < " 3 params: Array{Float64, length(model.param), 1} Vector of parameter values in declaration order " < < endl
< < endl
< < " ## Output " < < endl
< < " 4 residual: Array(Float64, model.eq_nbr, 1) Vector of residuals of the static model equations " < < endl
< < " in order of declaration of the equations. " < < endl
< < " Dynare may prepend auxiliary equations, see model.aux_vars " < < endl ;
2015-08-19 14:38:35 +02:00
StaticOutput < < " function static!(y::Vector{Float64}, x::Vector{Float64}, "
2015-08-19 12:11:08 +02:00
< < " params::Vector{Float64}, " < < endl
< < " residual::Vector{Float64}) " < < endl
2015-08-21 12:00:47 +02:00
< < " #= " < < endl < < comments . str ( ) < < " =# " < < endl
2016-05-20 02:33:31 +02:00
< < " @assert length(y) == " < < symbol_table . endo_nbr ( ) < < endl
< < " @assert length(x) == " < < symbol_table . exo_nbr ( ) < < endl
< < " @assert length(params) == " < < symbol_table . param_nbr ( ) < < endl
< < " @assert length(residual) == " < < equations . size ( ) < < endl
2015-08-21 16:54:11 +02:00
< < " # " < < endl
< < " # Model equations " < < endl
< < " # " < < endl
2015-09-03 13:50:02 +02:00
< < model_local_vars_output . str ( )
2015-08-18 15:41:27 +02:00
< < model_output . str ( )
< < " if ~isreal(residual) " < < endl
< < " residual = real(residual)+imag(residual).^2; " < < endl
< < " end " < < endl
< < " end " < < endl < < endl
2015-08-19 14:38:35 +02:00
< < " function static!(y::Vector{Float64}, x::Vector{Float64}, "
2015-08-19 12:11:08 +02:00
< < " params::Vector{Float64}, " < < endl
2015-08-21 12:00:47 +02:00
< < " residual::Vector{Float64}, g1::Matrix{Float64}) " < < endl ;
comments < < " 5 g1: Array(Float64, model.eq_nbr, length(model.endo)) Jacobian matrix of the static model equations; " < < endl
< < " columns: variables in declaration order " < < endl
< < " rows: equations in order of declaration " < < endl ;
StaticOutput < < " #= " < < endl < < comments . str ( ) < < " =# " < < endl
2015-08-20 15:08:24 +02:00
< < " @assert size(g1) == ( " < < equations . size ( ) < < " , " < < symbol_table . endo_nbr ( )
< < " ) " < < endl
2015-08-21 12:00:47 +02:00
< < " fill!(g1, 0.0) " < < endl
2015-08-18 15:41:27 +02:00
< < " static!(y, x, params, residual) " < < endl
2015-09-03 13:50:02 +02:00
< < model_local_vars_output . str ( )
2015-08-18 15:41:27 +02:00
< < " # " < < endl
< < " # Jacobian matrix " < < endl
2015-08-21 16:28:43 +02:00
< < " # " < < endl
2015-08-18 15:41:27 +02:00
< < jacobian_output . str ( )
< < " if ~isreal(g1) " < < endl
< < " g1 = real(g1)+2*imag(g1); " < < endl
< < " end " < < endl
< < " end " < < endl < < endl
2015-08-19 14:38:35 +02:00
< < " function static!(y::Vector{Float64}, x::Vector{Float64}, "
2015-08-19 12:11:08 +02:00
< < " params::Vector{Float64}, " < < endl
< < " residual::Vector{Float64}, g1::Matrix{Float64}, "
2015-08-21 12:00:47 +02:00
< < " g2::Matrix{Float64}) " < < endl ;
comments < < " 6 g2: spzeros(model.eq_nbr, length(model.endo)^2) Hessian matrix of the static model equations; " < < endl
< < " columns: variables in declaration order " < < endl
< < " rows: equations in order of declaration " < < endl ;
StaticOutput < < " #= " < < endl < < comments . str ( ) < < " =# " < < endl
2015-08-20 15:08:24 +02:00
< < " @assert size(g2) == ( " < < equations . size ( ) < < " , " < < g2ncols < < " ) " < < endl
2015-08-20 14:41:15 +02:00
< < " static!(y, x, params, residual, g1) " < < endl ;
2015-08-18 15:41:27 +02:00
if ( second_derivatives . size ( ) )
2015-09-03 13:50:02 +02:00
StaticOutput < < model_local_vars_output . str ( )
2015-08-20 14:41:15 +02:00
< < " # " < < endl
< < " # Hessian matrix " < < endl
< < " # " < < endl
< < hessian_output . str ( ) ;
2015-08-18 15:41:27 +02:00
// Initialize g3 matrix
2015-08-20 11:32:26 +02:00
int ncols = hessianColsNbr * JacobianColsNbr ;
2015-08-18 15:41:27 +02:00
StaticOutput < < " end " < < endl < < endl
2015-08-19 14:38:35 +02:00
< < " function static!(y::Vector{Float64}, x::Vector{Float64}, "
2015-08-19 12:11:08 +02:00
< < " params::Vector{Float64}, " < < endl
< < " residual::Vector{Float64}, g1::Matrix{Float64}, "
< < " g2::Matrix{Float64}, " < < endl
2015-08-21 12:00:47 +02:00
< < " g3::Matrix{Float64}) " < < endl ;
comments < < " 7 g3: spzeros(model.eq_nbr, length(model.endo)^3) Third derivatives matrix of the static model equations; " < < endl
< < " columns: variables in declaration order " < < endl
< < " rows: equations in order of declaration " < < endl ;
StaticOutput < < " #= " < < endl < < comments . str ( ) < < " =# " < < endl
2015-08-20 15:08:24 +02:00
< < " @assert size(g3) == ( " < < nrows < < " , " < < ncols < < " ) " < < endl
2015-08-20 14:41:15 +02:00
< < " static!(y, x, params, residual, g1, g2) " < < endl ;
2015-08-18 15:41:27 +02:00
if ( third_derivatives . size ( ) )
2015-09-03 13:50:02 +02:00
StaticOutput < < model_local_vars_output . str ( )
2015-08-20 14:41:15 +02:00
< < " # " < < endl
< < " # Third order derivatives " < < endl
< < " # " < < endl
< < third_derivatives_output . str ( ) ;
2015-08-18 15:41:27 +02:00
StaticOutput < < " end " < < endl ;
}
2011-08-18 12:44:11 +02:00
}
void
StaticModel : : writeStaticCFile ( const string & func_name ) const
{
// Writing comments and function definition command
string filename = func_name + " _static.c " ;
2011-12-22 14:55:57 +01:00
string filename_mex = func_name + " _static_mex.c " ;
2011-08-18 12:44:11 +02:00
ofstream output ;
output . open ( filename . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
output < < " /* " < < endl
< < " * " < < filename < < " : Computes static model for Dynare " < < endl
< < " * " < < endl
< < " * Warning : this file is generated automatically by Dynare " < < endl
< < " * from model file (.mod) " < < endl < < endl
< < " */ " < < endl
2016-10-17 14:59:07 +02:00
# if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
< < " #ifdef _MSC_VER " < < endl
< < " #define _USE_MATH_DEFINES " < < endl
< < " #endif " < < endl
# endif
2011-12-22 14:55:57 +01:00
< < " #include <math.h> " < < endl ;
if ( external_functions_table . get_total_number_of_unique_model_block_external_functions ( ) )
// External Matlab function, implies Static function will call mex
output < < " #include \" mex.h \" " < < endl ;
else
output < < " #include <stdlib.h> " < < endl ;
output < < " #define max(a, b) (((a) > (b)) ? (a) : (b)) " < < endl
2011-08-18 12:44:11 +02:00
< < " #define min(a, b) (((a) > (b)) ? (b) : (a)) " < < endl ;
2009-07-07 16:20:48 +02:00
2011-08-18 12:44:11 +02:00
// Write function definition if oPowerDeriv is used
writePowerDerivCHeader ( output ) ;
2016-10-14 16:11:44 +02:00
writeNormcdfCHeader ( output ) ;
2011-08-18 12:44:11 +02:00
// Writing the function body
2015-07-27 15:33:38 +02:00
writeStaticModel ( output , true , false ) ;
2011-08-18 12:44:11 +02:00
output < < " } " < < endl < < endl ;
2016-10-14 14:30:14 +02:00
writePowerDeriv ( output ) ;
2016-10-14 16:11:44 +02:00
writeNormcdf ( output ) ;
2011-12-22 14:55:57 +01:00
output . close ( ) ;
output . open ( filename_mex . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename_mex < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2011-08-18 12:44:11 +02:00
// Writing the gateway routine
2011-12-22 14:55:57 +01:00
output < < " /* " < < endl
< < " * " < < filename_mex < < " : The gateway routine used to call the Static function "
< < " located in " < < filename < < endl
< < " * " < < endl
< < " * Warning : this file is generated automatically by Dynare " < < endl
< < " * from model file (.mod) " < < endl < < endl
< < " */ " < < endl < < endl
< < " #include \" mex.h \" " < < endl < < endl
< < " void Static(double *y, double *x, int nb_row_x, double *params, double *residual, double *g1, double *v2); " < < endl
2011-08-18 12:44:11 +02:00
< < " void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) " < < endl
< < " { " < < endl
< < " double *y, *x, *params; " < < endl
< < " double *residual, *g1, *v2; " < < endl
< < " int nb_row_x; " < < endl
< < endl
< < " /* Create a pointer to the input matrix y. */ " < < endl
< < " y = mxGetPr(prhs[0]); " < < endl
< < endl
< < " /* Create a pointer to the input matrix x. */ " < < endl
< < " x = mxGetPr(prhs[1]); " < < endl
< < endl
< < " /* Create a pointer to the input matrix params. */ " < < endl
< < " params = mxGetPr(prhs[2]); " < < endl
< < endl
< < " /* Gets number of rows of matrix x. */ " < < endl
< < " nb_row_x = mxGetM(prhs[1]); " < < endl
< < endl
< < " residual = NULL; " < < endl
< < " if (nlhs >= 1) " < < endl
< < " { " < < endl
< < " /* Set the output pointer to the output matrix residual. */ " < < endl
< < " plhs[0] = mxCreateDoubleMatrix( " < < equations . size ( ) < < " ,1, mxREAL); " < < endl
< < " /* Create a C pointer to a copy of the output matrix residual. */ " < < endl
< < " residual = mxGetPr(plhs[0]); " < < endl
< < " } " < < endl
< < endl
< < " g1 = NULL; " < < endl
< < " if (nlhs >= 2) " < < endl
< < " { " < < endl
< < " /* Set the output pointer to the output matrix g1. */ " < < endl
< < " plhs[1] = mxCreateDoubleMatrix( " < < equations . size ( ) < < " , " < < symbol_table . endo_nbr ( ) < < " , mxREAL); " < < endl
< < " /* Create a C pointer to a copy of the output matrix g1. */ " < < endl
< < " g1 = mxGetPr(plhs[1]); " < < endl
< < " } " < < endl
< < endl
< < " v2 = NULL; " < < endl
< < " if (nlhs >= 3) " < < endl
< < " { " < < endl
< < " /* Set the output pointer to the output matrix v2. */ " < < endl
< < " plhs[2] = mxCreateDoubleMatrix( " < < NNZDerivatives [ 1 ] < < " , " < < 3
< < " , mxREAL); " < < endl
< < " v2 = mxGetPr(plhs[2]); " < < endl
< < " } " < < endl
< < endl
< < " /* Call the C subroutines. */ " < < endl
< < " Static(y, x, nb_row_x, params, residual, g1, v2); " < < endl
< < " } " < < endl < < endl ;
2009-07-07 16:20:48 +02:00
output . close ( ) ;
2009-04-14 16:39:53 +02:00
}
2009-12-16 14:21:31 +01:00
void
2015-07-27 15:59:13 +02:00
StaticModel : : writeStaticJuliaFile ( const string & basename ) const
2015-07-27 15:33:38 +02:00
{
2015-07-27 15:59:13 +02:00
string filename = basename + " Static.jl " ;
ofstream output ;
output . open ( filename . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2015-08-21 12:00:47 +02:00
output < < " module " < < basename < < " Static " < < endl
< < " # " < < endl
< < " # NB: this file was automatically generated by Dynare " < < endl
< < " # from " < < basename < < " .mod " < < endl
< < " # " < < endl
2015-08-18 15:45:23 +02:00
< < " using Utils " < < endl < < endl
2015-08-19 11:28:25 +02:00
< < " export static! " < < endl < < endl ;
2015-07-27 15:59:13 +02:00
writeStaticModel ( output , false , true ) ;
2015-08-18 15:41:27 +02:00
output < < " end " < < endl ;
2015-07-27 15:33:38 +02:00
}
void
2015-07-27 15:59:13 +02:00
StaticModel : : writeStaticFile ( const string & basename , bool block , bool bytecode , bool use_dll , bool julia ) const
2009-12-16 18:13:23 +01:00
{
int r ;
2009-04-14 16:39:53 +02:00
2009-12-16 18:13:23 +01:00
//assert(block);
2009-04-17 18:26:23 +02:00
2009-12-16 14:21:31 +01:00
# ifdef _WIN32
2009-12-16 18:13:23 +01:00
r = mkdir ( basename . c_str ( ) ) ;
2009-12-16 14:21:31 +01:00
# else
2009-12-16 18:13:23 +01:00
r = mkdir ( basename . c_str ( ) , 0777 ) ;
2009-12-16 14:21:31 +01:00
# endif
2009-12-16 18:13:23 +01:00
if ( r < 0 & & errno ! = EEXIST )
{
perror ( " ERROR " ) ;
exit ( EXIT_FAILURE ) ;
}
if ( block & & bytecode )
2010-10-11 19:21:32 +02:00
writeModelEquationsCode_Block ( basename + " _static " , basename , map_idx , map_idx2 ) ;
2010-01-22 11:03:29 +01:00
else if ( ! block & & bytecode )
writeModelEquationsCode ( basename + " _static " , basename , map_idx ) ;
2009-12-16 18:13:23 +01:00
else if ( block & & ! bytecode )
{
chdir ( basename . c_str ( ) ) ;
writeModelEquationsOrdered_M ( basename + " _static " ) ;
chdir ( " .. " ) ;
writeStaticBlockMFSFile ( basename ) ;
}
2011-08-18 12:44:11 +02:00
else if ( use_dll )
writeStaticCFile ( basename ) ;
2015-07-27 15:59:13 +02:00
else if ( julia )
writeStaticJuliaFile ( basename ) ;
2009-12-16 18:13:23 +01:00
else
writeStaticMFile ( basename ) ;
2016-04-04 17:11:03 +02:00
writeSetAuxiliaryVariables ( basename , julia ) ;
2009-12-16 18:13:23 +01:00
}
2009-04-27 19:15:14 +02:00
2016-08-12 11:50:57 +02:00
bool
StaticModel : : exoPresentInEqs ( ) const
{
for ( int i = 0 ; i < ( int ) equations . size ( ) ; i + + )
if ( equations [ i ] - > containsExogenous ( ) )
return true ;
return false ;
}
2009-04-27 19:15:14 +02:00
void
2009-12-16 14:21:31 +01:00
StaticModel : : writeStaticBlockMFSFile ( const string & basename ) const
2009-04-27 19:15:14 +02:00
{
2009-12-16 18:13:23 +01:00
string filename = basename + " _static.m " ;
2009-04-27 19:15:14 +02:00
2009-12-16 14:21:31 +01:00
ofstream output ;
output . open ( filename . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! output . is_open ( ) )
2009-04-27 19:15:14 +02:00
{
2009-12-16 14:21:31 +01:00
cerr < < " ERROR: Can't open file " < < filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
2009-04-27 19:15:14 +02:00
}
2009-12-16 14:21:31 +01:00
string func_name = basename + " _static " ;
2009-06-30 17:07:09 +02:00
2010-09-17 16:53:27 +02:00
output < < " function [residual, g1, y, var_index] = " < < func_name < < " (nblock, y, x, params) " < < endl
2009-12-16 14:21:31 +01:00
< < " residual = []; " < < endl
< < " g1 = []; " < < endl
2010-09-17 16:53:27 +02:00
< < " var_index = []; \n " < < endl
2009-12-16 14:21:31 +01:00
< < " switch nblock " < < endl ;
2009-06-30 17:07:09 +02:00
2009-12-16 14:21:31 +01:00
unsigned int nb_blocks = getNbBlocks ( ) ;
2009-06-30 17:07:09 +02:00
2009-12-16 18:13:23 +01:00
for ( int b = 0 ; b < ( int ) nb_blocks ; b + + )
2009-06-30 17:07:09 +02:00
{
2009-12-16 14:21:31 +01:00
set < int > local_var ;
2009-04-27 19:15:14 +02:00
2009-12-16 14:21:31 +01:00
output < < " case " < < b + 1 < < endl ;
2009-04-27 19:15:14 +02:00
2009-12-16 14:21:31 +01:00
BlockSimulationType simulation_type = getBlockSimulationType ( b ) ;
2009-12-16 18:13:23 +01:00
if ( simulation_type = = EVALUATE_BACKWARD | | simulation_type = = EVALUATE_FORWARD )
2010-09-17 16:53:27 +02:00
{
output < < " y_tmp = " < < func_name < < " _ " < < b + 1 < < " (y, x, params); \n " ;
ostringstream tmp ;
2011-02-04 16:25:38 +01:00
for ( int i = 0 ; i < ( int ) getBlockSize ( b ) ; i + + )
2010-09-17 16:53:27 +02:00
tmp < < " " < < getBlockVariableID ( b , i ) + 1 ;
output < < " var_index = [ " < < tmp . str ( ) < < " ]; \n " ;
output < < " residual = y(var_index) - y_tmp(var_index); \n " ;
output < < " y = y_tmp; \n " ;
}
2009-12-16 18:13:23 +01:00
else
output < < " [residual, y, g1] = " < < func_name < < " _ " < < b + 1 < < " (y, x, params); \n " ;
2010-09-17 16:53:27 +02:00
2009-04-28 19:11:48 +02:00
}
2009-12-16 14:21:31 +01:00
output < < " end " < < endl
< < " end " < < endl ;
output . close ( ) ;
}
2009-04-28 19:11:48 +02:00
2009-12-16 14:21:31 +01:00
void
StaticModel : : writeOutput ( ostream & output , bool block ) const
{
if ( ! block )
return ;
2009-06-30 17:07:09 +02:00
2009-12-16 14:21:31 +01:00
unsigned int nb_blocks = getNbBlocks ( ) ;
2012-06-06 16:36:56 +02:00
for ( int b = 0 ; b < ( int ) nb_blocks ; b + + )
{
BlockSimulationType simulation_type = getBlockSimulationType ( b ) ;
unsigned int block_size = getBlockSize ( b ) ;
ostringstream tmp_s , tmp_s_eq ;
tmp_s . str ( " " ) ;
tmp_s_eq . str ( " " ) ;
2012-06-08 19:00:57 +02:00
for ( unsigned int i = 0 ; i < block_size ; i + + )
2012-06-06 16:36:56 +02:00
{
tmp_s < < " " < < getBlockVariableID ( b , i ) + 1 ;
tmp_s_eq < < " " < < getBlockEquationID ( b , i ) + 1 ;
}
output < < " block_structure_stat.block( " < < b + 1 < < " ).Simulation_Type = " < < simulation_type < < " ; \n " ;
output < < " block_structure_stat.block( " < < b + 1 < < " ).endo_nbr = " < < block_size < < " ; \n " ;
output < < " block_structure_stat.block( " < < b + 1 < < " ).mfs = " < < getBlockMfs ( block ) < < " ; \n " ;
output < < " block_structure_stat.block( " < < b + 1 < < " ).equation = [ " < < tmp_s_eq . str ( ) < < " ]; \n " ;
output < < " block_structure_stat.block( " < < b + 1 < < " ).variable = [ " < < tmp_s . str ( ) < < " ]; \n " ;
}
output < < " M_.block_structure_stat.block = block_structure_stat.block; \n " ;
string cst_s ;
int nb_endo = symbol_table . endo_nbr ( ) ;
output < < " M_.block_structure_stat.variable_reordered = [ " ;
for ( int i = 0 ; i < nb_endo ; i + + )
output < < " " < < variable_reordered [ i ] + 1 ;
output < < " ]; \n " ;
output < < " M_.block_structure_stat.equation_reordered = [ " ;
for ( int i = 0 ; i < nb_endo ; i + + )
output < < " " < < equation_reordered [ i ] + 1 ;
output < < " ]; \n " ;
map < pair < int , int > , int > row_incidence ;
for ( first_derivatives_t : : const_iterator it = first_derivatives . begin ( ) ;
it ! = first_derivatives . end ( ) ; it + + )
{
int deriv_id = it - > first . second ;
if ( getTypeByDerivID ( deriv_id ) = = eEndogenous )
{
int eq = it - > first . first ;
int symb = getSymbIDByDerivID ( deriv_id ) ;
int var = symbol_table . getTypeSpecificID ( symb ) ;
//int lag = getLagByDerivID(deriv_id);
row_incidence [ make_pair ( eq , var ) ] = 1 ;
}
}
output < < " M_.block_structure_stat.incidence.sparse_IM = [ " ;
for ( map < pair < int , int > , int > : : const_iterator it = row_incidence . begin ( ) ; it ! = row_incidence . end ( ) ; it + + )
{
output < < it - > first . first + 1 < < " " < < it - > first . second + 1 < < " ; \n " ;
}
output < < " ]; \n " ;
2009-04-28 19:11:48 +02:00
}
2009-12-16 14:21:31 +01:00
SymbolType
StaticModel : : getTypeByDerivID ( int deriv_id ) const throw ( UnknownDerivIDException )
{
2012-11-29 18:07:48 +01:00
if ( deriv_id < symbol_table . endo_nbr ( ) )
return eEndogenous ;
else if ( deriv_id < symbol_table . endo_nbr ( ) + symbol_table . param_nbr ( ) )
return eParameter ;
else
throw UnknownDerivIDException ( ) ;
2009-04-27 19:15:14 +02:00
}
2009-04-30 15:14:33 +02:00
2009-12-16 14:21:31 +01:00
int
StaticModel : : getLagByDerivID ( int deriv_id ) const throw ( UnknownDerivIDException )
2009-04-30 15:14:33 +02:00
{
2009-12-16 14:21:31 +01:00
return 0 ;
2009-04-30 15:14:33 +02:00
}
2009-06-30 17:07:09 +02:00
2009-12-16 14:21:31 +01:00
int
StaticModel : : getSymbIDByDerivID ( int deriv_id ) const throw ( UnknownDerivIDException )
2009-06-30 17:07:09 +02:00
{
2012-11-29 18:07:48 +01:00
if ( deriv_id < symbol_table . endo_nbr ( ) )
return symbol_table . getID ( eEndogenous , deriv_id ) ;
else if ( deriv_id < symbol_table . endo_nbr ( ) + symbol_table . param_nbr ( ) )
return symbol_table . getID ( eParameter , deriv_id - symbol_table . endo_nbr ( ) ) ;
else
throw UnknownDerivIDException ( ) ;
2009-06-30 17:07:09 +02:00
}
2009-12-16 14:21:31 +01:00
int
StaticModel : : getDerivID ( int symb_id , int lag ) const throw ( UnknownDerivIDException )
2009-06-30 17:07:09 +02:00
{
2009-12-16 14:21:31 +01:00
if ( symbol_table . getType ( symb_id ) = = eEndogenous )
2012-11-29 18:07:48 +01:00
return symbol_table . getTypeSpecificID ( symb_id ) ;
else if ( symbol_table . getType ( symb_id ) = = eParameter )
return symbol_table . getTypeSpecificID ( symb_id ) + symbol_table . endo_nbr ( ) ;
2009-12-16 14:21:31 +01:00
else
return - 1 ;
}
2009-06-30 17:07:09 +02:00
2012-11-29 18:07:48 +01:00
void
StaticModel : : addAllParamDerivId ( set < int > & deriv_id_set )
{
for ( int i = 0 ; i < symbol_table . param_nbr ( ) ; i + + )
deriv_id_set . insert ( i + symbol_table . endo_nbr ( ) ) ;
}
2009-12-16 14:21:31 +01:00
map < pair < pair < int , pair < int , int > > , pair < int , int > > , int >
StaticModel : : get_Derivatives ( int block )
{
map < pair < pair < int , pair < int , int > > , pair < int , int > > , int > Derivatives ;
Derivatives . clear ( ) ;
int block_size = getBlockSize ( block ) ;
int block_nb_recursive = block_size - getBlockMfs ( block ) ;
int lag = 0 ;
2009-12-16 18:13:23 +01:00
for ( int eq = 0 ; eq < block_size ; eq + + )
2009-06-30 17:07:09 +02:00
{
2009-12-16 14:21:31 +01:00
int eqr = getBlockEquationID ( block , eq ) ;
2009-12-16 18:13:23 +01:00
for ( int var = 0 ; var < block_size ; var + + )
2009-06-30 17:07:09 +02:00
{
2009-12-16 14:21:31 +01:00
int varr = getBlockVariableID ( block , var ) ;
2009-12-16 18:13:23 +01:00
if ( dynamic_jacobian . find ( make_pair ( lag , make_pair ( eqr , varr ) ) ) ! = dynamic_jacobian . end ( ) )
2009-12-16 14:21:31 +01:00
{
bool OK = true ;
map < pair < pair < int , pair < int , int > > , pair < int , int > > , int > : : const_iterator its = Derivatives . find ( make_pair ( make_pair ( lag , make_pair ( eq , var ) ) , make_pair ( eqr , varr ) ) ) ;
2009-12-16 18:13:23 +01:00
if ( its ! = Derivatives . end ( ) )
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
if ( its - > second = = 2 )
OK = false ;
2009-12-16 14:21:31 +01:00
}
2009-06-30 17:07:09 +02:00
2009-12-16 18:13:23 +01:00
if ( OK )
2009-12-16 14:21:31 +01:00
{
2010-04-28 16:03:32 +02:00
if ( getBlockEquationType ( block , eq ) = = E_EVALUATE_S & & eq < block_nb_recursive )
2009-12-16 14:21:31 +01:00
//It's a normalized equation, we have to recompute the derivative using chain rule derivative function
Derivatives [ make_pair ( make_pair ( lag , make_pair ( eq , var ) ) , make_pair ( eqr , varr ) ) ] = 1 ;
else
//It's a feedback equation we can use the derivatives
Derivatives [ make_pair ( make_pair ( lag , make_pair ( eq , var ) ) , make_pair ( eqr , varr ) ) ] = 0 ;
}
2009-12-16 18:13:23 +01:00
if ( var < block_nb_recursive )
2009-12-16 14:21:31 +01:00
{
int eqs = getBlockEquationID ( block , var ) ;
2009-12-16 18:13:23 +01:00
for ( int vars = block_nb_recursive ; vars < block_size ; vars + + )
2009-12-16 14:21:31 +01:00
{
int varrs = getBlockVariableID ( block , vars ) ;
//A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation)
2009-12-16 18:13:23 +01:00
if ( Derivatives . find ( make_pair ( make_pair ( lag , make_pair ( var , vars ) ) , make_pair ( eqs , varrs ) ) ) ! = Derivatives . end ( ) )
2009-12-16 14:21:31 +01:00
Derivatives [ make_pair ( make_pair ( lag , make_pair ( eq , vars ) ) , make_pair ( eqr , varrs ) ) ] = 2 ;
}
}
}
2009-06-30 17:07:09 +02:00
}
}
2009-12-16 14:21:31 +01:00
2009-12-16 18:13:23 +01:00
return ( Derivatives ) ;
2009-06-30 17:07:09 +02:00
}
2009-07-07 16:20:48 +02:00
void
2010-09-16 19:00:48 +02:00
StaticModel : : computeChainRuleJacobian ( blocks_derivatives_t & blocks_derivatives )
2009-07-07 16:20:48 +02:00
{
2010-09-16 19:18:45 +02:00
map < int , expr_t > recursive_variables ;
2009-12-16 14:21:31 +01:00
unsigned int nb_blocks = getNbBlocks ( ) ;
2010-09-16 19:00:48 +02:00
blocks_derivatives = blocks_derivatives_t ( nb_blocks ) ;
2009-12-16 18:13:23 +01:00
for ( unsigned int block = 0 ; block < nb_blocks ; block + + )
2009-07-07 16:20:48 +02:00
{
2010-09-16 19:00:48 +02:00
block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives ;
2009-12-16 14:21:31 +01:00
recursive_variables . clear ( ) ;
BlockSimulationType simulation_type = getBlockSimulationType ( block ) ;
int block_size = getBlockSize ( block ) ;
int block_nb_mfs = getBlockMfs ( block ) ;
int block_nb_recursives = block_size - block_nb_mfs ;
2010-04-28 16:03:32 +02:00
if ( simulation_type = = SOLVE_TWO_BOUNDARIES_COMPLETE | | simulation_type = = SOLVE_TWO_BOUNDARIES_SIMPLE )
2009-07-10 16:22:40 +02:00
{
2010-09-16 19:00:48 +02:00
blocks_derivatives . push_back ( block_derivatives_equation_variable_laglead_nodeid_t ( 0 ) ) ;
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < block_nb_recursives ; i + + )
2009-12-16 14:21:31 +01:00
{
if ( getBlockEquationType ( block , i ) = = E_EVALUATE_S )
2010-09-16 19:18:45 +02:00
recursive_variables [ getDerivID ( symbol_table . getID ( eEndogenous , getBlockVariableID ( block , i ) ) , 0 ) ] = getBlockEquationRenormalizedExpr ( block , i ) ;
2009-12-16 14:21:31 +01:00
else
2010-09-16 19:18:45 +02:00
recursive_variables [ getDerivID ( symbol_table . getID ( eEndogenous , getBlockVariableID ( block , i ) ) , 0 ) ] = getBlockEquationExpr ( block , i ) ;
2009-12-16 14:21:31 +01:00
}
map < pair < pair < int , pair < int , int > > , pair < int , int > > , int > Derivatives = get_Derivatives ( block ) ;
map < pair < pair < int , pair < int , int > > , pair < int , int > > , int > : : const_iterator it = Derivatives . begin ( ) ;
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < ( int ) Derivatives . size ( ) ; i + + )
2009-07-10 16:22:40 +02:00
{
2009-12-16 14:21:31 +01:00
int Deriv_type = it - > second ;
pair < pair < int , pair < int , int > > , pair < int , int > > it_l ( it - > first ) ;
it + + ;
int lag = it_l . first . first ;
int eq = it_l . first . second . first ;
int var = it_l . first . second . second ;
int eqr = it_l . second . first ;
int varr = it_l . second . second ;
2009-12-16 18:13:23 +01:00
if ( Deriv_type = = 0 )
2009-12-16 14:21:31 +01:00
first_chain_rule_derivatives [ make_pair ( eqr , make_pair ( varr , lag ) ) ] = first_derivatives [ make_pair ( eqr , getDerivID ( symbol_table . getID ( eEndogenous , varr ) , lag ) ) ] ;
else if ( Deriv_type = = 1 )
first_chain_rule_derivatives [ make_pair ( eqr , make_pair ( varr , lag ) ) ] = ( equation_type_and_normalized_equation [ eqr ] . second ) - > getChainRuleDerivative ( getDerivID ( symbol_table . getID ( eEndogenous , varr ) , lag ) , recursive_variables ) ;
else if ( Deriv_type = = 2 )
2009-07-10 16:22:40 +02:00
{
2009-12-16 18:13:23 +01:00
if ( getBlockEquationType ( block , eq ) = = E_EVALUATE_S & & eq < block_nb_recursives )
2009-12-16 14:21:31 +01:00
first_chain_rule_derivatives [ make_pair ( eqr , make_pair ( varr , lag ) ) ] = ( equation_type_and_normalized_equation [ eqr ] . second ) - > getChainRuleDerivative ( getDerivID ( symbol_table . getID ( eEndogenous , varr ) , lag ) , recursive_variables ) ;
else
first_chain_rule_derivatives [ make_pair ( eqr , make_pair ( varr , lag ) ) ] = equations [ eqr ] - > getChainRuleDerivative ( getDerivID ( symbol_table . getID ( eEndogenous , varr ) , lag ) , recursive_variables ) ;
2009-07-10 16:22:40 +02:00
}
2009-12-16 18:13:23 +01:00
tmp_derivatives . push_back ( make_pair ( make_pair ( eq , var ) , make_pair ( lag , first_chain_rule_derivatives [ make_pair ( eqr , make_pair ( varr , lag ) ) ] ) ) ) ;
2009-07-10 16:22:40 +02:00
}
}
2010-10-11 19:21:32 +02:00
else
2009-07-07 16:20:48 +02:00
{
2010-09-16 19:00:48 +02:00
blocks_derivatives . push_back ( block_derivatives_equation_variable_laglead_nodeid_t ( 0 ) ) ;
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < block_nb_recursives ; i + + )
2009-07-07 17:58:17 +02:00
{
2009-12-16 14:21:31 +01:00
if ( getBlockEquationType ( block , i ) = = E_EVALUATE_S )
2010-09-16 19:18:45 +02:00
recursive_variables [ getDerivID ( symbol_table . getID ( eEndogenous , getBlockVariableID ( block , i ) ) , 0 ) ] = getBlockEquationRenormalizedExpr ( block , i ) ;
2009-12-16 14:21:31 +01:00
else
2010-09-16 19:18:45 +02:00
recursive_variables [ getDerivID ( symbol_table . getID ( eEndogenous , getBlockVariableID ( block , i ) ) , 0 ) ] = getBlockEquationExpr ( block , i ) ;
2009-12-16 14:21:31 +01:00
}
2009-12-16 18:13:23 +01:00
for ( int eq = block_nb_recursives ; eq < block_size ; eq + + )
2009-12-16 14:21:31 +01:00
{
int eqr = getBlockEquationID ( block , eq ) ;
2009-12-16 18:13:23 +01:00
for ( int var = block_nb_recursives ; var < block_size ; var + + )
2009-12-16 14:21:31 +01:00
{
int varr = getBlockVariableID ( block , var ) ;
2010-09-16 19:18:45 +02:00
expr_t d1 = equations [ eqr ] - > getChainRuleDerivative ( getDerivID ( symbol_table . getID ( eEndogenous , varr ) , 0 ) , recursive_variables ) ;
2009-12-16 14:21:31 +01:00
if ( d1 = = Zero )
continue ;
first_chain_rule_derivatives [ make_pair ( eqr , make_pair ( varr , 0 ) ) ] = d1 ;
tmp_derivatives . push_back (
2009-12-16 18:13:23 +01:00
make_pair ( make_pair ( eq , var ) , make_pair ( 0 , first_chain_rule_derivatives [ make_pair ( eqr , make_pair ( varr , 0 ) ) ] ) ) ) ;
2009-12-16 14:21:31 +01:00
}
2009-07-07 17:58:17 +02:00
}
2009-07-07 16:20:48 +02:00
}
2009-12-16 14:21:31 +01:00
blocks_derivatives [ block ] = tmp_derivatives ;
2009-07-07 16:20:48 +02:00
}
}
void
2009-12-16 14:21:31 +01:00
StaticModel : : collect_block_first_order_derivatives ( )
2009-07-07 16:20:48 +02:00
{
2009-12-16 14:21:31 +01:00
//! vector for an equation or a variable indicates the block number
vector < int > equation_2_block , variable_2_block ;
unsigned int nb_blocks = getNbBlocks ( ) ;
equation_2_block = vector < int > ( equation_reordered . size ( ) ) ;
variable_2_block = vector < int > ( variable_reordered . size ( ) ) ;
2009-12-16 18:13:23 +01:00
for ( unsigned int block = 0 ; block < nb_blocks ; block + + )
2009-07-07 16:20:48 +02:00
{
2009-12-16 14:21:31 +01:00
unsigned int block_size = getBlockSize ( block ) ;
2009-12-16 18:13:23 +01:00
for ( unsigned int i = 0 ; i < block_size ; i + + )
2009-12-16 14:21:31 +01:00
{
equation_2_block [ getBlockEquationID ( block , i ) ] = block ;
variable_2_block [ getBlockVariableID ( block , i ) ] = block ;
}
2009-07-07 16:20:48 +02:00
}
2010-09-16 19:00:48 +02:00
derivative_endo = vector < derivative_t > ( nb_blocks ) ;
2009-12-16 18:13:23 +01:00
endo_max_leadlag_block = vector < pair < int , int > > ( nb_blocks , make_pair ( 0 , 0 ) ) ;
max_leadlag_block = vector < pair < int , int > > ( nb_blocks , make_pair ( 0 , 0 ) ) ;
2010-09-16 19:00:48 +02:00
for ( first_derivatives_t : : iterator it2 = first_derivatives . begin ( ) ;
2009-12-16 14:21:31 +01:00
it2 ! = first_derivatives . end ( ) ; it2 + + )
2009-07-13 15:44:46 +02:00
{
2009-12-16 14:21:31 +01:00
int eq = it2 - > first . first ;
int var = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( it2 - > first . second ) ) ;
int lag = 0 ;
int block_eq = equation_2_block [ eq ] ;
int block_var = variable_2_block [ var ] ;
max_leadlag_block [ block_eq ] = make_pair ( 0 , 0 ) ;
max_leadlag_block [ block_eq ] = make_pair ( 0 , 0 ) ;
endo_max_leadlag_block [ block_eq ] = make_pair ( 0 , 0 ) ;
endo_max_leadlag_block [ block_eq ] = make_pair ( 0 , 0 ) ;
2010-09-16 19:00:48 +02:00
derivative_t tmp_derivative ;
lag_var_t lag_var ;
2009-12-16 14:21:31 +01:00
if ( getTypeByDerivID ( it2 - > first . second ) = = eEndogenous & & block_eq = = block_var )
{
tmp_derivative = derivative_endo [ block_eq ] ;
tmp_derivative [ make_pair ( lag , make_pair ( eq , var ) ) ] = first_derivatives [ make_pair ( eq , getDerivID ( symbol_table . getID ( eEndogenous , var ) , lag ) ) ] ;
derivative_endo [ block_eq ] = tmp_derivative ;
}
2009-07-13 15:44:46 +02:00
}
}
2009-07-07 17:58:17 +02:00
void
2009-12-16 14:21:31 +01:00
StaticModel : : writeChainRuleDerivative ( ostream & output , int eqr , int varr , int lag ,
2009-12-16 18:13:23 +01:00
ExprNodeOutputType output_type ,
2010-09-16 19:00:48 +02:00
const temporary_terms_t & temporary_terms ) const
2009-07-07 17:58:17 +02:00
{
2010-09-16 19:18:45 +02:00
map < pair < int , pair < int , int > > , expr_t > : : const_iterator it = first_chain_rule_derivatives . find ( make_pair ( eqr , make_pair ( varr , lag ) ) ) ;
2009-12-16 14:21:31 +01:00
if ( it ! = first_chain_rule_derivatives . end ( ) )
( it - > second ) - > writeOutput ( output , output_type , temporary_terms ) ;
else
output < < 0 ;
}
2009-07-07 17:58:17 +02:00
2009-12-16 14:21:31 +01:00
void
StaticModel : : writeLatexFile ( const string & basename ) const
2009-12-16 18:13:23 +01:00
{
2015-07-15 08:58:15 +02:00
writeLatexModelFile ( basename + " _static " , oLatexStaticModel ) ;
2009-12-16 18:13:23 +01:00
}
2009-07-07 17:58:17 +02:00
2009-09-30 17:10:31 +02:00
void
2010-04-27 17:04:52 +02:00
StaticModel : : writeAuxVarInitval ( ostream & output , ExprNodeOutputType output_type ) const
2009-09-30 17:10:31 +02:00
{
2011-09-28 22:08:45 +02:00
for ( int i = 0 ; i < ( int ) aux_equations . size ( ) ; i + + )
2009-09-30 17:10:31 +02:00
{
2010-04-27 17:04:52 +02:00
dynamic_cast < ExprNode * > ( aux_equations [ i ] ) - > writeOutput ( output , output_type ) ;
2009-09-30 17:10:31 +02:00
output < < " ; " < < endl ;
}
}
2011-07-24 20:52:03 +02:00
2016-04-04 17:11:03 +02:00
void StaticModel : : writeSetAuxiliaryVariables ( const string & basename , const bool julia ) const
2011-07-24 20:52:03 +02:00
{
2016-04-04 17:11:03 +02:00
2011-07-24 20:52:03 +02:00
string func_name = basename + " _set_auxiliary_variables " ;
2015-08-17 15:36:18 +02:00
string filename = julia ? func_name + " .jl " : func_name + " .m " ;
string comment = julia ? " # " : " % " ;
2011-07-24 20:52:03 +02:00
ofstream output ;
output . open ( filename . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
output < < " function y = " < < func_name + " (y, x, params) " < < endl
2015-08-17 15:36:18 +02:00
< < comment < < endl
< < comment < < " Status : Computes static model for Dynare " < < endl
< < comment < < endl
< < comment < < " Warning : this file is generated automatically by Dynare " < < endl
< < comment < < " from model file (.mod) " < < endl
2011-07-24 20:52:03 +02:00
< < endl ;
2016-04-04 17:11:03 +02:00
writeAuxVarRecursiveDefinitions ( output , oMatlabStaticModel ) ;
}
2015-05-12 16:36:03 +02:00
2016-04-04 17:11:03 +02:00
void
StaticModel : : writeAuxVarRecursiveDefinitions ( ostream & output , ExprNodeOutputType output_type ) const
{
deriv_node_temp_terms_t tef_terms ;
2016-04-08 21:45:44 +02:00
temporary_terms_t temporary_terms ;
for ( int i = 0 ; i < ( int ) aux_equations . size ( ) ; i + + )
if ( dynamic_cast < ExprNode * > ( aux_equations [ i ] ) - > containsExternalFunction ( ) )
dynamic_cast < ExprNode * > ( aux_equations [ i ] ) - > writeExternalFunctionOutput ( output , oMatlabStaticModel ,
temporary_terms , tef_terms ) ;
2011-09-28 22:08:45 +02:00
for ( int i = 0 ; i < ( int ) aux_equations . size ( ) ; i + + )
2011-07-24 20:52:03 +02:00
{
2016-04-08 21:45:44 +02:00
dynamic_cast < ExprNode * > ( aux_equations [ i ] - > substituteStaticAuxiliaryDefinition ( ) ) - > writeOutput ( output , output_type , temporary_terms , tef_terms ) ;
2011-07-24 20:52:03 +02:00
output < < " ; " < < endl ;
}
}
2012-11-29 18:07:48 +01:00
void
2015-07-28 12:29:42 +02:00
StaticModel : : writeParamsDerivativesFile ( const string & basename , bool julia ) const
2012-11-29 18:07:48 +01:00
{
if ( ! residuals_params_derivatives . size ( )
& & ! residuals_params_second_derivatives . size ( )
& & ! jacobian_params_derivatives . size ( )
& & ! jacobian_params_second_derivatives . size ( )
& & ! hessian_params_derivatives . size ( ) )
return ;
2015-07-28 12:29:42 +02:00
ExprNodeOutputType output_type = ( julia ? oJuliaStaticModel : oMatlabStaticModel ) ;
2017-01-05 17:31:48 +01:00
ostringstream model_local_vars_output ; // Used for storing model local vars
ostringstream model_output ; // Used for storing model
ostringstream jacobian_output ; // Used for storing jacobian equations
ostringstream hessian_output ; // Used for storing Hessian equations
ostringstream hessian1_output ; // Used for storing Hessian equations
ostringstream third_derivs_output ; // Used for storing third order derivatives equations
ostringstream third_derivs1_output ; // Used for storing third order derivatives equations
2012-11-29 18:07:48 +01:00
deriv_node_temp_terms_t tef_terms ;
2017-01-05 17:31:48 +01:00
writeModelLocalVariables ( model_local_vars_output , output_type , tef_terms ) ;
2012-11-29 18:07:48 +01:00
2016-06-01 15:28:11 +02:00
temporary_terms_t temp_terms_empty ;
2017-01-05 17:31:48 +01:00
writeTemporaryTerms ( params_derivs_temporary_terms , temp_terms_empty , model_output , output_type , tef_terms ) ;
2012-11-29 18:07:48 +01:00
for ( first_derivatives_t : : const_iterator it = residuals_params_derivatives . begin ( ) ;
it ! = residuals_params_derivatives . end ( ) ; it + + )
{
int eq = it - > first . first ;
int param = it - > first . second ;
expr_t d1 = it - > second ;
int param_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( param ) ) + 1 ;
2017-01-05 17:31:48 +01:00
jacobian_output < < " rp " < < LEFT_ARRAY_SUBSCRIPT ( output_type )
< < eq + 1 < < " , " < < param_col
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " ;
d1 - > writeOutput ( jacobian_output , output_type , params_derivs_temporary_terms , tef_terms ) ;
jacobian_output < < " ; " < < endl ;
2012-11-29 18:07:48 +01:00
}
for ( second_derivatives_t : : const_iterator it = jacobian_params_derivatives . begin ( ) ;
it ! = jacobian_params_derivatives . end ( ) ; it + + )
{
int eq = it - > first . first ;
int var = it - > first . second . first ;
int param = it - > first . second . second ;
expr_t d2 = it - > second ;
int var_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( var ) ) + 1 ;
int param_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( param ) ) + 1 ;
2017-01-05 17:31:48 +01:00
hessian_output < < " gp " < < LEFT_ARRAY_SUBSCRIPT ( output_type )
< < eq + 1 < < " , " < < var_col < < " , " < < param_col
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " ;
d2 - > writeOutput ( hessian_output , output_type , params_derivs_temporary_terms , tef_terms ) ;
hessian_output < < " ; " < < endl ;
2012-11-29 18:07:48 +01:00
}
int i = 1 ;
for ( second_derivatives_t : : const_iterator it = residuals_params_second_derivatives . begin ( ) ;
it ! = residuals_params_second_derivatives . end ( ) ; + + it , i + + )
{
int eq = it - > first . first ;
int param1 = it - > first . second . first ;
int param2 = it - > first . second . second ;
expr_t d2 = it - > second ;
int param1_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( param1 ) ) + 1 ;
int param2_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( param2 ) ) + 1 ;
2017-01-05 17:31:48 +01:00
hessian1_output < < " rpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,1 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < eq + 1 < < " ; " < < endl
< < " rpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,2 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < param1_col < < " ; " < < endl
< < " rpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,3 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < param2_col < < " ; " < < endl
< < " rpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,4 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " ;
d2 - > writeOutput ( hessian1_output , output_type , params_derivs_temporary_terms , tef_terms ) ;
hessian1_output < < " ; " < < endl ;
2012-11-29 18:07:48 +01:00
}
i = 1 ;
for ( third_derivatives_t : : const_iterator it = jacobian_params_second_derivatives . begin ( ) ;
it ! = jacobian_params_second_derivatives . end ( ) ; + + it , i + + )
{
int eq = it - > first . first ;
int var = it - > first . second . first ;
int param1 = it - > first . second . second . first ;
int param2 = it - > first . second . second . second ;
expr_t d2 = it - > second ;
int var_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( var ) ) + 1 ;
int param1_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( param1 ) ) + 1 ;
int param2_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( param2 ) ) + 1 ;
2017-01-05 17:31:48 +01:00
third_derivs_output < < " gpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,1 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < eq + 1 < < " ; " < < endl
< < " gpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,2 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < var_col < < " ; " < < endl
< < " gpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,3 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < param1_col < < " ; " < < endl
< < " gpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,4 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < param2_col < < " ; " < < endl
< < " gpp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,5 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " ;
d2 - > writeOutput ( third_derivs_output , output_type , params_derivs_temporary_terms , tef_terms ) ;
third_derivs_output < < " ; " < < endl ;
2012-11-29 18:07:48 +01:00
}
i = 1 ;
for ( third_derivatives_t : : const_iterator it = hessian_params_derivatives . begin ( ) ;
it ! = hessian_params_derivatives . end ( ) ; + + it , i + + )
{
int eq = it - > first . first ;
int var1 = it - > first . second . first ;
int var2 = it - > first . second . second . first ;
int param = it - > first . second . second . second ;
expr_t d2 = it - > second ;
int var1_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( var1 ) ) + 1 ;
int var2_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( var2 ) ) + 1 ;
int param_col = symbol_table . getTypeSpecificID ( getSymbIDByDerivID ( param ) ) + 1 ;
2017-01-05 17:31:48 +01:00
third_derivs1_output < < " hp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,1 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < eq + 1 < < " ; " < < endl
< < " hp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,2 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < var1_col < < " ; " < < endl
< < " hp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,3 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < var2_col < < " ; " < < endl
< < " hp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,4 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " < < param_col < < " ; " < < endl
< < " hp " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) < < i < < " ,5 "
< < RIGHT_ARRAY_SUBSCRIPT ( output_type ) < < " = " ;
d2 - > writeOutput ( third_derivs1_output , output_type , params_derivs_temporary_terms , tef_terms ) ;
third_derivs1_output < < " ; " < < endl ;
2012-11-29 18:07:48 +01:00
}
2017-01-05 17:31:48 +01:00
ofstream paramsDerivsFile ;
string filename = julia ? basename + " StaticParamsDerivs.jl " : basename + " _static_params_derivs.m " ;
paramsDerivsFile . open ( filename . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! paramsDerivsFile . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
if ( ! julia )
{
// Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201
map < string , string > tmp_paren_vars ;
bool message_printed = false ;
fixNestedParenthesis ( model_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( model_local_vars_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( jacobian_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( hessian_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( hessian1_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( third_derivs_output , tmp_paren_vars , message_printed ) ;
fixNestedParenthesis ( third_derivs1_output , tmp_paren_vars , message_printed ) ;
paramsDerivsFile < < " function [rp, gp, rpp, gpp, hp] = " < < basename < < " _static_params_derivs(y, x, params) " < < endl
< < " % " < < endl
< < " % Status : Computes derivatives of the static model with respect to the parameters " < < endl
< < " % " < < endl
< < " % Inputs : " < < endl
< < " % y [M_.endo_nbr by 1] double vector of endogenous variables in declaration order " < < endl
< < " % x [M_.exo_nbr by 1] double vector of exogenous variables in declaration order " < < endl
< < " % params [M_.param_nbr by 1] double vector of parameter values in declaration order " < < endl
< < " % " < < endl
< < " % Outputs: " < < endl
< < " % rp [M_.eq_nbr by #params] double Jacobian matrix of static model equations with respect to parameters " < < endl
< < " % Dynare may prepend or append auxiliary equations, see M_.aux_vars " < < endl
< < " % gp [M_.endo_nbr by M_.endo_nbr by #params] double Derivative of the Jacobian matrix of the static model equations with respect to the parameters " < < endl
< < " % rows: variables in declaration order " < < endl
< < " % rows: equations in order of declaration " < < endl
< < " % rpp [#second_order_residual_terms by 4] double Hessian matrix of second derivatives of residuals with respect to parameters; " < < endl
< < " % rows: respective derivative term " < < endl
< < " % 1st column: equation number of the term appearing " < < endl
< < " % 2nd column: number of the first parameter in derivative " < < endl
< < " % 3rd column: number of the second parameter in derivative " < < endl
< < " % 4th column: value of the Hessian term " < < endl
< < " % gpp [#second_order_Jacobian_terms by 5] double Hessian matrix of second derivatives of the Jacobian with respect to the parameters; " < < endl
< < " % rows: respective derivative term " < < endl
< < " % 1st column: equation number of the term appearing " < < endl
< < " % 2nd column: column number of variable in Jacobian of the static model " < < endl
< < " % 3rd column: number of the first parameter in derivative " < < endl
< < " % 4th column: number of the second parameter in derivative " < < endl
< < " % 5th column: value of the Hessian term " < < endl
< < " % " < < endl
< < " % " < < endl
< < " % Warning : this file is generated automatically by Dynare " < < endl
< < " % from model file (.mod) " < < endl < < endl
< < model_local_vars_output . str ( )
< < model_output . str ( )
2017-02-24 11:20:54 +01:00
< < " rp = zeros( " < < equations . size ( ) < < " , "
2017-01-05 17:31:48 +01:00
< < symbol_table . param_nbr ( ) < < " ); " < < endl
< < jacobian_output . str ( )
2017-02-24 11:20:54 +01:00
< < " gp = zeros( " < < equations . size ( ) < < " , " < < symbol_table . endo_nbr ( ) < < " , "
2017-01-05 17:31:48 +01:00
< < symbol_table . param_nbr ( ) < < " ); " < < endl
< < hessian_output . str ( )
< < " if nargout >= 3 " < < endl
< < " rpp = zeros( " < < residuals_params_second_derivatives . size ( ) < < " ,4); " < < endl
< < hessian1_output . str ( )
< < " gpp = zeros( " < < jacobian_params_second_derivatives . size ( ) < < " ,5); " < < endl
< < third_derivs_output . str ( )
< < " end " < < endl
< < " if nargout >= 5 " < < endl
< < " hp = zeros( " < < hessian_params_derivatives . size ( ) < < " ,5); " < < endl
< < third_derivs1_output . str ( )
< < " end " < < endl
< < " end " < < endl ;
}
else
paramsDerivsFile < < " module " < < basename < < " StaticParamsDerivs " < < endl
< < " # " < < endl
< < " # NB: this file was automatically generated by Dynare " < < endl
< < " # from " < < basename < < " .mod " < < endl
< < " # " < < endl
< < " export params_derivs " < < endl < < endl
< < " function params_derivs(y, x, params) " < < endl
< < model_local_vars_output . str ( )
< < model_output . str ( )
2017-02-24 11:20:54 +01:00
< < " rp = zeros( " < < equations . size ( ) < < " , "
2017-01-05 17:31:48 +01:00
< < symbol_table . param_nbr ( ) < < " ); " < < endl
< < jacobian_output . str ( )
2017-02-24 11:20:54 +01:00
< < " gp = zeros( " < < equations . size ( ) < < " , " < < symbol_table . endo_nbr ( ) < < " , "
2017-01-05 17:31:48 +01:00
< < symbol_table . param_nbr ( ) < < " ); " < < endl
< < hessian_output . str ( )
< < " rpp = zeros( " < < residuals_params_second_derivatives . size ( ) < < " ,4); " < < endl
< < hessian1_output . str ( )
< < " gpp = zeros( " < < jacobian_params_second_derivatives . size ( ) < < " ,5); " < < endl
< < third_derivs_output . str ( )
< < " hp = zeros( " < < hessian_params_derivatives . size ( ) < < " ,5); " < < endl
< < third_derivs1_output . str ( )
< < " (rp, gp, rpp, gpp, hp) " < < endl
< < " end " < < endl
< < " end " < < endl ;
2012-11-29 18:07:48 +01:00
paramsDerivsFile . close ( ) ;
}