2009-04-14 16:39:53 +02:00
/*
2023-01-05 16:40:04 +01:00
* Copyright © 2003 - 2023 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
2021-06-09 16:52:20 +02:00
* along with Dynare . If not , see < https : //www.gnu.org/licenses/>.
2009-04-14 16:39:53 +02:00
*/
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-04-28 19:11:48 +02:00
# include <algorithm>
2021-04-23 17:30:38 +02:00
# include <sstream>
2022-09-13 15:43:13 +02:00
# include <numeric>
2009-12-16 14:21:31 +01:00
2018-06-27 15:01:31 +02:00
# include "StaticModel.hh"
2018-10-10 17:08:54 +02:00
# include "DynamicModel.hh"
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 ,
2018-09-14 17:04:06 +02:00
ExternalFunctionsTable & external_functions_table_arg ) :
2018-10-04 17:18:27 +02:00
ModelTree { symbol_table_arg , num_constants_arg , external_functions_table_arg }
2009-12-16 14:21:31 +01:00
{
}
2009-04-14 16:47:57 +02:00
2018-10-09 18:27:19 +02:00
StaticModel : : StaticModel ( const StaticModel & m ) :
2020-05-06 17:13:47 +02:00
ModelTree { m }
2018-10-09 18:27:19 +02:00
{
}
StaticModel &
StaticModel : : operator = ( const StaticModel & m )
{
ModelTree : : operator = ( m ) ;
return * this ;
}
2018-10-10 17:08:54 +02:00
StaticModel : : StaticModel ( const DynamicModel & m ) :
2019-12-20 16:59:30 +01:00
ModelTree { m . symbol_table , m . num_constants , m . external_functions_table }
2018-10-10 17:08:54 +02:00
{
// Convert model local variables (need to be done first)
for ( int it : m . local_variables_vector )
2022-09-27 12:51:22 +02:00
AddLocalVariable ( it , m . local_variables_table . at ( it ) - > toStatic ( * this ) ) ;
2018-10-10 17:08:54 +02:00
// Convert equations
int static_only_index = 0 ;
2020-02-20 15:29:10 +01:00
set < int > dynamic_equations = m . equation_tags . getDynamicEqns ( ) ;
2018-10-10 17:08:54 +02:00
for ( int i = 0 ; i < static_cast < int > ( m . equations . size ( ) ) ; i + + )
2020-02-20 15:29:10 +01:00
try
{
// If equation is dynamic, replace it by an equation marked [static]
2022-05-04 16:01:34 +02:00
if ( dynamic_equations . contains ( i ) )
2018-10-10 17:08:54 +02:00
{
2020-02-20 15:29:10 +01:00
auto [ static_only_equations ,
static_only_equations_lineno ,
static_only_equations_equation_tags ] = m . getStaticOnlyEquationsInfo ( ) ;
addEquation ( static_only_equations [ static_only_index ] - > toStatic ( * this ) ,
static_only_equations_lineno [ static_only_index ] ,
static_only_equations_equation_tags . getTagsByEqn ( static_only_index ) ) ;
static_only_index + + ;
2018-10-10 17:08:54 +02:00
}
2020-02-20 15:29:10 +01:00
else
addEquation ( m . equations [ i ] - > toStatic ( * this ) ,
m . equations_lineno [ i ] ,
m . equation_tags . getTagsByEqn ( i ) ) ;
}
catch ( DataTree : : DivisionByZeroException )
{
2021-12-10 12:44:42 +01:00
cerr < < " ...division by zero error encountered when converting equation " < < i < < " to static " < < endl ;
2020-02-20 15:29:10 +01:00
exit ( EXIT_FAILURE ) ;
}
2018-10-10 17:08:54 +02:00
// Convert auxiliary equations
for ( auto aux_eq : m . aux_equations )
addAuxEquation ( aux_eq - > toStatic ( * this ) ) ;
2019-12-02 19:21:14 +01:00
2023-01-17 16:38:03 +01:00
cutoff = m . cutoff ;
mfs = m . mfs ;
2019-12-02 19:21:14 +01:00
user_set_add_flags = m . user_set_add_flags ;
user_set_subst_flags = m . user_set_subst_flags ;
user_set_add_libs = m . user_set_add_libs ;
user_set_subst_libs = m . user_set_subst_libs ;
user_set_compiler = m . user_set_compiler ;
2018-10-10 17:08:54 +02:00
}
2009-12-16 14:21:31 +01:00
void
2020-05-19 17:43:35 +02:00
StaticModel : : writeStaticBytecode ( const string & basename ) const
2010-01-22 11:03:29 +01:00
{
2022-07-13 13:04:10 +02:00
// First write the .bin file
int u_count_int { writeBytecodeBinFile ( basename + " /model/bytecode/static.bin " , false ) } ;
BytecodeWriter code_file { basename + " /model/bytecode/static.cod " } ;
2022-09-13 15:43:13 +02:00
vector < int > eq_idx ( equations . size ( ) ) ;
iota ( eq_idx . begin ( ) , eq_idx . end ( ) , 0 ) ;
vector < int > endo_idx ( symbol_table . endo_nbr ( ) ) ;
iota ( endo_idx . begin ( ) , endo_idx . end ( ) , 0 ) ;
2022-07-13 13:04:10 +02:00
// Declare temporary terms and the (single) block
code_file < < FDIMST_ { static_cast < int > ( temporary_terms_derivatives [ 0 ] . size ( )
+ temporary_terms_derivatives [ 1 ] . size ( ) ) }
< < FBEGINBLOCK_ { symbol_table . endo_nbr ( ) ,
BlockSimulationType : : solveForwardComplete ,
0 ,
symbol_table . endo_nbr ( ) ,
2022-09-13 15:43:13 +02:00
endo_idx ,
eq_idx ,
2022-07-13 13:04:10 +02:00
false ,
symbol_table . endo_nbr ( ) ,
u_count_int ,
symbol_table . endo_nbr ( ) } ;
writeBytecodeHelper < false > ( code_file ) ;
2010-01-22 11:03:29 +01:00
}
void
2020-05-19 17:43:35 +02:00
StaticModel : : writeStaticBlockBytecode ( const string & basename ) const
2009-12-16 18:13:23 +01:00
{
2023-01-09 13:35:49 +01:00
BytecodeWriter code_file { basename + " /model/bytecode/block/static.cod " } ;
2009-12-16 14:21:31 +01:00
2023-01-09 13:35:49 +01:00
const filesystem : : path bin_filename { basename + " /model/bytecode/block/static.bin " } ;
2022-07-19 18:24:36 +02:00
ofstream bin_file { bin_filename , ios : : out | ios : : binary } ;
if ( ! bin_file . is_open ( ) )
2009-12-16 18:13:23 +01:00
{
2023-01-05 16:40:04 +01:00
cerr < < R " (Error : Can't open file " ) " << bin_filename.string() << R " ( " for writing) " < < endl ;
2009-12-16 18:13:23 +01:00
exit ( EXIT_FAILURE ) ;
}
2022-07-19 18:24:36 +02:00
// Temporary variables declaration
code_file < < FDIMST_ { static_cast < int > ( blocks_temporary_terms_idxs . size ( ) ) } ;
for ( int block { 0 } ; block < static_cast < int > ( blocks . size ( ) ) ; block + + )
2009-12-16 18:13:23 +01:00
{
2022-07-19 18:24:36 +02:00
const BlockSimulationType simulation_type { blocks [ block ] . simulation_type } ;
const int block_size { blocks [ block ] . size } ;
const int u_count { simulation_type = = BlockSimulationType : : solveBackwardComplete
| | simulation_type = = BlockSimulationType : : solveForwardComplete
? writeBlockBytecodeBinFile ( bin_file , block )
: 0 } ;
code_file < < FBEGINBLOCK_ { blocks [ block ] . mfs_size ,
simulation_type ,
blocks [ block ] . first_equation ,
block_size ,
endo_idx_block2orig ,
eq_idx_block2orig ,
blocks [ block ] . linear ,
symbol_table . endo_nbr ( ) ,
u_count ,
block_size } ;
writeBlockBytecodeHelper < false > ( code_file , block ) ;
2009-12-16 18:13:23 +01:00
}
2022-07-19 18:24:36 +02:00
code_file < < FEND_ { } ;
2009-12-16 18:13:23 +01:00
}
2009-12-16 14:21:31 +01:00
void
2023-03-27 17:20:34 +02:00
StaticModel : : computingPass ( int derivsOrder , int paramsDerivsOrder , const eval_context_t & eval_context , bool no_tmp_terms , bool block , bool use_dll )
2009-12-16 14:21:31 +01:00
{
2011-06-22 11:34:25 +02:00
initializeVariablesAndEquations ( ) ;
2022-01-21 14:31:29 +01:00
/* In both MATLAB and Julia, tensors for higher-order derivatives are stored
in matrices whose columns correspond to variable multi - indices . Since we
currently are limited to 32 - bit signed integers ( hence 31 bits ) for matrix
indices , check that we will not overflow ( see # 89 ) . Note that such a check
is not needed for parameter derivatives , since tensors for those are not
stored as matrices . This check is implemented at this place for symmetry
with DynamicModel : : computingPass ( ) . */
if ( log2 ( symbol_table . endo_nbr ( ) ) * derivsOrder > = numeric_limits < int > : : digits )
{
2022-09-14 17:48:33 +02:00
cerr < < " ERROR: The derivatives matrix of the " < < modelClassName ( ) < < " is too large. Please decrease the approximation order. " < < endl ;
2022-01-21 14:31:29 +01:00
exit ( EXIT_FAILURE ) ;
}
2018-11-22 14:32:40 +01:00
// Compute derivatives w.r. to all endogenous
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
{
2018-07-17 18:34:07 +02:00
int id = symbol_table . getID ( SymbolType : : endogenous , i ) ;
2016-03-21 20:42:35 +01:00
vars . insert ( getDerivID ( id , 0 ) ) ;
2017-06-14 07:01:31 +02:00
}
2009-12-16 14:21:31 +01:00
// Launch computations
2022-09-14 17:48:33 +02:00
cout < < " Computing " < < modelClassName ( ) < < " derivatives (order " < < derivsOrder < < " ). " < < endl ;
2009-12-16 14:21:31 +01:00
2018-11-22 14:32:40 +01:00
computeDerivatives ( derivsOrder , vars ) ;
2013-11-22 21:09:04 +01:00
2016-05-18 12:26:19 +02:00
if ( paramsDerivsOrder > 0 )
2012-11-29 18:07:48 +01:00
{
2022-09-14 17:48:33 +02:00
cout < < " Computing " < < modelClassName ( ) < < " derivatives w.r.t. parameters (order " < < paramsDerivsOrder < < " ). " < < endl ;
2016-05-18 12:26:19 +02:00
computeParamsDerivatives ( paramsDerivsOrder ) ;
2012-11-29 18:07:48 +01:00
}
2023-03-27 17:20:34 +02:00
computeTemporaryTerms ( ! use_dll , no_tmp_terms ) ;
2009-12-16 14:21:31 +01:00
2022-09-13 15:43:13 +02:00
if ( paramsDerivsOrder > 0 & & ! no_tmp_terms )
computeParamsDerivativesTemporaryTerms ( ) ;
2009-12-16 14:21:31 +01:00
2022-09-28 16:31:51 +02:00
computingPassBlock ( eval_context , no_tmp_terms ) ;
if ( ! block_decomposed & & block )
2010-01-22 11:03:29 +01:00
{
2022-09-13 15:43:13 +02:00
cerr < < " ERROR: Block decomposition requested but failed. If your model does not have a steady state, you may want to try the 'no_static' option of the 'model' block. " < < endl ;
exit ( EXIT_FAILURE ) ;
2010-01-22 11:03:29 +01:00
}
2009-04-14 16:39:53 +02:00
}
void
2018-06-27 15:01:31 +02:00
StaticModel : : writeStaticMFile ( const string & basename ) const
2009-04-14 16:39:53 +02:00
{
2022-07-11 17:33:09 +02:00
auto [ d_output , tt_output ] = writeModelFileHelper < ExprNodeOutputType : : matlabStaticModel > ( ) ;
ostringstream init_output , end_output ;
init_output < < " residual = zeros( " < < equations . size ( ) < < " , 1); " ;
end_output < < " if ~isreal(residual) " < < endl
< < " residual = real(residual)+imag(residual).^2; " < < endl
< < " end " ;
2022-07-12 17:04:41 +02:00
writeStaticMFileHelper ( basename , " static_resid " , " residual " , " static_resid_tt " ,
2022-09-26 14:53:36 +02:00
temporary_terms_derivatives [ 0 ] . size ( ) ,
2022-07-11 17:33:09 +02:00
" " , init_output , end_output ,
d_output [ 0 ] , tt_output [ 0 ] ) ;
init_output . str ( " " ) ;
end_output . str ( " " ) ;
init_output < < " g1 = zeros( " < < equations . size ( ) < < " , " < < symbol_table . endo_nbr ( ) < < " ); " ;
end_output < < " if ~isreal(g1) " < < endl
< < " g1 = real(g1)+2*imag(g1); " < < endl
< < " end " ;
2022-07-12 17:04:41 +02:00
writeStaticMFileHelper ( basename , " static_g1 " , " g1 " , " static_g1_tt " ,
2022-09-26 14:53:36 +02:00
temporary_terms_derivatives [ 0 ] . size ( ) + temporary_terms_derivatives [ 1 ] . size ( ) ,
2022-07-11 17:33:09 +02:00
" static_resid_tt " ,
init_output , end_output ,
d_output [ 1 ] , tt_output [ 1 ] ) ;
2022-07-12 17:04:41 +02:00
writeStaticMWrapperFunction ( basename , " g1 " ) ;
2022-07-11 17:33:09 +02:00
// For order ≥ 2
int ncols { symbol_table . endo_nbr ( ) } ;
2022-09-26 14:53:36 +02:00
int ntt { static_cast < int > ( temporary_terms_derivatives [ 0 ] . size ( ) + temporary_terms_derivatives [ 1 ] . size ( ) ) } ;
2022-07-11 17:33:09 +02:00
for ( size_t i { 2 } ; i < derivatives . size ( ) ; i + + )
{
ncols * = symbol_table . endo_nbr ( ) ;
ntt + = temporary_terms_derivatives [ i ] . size ( ) ;
string gname { " g " + to_string ( i ) } ;
string gprevname { " g " + to_string ( i - 1 ) } ;
init_output . str ( " " ) ;
end_output . str ( " " ) ;
if ( derivatives [ i ] . size ( ) )
{
init_output < < gname < < " _i = zeros( " < < NNZDerivatives [ i ] < < " ,1); " < < endl
< < gname < < " _j = zeros( " < < NNZDerivatives [ i ] < < " ,1); " < < endl
< < gname < < " _v = zeros( " < < NNZDerivatives [ i ] < < " ,1); " < < endl ;
end_output < < gname < < " = sparse( "
< < gname < < " _i, " < < gname < < " _j, " < < gname < < " _v, "
< < equations . size ( ) < < " , " < < ncols < < " ); " ;
}
else
init_output < < gname < < " = sparse([],[],[], " < < equations . size ( ) < < " , " < < ncols < < " ); " ;
2022-07-12 17:04:41 +02:00
writeStaticMFileHelper ( basename , " static_ " + gname , gname ,
2022-07-11 17:33:09 +02:00
" static_ " + gname + " _tt " ,
ntt ,
" static_ " + gprevname + " _tt " ,
init_output , end_output ,
d_output [ i ] , tt_output [ i ] ) ;
if ( i < = 3 )
2022-07-12 17:04:41 +02:00
writeStaticMWrapperFunction ( basename , gname ) ;
2022-07-11 17:33:09 +02:00
}
2022-07-12 17:04:41 +02:00
writeStaticMCompatFile ( basename ) ;
2018-03-27 17:14:30 +02:00
}
2009-12-16 14:21:31 +01:00
2018-03-27 17:14:30 +02:00
void
2022-07-12 17:04:41 +02:00
StaticModel : : writeStaticMWrapperFunction ( const string & basename , const string & ending ) const
2018-03-27 17:14:30 +02:00
{
string name ;
if ( ending = = " g1 " )
2018-06-27 15:01:31 +02:00
name = " static_resid_g1 " ;
2018-03-27 17:14:30 +02:00
else if ( ending = = " g2 " )
2018-06-27 15:01:31 +02:00
name = " static_resid_g1_g2 " ;
2018-03-27 17:14:30 +02:00
else if ( ending = = " g3 " )
2018-06-27 15:01:31 +02:00
name = " static_resid_g1_g2_g3 " ;
2018-03-27 17:14:30 +02:00
2022-10-11 15:59:56 +02:00
filesystem : : path filename { packageDir ( basename ) / ( name + " .m " ) } ;
2022-07-11 16:09:07 +02:00
ofstream output { filename , ios : : out | ios : : binary } ;
2009-12-16 14:21:31 +01:00
if ( ! output . is_open ( ) )
{
2022-10-11 15:59:56 +02:00
cerr < < " ERROR: Can't open file " < < filename . string ( ) < < " for writing " < < endl ;
2018-03-27 17:14:30 +02:00
exit ( EXIT_FAILURE ) ;
}
if ( ending = = " g1 " )
2018-05-23 17:32:50 +02:00
output < < " function [residual, g1] = " < < name < < " (T, y, x, params, T_flag) " < < endl
< < " % function [residual, g1] = " < < name < < " (T, y, x, params, T_flag) " < < endl ;
2018-03-27 17:14:30 +02:00
else if ( ending = = " g2 " )
2018-05-23 17:32:50 +02:00
output < < " function [residual, g1, g2] = " < < name < < " (T, y, x, params, T_flag) " < < endl
< < " % function [residual, g1, g2] = " < < name < < " (T, y, x, params, T_flag) " < < endl ;
2018-03-27 17:14:30 +02:00
else if ( ending = = " g3 " )
2018-05-23 17:32:50 +02:00
output < < " function [residual, g1, g2, g3] = " < < name < < " (T, y, x, params, T_flag) " < < endl
< < " % function [residual, g1, g2, g3] = " < < name < < " (T, y, x, params, T_flag) " < < endl ;
2018-03-27 17:14:30 +02:00
output < < " % " < < endl
< < " % Wrapper function automatically created by Dynare " < < endl
2018-05-23 17:32:50 +02:00
< < " % " < < endl
< < endl
< < " if T_flag " < < endl
2018-06-27 15:01:31 +02:00
< < " T = " < < basename < < " .static_ " < < ending < < " _tt(T, y, x, params); " < < endl
2018-05-23 17:32:50 +02:00
< < " end " < < endl ;
2018-03-27 17:14:30 +02:00
if ( ending = = " g1 " )
2018-06-27 15:01:31 +02:00
output < < " residual = " < < basename < < " .static_resid(T, y, x, params, false); " < < endl
< < " g1 = " < < basename < < " .static_g1(T, y, x, params, false); " < < endl ;
2018-03-27 17:14:30 +02:00
else if ( ending = = " g2 " )
2018-06-27 15:01:31 +02:00
output < < " [residual, g1] = " < < basename < < " .static_resid_g1(T, y, x, params, false); " < < endl
< < " g2 = " < < basename < < " .static_g2(T, y, x, params, false); " < < endl ;
2018-03-27 17:14:30 +02:00
else if ( ending = = " g3 " )
2018-06-27 15:01:31 +02:00
output < < " [residual, g1, g2] = " < < basename < < " .static_resid_g1_g2(T, y, x, params, false); " < < endl
< < " g3 = " < < basename < < " .static_g3(T, y, x, params, false); " < < endl ;
2018-03-27 17:14:30 +02:00
output < < endl < < " end " < < endl ;
output . close ( ) ;
}
void
2022-07-12 17:04:41 +02:00
StaticModel : : writeStaticMFileHelper ( const string & basename ,
2018-06-27 15:01:31 +02:00
const string & name , const string & retvalname ,
2018-03-27 17:14:30 +02:00
const string & name_tt , size_t ttlen ,
const string & previous_tt_name ,
const ostringstream & init_s , const ostringstream & end_s ,
const ostringstream & s , const ostringstream & s_tt ) const
{
2022-10-11 15:59:56 +02:00
filesystem : : path filename { packageDir ( basename ) / ( name_tt + " .m " ) } ;
2022-07-11 16:09:07 +02:00
ofstream output { filename , ios : : out | ios : : binary } ;
2018-03-27 17:14:30 +02:00
if ( ! output . is_open ( ) )
{
2022-10-11 15:59:56 +02:00
cerr < < " ERROR: Can't open file " < < filename . string ( ) < < " for writing " < < endl ;
2009-12-16 14:21:31 +01:00
exit ( EXIT_FAILURE ) ;
}
2018-03-27 17:14:30 +02:00
output < < " function T = " < < name_tt < < " (T, y, x, params) " < < endl
< < " % function T = " < < name_tt < < " (T, y, x, params) " < < endl
2009-07-07 16:20:48 +02:00
< < " % " < < endl
2018-03-27 17:14:30 +02:00
< < " % File created by Dynare Preprocessor from .mod file " < < endl
2009-07-07 16:20:48 +02:00
< < " % " < < endl
2018-03-27 17:14:30 +02:00
< < " % Inputs: " < < endl
< < " % T [#temp variables by 1] double vector of temporary terms to be filled by function " < < 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
2013-07-28 11:32:14 +02:00
< < " % " < < endl
2018-03-27 17:14:30 +02:00
< < " % Output: " < < endl
< < " % T [#temp variables by 1] double vector of temporary terms " < < endl
< < " % " < < endl < < endl
< < " assert(length(T) >= " < < ttlen < < " ); " < < endl
< < endl ;
if ( ! previous_tt_name . empty ( ) )
2018-06-27 15:01:31 +02:00
output < < " T = " < < basename < < " . " < < previous_tt_name < < " (T, y, x, params); " < < endl < < endl ;
2018-03-27 17:14:30 +02:00
output < < s_tt . str ( ) < < endl
< < " end " < < endl ;
output . close ( ) ;
2022-10-11 15:59:56 +02:00
filename = packageDir ( basename ) / ( name + " .m " ) ;
2018-06-27 15:12:12 +02:00
output . open ( filename , ios : : out | ios : : binary ) ;
2018-03-27 17:14:30 +02:00
if ( ! output . is_open ( ) )
{
2022-10-11 15:59:56 +02:00
cerr < < " ERROR: Can't open file " < < filename . string ( ) < < " for writing " < < endl ;
2018-03-27 17:14:30 +02:00
exit ( EXIT_FAILURE ) ;
}
output < < " function " < < retvalname < < " = " < < name < < " (T, y, x, params, T_flag) " < < endl
< < " % function " < < retvalname < < " = " < < name < < " (T, y, x, params, T_flag) " < < endl
2013-07-28 11:32:14 +02:00
< < " % " < < endl
2018-03-27 17:14:30 +02:00
< < " % File created by Dynare Preprocessor from .mod file " < < endl
2017-06-14 07:01:31 +02:00
< < " % " < < endl
2018-03-27 17:14:30 +02:00
< < " % Inputs: " < < endl
< < " % T [#temp variables by 1] double vector of temporary terms to be filled by function " < < 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
< < " % to evaluate the model " < < endl
< < " % T_flag boolean boolean flag saying whether or not to calculate temporary terms " < < endl
< < " % " < < endl
< < " % Output: " < < endl
< < " % " < < retvalname < < endl
< < " % " < < endl < < endl ;
if ( ! name_tt . empty ( ) )
output < < " if T_flag " < < endl
2018-06-27 15:01:31 +02:00
< < " T = " < < basename < < " . " < < name_tt < < " (T, y, x, params); " < < endl
2018-03-27 17:14:30 +02:00
< < " end " < < endl ;
output < < init_s . str ( ) < < endl
< < s . str ( )
< < end_s . str ( ) < < endl
< < " end " < < endl ;
2011-08-18 12:44:11 +02:00
output . close ( ) ;
}
2018-05-23 16:10:26 +02:00
void
2022-07-12 17:04:41 +02:00
StaticModel : : writeStaticMCompatFile ( const string & basename ) const
2018-05-23 16:10:26 +02:00
{
2022-10-11 15:59:56 +02:00
filesystem : : path filename { packageDir ( basename ) / " static.m " } ;
2022-07-11 16:09:07 +02:00
ofstream output { filename , ios : : out | ios : : binary } ;
2018-05-23 16:10:26 +02:00
if ( ! output . is_open ( ) )
{
2022-10-11 15:59:56 +02:00
cerr < < " ERROR: Can't open file " < < filename . string ( ) < < " for writing " < < endl ;
2018-05-23 16:10:26 +02:00
exit ( EXIT_FAILURE ) ;
}
2022-09-26 14:53:36 +02:00
int ntt { static_cast < int > ( temporary_terms_derivatives [ 0 ] . size ( ) + temporary_terms_derivatives [ 1 ] . size ( ) + temporary_terms_derivatives [ 2 ] . size ( ) + temporary_terms_derivatives [ 3 ] . size ( ) ) } ;
2018-05-23 16:10:26 +02:00
2018-06-27 15:01:31 +02:00
output < < " function [residual, g1, g2, g3] = static(y, x, params) " < < endl
2018-05-23 16:10:26 +02:00
< < " T = NaN( " < < ntt < < " , 1); " < < endl
< < " if nargout <= 1 " < < endl
2018-06-27 15:01:31 +02:00
< < " residual = " < < basename < < " .static_resid(T, y, x, params, true); " < < endl
2018-05-23 16:10:26 +02:00
< < " elseif nargout == 2 " < < endl
2018-06-27 15:01:31 +02:00
< < " [residual, g1] = " < < basename < < " .static_resid_g1(T, y, x, params, true); " < < endl
2018-05-23 16:10:26 +02:00
< < " elseif nargout == 3 " < < endl
2018-06-27 15:01:31 +02:00
< < " [residual, g1, g2] = " < < basename < < " .static_resid_g1_g2(T, y, x, params, true); " < < endl
2018-05-23 16:10:26 +02:00
< < " else " < < endl
2018-06-27 15:01:31 +02:00
< < " [residual, g1, g2, g3] = " < < basename < < " .static_resid_g1_g2_g3(T, y, x, params, true); " < < endl
2018-05-23 16:10:26 +02:00
< < " end " < < endl
< < " end " < < endl ;
output . close ( ) ;
}
2015-07-27 15:33:38 +02:00
void
2023-03-20 17:54:32 +01:00
StaticModel : : writeStaticFile ( const string & basename , bool use_dll , const string & mexext , const filesystem : : path & matlabroot , bool julia ) const
2009-12-16 18:13:23 +01:00
{
2020-06-23 15:13:04 +02:00
filesystem : : path model_dir { basename } ;
model_dir / = " model " ;
if ( use_dll )
2022-10-05 17:45:18 +02:00
{
create_directories ( model_dir / " src " / " sparse " ) ;
if ( block_decomposed )
create_directories ( model_dir / " src " / " sparse " / " block " ) ;
}
2022-09-20 18:28:30 +02:00
if ( julia )
create_directories ( model_dir / " julia " ) ;
else
2022-10-11 15:59:56 +02:00
{
auto plusfolder { packageDir ( basename ) } ;
/* The following is not a duplicate of the same call from
ModFile : : writeMOutput ( ) , because of planner_objective which needs its
+ objective subdirectory */
create_directories ( plusfolder ) ;
2022-09-30 17:26:43 +02:00
auto sparsefolder { plusfolder / " +sparse " } ;
create_directories ( sparsefolder ) ;
if ( block_decomposed )
create_directories ( sparsefolder / " +block " ) ;
2022-11-09 14:46:02 +01:00
create_directories ( plusfolder / " +debug " ) ;
2022-10-11 15:59:56 +02:00
}
2023-01-09 13:35:49 +01:00
create_directories ( model_dir / " bytecode " / " block " ) ;
2020-06-23 15:13:04 +02:00
2022-09-20 18:28:30 +02:00
// Legacy representation
2023-01-13 12:05:04 +01:00
if ( use_dll )
2023-03-20 17:54:32 +01:00
writeModelCFile < false > ( basename , mexext , matlabroot ) ;
2023-01-13 12:05:04 +01:00
else if ( ! julia ) // M-files
writeStaticMFile ( basename ) ;
// The legacy representation is no longer produced for Julia
2023-01-09 13:35:49 +01:00
writeStaticBytecode ( basename ) ;
if ( block_decomposed )
writeStaticBlockBytecode ( basename ) ;
2020-06-23 15:13:04 +02:00
2022-09-20 18:28:30 +02:00
// Sparse representation
2022-09-30 17:26:43 +02:00
if ( use_dll )
2023-03-20 17:54:32 +01:00
writeSparseModelCFiles < false > ( basename , mexext , matlabroot ) ;
2022-09-30 17:26:43 +02:00
else if ( julia )
2022-09-20 18:28:30 +02:00
writeSparseModelJuliaFiles < false > ( basename ) ;
2022-09-30 17:26:43 +02:00
else // MATLAB/Octave
writeSparseModelMFiles < false > ( basename ) ;
2022-09-20 18:28:30 +02:00
2023-03-28 16:37:05 +02:00
writeSetAuxiliaryVariablesFile < false > ( basename , julia ) ;
2022-11-09 14:46:02 +01:00
// Support for model debugging
if ( ! julia )
writeDebugModelMFiles < false > ( basename ) ;
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
{
2018-06-04 12:26:16 +02:00
for ( auto equation : equations )
2021-07-20 12:10:58 +02:00
if ( equation - > hasExogenous ( ) )
return true ;
2016-08-12 11:50:57 +02:00
return false ;
}
2009-12-16 14:21:31 +01:00
void
2022-09-28 16:54:03 +02:00
StaticModel : : writeDriverOutput ( ostream & output ) const
2009-12-16 14:21:31 +01:00
{
2019-04-12 15:41:52 +02:00
output < < " M_.static_tmp_nbr = [ " ;
2019-12-20 16:59:30 +01:00
for ( const auto & temporary_terms_derivative : temporary_terms_derivatives )
2019-12-16 19:42:59 +01:00
output < < temporary_terms_derivative . size ( ) < < " ; " ;
2019-04-12 15:41:52 +02:00
output < < " ]; " < < endl ;
2018-03-28 18:53:02 +02:00
2022-09-28 16:54:03 +02:00
if ( block_decomposed )
2022-09-21 18:47:14 +02:00
writeBlockDriverOutput ( output ) ;
2022-09-14 17:07:08 +02:00
writeDriverSparseIndicesHelper < false > ( output ) ;
2022-09-21 18:47:14 +02:00
}
2009-06-30 17:07:09 +02:00
2022-09-21 18:47:14 +02:00
void
StaticModel : : writeBlockDriverOutput ( ostream & output ) const
{
2020-05-20 11:35:14 +02:00
for ( int blk = 0 ; blk < static_cast < int > ( blocks . size ( ) ) ; blk + + )
2012-06-06 16:36:56 +02:00
{
2022-12-12 14:08:03 +01:00
output < < " M_.block_structure_stat.block( " < < blk + 1 < < " ).Simulation_Type = " < < static_cast < int > ( blocks [ blk ] . simulation_type ) < < " ; " < < endl
< < " M_.block_structure_stat.block( " < < blk + 1 < < " ).endo_nbr = " < < blocks [ blk ] . size < < " ; " < < endl
< < " M_.block_structure_stat.block( " < < blk + 1 < < " ).mfs = " < < blocks [ blk ] . mfs_size < < " ; " < < endl
< < " M_.block_structure_stat.block( " < < blk + 1 < < " ).equation = [ " ;
2020-05-20 11:35:14 +02:00
for ( int eq = 0 ; eq < blocks [ blk ] . size ; eq + + )
output < < " " < < getBlockEquationID ( blk , eq ) + 1 ;
2020-05-05 17:49:58 +02:00
output < < " ]; " < < endl
2022-12-12 14:08:03 +01:00
< < " M_.block_structure_stat.block( " < < blk + 1 < < " ).variable = [ " ;
2020-05-20 11:35:14 +02:00
for ( int var = 0 ; var < blocks [ blk ] . size ; var + + )
output < < " " < < getBlockVariableID ( blk , var ) + 1 ;
2020-05-05 17:49:58 +02:00
output < < " ]; " < < endl ;
2012-06-06 16:36:56 +02:00
}
2022-12-12 14:08:03 +01:00
output < < " M_.block_structure_stat.variable_reordered = [ " ;
2020-05-05 17:49:58 +02:00
for ( int i = 0 ; i < symbol_table . endo_nbr ( ) ; i + + )
2020-04-17 14:55:55 +02:00
output < < " " < < endo_idx_block2orig [ i ] + 1 ;
2018-11-23 17:19:59 +01:00
output < < " ]; " < < endl
< < " M_.block_structure_stat.equation_reordered = [ " ;
2020-05-05 17:49:58 +02:00
for ( int i = 0 ; i < symbol_table . endo_nbr ( ) ; i + + )
2020-04-17 14:55:55 +02:00
output < < " " < < eq_idx_block2orig [ i ] + 1 ;
2018-11-23 17:19:59 +01:00
output < < " ]; " < < endl ;
2017-06-14 07:01:31 +02:00
2020-05-05 17:49:58 +02:00
set < pair < int , int > > row_incidence ;
for ( const auto & [ indices , d1 ] : derivatives [ 1 ] )
if ( int deriv_id = indices [ 1 ] ;
getTypeByDerivID ( deriv_id ) = = SymbolType : : endogenous )
{
int eq = indices [ 0 ] ;
2022-07-12 15:28:52 +02:00
int var { getTypeSpecificIDByDerivID ( deriv_id ) } ;
2020-05-05 17:49:58 +02:00
row_incidence . emplace ( eq , var ) ;
}
2020-05-20 11:35:14 +02:00
output < < " M_.block_structure_stat.incidence.sparse_IM = [ " < < endl ;
2020-05-05 17:49:58 +02:00
for ( auto [ eq , var ] : row_incidence )
2020-05-20 11:35:14 +02:00
output < < " " < < eq + 1 < < " " < < var + 1 < < " ; " < < endl ;
2020-05-25 18:35:36 +02:00
output < < " ]; " < < endl
< < " M_.block_structure_stat.tmp_nbr = " < < blocks_temporary_terms_idxs . size ( )
< < " ; " < < endl ;
2022-09-28 19:18:15 +02:00
writeBlockDriverSparseIndicesHelper < false > ( output ) ;
2009-04-28 19:11:48 +02:00
}
2009-12-16 14:21:31 +01:00
SymbolType
2018-06-04 12:50:53 +02:00
StaticModel : : getTypeByDerivID ( int deriv_id ) const noexcept ( false )
2009-12-16 14:21:31 +01:00
{
2012-11-29 18:07:48 +01:00
if ( deriv_id < symbol_table . endo_nbr ( ) )
2018-07-17 18:34:07 +02:00
return SymbolType : : endogenous ;
2012-11-29 18:07:48 +01:00
else if ( deriv_id < symbol_table . endo_nbr ( ) + symbol_table . param_nbr ( ) )
2018-07-17 18:34:07 +02:00
return SymbolType : : parameter ;
2012-11-29 18:07:48 +01:00
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
2022-06-24 17:10:12 +02:00
StaticModel : : getLagByDerivID ( [[maybe_unused]] int deriv_id ) const noexcept ( false )
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
2018-06-04 12:50:53 +02:00
StaticModel : : getSymbIDByDerivID ( int deriv_id ) const noexcept ( false )
2009-06-30 17:07:09 +02:00
{
2012-11-29 18:07:48 +01:00
if ( deriv_id < symbol_table . endo_nbr ( ) )
2018-07-17 18:34:07 +02:00
return symbol_table . getID ( SymbolType : : endogenous , deriv_id ) ;
2012-11-29 18:07:48 +01:00
else if ( deriv_id < symbol_table . endo_nbr ( ) + symbol_table . param_nbr ( ) )
2018-07-17 18:34:07 +02:00
return symbol_table . getID ( SymbolType : : parameter , deriv_id - symbol_table . endo_nbr ( ) ) ;
2012-11-29 18:07:48 +01:00
else
throw UnknownDerivIDException ( ) ;
2009-06-30 17:07:09 +02:00
}
2022-07-12 15:28:52 +02:00
int
StaticModel : : getTypeSpecificIDByDerivID ( int deriv_id ) const
{
if ( deriv_id < symbol_table . endo_nbr ( ) )
return deriv_id ;
else if ( deriv_id < symbol_table . endo_nbr ( ) + symbol_table . param_nbr ( ) )
return deriv_id - symbol_table . endo_nbr ( ) ;
else
throw UnknownDerivIDException ( ) ;
}
2009-12-16 14:21:31 +01:00
int
2022-06-24 17:10:12 +02:00
StaticModel : : getDerivID ( int symb_id , [[maybe_unused]] int lag ) const noexcept ( false )
2009-06-30 17:07:09 +02:00
{
2018-07-17 18:34:07 +02:00
if ( symbol_table . getType ( symb_id ) = = SymbolType : : endogenous )
2012-11-29 18:07:48 +01:00
return symbol_table . getTypeSpecificID ( symb_id ) ;
2018-07-17 18:34:07 +02:00
else if ( symbol_table . getType ( symb_id ) = = SymbolType : : parameter )
2012-11-29 18:07:48 +01:00
return symbol_table . getTypeSpecificID ( symb_id ) + symbol_table . endo_nbr ( ) ;
2009-12-16 14:21:31 +01:00
else
2022-06-08 12:45:33 +02:00
/* See the special treatment in VariableNode::prepareForDerivation(),
VariableNode : : computeDerivative ( ) and VariableNode : : getChainRuleDerivative ( ) */
throw UnknownDerivIDException { } ;
2009-12-16 14:21:31 +01:00
}
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-07-07 16:20:48 +02:00
void
2020-03-19 17:46:10 +01:00
StaticModel : : computeChainRuleJacobian ( )
2009-07-07 16:20:48 +02:00
{
2020-04-23 16:04:02 +02:00
int nb_blocks = blocks . size ( ) ;
2022-09-28 19:18:15 +02:00
2020-03-19 17:46:10 +01:00
blocks_derivatives . resize ( nb_blocks ) ;
2022-09-28 19:18:15 +02:00
blocks_jacobian_sparse_column_major_order . resize ( nb_blocks ) ;
blocks_jacobian_sparse_colptr . resize ( nb_blocks ) ;
2020-04-24 12:29:02 +02:00
for ( int blk = 0 ; blk < nb_blocks ; blk + + )
2009-07-07 16:20:48 +02:00
{
2020-04-24 12:29:02 +02:00
int nb_recursives = blocks [ blk ] . getRecursiveSize ( ) ;
2022-09-28 19:18:15 +02:00
BlockSimulationType simulation_type { blocks [ blk ] . simulation_type } ;
2020-04-24 12:29:02 +02:00
2020-10-02 18:31:55 +02:00
map < int , BinaryOpNode * > recursive_vars ;
2020-04-24 12:29:02 +02:00
for ( int i = 0 ; i < nb_recursives ; i + + )
2009-07-10 16:22:40 +02:00
{
2020-04-24 12:29:02 +02:00
int deriv_id = getDerivID ( symbol_table . getID ( SymbolType : : endogenous , getBlockVariableID ( blk , i ) ) , 0 ) ;
2020-05-20 11:49:32 +02:00
if ( getBlockEquationType ( blk , i ) = = EquationType : : evaluateRenormalized )
2020-04-24 12:29:02 +02:00
recursive_vars [ deriv_id ] = getBlockEquationRenormalizedExpr ( blk , i ) ;
else
recursive_vars [ deriv_id ] = getBlockEquationExpr ( blk , i ) ;
2009-07-10 16:22:40 +02:00
}
2020-04-24 12:29:02 +02:00
2022-09-28 19:18:15 +02:00
assert ( simulation_type ! = BlockSimulationType : : solveTwoBoundariesSimple
& & simulation_type ! = BlockSimulationType : : solveTwoBoundariesComplete ) ;
2020-04-24 12:29:02 +02:00
int size = blocks [ blk ] . size ;
2023-03-02 17:49:16 +01:00
map < expr_t , set < int > > non_null_chain_rule_derivatives ;
2022-11-08 12:28:48 +01:00
map < pair < expr_t , int > , expr_t > chain_rule_deriv_cache ;
2020-04-24 12:29:02 +02:00
for ( int eq = nb_recursives ; eq < size ; eq + + )
2009-07-07 16:20:48 +02:00
{
2020-04-24 12:29:02 +02:00
int eq_orig = getBlockEquationID ( blk , eq ) ;
for ( int var = nb_recursives ; var < size ; var + + )
2009-12-16 14:21:31 +01:00
{
2020-04-24 12:29:02 +02:00
int var_orig = getBlockVariableID ( blk , var ) ;
2023-03-02 17:49:16 +01:00
expr_t d1 = equations [ eq_orig ] - > getChainRuleDerivative ( getDerivID ( symbol_table . getID ( SymbolType : : endogenous , var_orig ) , 0 ) , recursive_vars , non_null_chain_rule_derivatives , chain_rule_deriv_cache ) ;
2020-05-06 14:21:33 +02:00
if ( d1 ! = Zero )
blocks_derivatives [ blk ] [ { eq , var , 0 } ] = d1 ;
2009-07-07 17:58:17 +02:00
}
2009-07-07 16:20:48 +02:00
}
2022-09-28 19:18:15 +02:00
// Compute the sparse representation of the Jacobian
if ( simulation_type ! = BlockSimulationType : : evaluateForward
& & simulation_type ! = BlockSimulationType : : evaluateBackward )
{
for ( const auto & [ indices , d1 ] : blocks_derivatives [ blk ] )
{
auto & [ eq , var , lag ] { indices } ;
assert ( lag = = 0 ) ;
2023-01-17 19:04:07 +01:00
if ( eq > = nb_recursives & & var > = nb_recursives )
2023-03-23 12:58:42 +01:00
blocks_jacobian_sparse_column_major_order [ blk ] . try_emplace ( { eq - nb_recursives , var - nb_recursives } , d1 ) ;
2022-09-28 19:18:15 +02:00
}
2023-01-19 10:03:37 +01:00
blocks_jacobian_sparse_colptr [ blk ] = computeCSCColPtr ( blocks_jacobian_sparse_column_major_order [ blk ] , blocks [ blk ] . mfs_size ) ;
2022-09-28 19:18:15 +02:00
}
2009-07-07 16:20:48 +02:00
}
}
2009-12-16 14:21:31 +01:00
void
2017-08-24 15:35:10 +02:00
StaticModel : : writeLatexFile ( const string & basename , bool write_equation_tags ) const
2009-12-16 18:13:23 +01:00
{
2019-07-11 17:33:53 +02:00
writeLatexModelFile ( basename , " static " , ExprNodeOutputType : : latexStaticModel , write_equation_tags ) ;
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
{
2018-06-04 12:26:16 +02:00
for ( auto aux_equation : aux_equations )
2009-09-30 17:10:31 +02:00
{
2018-06-04 12:26:16 +02:00
dynamic_cast < ExprNode * > ( aux_equation ) - > writeOutput ( output , output_type ) ;
2009-09-30 17:10:31 +02:00
output < < " ; " < < endl ;
}
}
2011-07-24 20:52:03 +02:00
2017-08-30 11:32:01 +02:00
void
StaticModel : : writeLatexAuxVarRecursiveDefinitions ( ostream & output ) const
{
deriv_node_temp_terms_t tef_terms ;
temporary_terms_t temporary_terms ;
2019-12-20 16:59:30 +01:00
temporary_terms_idxs_t temporary_terms_idxs ;
2018-06-04 12:26:16 +02:00
for ( auto aux_equation : aux_equations )
if ( dynamic_cast < ExprNode * > ( aux_equation ) - > containsExternalFunction ( ) )
2018-09-05 18:27:13 +02:00
dynamic_cast < ExprNode * > ( aux_equation ) - > writeExternalFunctionOutput ( output , ExprNodeOutputType : : latexStaticModel ,
2019-12-20 16:59:30 +01:00
temporary_terms , temporary_terms_idxs , tef_terms ) ;
2018-06-04 12:26:16 +02:00
for ( auto aux_equation : aux_equations )
2017-08-30 11:32:01 +02:00
{
2019-04-03 16:32:52 +02:00
output < < R " ( \b egin{dmath}) " < < endl ;
2023-03-23 18:28:13 +01:00
dynamic_cast < ExprNode * > ( aux_equation ) - > writeOutput ( output , ExprNodeOutputType : : latexStaticModel ) ;
2019-04-03 16:32:52 +02:00
output < < endl < < R " ( \ end{dmath}) " < < endl ;
2017-08-30 11:32:01 +02:00
}
}
2017-10-16 17:24:55 +02:00
void
StaticModel : : writeJsonAuxVarRecursiveDefinitions ( ostream & output ) const
{
deriv_node_temp_terms_t tef_terms ;
temporary_terms_t temporary_terms ;
2018-06-04 12:26:16 +02:00
for ( auto aux_equation : aux_equations )
if ( dynamic_cast < ExprNode * > ( aux_equation ) - > containsExternalFunction ( ) )
2017-10-16 17:24:55 +02:00
{
vector < string > efout ;
2018-06-04 12:26:16 +02:00
dynamic_cast < ExprNode * > ( aux_equation ) - > writeJsonExternalFunctionOutput ( efout ,
2019-12-16 19:42:59 +01:00
temporary_terms ,
tef_terms ,
false ) ;
2022-06-03 16:24:26 +02:00
for ( bool printed_something { false } ;
const auto & it : efout )
2017-10-16 17:24:55 +02:00
{
2022-06-03 16:24:26 +02:00
if ( exchange ( printed_something , true ) )
2017-10-16 17:24:55 +02:00
output < < " , " ;
2022-06-03 16:24:26 +02:00
output < < it ;
2017-10-16 17:24:55 +02:00
}
}
2018-06-04 12:26:16 +02:00
for ( auto aux_equation : aux_equations )
2017-10-16 17:24:55 +02:00
{
2019-04-03 16:32:52 +02:00
output < < R " (, { " lhs " : " ) " ;
2018-11-28 14:32:26 +01:00
aux_equation - > arg1 - > writeJsonOutput ( output , temporary_terms , tef_terms , false ) ;
2019-04-03 16:32:52 +02:00
output < < R " ( " , " rhs " : " ) " ;
2023-03-23 18:28:13 +01:00
aux_equation - > arg2 - > writeJsonOutput ( output , temporary_terms , tef_terms , false ) ;
2019-04-03 16:32:52 +02:00
output < < R " ( " } ) " ;
2017-10-16 17:24:55 +02:00
}
}
2017-02-27 15:40:34 +01:00
void
StaticModel : : writeJsonOutput ( ostream & output ) const
{
2022-09-26 13:07:18 +02:00
output < < R " ( " static_tmp_nbr " : [) " ;
for ( bool printed_something { false } ;
const auto & tts : temporary_terms_derivatives )
{
if ( exchange ( printed_something , true ) )
output < < " , " ;
output < < tts . size ( ) ;
}
output < < " ], " ;
2022-09-14 17:07:08 +02:00
writeJsonSparseIndicesHelper < false > ( output ) ;
2017-02-27 15:40:34 +01:00
}
2017-02-20 12:18:11 +01:00
void
2017-03-02 18:34:18 +01:00
StaticModel : : writeJsonComputingPassOutput ( ostream & output , bool writeDetails ) const
2017-02-20 12:18:11 +01:00
{
2022-07-12 17:47:02 +02:00
auto [ mlv_output , d_output ] { writeJsonComputingPassOutputHelper < false > ( writeDetails ) } ;
2017-02-20 12:18:11 +01:00
2017-03-02 18:34:18 +01:00
if ( writeDetails )
2019-04-03 16:32:52 +02:00
output < < R " ( " static_model " : {) " ;
2017-03-02 18:34:18 +01:00
else
2019-04-03 16:32:52 +02:00
output < < R " ( " static_model_simple " : {) " ;
2022-07-12 17:47:02 +02:00
output < < mlv_output . str ( ) ;
2019-04-18 17:07:55 +02:00
for ( const auto & it : d_output )
output < < " , " < < it . str ( ) ;
output < < " } " ;
2017-02-20 12:18:11 +01:00
}
void
2022-07-12 17:47:02 +02:00
StaticModel : : writeJsonParamsDerivatives ( ostream & output , bool writeDetails ) const
2017-02-20 12:18:11 +01:00
{
2018-11-15 16:39:53 +01:00
if ( ! params_derivatives . size ( ) )
2017-02-20 12:18:11 +01:00
return ;
2022-07-12 17:47:02 +02:00
auto [ mlv_output , tt_output , rp_output , gp_output , rpp_output , gpp_output , hp_output , g3p_output ]
{ writeJsonParamsDerivativesHelper < false > ( writeDetails ) } ;
// g3p_output is ignored
2017-02-20 12:18:11 +01:00
2017-03-02 18:34:18 +01:00
if ( writeDetails )
2019-04-03 16:32:52 +02:00
output < < R " ( " static_model_params_derivative " : {) " ;
2017-03-02 18:34:18 +01:00
else
2019-04-03 16:32:52 +02:00
output < < R " ( " static_model_params_derivatives_simple " : {) " ;
2022-07-12 17:47:02 +02:00
output < < mlv_output . str ( )
< < " , " < < tt_output . str ( )
< < " , " < < rp_output . str ( )
< < " , " < < gp_output . str ( )
< < " , " < < rpp_output . str ( )
< < " , " < < gpp_output . str ( )
< < " , " < < hp_output . str ( )
2017-02-20 12:18:11 +01:00
< < " } " ;
}
Ramsey: write derivatives of static model w.r.t. Lagrange multipliers in a new file
The computing of the Ramsey steady state relies on the fact that Lagrange
multipliers appear linearly in the system to be solved. Instead of directly
solving for the Lagrange multipliers along with the other variables,
dyn_ramsey_static.m reduces the size of the problem by always computing the
value of the multipliers that minimizes the residuals, given the other
variables (using a minimum norm solution, easy to compute because of the
linearity property). That function thus needs the derivatives of the optimality
FOCs with respect to the multipliers. The problem is that, when multipliers
appear in an auxiliary variable related to a lead/lag, then those derivatives
need to be retrieved by a chain rule derivation, which cannot be easily done
with the regular static file.
This commit implements the creation of a new file,
ramsey_multipliers_static_g1.{m,mex}, that provides exactly the needed
derivatives w.r.t. Lagrange multipliers through chain rule derivation.
Ref. dynare#633, dynare#1119, dynare#1133
2023-03-24 18:58:12 +01:00
void
StaticModel : : computeRamseyMultipliersDerivatives ( int ramsey_orig_endo_nbr , bool is_matlab ,
bool no_tmp_terms )
{
// Compute derivation IDs of Lagrange multipliers
set < int > mult_symb_ids { symbol_table . getLagrangeMultipliers ( ) } ;
vector < int > mult_deriv_ids ;
for ( int symb_id : mult_symb_ids )
mult_deriv_ids . push_back ( getDerivID ( symb_id , 0 ) ) ;
// Compute the list of aux vars for which to apply the chain rule derivation
map < int , BinaryOpNode * > recursive_variables ;
for ( auto aux_eq : aux_equations )
{
auto varexpr { dynamic_cast < VariableNode * > ( aux_eq - > arg1 ) } ;
assert ( varexpr & & symbol_table . isAuxiliaryVariable ( varexpr - > symb_id ) ) ;
/* Determine whether the auxiliary variable has been added after the last
Lagrange multiplier . We use the guarantee given by SymbolTable that
symbol IDs are increasing . */
if ( varexpr - > symb_id > * mult_symb_ids . crbegin ( ) )
recursive_variables . emplace ( getDerivID ( varexpr - > symb_id , 0 ) , aux_eq ) ;
}
// Compute the chain rule derivatives w.r.t. multipliers
map < expr_t , set < int > > non_null_chain_rule_derivatives ;
map < pair < expr_t , int > , expr_t > cache ;
for ( int eq { 0 } ; eq < ramsey_orig_endo_nbr ; eq + + )
for ( int mult { 0 } ; mult < static_cast < int > ( mult_deriv_ids . size ( ) ) ; mult + + )
if ( expr_t d { equations [ eq ] - > getChainRuleDerivative ( mult_deriv_ids [ mult ] , recursive_variables ,
non_null_chain_rule_derivatives , cache ) } ;
d ! = Zero )
ramsey_multipliers_derivatives . try_emplace ( { eq , mult } , d ) ;
// Compute the temporary terms
map < pair < int , int > , temporary_terms_t > temp_terms_map ;
map < expr_t , pair < int , pair < int , int > > > reference_count ;
for ( const auto & [ row_col , d ] : ramsey_multipliers_derivatives )
d - > computeTemporaryTerms ( { 1 , 0 } , temp_terms_map , reference_count , is_matlab ) ;
/* If the user has specified the notmpterms option, clear all temporary
terms , except those that correspond to external functions ( since they are
not optional ) */
if ( no_tmp_terms )
for ( auto & it : temp_terms_map )
erase_if ( it . second ,
[ ] ( expr_t e ) { return ! dynamic_cast < AbstractExternalFunctionNode * > ( e ) ; } ) ;
ramsey_multipliers_derivatives_temporary_terms = move ( temp_terms_map [ { 1 , 0 } ] ) ;
for ( int idx { 0 } ;
auto it : ramsey_multipliers_derivatives_temporary_terms )
ramsey_multipliers_derivatives_temporary_terms_idxs [ it ] = idx + + ;
// Compute the CSC format
ramsey_multipliers_derivatives_sparse_colptr = computeCSCColPtr ( ramsey_multipliers_derivatives ,
mult_deriv_ids . size ( ) ) ;
}
void
StaticModel : : writeDriverRamseyMultipliersDerivativesSparseIndices ( ostream & output ) const
{
output < < " M_.ramsey_multipliers_static_g1_sparse_rowval = int32([ " ;
for ( auto & [ row_col , d ] : ramsey_multipliers_derivatives )
output < < row_col . first + 1 < < ' ' ;
output < < " ]); " < < endl
< < " M_.ramsey_multipliers_static_g1_sparse_colval = int32([ " ;
for ( auto & [ row_col , d ] : ramsey_multipliers_derivatives )
output < < row_col . second + 1 < < ' ' ;
output < < " ]); " < < endl
< < " M_.ramsey_multipliers_static_g1_sparse_colptr = int32([ " ;
for ( int it : ramsey_multipliers_derivatives_sparse_colptr )
output < < it + 1 < < ' ' ;
output < < " ]); " < < endl ;
}
void
StaticModel : : writeRamseyMultipliersDerivativesMFile ( const string & basename , int ramsey_orig_endo_nbr ) const
{
constexpr auto output_type { ExprNodeOutputType : : matlabStaticModel } ;
filesystem : : path filename { packageDir ( basename ) / " ramsey_multipliers_static_g1.m " } ;
ofstream output_file { filename , ios : : out | ios : : binary } ;
if ( ! output_file . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename . string ( ) < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
output_file < < " function g1m = ramsey_multipliers_static_g1(y, x, params, sparse_rowval, sparse_colval, sparse_colptr) " < < endl
< < " g1m_v=NaN( " < < ramsey_multipliers_derivatives . size ( ) < < " ,1); " < < endl ;
writeRamseyMultipliersDerivativesHelper < output_type > ( output_file ) ;
// On MATLAB < R2020a, sparse() does not accept int32 indices
output_file < < " if ~isoctave && matlab_ver_less_than('9.8') " < < endl
< < " sparse_rowval = double(sparse_rowval); " < < endl
< < " sparse_colval = double(sparse_colval); " < < endl
< < " end " < < endl
< < " g1m = sparse(sparse_rowval, sparse_colval, g1m_v, " < < ramsey_orig_endo_nbr < < " , " < < symbol_table . getLagrangeMultipliers ( ) . size ( ) < < " ); " < < endl
< < " end " < < endl ;
output_file . close ( ) ;
}
void
StaticModel : : writeRamseyMultipliersDerivativesCFile ( const string & basename , const string & mexext , const filesystem : : path & matlabroot , int ramsey_orig_endo_nbr ) const
{
constexpr auto output_type { ExprNodeOutputType : : CStaticModel } ;
const filesystem : : path model_src_dir { filesystem : : path { basename } / " model " / " src " } ;
const int xlen { symbol_table . exo_nbr ( ) + symbol_table . exo_det_nbr ( ) } ;
const int nzval { static_cast < int > ( ramsey_multipliers_derivatives . size ( ) ) } ;
const int ncols { static_cast < int > ( symbol_table . getLagrangeMultipliers ( ) . size ( ) ) } ;
const filesystem : : path p { model_src_dir / " ramsey_multipliers_static_g1.c " } ;
ofstream output { p , ios : : out | ios : : binary } ;
if ( ! output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < p . string ( ) < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
output < < " #include <math.h> " < < endl < < endl
< < R " (#include " mex . h " ) " < < endl // Needed for calls to external functions
< < endl ;
writePowerDeriv ( output ) ;
output < < endl
< < " void ramsey_multipliers_static_g1(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T, double *restrict g1m_v) " < < endl
< < " { " < < endl ;
writeRamseyMultipliersDerivativesHelper < output_type > ( output ) ;
output < < " } " < < endl
< < endl
< < " void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) " < < endl
< < " { " < < endl
< < " if (nrhs != 6) " < < endl
< < R " ( mexErrMsgTxt( " Accepts exactly 6 input arguments " );) " < < endl
< < " if (nlhs != 1) " < < endl
< < R " ( mexErrMsgTxt( " Accepts exactly 1 output argument " );) " < < endl
< < " if (!(mxIsDouble(prhs[0]) && !mxIsComplex(prhs[0]) && !mxIsSparse(prhs[0]) && mxGetNumberOfElements(prhs[0]) == " < < symbol_table . endo_nbr ( ) < < " )) " < < endl
< < R " ( mexErrMsgTxt( " y must be a real dense numeric array with ) " << symbol_table.endo_nbr() << R " ( elements " );) " < < endl
< < " const double *restrict y = mxGetPr(prhs[0]); " < < endl
< < " if (!(mxIsDouble(prhs[1]) && !mxIsComplex(prhs[1]) && !mxIsSparse(prhs[1]) && mxGetNumberOfElements(prhs[1]) == " < < xlen < < " )) " < < endl
< < R " ( mexErrMsgTxt( " x must be a real dense numeric array with ) " << xlen << R " ( elements " );) " < < endl
< < " const double *restrict x = mxGetPr(prhs[1]); " < < endl
< < " if (!(mxIsDouble(prhs[2]) && !mxIsComplex(prhs[2]) && !mxIsSparse(prhs[2]) && mxGetNumberOfElements(prhs[2]) == " < < symbol_table . param_nbr ( ) < < " )) " < < endl
< < R " ( mexErrMsgTxt( " params must be a real dense numeric array with ) " << symbol_table.param_nbr() << R " ( elements " );) " < < endl
< < " const double *restrict params = mxGetPr(prhs[2]); " < < endl
< < " if (!(mxIsInt32(prhs[3]) && mxGetNumberOfElements(prhs[3]) == " < < nzval < < " )) " < < endl
< < R " ( mexErrMsgTxt( " sparse_rowval must be an int32 array with ) " << nzval << R " ( elements " );) " < < endl
< < " if (!(mxIsInt32(prhs[5]) && mxGetNumberOfElements(prhs[5]) == " < < ncols + 1 < < " )) " < < endl
< < R " ( mexErrMsgTxt( " sparse_colptr must be an int32 array with ) " << ncols+1 << R " ( elements " );) " < < endl
< < " #if MX_HAS_INTERLEAVED_COMPLEX " < < endl
< < " const int32_T *restrict sparse_rowval = mxGetInt32s(prhs[3]); " < < endl
< < " const int32_T *restrict sparse_colptr = mxGetInt32s(prhs[5]); " < < endl
< < " #else " < < endl
< < " const int32_T *restrict sparse_rowval = (int32_T *) mxGetData(prhs[3]); " < < endl
< < " const int32_T *restrict sparse_colptr = (int32_T *) mxGetData(prhs[5]); " < < endl
< < " #endif " < < endl
< < " plhs[0] = mxCreateSparse( " < < ramsey_orig_endo_nbr < < " , " < < ncols < < " , " < < nzval < < " , mxREAL); " < < endl
< < " mwIndex *restrict ir = mxGetIr(plhs[0]), *restrict jc = mxGetJc(plhs[0]); " < < endl
< < " for (mwSize i = 0; i < " < < nzval < < " ; i++) " < < endl
< < " *ir++ = *sparse_rowval++ - 1; " < < endl
< < " for (mwSize i = 0; i < " < < ncols + 1 < < " ; i++) " < < endl
< < " *jc++ = *sparse_colptr++ - 1; " < < endl
< < " mxArray *T_mx = mxCreateDoubleMatrix( " < < ramsey_multipliers_derivatives_temporary_terms . size ( ) < < " , 1, mxREAL); " < < endl
< < " ramsey_multipliers_static_g1(y, x, params, mxGetPr(T_mx), mxGetPr(plhs[0])); " < < endl
< < " mxDestroyArray(T_mx); " < < endl
< < " } " < < endl ;
output . close ( ) ;
compileMEX ( packageDir ( basename ) , " ramsey_multipliers_static_g1 " , mexext , { p } , matlabroot ) ;
}