2008-02-03 11:28:36 +01:00
/*
2020-01-07 12:26:40 +01:00
* Copyright © 2003 - 2020 Dynare Team
2008-02-03 11:28:36 +01: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/>.
*/
# include "ModelTree.hh"
2009-12-16 14:21:31 +01:00
# include "MinimumFeedbackSet.hh"
2019-04-23 14:51:14 +02:00
# pragma GCC diagnostic push
# pragma GCC diagnostic ignored "-Wold-style-cast"
2019-09-11 16:24:09 +02:00
# pragma GCC diagnostic ignored "-Wsign-compare"
# pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
2009-12-16 14:21:31 +01:00
# include <boost/graph/adjacency_list.hpp>
# include <boost/graph/max_cardinality_matching.hpp>
# include <boost/graph/strong_components.hpp>
# include <boost/graph/topological_sort.hpp>
2019-04-23 14:51:14 +02:00
# pragma GCC diagnostic pop
2009-12-16 14:21:31 +01:00
2019-11-04 15:50:26 +01:00
# ifdef __APPLE__
2019-12-20 16:59:30 +01:00
# include <mach-o / dyld.h>
2019-11-04 15:50:26 +01:00
# endif
2020-01-07 12:41:27 +01:00
# include <regex>
2009-12-16 14:21:31 +01:00
using namespace MFS ;
2018-10-09 18:27:19 +02:00
void
ModelTree : : copyHelper ( const ModelTree & m )
{
2018-10-10 13:07:25 +02:00
auto f = [ this ] ( expr_t e ) { return e - > clone ( * this ) ; } ;
2018-10-09 18:27:19 +02:00
// Equations
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . equations )
2018-10-09 18:27:19 +02:00
equations . push_back ( dynamic_cast < BinaryOpNode * > ( f ( it ) ) ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . aux_equations )
2018-10-09 18:27:19 +02:00
aux_equations . push_back ( dynamic_cast < BinaryOpNode * > ( f ( it ) ) ) ;
2018-11-15 16:39:53 +01:00
auto convert_deriv_map = [ f ] ( map < vector < int > , expr_t > dm )
2019-12-20 16:59:30 +01:00
{
map < vector < int > , expr_t > dm2 ;
for ( const auto & it : dm )
dm2 . emplace ( it . first , f ( it . second ) ) ;
return dm2 ;
} ;
2018-11-15 16:39:53 +01:00
2018-10-09 18:27:19 +02:00
// Derivatives
2018-11-15 16:39:53 +01:00
for ( const auto & it : m . derivatives )
derivatives . push_back ( convert_deriv_map ( it ) ) ;
for ( const auto & it : m . params_derivatives )
params_derivatives [ it . first ] = convert_deriv_map ( it . second ) ;
auto convert_temporary_terms_t = [ f ] ( temporary_terms_t tt )
2019-12-20 16:59:30 +01:00
{
temporary_terms_t tt2 ;
for ( const auto & it : tt )
tt2 . insert ( f ( it ) ) ;
return tt2 ;
} ;
2018-10-09 18:27:19 +02:00
// Temporary terms
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . temporary_terms )
2018-10-09 18:27:19 +02:00
temporary_terms . insert ( f ( it ) ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . temporary_terms_mlv )
2018-10-09 18:27:19 +02:00
temporary_terms_mlv [ f ( it . first ) ] = f ( it . second ) ;
2018-11-15 16:39:53 +01:00
for ( const auto & it : m . temporary_terms_derivatives )
temporary_terms_derivatives . push_back ( convert_temporary_terms_t ( it ) ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . temporary_terms_idxs )
2018-10-09 18:27:19 +02:00
temporary_terms_idxs [ f ( it . first ) ] = it . second ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . params_derivs_temporary_terms )
2018-11-16 18:24:06 +01:00
params_derivs_temporary_terms [ it . first ] = convert_temporary_terms_t ( it . second ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . params_derivs_temporary_terms_idxs )
2018-10-09 18:27:19 +02:00
params_derivs_temporary_terms_idxs [ f ( it . first ) ] = it . second ;
// Other stuff
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . trend_symbols_map )
2018-10-09 18:27:19 +02:00
trend_symbols_map [ it . first ] = f ( it . second ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : m . nonstationary_symbols_map )
2019-02-21 10:30:25 +01:00
nonstationary_symbols_map [ it . first ] = { it . second . first , f ( it . second . second ) } ;
2018-10-09 18:27:19 +02:00
}
2018-11-22 14:36:03 +01:00
ModelTree : : ModelTree ( SymbolTable & symbol_table_arg ,
NumericalConstants & num_constants_arg ,
ExternalFunctionsTable & external_functions_table_arg ,
bool is_dynamic_arg ) :
2019-12-20 16:59:30 +01:00
DataTree { symbol_table_arg , num_constants_arg , external_functions_table_arg , is_dynamic_arg } ,
2018-11-22 14:36:03 +01:00
derivatives ( 4 ) ,
NNZDerivatives ( 4 , 0 ) ,
temporary_terms_derivatives ( 4 )
{
}
2018-10-09 18:27:19 +02:00
ModelTree : : ModelTree ( const ModelTree & m ) :
2019-12-20 16:59:30 +01:00
DataTree { m } ,
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 } ,
equations_lineno { m . equations_lineno } ,
equation_tags { m . equation_tags } ,
equation_tags_xref { m . equation_tags_xref } ,
NNZDerivatives { m . NNZDerivatives } ,
equation_reordered { m . equation_reordered } ,
variable_reordered { m . variable_reordered } ,
inv_equation_reordered { m . inv_equation_reordered } ,
inv_variable_reordered { m . inv_variable_reordered } ,
is_equation_linear { m . is_equation_linear } ,
endo2eq { m . endo2eq } ,
epilogue { m . epilogue } ,
prologue { m . prologue } ,
block_lag_lead { m . block_lag_lead } ,
cutoff { m . cutoff } ,
mfs { m . mfs }
2018-10-09 18:27:19 +02:00
{
copyHelper ( m ) ;
}
ModelTree &
ModelTree : : operator = ( const ModelTree & m )
{
DataTree : : operator = ( m ) ;
equations . clear ( ) ;
equations_lineno = m . equations_lineno ;
aux_equations . clear ( ) ;
equation_tags = m . equation_tags ;
2019-10-29 12:56:53 +01:00
equation_tags_xref = m . equation_tags_xref ;
2018-10-09 18:27:19 +02:00
NNZDerivatives = m . NNZDerivatives ;
2018-11-15 16:39:53 +01:00
derivatives . clear ( ) ;
params_derivatives . clear ( ) ;
2018-10-09 18:27:19 +02:00
temporary_terms . clear ( ) ;
temporary_terms_mlv . clear ( ) ;
2018-11-15 16:39:53 +01:00
temporary_terms_derivatives . clear ( ) ;
2018-10-09 18:27:19 +02:00
params_derivs_temporary_terms . clear ( ) ;
params_derivs_temporary_terms_idxs . clear ( ) ;
trend_symbols_map . clear ( ) ;
nonstationary_symbols_map . clear ( ) ;
equation_reordered = m . equation_reordered ;
variable_reordered = m . variable_reordered ;
inv_equation_reordered = m . inv_equation_reordered ;
inv_variable_reordered = m . inv_variable_reordered ;
is_equation_linear = m . is_equation_linear ;
endo2eq = m . endo2eq ;
epilogue = m . epilogue ;
prologue = m . prologue ;
block_lag_lead = m . block_lag_lead ;
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-09 18:27:19 +02:00
copyHelper ( m ) ;
return * this ;
}
2009-12-21 11:29:21 +01:00
bool
2010-09-16 19:00:48 +02:00
ModelTree : : computeNormalization ( const jacob_map_t & contemporaneous_jacobian , bool verbose )
2009-12-16 14:21:31 +01:00
{
2017-02-24 11:20:54 +01:00
const int n = equations . size ( ) ;
2009-12-16 14:21:31 +01:00
assert ( n = = symbol_table . endo_nbr ( ) ) ;
2018-06-05 14:56:34 +02:00
using BipartiteGraph = boost : : adjacency_list < boost : : vecS , boost : : vecS , boost : : undirectedS > ;
2009-12-16 14:21:31 +01:00
/*
Vertices 0 to n - 1 are for endogenous ( using type specific ID )
Vertices n to 2 * n - 1 are for equations ( using equation no . )
*/
BipartiteGraph g ( 2 * n ) ;
// Fill in the graph
2018-06-04 14:17:36 +02:00
set < pair < int , int > > endo ;
2009-12-16 14:21:31 +01:00
2019-12-20 16:59:30 +01:00
for ( const auto & it : contemporaneous_jacobian )
2018-06-04 12:26:16 +02:00
add_edge ( it . first . first + n , it . first . second , g ) ;
2009-12-16 14:21:31 +01:00
// Compute maximum cardinality matching
vector < int > mate_map ( 2 * n ) ;
# if 1
bool check = checked_edmonds_maximum_cardinality_matching ( g , & mate_map [ 0 ] ) ;
# else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
2018-11-23 17:19:59 +01:00
fill ( mate_map . begin ( ) , mate_map . end ( ) , boost : : graph_traits < BipartiteGraph > : : null_vertex ( ) ) ;
2009-12-16 14:21:31 +01:00
2019-12-16 19:42:59 +01:00
auto natural_endo2eqs = computeNormalizedEquations ( ) ;
2009-12-16 14:21:31 +01:00
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < symbol_table . endo_nbr ( ) ; i + + )
2009-12-16 14:21:31 +01:00
{
if ( natural_endo2eqs . count ( i ) = = 0 )
continue ;
int j = natural_endo2eqs . find ( i ) - > second ;
put ( & mate_map [ 0 ] , i , n + j ) ;
put ( & mate_map [ 0 ] , n + j , i ) ;
}
2018-11-23 17:19:59 +01:00
boost : : edmonds_augmenting_path_finder < BipartiteGraph , int * , boost : : property_map < BipartiteGraph , boost : : vertex_index_t > : : type > augmentor ( g , & mate_map [ 0 ] , get ( boost : : vertex_index , g ) ) ;
while ( augmentor . augment_matching ( ) )
2009-12-16 14:21:31 +01:00
{
2018-11-23 17:19:59 +01:00
} ;
2009-12-16 14:21:31 +01:00
augmentor . get_current_matching ( & mate_map [ 0 ] ) ;
2018-11-23 17:19:59 +01:00
bool check = boost : : maximum_cardinality_matching_verifier < BipartiteGraph , int * , boost : : property_map < BipartiteGraph , boost : : vertex_index_t > : : type > : : verify_matching ( g , & mate_map [ 0 ] , get ( boost : : vertex_index , g ) ) ;
2009-12-16 14:21:31 +01:00
# endif
assert ( check ) ;
# ifdef DEBUG
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < n ; i + + )
2009-12-16 14:21:31 +01:00
cout < < " Endogenous " < < symbol_table . getName ( symbol_table . getID ( eEndogenous , i ) )
< < " matched with equation " < < ( mate_map [ i ] - n + 1 ) < < endl ;
# endif
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
2017-02-24 11:20:54 +01:00
endo2eq . resize ( equations . size ( ) ) ;
2018-11-23 17:19:59 +01:00
transform ( mate_map . begin ( ) , mate_map . begin ( ) + n , endo2eq . begin ( ) , [ = ] ( int i ) { return i - n ; } ) ;
2009-12-16 14:21:31 +01:00
# ifdef DEBUG
2019-12-16 19:42:59 +01:00
auto natural_endo2eqs = computeNormalizedEquations ( natural_endo2eqs ) ;
2009-12-16 14:21:31 +01:00
int n1 = 0 , n2 = 0 ;
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < symbol_table . endo_nbr ( ) ; i + + )
2009-12-16 14:21:31 +01:00
{
if ( natural_endo2eqs . count ( i ) = = 0 )
continue ;
n1 + + ;
2018-11-23 17:19:59 +01:00
auto x = natural_endo2eqs . equal_range ( i ) ;
if ( find_if ( x . first , x . second , [ = ] ( auto y ) { return y . second = = endo2eq [ i ] ; } ) = = x . second )
cout < < " Natural normalization of variable " < < symbol_table . getName ( symbol_table . getID ( SymbolType : : endogenous , i ) )
2009-12-16 14:21:31 +01:00
< < " not used. " < < endl ;
else
n2 + + ;
}
cout < < " Used " < < n2 < < " natural normalizations out of " < < n1 < < " , for a total of " < < n < < " equations. " < < endl ;
# endif
// Check if all variables are normalized
2019-12-16 19:42:59 +01:00
if ( auto it = find ( mate_map . begin ( ) , mate_map . begin ( ) + n , boost : : graph_traits < BipartiteGraph > : : null_vertex ( ) ) ;
it ! = mate_map . begin ( ) + n )
2009-12-21 11:29:21 +01:00
{
if ( verbose )
cerr < < " ERROR: Could not normalize the model. Variable "
2018-07-17 18:34:07 +02:00
< < symbol_table . getName ( symbol_table . getID ( SymbolType : : endogenous , it - mate_map . begin ( ) ) )
2009-12-21 11:29:21 +01:00
< < " is not in the maximum cardinality matching. " < < endl ;
check = false ;
}
return check ;
2009-12-16 14:21:31 +01:00
}
void
2010-09-16 19:00:48 +02:00
ModelTree : : computeNonSingularNormalization ( jacob_map_t & contemporaneous_jacobian , double cutoff , jacob_map_t & static_jacobian , dynamic_jacob_map_t & dynamic_jacobian )
2009-12-16 14:21:31 +01:00
{
2009-12-21 11:29:21 +01:00
bool check = false ;
2009-12-16 14:21:31 +01:00
cout < < " Normalizing the model... " < < endl ;
2017-02-24 11:20:54 +01:00
int n = equations . size ( ) ;
2009-12-16 14:21:31 +01:00
2009-12-21 11:29:21 +01:00
// compute the maximum value of each row of the contemporaneous Jacobian matrix
//jacob_map normalized_contemporaneous_jacobian;
2010-09-16 19:00:48 +02:00
jacob_map_t normalized_contemporaneous_jacobian ( contemporaneous_jacobian ) ;
2009-12-21 11:29:21 +01:00
vector < double > max_val ( n , 0.0 ) ;
2018-11-23 17:19:59 +01:00
for ( const auto & it : contemporaneous_jacobian )
if ( fabs ( it . second ) > max_val [ it . first . first ] )
max_val [ it . first . first ] = fabs ( it . second ) ;
2009-12-16 14:21:31 +01:00
2019-12-20 16:59:30 +01:00
for ( auto & iter : normalized_contemporaneous_jacobian )
2018-06-04 12:26:16 +02:00
iter . second / = max_val [ iter . first . first ] ;
2009-12-21 11:29:21 +01:00
//We start with the highest value of the cutoff and try to normalize the model
double current_cutoff = 0.99999999 ;
int suppressed = 0 ;
while ( ! check & & current_cutoff > 1e-19 )
2009-12-16 14:21:31 +01:00
{
2010-09-16 19:00:48 +02:00
jacob_map_t tmp_normalized_contemporaneous_jacobian ;
2009-12-21 11:29:21 +01:00
int suppress = 0 ;
2019-12-20 16:59:30 +01:00
for ( auto & iter : normalized_contemporaneous_jacobian )
2018-06-04 12:26:16 +02:00
if ( fabs ( iter . second ) > max ( current_cutoff , cutoff ) )
2018-06-04 16:36:46 +02:00
tmp_normalized_contemporaneous_jacobian [ { iter . first . first , iter . first . second } ] = iter . second ;
2009-12-21 11:29:21 +01:00
else
suppress + + ;
if ( suppress ! = suppressed )
check = computeNormalization ( tmp_normalized_contemporaneous_jacobian , false ) ;
suppressed = suppress ;
if ( ! check )
2009-12-16 14:21:31 +01:00
{
2009-12-21 11:29:21 +01:00
current_cutoff / = 2 ;
// In this last case try to normalize with the complete jacobian
if ( current_cutoff < = 1e-19 )
check = computeNormalization ( normalized_contemporaneous_jacobian , false ) ;
2009-12-16 14:21:31 +01:00
}
}
2009-12-21 11:29:21 +01:00
if ( ! check )
2009-12-16 14:21:31 +01:00
{
2009-12-21 11:29:21 +01:00
cout < < " Normalization failed with cutoff, trying symbolic normalization... " < < endl ;
//if no non-singular normalization can be found, try to find a normalization even with a potential singularity
2010-09-16 19:00:48 +02:00
jacob_map_t tmp_normalized_contemporaneous_jacobian ;
2018-06-04 14:17:36 +02:00
set < pair < int , int > > endo ;
2009-12-21 11:29:21 +01:00
for ( int i = 0 ; i < n ; i + + )
2009-12-16 14:21:31 +01:00
{
endo . clear ( ) ;
equations [ i ] - > collectEndogenous ( endo ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : endo )
2018-06-04 16:36:46 +02:00
tmp_normalized_contemporaneous_jacobian [ { i , it . first } ] = 1 ;
2009-12-16 14:21:31 +01:00
}
2009-12-21 11:29:21 +01:00
check = computeNormalization ( tmp_normalized_contemporaneous_jacobian , true ) ;
if ( check )
2009-12-16 14:21:31 +01:00
{
2009-12-21 11:29:21 +01:00
// Update the jacobian matrix
2019-12-16 19:42:59 +01:00
for ( const auto & [ key , ignore ] : tmp_normalized_contemporaneous_jacobian )
2009-12-21 11:29:21 +01:00
{
2019-12-16 19:42:59 +01:00
if ( static_jacobian . find ( { key . first , key . second } ) = = static_jacobian . end ( ) )
static_jacobian [ { key . first , key . second } ] = 0 ;
if ( dynamic_jacobian . find ( { 0 , key . first , key . second } ) = = dynamic_jacobian . end ( ) )
dynamic_jacobian [ { 0 , key . first , key . second } ] = nullptr ;
if ( contemporaneous_jacobian . find ( { key . first , key . second } ) = = contemporaneous_jacobian . end ( ) )
contemporaneous_jacobian [ { key . first , key . second } ] = 0 ;
2014-07-15 17:32:12 +02:00
try
{
2019-12-16 19:42:59 +01:00
if ( derivatives [ 1 ] . find ( { key . first , getDerivID ( symbol_table . getID ( SymbolType : : endogenous , key . second ) , 0 ) } ) = = derivatives [ 1 ] . end ( ) )
derivatives [ 1 ] [ { key . first , getDerivID ( symbol_table . getID ( SymbolType : : endogenous , key . second ) , 0 ) } ] = Zero ;
2014-07-15 17:32:12 +02:00
}
2017-06-14 07:01:31 +02:00
catch ( DataTree : : UnknownDerivIDException & e )
2014-07-15 17:32:12 +02:00
{
2019-12-16 19:42:59 +01:00
cerr < < " The variable " < < symbol_table . getName ( symbol_table . getID ( SymbolType : : endogenous , key . second ) )
2014-07-15 17:32:12 +02:00
< < " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block. " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2009-12-21 11:29:21 +01:00
}
2009-12-16 14:21:31 +01:00
}
}
2009-12-21 11:29:21 +01:00
if ( ! check )
{
cerr < < " No normalization could be computed. Aborting. " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2009-12-16 14:21:31 +01:00
}
2019-12-16 19:42:59 +01:00
multimap < int , int >
ModelTree : : computeNormalizedEquations ( ) const
2009-12-16 14:21:31 +01:00
{
2019-12-16 19:42:59 +01:00
multimap < int , int > endo2eqs ;
2017-06-14 23:49:10 +02:00
for ( size_t i = 0 ; i < equations . size ( ) ; i + + )
2009-12-16 14:21:31 +01:00
{
2019-12-16 19:42:59 +01:00
auto lhs = dynamic_cast < VariableNode * > ( equations [ i ] - > arg1 ) ;
if ( ! lhs )
2009-12-16 14:21:31 +01:00
continue ;
2018-11-28 14:32:26 +01:00
int symb_id = lhs - > symb_id ;
2018-07-17 18:34:07 +02:00
if ( symbol_table . getType ( symb_id ) ! = SymbolType : : endogenous )
2009-12-16 14:21:31 +01:00
continue ;
2018-06-04 14:17:36 +02:00
set < pair < int , int > > endo ;
2018-11-28 14:32:26 +01:00
equations [ i ] - > arg2 - > collectEndogenous ( endo ) ;
2018-06-04 16:36:46 +02:00
if ( endo . find ( { symbol_table . getTypeSpecificID ( symb_id ) , 0 } ) ! = endo . end ( ) )
2009-12-16 14:21:31 +01:00
continue ;
2019-04-23 11:07:32 +02:00
endo2eqs . emplace ( symbol_table . getTypeSpecificID ( symb_id ) , static_cast < int > ( i ) ) ;
2019-12-16 19:42:59 +01:00
cout < < " Endogenous " < < symbol_table . getName ( symb_id ) < < " normalized in equation " < < i + 1 < < endl ;
2009-12-16 14:21:31 +01:00
}
2019-12-16 19:42:59 +01:00
return endo2eqs ;
2009-12-16 14:21:31 +01:00
}
void
2010-09-16 19:00:48 +02:00
ModelTree : : evaluateAndReduceJacobian ( const eval_context_t & eval_context , jacob_map_t & contemporaneous_jacobian , jacob_map_t & static_jacobian , dynamic_jacob_map_t & dynamic_jacobian , double cutoff , bool verbose )
2009-12-16 14:21:31 +01:00
{
int nb_elements_contemparenous_Jacobian = 0 ;
2018-11-15 16:39:53 +01:00
set < vector < int > > jacobian_elements_to_delete ;
2019-12-16 19:42:59 +01:00
for ( const auto & [ indices , d1 ] : derivatives [ 1 ] )
2009-12-16 14:21:31 +01:00
{
2019-12-16 19:42:59 +01:00
int deriv_id = indices [ 1 ] ;
2018-07-17 18:34:07 +02:00
if ( getTypeByDerivID ( deriv_id ) = = SymbolType : : endogenous )
2009-12-16 14:21:31 +01:00
{
2019-12-16 19:42:59 +01:00
int eq = indices [ 0 ] ;
2009-12-16 14:21:31 +01:00
int symb = getSymbIDByDerivID ( deriv_id ) ;
int var = symbol_table . getTypeSpecificID ( symb ) ;
int lag = getLagByDerivID ( deriv_id ) ;
double val = 0 ;
try
{
2019-12-16 19:42:59 +01:00
val = d1 - > eval ( eval_context ) ;
2009-12-16 14:21:31 +01:00
}
2010-12-10 11:50:27 +01:00
catch ( ExprNode : : EvalExternalFunctionException & e )
{
val = 1 ;
}
2009-12-16 14:21:31 +01:00
catch ( ExprNode : : EvalException & e )
{
2014-01-27 16:41:43 +01:00
cerr < < " ERROR: evaluation of Jacobian failed for equation " < < eq + 1 < < " (line " < < equations_lineno [ eq ] < < " ) and variable " < < symbol_table . getName ( symb ) < < " ( " < < lag < < " ) [ " < < symb < < " ] ! " < < endl ;
2019-12-16 19:42:59 +01:00
d1 - > writeOutput ( cerr , ExprNodeOutputType : : matlabDynamicModelSparse , temporary_terms , { } ) ;
2009-12-16 14:21:31 +01:00
cerr < < endl ;
exit ( EXIT_FAILURE ) ;
}
if ( fabs ( val ) < cutoff )
{
if ( verbose )
cout < < " the coefficient related to variable " < < var < < " with lag " < < lag < < " in equation " < < eq < < " is equal to " < < val < < " and is set to 0 in the incidence matrix (size= " < < symbol_table . endo_nbr ( ) < < " ) " < < endl ;
2019-03-08 18:57:17 +01:00
jacobian_elements_to_delete . insert ( { eq , deriv_id } ) ;
2009-12-16 14:21:31 +01:00
}
else
{
if ( lag = = 0 )
{
nb_elements_contemparenous_Jacobian + + ;
2018-06-04 16:36:46 +02:00
contemporaneous_jacobian [ { eq , var } ] = val ;
2009-12-16 14:21:31 +01:00
}
2018-06-04 16:36:46 +02:00
if ( static_jacobian . find ( { eq , var } ) ! = static_jacobian . end ( ) )
static_jacobian [ { eq , var } ] + = val ;
2009-12-16 14:21:31 +01:00
else
2018-06-04 16:36:46 +02:00
static_jacobian [ { eq , var } ] = val ;
2019-12-16 19:42:59 +01:00
dynamic_jacobian [ { lag , eq , var } ] = d1 ;
2009-12-16 14:21:31 +01:00
}
}
}
// Get rid of the elements of the Jacobian matrix below the cutoff
2019-12-20 16:59:30 +01:00
for ( const auto & it : jacobian_elements_to_delete )
2018-11-15 16:39:53 +01:00
derivatives [ 1 ] . erase ( it ) ;
2009-12-16 14:21:31 +01:00
2009-12-16 18:13:23 +01:00
if ( jacobian_elements_to_delete . size ( ) > 0 )
2009-12-16 14:21:31 +01:00
{
2018-11-15 16:39:53 +01:00
cout < < jacobian_elements_to_delete . size ( ) < < " elements among " < < derivatives [ 1 ] . size ( ) < < " in the incidence matrices are below the cutoff ( " < < cutoff < < " ) and are discarded " < < endl
2009-12-16 14:21:31 +01:00
< < " The contemporaneous incidence matrix has " < < nb_elements_contemparenous_Jacobian < < " elements " < < endl ;
}
}
2019-12-20 16:59:30 +01:00
vector < pair < int , int > >
2018-09-21 17:13:19 +02:00
ModelTree : : select_non_linear_equations_and_variables ( vector < bool > is_equation_linear , const dynamic_jacob_map_t & dynamic_jacobian , vector < int > & equation_reordered , vector < int > & variable_reordered ,
vector < int > & inv_equation_reordered , vector < int > & inv_variable_reordered ,
lag_lead_vector_t & equation_lag_lead , lag_lead_vector_t & variable_lag_lead ,
vector < unsigned int > & n_static , vector < unsigned int > & n_forward , vector < unsigned int > & n_backward , vector < unsigned int > & n_mixed )
{
vector < int > eq2endo ( equations . size ( ) , 0 ) ;
/*equation_reordered.resize(equations.size());
2019-12-20 16:59:30 +01:00
variable_reordered . resize ( equations . size ( ) ) ; */
2018-09-21 17:13:19 +02:00
unsigned int num = 0 ;
2018-11-23 17:19:59 +01:00
for ( auto it : endo2eq )
if ( ! is_equation_linear [ it ] )
2018-09-21 17:13:19 +02:00
num + + ;
vector < int > endo2block = vector < int > ( endo2eq . size ( ) , 1 ) ;
2019-12-16 19:42:59 +01:00
vector < pair < set < int > , pair < set < int > , vector < int > > > > components_set ( num ) ;
2018-09-21 17:13:19 +02:00
int i = 0 , j = 0 ;
2018-11-23 17:19:59 +01:00
for ( auto it : endo2eq )
2019-12-16 19:42:59 +01:00
if ( ! is_equation_linear [ it ] )
{
equation_reordered [ i ] = it ;
variable_reordered [ i ] = j ;
endo2block [ j ] = 0 ;
components_set [ endo2block [ j ] ] . first . insert ( i ) ;
i + + ;
j + + ;
}
2018-09-21 17:13:19 +02:00
getVariableLeadLagByBlock ( dynamic_jacobian , endo2block , endo2block . size ( ) , equation_lag_lead , variable_lag_lead , equation_reordered , variable_reordered ) ;
n_static = vector < unsigned int > ( endo2eq . size ( ) , 0 ) ;
n_forward = vector < unsigned int > ( endo2eq . size ( ) , 0 ) ;
n_backward = vector < unsigned int > ( endo2eq . size ( ) , 0 ) ;
n_mixed = vector < unsigned int > ( endo2eq . size ( ) , 0 ) ;
for ( unsigned int i = 0 ; i < endo2eq . size ( ) ; i + + )
{
2019-12-16 19:42:59 +01:00
if ( variable_lag_lead [ variable_reordered [ i ] ] . first ! = 0 & & variable_lag_lead [ variable_reordered [ i ] ] . second ! = 0 )
2018-09-21 17:13:19 +02:00
n_mixed [ i ] + + ;
else if ( variable_lag_lead [ variable_reordered [ i ] ] . first = = 0 & & variable_lag_lead [ variable_reordered [ i ] ] . second ! = 0 )
n_forward [ i ] + + ;
else if ( variable_lag_lead [ variable_reordered [ i ] ] . first ! = 0 & & variable_lag_lead [ variable_reordered [ i ] ] . second = = 0 )
n_backward [ i ] + + ;
else if ( variable_lag_lead [ variable_reordered [ i ] ] . first = = 0 & & variable_lag_lead [ variable_reordered [ i ] ] . second = = 0 )
n_static [ i ] + + ;
}
cout . flush ( ) ;
int nb_endo = is_equation_linear . size ( ) ;
2019-02-21 10:30:25 +01:00
vector < pair < int , int > > blocks ( 1 , { i , i } ) ;
2018-11-23 17:19:59 +01:00
inv_equation_reordered . resize ( nb_endo ) ;
inv_variable_reordered . resize ( nb_endo ) ;
2018-09-21 17:13:19 +02:00
for ( int i = 0 ; i < nb_endo ; i + + )
{
inv_variable_reordered [ variable_reordered [ i ] ] = i ;
inv_equation_reordered [ equation_reordered [ i ] ] = i ;
}
return blocks ;
}
bool
ModelTree : : computeNaturalNormalization ( )
{
bool bool_result = true ;
2019-12-20 16:59:30 +01:00
set < pair < int , int > > result ;
2018-09-21 17:13:19 +02:00
endo2eq . resize ( equations . size ( ) ) ;
2019-04-23 11:07:32 +02:00
for ( int eq = 0 ; eq < static_cast < int > ( equations . size ( ) ) ; eq + + )
2018-09-21 17:13:19 +02:00
if ( ! is_equation_linear [ eq ] )
{
BinaryOpNode * eq_node = equations [ eq ] ;
2018-11-28 14:32:26 +01:00
expr_t lhs = eq_node - > arg1 ;
2018-09-21 17:13:19 +02:00
result . clear ( ) ;
lhs - > collectDynamicVariables ( SymbolType : : endogenous , result ) ;
if ( result . size ( ) = = 1 & & result . begin ( ) - > second = = 0 )
{
//check if the endogenous variable has not been already used in an other match !
2019-12-16 19:42:59 +01:00
if ( find ( endo2eq . begin ( ) , endo2eq . end ( ) , result . begin ( ) - > first ) = = endo2eq . end ( ) )
endo2eq [ result . begin ( ) - > first ] = eq ;
else
2018-09-21 17:13:19 +02:00
{
bool_result = false ;
break ;
}
}
}
return bool_result ;
}
2009-12-16 14:21:31 +01:00
void
2010-09-16 19:00:48 +02:00
ModelTree : : computePrologueAndEpilogue ( const jacob_map_t & static_jacobian_arg , vector < int > & equation_reordered , vector < int > & variable_reordered )
2009-12-16 14:21:31 +01:00
{
2017-02-24 11:20:54 +01:00
vector < int > eq2endo ( equations . size ( ) , 0 ) ;
equation_reordered . resize ( equations . size ( ) ) ;
variable_reordered . resize ( equations . size ( ) ) ;
int n = equations . size ( ) ;
2018-11-23 17:19:59 +01:00
vector < bool > IM ( n * n ) ;
2009-12-16 14:21:31 +01:00
int i = 0 ;
2018-11-23 17:19:59 +01:00
for ( auto it : endo2eq )
2009-12-16 14:21:31 +01:00
{
2018-11-23 17:19:59 +01:00
eq2endo [ it ] = i ;
2009-12-16 14:21:31 +01:00
equation_reordered [ i ] = i ;
2018-11-23 17:19:59 +01:00
variable_reordered [ it ] = i ;
i + + ;
2017-06-14 07:01:31 +02:00
}
if ( cutoff = = 0 )
{
2018-06-04 14:17:36 +02:00
set < pair < int , int > > endo ;
2016-03-25 15:38:49 +01:00
for ( int i = 0 ; i < n ; i + + )
2017-06-14 07:01:31 +02:00
{
2016-03-25 15:38:49 +01:00
endo . clear ( ) ;
2017-06-14 07:01:31 +02:00
equations [ i ] - > collectEndogenous ( endo ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : endo )
2018-06-04 12:26:16 +02:00
IM [ i * n + endo2eq [ it . first ] ] = true ;
2017-06-14 07:01:31 +02:00
}
}
else
2019-12-20 16:59:30 +01:00
for ( const auto & it : static_jacobian_arg )
2018-06-04 12:26:16 +02:00
IM [ it . first . first * n + endo2eq [ it . first . second ] ] = true ;
2009-12-16 14:21:31 +01:00
bool something_has_been_done = true ;
prologue = 0 ;
int k = 0 ;
// Find the prologue equations and place first the AR(1) shock equations first
while ( something_has_been_done )
{
int tmp_prologue = prologue ;
something_has_been_done = false ;
2009-12-16 18:13:23 +01:00
for ( int i = prologue ; i < n ; i + + )
2009-12-16 14:21:31 +01:00
{
int nze = 0 ;
2009-12-16 18:13:23 +01:00
for ( int j = tmp_prologue ; j < n ; j + + )
if ( IM [ i * n + j ] )
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
nze + + ;
2009-12-16 14:21:31 +01:00
k = j ;
}
2009-12-16 18:13:23 +01:00
if ( nze = = 1 )
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
for ( int j = 0 ; j < n ; j + + )
2009-12-16 14:21:31 +01:00
{
bool tmp_bool = IM [ tmp_prologue * n + j ] ;
IM [ tmp_prologue * n + j ] = IM [ i * n + j ] ;
IM [ i * n + j ] = tmp_bool ;
}
int tmp = equation_reordered [ tmp_prologue ] ;
equation_reordered [ tmp_prologue ] = equation_reordered [ i ] ;
equation_reordered [ i ] = tmp ;
2009-12-16 18:13:23 +01:00
for ( int j = 0 ; j < n ; j + + )
2009-12-16 14:21:31 +01:00
{
bool tmp_bool = IM [ j * n + tmp_prologue ] ;
IM [ j * n + tmp_prologue ] = IM [ j * n + k ] ;
IM [ j * n + k ] = tmp_bool ;
}
tmp = variable_reordered [ tmp_prologue ] ;
variable_reordered [ tmp_prologue ] = variable_reordered [ k ] ;
variable_reordered [ k ] = tmp ;
tmp_prologue + + ;
something_has_been_done = true ;
}
}
prologue = tmp_prologue ;
}
something_has_been_done = true ;
epilogue = 0 ;
// Find the epilogue equations
while ( something_has_been_done )
{
int tmp_epilogue = epilogue ;
something_has_been_done = false ;
2019-04-23 11:07:32 +02:00
for ( int i = prologue ; i < n - static_cast < int > ( epilogue ) ; i + + )
2009-12-16 14:21:31 +01:00
{
int nze = 0 ;
2009-12-16 18:13:23 +01:00
for ( int j = prologue ; j < n - tmp_epilogue ; j + + )
if ( IM [ j * n + i ] )
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
nze + + ;
2009-12-16 14:21:31 +01:00
k = j ;
}
2009-12-16 18:13:23 +01:00
if ( nze = = 1 )
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
for ( int j = 0 ; j < n ; j + + )
2009-12-16 14:21:31 +01:00
{
bool tmp_bool = IM [ ( n - 1 - tmp_epilogue ) * n + j ] ;
IM [ ( n - 1 - tmp_epilogue ) * n + j ] = IM [ k * n + j ] ;
IM [ k * n + j ] = tmp_bool ;
}
int tmp = equation_reordered [ n - 1 - tmp_epilogue ] ;
equation_reordered [ n - 1 - tmp_epilogue ] = equation_reordered [ k ] ;
equation_reordered [ k ] = tmp ;
2009-12-16 18:13:23 +01:00
for ( int j = 0 ; j < n ; j + + )
2009-12-16 14:21:31 +01:00
{
bool tmp_bool = IM [ j * n + n - 1 - tmp_epilogue ] ;
IM [ j * n + n - 1 - tmp_epilogue ] = IM [ j * n + i ] ;
IM [ j * n + i ] = tmp_bool ;
}
tmp = variable_reordered [ n - 1 - tmp_epilogue ] ;
variable_reordered [ n - 1 - tmp_epilogue ] = variable_reordered [ i ] ;
variable_reordered [ i ] = tmp ;
tmp_epilogue + + ;
something_has_been_done = true ;
}
}
epilogue = tmp_epilogue ;
}
}
2010-09-16 19:00:48 +02:00
equation_type_and_normalized_equation_t
2018-11-23 17:19:59 +01:00
ModelTree : : equationTypeDetermination ( const map < tuple < int , int , int > , expr_t > & first_order_endo_derivatives , const vector < int > & Index_Var_IM , const vector < int > & Index_Equ_IM , int mfs ) const
2009-12-16 14:21:31 +01:00
{
2011-08-19 15:17:04 +02:00
expr_t lhs ;
2009-12-16 14:21:31 +01:00
BinaryOpNode * eq_node ;
EquationType Equation_Simulation_Type ;
2010-09-16 19:00:48 +02:00
equation_type_and_normalized_equation_t V_Equation_Simulation_Type ( equations . size ( ) ) ;
2009-12-16 14:21:31 +01:00
for ( unsigned int i = 0 ; i < equations . size ( ) ; i + + )
{
int eq = Index_Equ_IM [ i ] ;
int var = Index_Var_IM [ i ] ;
eq_node = equations [ eq ] ;
2018-11-28 14:32:26 +01:00
lhs = eq_node - > arg1 ;
2009-12-16 14:21:31 +01:00
Equation_Simulation_Type = E_SOLVE ;
2018-11-23 17:19:59 +01:00
auto derivative = first_order_endo_derivatives . find ( { eq , var , 0 } ) ;
2010-09-16 19:18:45 +02:00
pair < bool , expr_t > res ;
2009-12-16 18:13:23 +01:00
if ( derivative ! = first_order_endo_derivatives . end ( ) )
2009-12-16 14:21:31 +01:00
{
2018-06-04 14:17:36 +02:00
set < pair < int , int > > result ;
2009-12-16 14:21:31 +01:00
derivative - > second - > collectEndogenous ( result ) ;
2018-06-04 16:36:46 +02:00
auto d_endo_variable = result . find ( { var , 0 } ) ;
2009-12-16 14:21:31 +01:00
//Determine whether the equation could be evaluated rather than to be solved
2018-07-17 18:34:07 +02:00
if ( lhs - > isVariableNodeEqualTo ( SymbolType : : endogenous , Index_Var_IM [ i ] , 0 ) & & derivative - > second - > isNumConstNodeEqualTo ( 1 ) )
2018-11-23 17:19:59 +01:00
Equation_Simulation_Type = E_EVALUATE ;
2009-12-16 14:21:31 +01:00
else
{
2018-12-03 15:05:49 +01:00
vector < tuple < int , expr_t , expr_t > > List_of_Op_RHS ;
2019-12-20 16:59:30 +01:00
res = equations [ eq ] - > normalizeEquation ( var , List_of_Op_RHS ) ;
2009-12-16 18:13:23 +01:00
if ( mfs = = 2 )
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
if ( d_endo_variable = = result . end ( ) & & res . second )
2009-12-16 14:21:31 +01:00
Equation_Simulation_Type = E_EVALUATE_S ;
}
2009-12-16 18:13:23 +01:00
else if ( mfs = = 3 )
2009-12-16 14:21:31 +01:00
{
2009-12-16 18:13:23 +01:00
if ( res . second ) // The equation could be solved analytically
2009-12-16 14:21:31 +01:00
Equation_Simulation_Type = E_EVALUATE_S ;
}
}
}
2018-06-04 16:36:46 +02:00
V_Equation_Simulation_Type [ eq ] = { Equation_Simulation_Type , dynamic_cast < BinaryOpNode * > ( res . second ) } ;
2009-12-16 14:21:31 +01:00
}
2018-11-23 17:19:59 +01:00
return V_Equation_Simulation_Type ;
2009-12-16 14:21:31 +01:00
}
void
2010-09-16 19:00:48 +02:00
ModelTree : : getVariableLeadLagByBlock ( const dynamic_jacob_map_t & dynamic_jacobian , const vector < int > & components_set , int nb_blck_sim , lag_lead_vector_t & equation_lead_lag , lag_lead_vector_t & variable_lead_lag , const vector < int > & equation_reordered , const vector < int > & variable_reordered ) const
2009-12-16 14:21:31 +01:00
{
int nb_endo = symbol_table . endo_nbr ( ) ;
2018-06-04 16:36:46 +02:00
variable_lead_lag = lag_lead_vector_t ( nb_endo , { 0 , 0 } ) ;
equation_lead_lag = lag_lead_vector_t ( nb_endo , { 0 , 0 } ) ;
2009-12-16 14:21:31 +01:00
vector < int > variable_blck ( nb_endo ) , equation_blck ( nb_endo ) ;
for ( int i = 0 ; i < nb_endo ; i + + )
{
2019-04-23 11:07:32 +02:00
if ( i < static_cast < int > ( prologue ) )
2009-12-16 14:21:31 +01:00
{
variable_blck [ variable_reordered [ i ] ] = i ;
equation_blck [ equation_reordered [ i ] ] = i ;
}
2019-04-23 11:07:32 +02:00
else if ( i < static_cast < int > ( components_set . size ( ) + prologue ) )
2009-12-16 14:21:31 +01:00
{
variable_blck [ variable_reordered [ i ] ] = components_set [ i - prologue ] + prologue ;
equation_blck [ equation_reordered [ i ] ] = components_set [ i - prologue ] + prologue ;
}
else
{
variable_blck [ variable_reordered [ i ] ] = i - ( nb_endo - nb_blck_sim - prologue - epilogue ) ;
equation_blck [ equation_reordered [ i ] ] = i - ( nb_endo - nb_blck_sim - prologue - epilogue ) ;
}
}
2019-12-20 16:59:30 +01:00
for ( const auto & it : dynamic_jacobian )
2009-12-16 14:21:31 +01:00
{
2019-09-11 15:59:23 +02:00
auto [ lag , j_1 , i_1 ] = it . first ;
2009-12-16 14:21:31 +01:00
if ( variable_blck [ i_1 ] = = equation_blck [ j_1 ] )
{
if ( lag > variable_lead_lag [ i_1 ] . second )
2018-06-04 16:36:46 +02:00
variable_lead_lag [ i_1 ] = { variable_lead_lag [ i_1 ] . first , lag } ;
2009-12-16 14:21:31 +01:00
if ( lag < - variable_lead_lag [ i_1 ] . first )
2018-06-04 16:36:46 +02:00
variable_lead_lag [ i_1 ] = { - lag , variable_lead_lag [ i_1 ] . second } ;
2009-12-16 14:21:31 +01:00
if ( lag > equation_lead_lag [ j_1 ] . second )
2018-06-04 16:36:46 +02:00
equation_lead_lag [ j_1 ] = { equation_lead_lag [ j_1 ] . first , lag } ;
2009-12-16 14:21:31 +01:00
if ( lag < - equation_lead_lag [ j_1 ] . first )
2018-06-04 16:36:46 +02:00
equation_lead_lag [ j_1 ] = { - lag , equation_lead_lag [ j_1 ] . second } ;
2009-12-16 14:21:31 +01:00
}
}
}
void
2018-06-04 14:17:36 +02:00
ModelTree : : computeBlockDecompositionAndFeedbackVariablesForEachBlock ( const jacob_map_t & static_jacobian , const dynamic_jacob_map_t & dynamic_jacobian , vector < int > & equation_reordered , vector < int > & variable_reordered , vector < pair < int , int > > & blocks , const equation_type_and_normalized_equation_t & Equation_Type , bool verbose_ , bool select_feedback_variable , int mfs , vector < int > & inv_equation_reordered , vector < int > & inv_variable_reordered , lag_lead_vector_t & equation_lag_lead , lag_lead_vector_t & variable_lag_lead , vector < unsigned int > & n_static , vector < unsigned int > & n_forward , vector < unsigned int > & n_backward , vector < unsigned int > & n_mixed ) const
2009-12-16 14:21:31 +01:00
{
int nb_var = variable_reordered . size ( ) ;
int n = nb_var - prologue - epilogue ;
2010-12-13 14:23:04 +01:00
AdjacencyList_t G2 ( n ) ;
// It is necessary to manually initialize vertex_index property since this graph uses listS and not vecS as underlying vertex container
2018-06-05 14:56:34 +02:00
auto v_index = get ( boost : : vertex_index , G2 ) ;
2010-12-13 14:23:04 +01:00
for ( int i = 0 ; i < n ; i + + )
put ( v_index , vertex ( i , G2 ) , i ) ;
2009-12-16 14:21:31 +01:00
vector < int > reverse_equation_reordered ( nb_var ) , reverse_variable_reordered ( nb_var ) ;
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < nb_var ; i + + )
2009-12-16 14:21:31 +01:00
{
reverse_equation_reordered [ equation_reordered [ i ] ] = i ;
reverse_variable_reordered [ variable_reordered [ i ] ] = i ;
2017-06-14 07:01:31 +02:00
}
jacob_map_t tmp_normalized_contemporaneous_jacobian ;
if ( cutoff = = 0 )
{
2018-06-04 14:17:36 +02:00
set < pair < int , int > > endo ;
2016-03-25 15:38:49 +01:00
for ( int i = 0 ; i < nb_var ; i + + )
2017-06-14 07:01:31 +02:00
{
2016-03-25 15:38:49 +01:00
endo . clear ( ) ;
2017-06-14 07:01:31 +02:00
equations [ i ] - > collectEndogenous ( endo ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : endo )
2018-06-04 16:36:46 +02:00
tmp_normalized_contemporaneous_jacobian [ { i , it . first } ] = 1 ;
2017-06-14 07:01:31 +02:00
}
}
else
tmp_normalized_contemporaneous_jacobian = static_jacobian ;
2019-12-16 19:42:59 +01:00
for ( const auto & [ key , value ] : tmp_normalized_contemporaneous_jacobian )
if ( reverse_equation_reordered [ key . first ] > = static_cast < int > ( prologue ) & & reverse_equation_reordered [ key . first ] < static_cast < int > ( nb_var - epilogue )
& & reverse_variable_reordered [ key . second ] > = static_cast < int > ( prologue ) & & reverse_variable_reordered [ key . second ] < static_cast < int > ( nb_var - epilogue )
& & key . first ! = endo2eq [ key . second ] )
add_edge ( vertex ( reverse_equation_reordered [ endo2eq [ key . second ] ] - prologue , G2 ) ,
vertex ( reverse_equation_reordered [ key . first ] - prologue , G2 ) ,
2010-12-13 14:23:04 +01:00
G2 ) ;
2009-12-16 14:21:31 +01:00
vector < int > endo2block ( num_vertices ( G2 ) ) , discover_time ( num_vertices ( G2 ) ) ;
2018-06-05 14:56:34 +02:00
boost : : iterator_property_map < int * , boost : : property_map < AdjacencyList_t , boost : : vertex_index_t > : : type , int , int & > endo2block_map ( & endo2block [ 0 ] , get ( boost : : vertex_index , G2 ) ) ;
2009-12-16 14:21:31 +01:00
// Compute strongly connected components
2010-12-13 14:23:04 +01:00
int num = strong_components ( G2 , endo2block_map ) ;
2009-12-16 14:21:31 +01:00
2018-06-04 16:36:46 +02:00
blocks = vector < pair < int , int > > ( num , { 0 , 0 } ) ;
2009-12-16 14:21:31 +01:00
// Create directed acyclic graph associated to the strongly connected components
2018-06-05 14:56:34 +02:00
using DirectedGraph = boost : : adjacency_list < boost : : vecS , boost : : vecS , boost : : directedS > ;
2009-12-16 14:21:31 +01:00
DirectedGraph dag ( num ) ;
2009-12-16 18:13:23 +01:00
for ( unsigned int i = 0 ; i < num_vertices ( G2 ) ; i + + )
2009-12-16 14:21:31 +01:00
{
2010-12-13 14:23:04 +01:00
AdjacencyList_t : : out_edge_iterator it_out , out_end ;
AdjacencyList_t : : vertex_descriptor vi = vertex ( i , G2 ) ;
2009-12-16 14:21:31 +01:00
for ( tie ( it_out , out_end ) = out_edges ( vi , G2 ) ; it_out ! = out_end ; + + it_out )
{
2010-12-13 14:23:04 +01:00
int t_b = endo2block_map [ target ( * it_out , G2 ) ] ;
int s_b = endo2block_map [ source ( * it_out , G2 ) ] ;
2009-12-16 14:21:31 +01:00
if ( s_b ! = t_b )
add_edge ( s_b , t_b , dag ) ;
}
}
// Compute topological sort of DAG (ordered list of unordered SCC)
deque < int > ordered2unordered ;
topological_sort ( dag , front_inserter ( ordered2unordered ) ) ; // We use a front inserter because topological_sort returns the inverse order
// Construct mapping from unordered SCC to ordered SCC
vector < int > unordered2ordered ( num ) ;
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < num ; i + + )
2009-12-16 14:21:31 +01:00
unordered2ordered [ ordered2unordered [ i ] ] = i ;
//This vector contains for each block:
// - first set = equations belonging to the block,
// - second set = the feeback variables,
// - third vector = the reordered non-feedback variables.
2018-11-23 17:19:59 +01:00
vector < tuple < set < int > , set < int > , vector < int > > > components_set ( num ) ;
2009-12-16 14:21:31 +01:00
for ( unsigned int i = 0 ; i < endo2block . size ( ) ; i + + )
{
endo2block [ i ] = unordered2ordered [ endo2block [ i ] ] ;
blocks [ endo2block [ i ] ] . first + + ;
2018-11-23 17:19:59 +01:00
get < 0 > ( components_set [ endo2block [ i ] ] ) . insert ( i ) ;
2009-12-16 14:21:31 +01:00
}
2010-09-16 17:51:50 +02:00
getVariableLeadLagByBlock ( dynamic_jacobian , endo2block , num , equation_lag_lead , variable_lag_lead , equation_reordered , variable_reordered ) ;
2009-12-16 14:21:31 +01:00
vector < int > tmp_equation_reordered ( equation_reordered ) , tmp_variable_reordered ( variable_reordered ) ;
int order = prologue ;
//Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
2009-12-16 18:13:23 +01:00
if ( select_feedback_variable )
2009-12-16 14:21:31 +01:00
{
for ( int i = 0 ; i < n ; i + + )
if ( Equation_Type [ equation_reordered [ i + prologue ] ] . first = = E_SOLVE
2010-04-28 16:03:32 +02:00
| | variable_lag_lead [ variable_reordered [ i + prologue ] ] . second > 0
| | variable_lag_lead [ variable_reordered [ i + prologue ] ] . first > 0
| | equation_lag_lead [ equation_reordered [ i + prologue ] ] . second > 0
| | equation_lag_lead [ equation_reordered [ i + prologue ] ] . first > 0
| | mfs = = 0 )
2010-12-13 14:23:04 +01:00
add_edge ( vertex ( i , G2 ) , vertex ( i , G2 ) , G2 ) ;
2009-12-16 14:21:31 +01:00
}
else
2019-12-16 19:42:59 +01:00
for ( int i = 0 ; i < n ; i + + )
if ( Equation_Type [ equation_reordered [ i + prologue ] ] . first = = E_SOLVE | | mfs = = 0 )
add_edge ( vertex ( i , G2 ) , vertex ( i , G2 ) , G2 ) ;
2010-07-23 11:20:24 +02:00
//Determines the dynamic structure of each equation
n_static = vector < unsigned int > ( prologue + num + epilogue , 0 ) ;
n_forward = vector < unsigned int > ( prologue + num + epilogue , 0 ) ;
n_backward = vector < unsigned int > ( prologue + num + epilogue , 0 ) ;
n_mixed = vector < unsigned int > ( prologue + num + epilogue , 0 ) ;
2019-04-23 11:07:32 +02:00
for ( int i = 0 ; i < static_cast < int > ( prologue ) ; i + + )
2019-12-16 19:42:59 +01:00
if ( variable_lag_lead [ tmp_variable_reordered [ i ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ i ] ] . second ! = 0 )
n_mixed [ i ] + + ;
else if ( variable_lag_lead [ tmp_variable_reordered [ i ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ i ] ] . second ! = 0 )
n_forward [ i ] + + ;
else if ( variable_lag_lead [ tmp_variable_reordered [ i ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ i ] ] . second = = 0 )
n_backward [ i ] + + ;
else if ( variable_lag_lead [ tmp_variable_reordered [ i ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ i ] ] . second = = 0 )
n_static [ i ] + + ;
2009-12-16 14:21:31 +01:00
//For each block, the minimum set of feedback variable is computed
// and the non-feedback variables are reordered to get
// a sub-recursive block without feedback variables
for ( int i = 0 ; i < num ; i + + )
{
2018-11-23 17:19:59 +01:00
AdjacencyList_t G = extract_subgraph ( G2 , get < 0 > ( components_set [ i ] ) ) ;
2009-12-16 14:21:31 +01:00
set < int > feed_back_vertices ;
2010-09-16 19:00:48 +02:00
AdjacencyList_t G1 = Minimal_set_of_feedback_vertex ( feed_back_vertices , G ) ;
2018-06-05 14:56:34 +02:00
auto v_index = get ( boost : : vertex_index , G ) ;
2018-11-23 17:19:59 +01:00
get < 1 > ( components_set [ i ] ) = feed_back_vertices ;
2009-12-16 14:21:31 +01:00
blocks [ i ] . second = feed_back_vertices . size ( ) ;
vector < int > Reordered_Vertice ;
Reorder_the_recursive_variables ( G , feed_back_vertices , Reordered_Vertice ) ;
//First we have the recursive equations conditional on feedback variables
2010-07-23 11:20:24 +02:00
for ( int j = 0 ; j < 4 ; j + + )
2019-12-16 19:42:59 +01:00
for ( int its : Reordered_Vertice )
{
bool something_done = false ;
if ( j = = 2 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . second ! = 0 )
{
n_mixed [ prologue + i ] + + ;
something_done = true ;
}
else if ( j = = 3 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . second ! = 0 )
{
n_forward [ prologue + i ] + + ;
something_done = true ;
}
else if ( j = = 1 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . second = = 0 )
{
n_backward [ prologue + i ] + + ;
something_done = true ;
}
else if ( j = = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . second = = 0 )
{
n_static [ prologue + i ] + + ;
something_done = true ;
}
if ( something_done )
{
equation_reordered [ order ] = tmp_equation_reordered [ its + prologue ] ;
variable_reordered [ order ] = tmp_variable_reordered [ its + prologue ] ;
order + + ;
}
}
2018-11-23 17:19:59 +01:00
get < 2 > ( components_set [ i ] ) = Reordered_Vertice ;
2009-12-16 14:21:31 +01:00
//Second we have the equations related to the feedback variables
2010-07-23 11:20:24 +02:00
for ( int j = 0 ; j < 4 ; j + + )
2019-12-16 19:42:59 +01:00
for ( int feed_back_vertice : feed_back_vertices )
{
bool something_done = false ;
if ( j = = 2 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ] . second ! = 0 )
{
n_mixed [ prologue + i ] + + ;
something_done = true ;
}
else if ( j = = 3 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ] . second ! = 0 )
{
n_forward [ prologue + i ] + + ;
something_done = true ;
}
else if ( j = = 1 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ] . second = = 0 )
{
n_backward [ prologue + i ] + + ;
something_done = true ;
}
else if ( j = = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ] . second = = 0 )
{
n_static [ prologue + i ] + + ;
something_done = true ;
}
if ( something_done )
{
equation_reordered [ order ] = tmp_equation_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ;
variable_reordered [ order ] = tmp_variable_reordered [ v_index [ vertex ( feed_back_vertice , G ) ] + prologue ] ;
order + + ;
}
}
2009-12-16 14:21:31 +01:00
}
2010-07-23 11:20:24 +02:00
2019-04-23 11:07:32 +02:00
for ( int i = 0 ; i < static_cast < int > ( epilogue ) ; i + + )
2019-12-16 19:42:59 +01:00
if ( variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . second ! = 0 )
n_mixed [ prologue + num + i ] + + ;
else if ( variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . second ! = 0 )
n_forward [ prologue + num + i ] + + ;
else if ( variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . second = = 0 )
n_backward [ prologue + num + i ] + + ;
else if ( variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . second = = 0 )
n_static [ prologue + num + i ] + + ;
2010-07-23 11:20:24 +02:00
2018-11-23 17:19:59 +01:00
inv_equation_reordered . resize ( nb_var ) ;
inv_variable_reordered . resize ( nb_var ) ;
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < nb_var ; i + + )
2009-12-16 14:21:31 +01:00
{
inv_variable_reordered [ variable_reordered [ i ] ] = i ;
inv_equation_reordered [ equation_reordered [ i ] ] = i ;
}
}
2009-12-16 18:13:23 +01:00
void
2018-06-04 14:17:36 +02:00
ModelTree : : printBlockDecomposition ( const vector < pair < int , int > > & blocks ) const
2009-12-16 14:21:31 +01:00
{
2019-12-16 19:42:59 +01:00
int largest_block = 0 ,
Nb_SimulBlocks = 0 ,
Nb_feedback_variable = 0 ;
2009-12-16 14:21:31 +01:00
unsigned int Nb_TotalBlocks = getNbBlocks ( ) ;
for ( unsigned int block = 0 ; block < Nb_TotalBlocks ; block + + )
{
BlockSimulationType simulation_type = getBlockSimulationType ( block ) ;
if ( simulation_type = = SOLVE_FORWARD_COMPLETE | | simulation_type = = SOLVE_BACKWARD_COMPLETE | | simulation_type = = SOLVE_TWO_BOUNDARIES_COMPLETE )
{
Nb_SimulBlocks + + ;
int size = getBlockSize ( block ) ;
if ( size > largest_block )
{
largest_block = size ;
2011-09-20 15:02:27 +02:00
Nb_feedback_variable = getBlockMfs ( block ) ;
2009-12-16 14:21:31 +01:00
}
}
}
int Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks ;
cout < < Nb_TotalBlocks < < " block(s) found: " < < endl
< < " " < < Nb_RecursBlocks < < " recursive block(s) and " < < Nb_SimulBlocks < < " simultaneous block(s). " < < endl
< < " the largest simultaneous block has " < < largest_block < < " equation(s) " < < endl
< < " and " < < Nb_feedback_variable < < " feedback variable(s). " < < endl ;
}
2010-09-16 19:00:48 +02:00
block_type_firstequation_size_mfs_t
2018-11-23 17:19:59 +01:00
ModelTree : : reduceBlocksAndTypeDetermination ( const dynamic_jacob_map_t & dynamic_jacobian , vector < pair < int , int > > & blocks , const equation_type_and_normalized_equation_t & Equation_Type , const vector < int > & variable_reordered , const vector < int > & equation_reordered , vector < unsigned int > & n_static , vector < unsigned int > & n_forward , vector < unsigned int > & n_backward , vector < unsigned int > & n_mixed , vector < tuple < int , int , int , int > > & block_col_type , bool linear_decomposition )
2009-12-16 14:21:31 +01:00
{
int i = 0 ;
int count_equ = 0 , blck_count_simult = 0 ;
int Blck_Size , MFS_Size ;
int Lead , Lag ;
2010-09-16 19:00:48 +02:00
block_type_firstequation_size_mfs_t block_type_size_mfs ;
2009-12-16 14:21:31 +01:00
BlockSimulationType Simulation_Type , prev_Type = UNKNOWN ;
int eq = 0 ;
2019-12-16 19:42:59 +01:00
unsigned int l_n_static = 0 , l_n_forward = 0 , l_n_backward = 0 , l_n_mixed = 0 ;
2019-04-23 11:07:32 +02:00
for ( i = 0 ; i < static_cast < int > ( prologue + blocks . size ( ) + epilogue ) ; i + + )
2009-12-16 14:21:31 +01:00
{
int first_count_equ = count_equ ;
2019-04-23 11:07:32 +02:00
if ( i < static_cast < int > ( prologue ) )
2009-12-16 14:21:31 +01:00
{
Blck_Size = 1 ;
MFS_Size = 1 ;
}
2019-04-23 11:07:32 +02:00
else if ( i < static_cast < int > ( prologue + blocks . size ( ) ) )
2009-12-16 14:21:31 +01:00
{
Blck_Size = blocks [ blck_count_simult ] . first ;
MFS_Size = blocks [ blck_count_simult ] . second ;
blck_count_simult + + ;
}
2019-04-23 11:07:32 +02:00
else if ( i < static_cast < int > ( prologue + blocks . size ( ) + epilogue ) )
2009-12-16 14:21:31 +01:00
{
Blck_Size = 1 ;
MFS_Size = 1 ;
}
Lag = Lead = 0 ;
2018-06-04 14:17:36 +02:00
set < pair < int , int > > endo ;
2019-12-16 19:42:59 +01:00
for ( count_equ = first_count_equ ; count_equ < Blck_Size + first_count_equ ; count_equ + + )
2009-12-16 14:21:31 +01:00
{
endo . clear ( ) ;
equations [ equation_reordered [ count_equ ] ] - > collectEndogenous ( endo ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : endo )
2009-12-16 14:21:31 +01:00
{
2018-06-04 12:26:16 +02:00
int curr_variable = it . first ;
int curr_lag = it . second ;
2018-09-21 17:13:19 +02:00
if ( linear_decomposition )
{
2018-11-23 17:19:59 +01:00
if ( dynamic_jacobian . find ( { curr_lag , equation_reordered [ count_equ ] , curr_variable } ) ! = dynamic_jacobian . end ( ) )
2018-09-21 17:13:19 +02:00
{
if ( curr_lag > Lead )
Lead = curr_lag ;
else if ( - curr_lag > Lag )
Lag = - curr_lag ;
}
}
else
{
2019-12-16 19:42:59 +01:00
if ( find ( variable_reordered . begin ( ) + first_count_equ , variable_reordered . begin ( ) + ( first_count_equ + Blck_Size ) , curr_variable )
! = variable_reordered . begin ( ) + ( first_count_equ + Blck_Size )
& & dynamic_jacobian . find ( { curr_lag , equation_reordered [ count_equ ] , curr_variable } ) ! = dynamic_jacobian . end ( ) )
{
if ( curr_lag > Lead )
Lead = curr_lag ;
else if ( - curr_lag > Lag )
Lag = - curr_lag ;
}
2018-09-21 17:13:19 +02:00
}
2009-12-16 14:21:31 +01:00
}
}
2019-12-16 19:42:59 +01:00
if ( Lag > 0 & & Lead > 0 )
2009-12-16 14:21:31 +01:00
{
if ( Blck_Size = = 1 )
Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE ;
else
Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE ;
}
else if ( Blck_Size > 1 )
{
if ( Lead > 0 )
Simulation_Type = SOLVE_BACKWARD_COMPLETE ;
else
Simulation_Type = SOLVE_FORWARD_COMPLETE ;
}
else
{
if ( Lead > 0 )
Simulation_Type = SOLVE_BACKWARD_SIMPLE ;
else
Simulation_Type = SOLVE_FORWARD_SIMPLE ;
}
2010-07-23 11:20:24 +02:00
l_n_static = n_static [ i ] ;
l_n_forward = n_forward [ i ] ;
l_n_backward = n_backward [ i ] ;
l_n_mixed = n_mixed [ i ] ;
2009-12-16 14:21:31 +01:00
if ( Blck_Size = = 1 )
{
2010-04-28 16:03:32 +02:00
if ( Equation_Type [ equation_reordered [ eq ] ] . first = = E_EVALUATE | | Equation_Type [ equation_reordered [ eq ] ] . first = = E_EVALUATE_S )
2009-12-16 14:21:31 +01:00
{
if ( Simulation_Type = = SOLVE_BACKWARD_SIMPLE )
Simulation_Type = EVALUATE_BACKWARD ;
else if ( Simulation_Type = = SOLVE_FORWARD_SIMPLE )
Simulation_Type = EVALUATE_FORWARD ;
}
if ( i > 0 )
{
2011-10-12 14:58:29 +02:00
bool is_lead = false , is_lag = false ;
2018-11-23 17:19:59 +01:00
int c_Size = get < 2 > ( block_type_size_mfs [ block_type_size_mfs . size ( ) - 1 ] ) ;
int first_equation = get < 1 > ( block_type_size_mfs [ block_type_size_mfs . size ( ) - 1 ] ) ;
2019-12-20 16:59:30 +01:00
if ( c_Size > 0 & & ( ( prev_Type = = EVALUATE_FORWARD & & Simulation_Type = = EVALUATE_FORWARD & & ! is_lead )
| | ( prev_Type = = EVALUATE_BACKWARD & & Simulation_Type = = EVALUATE_BACKWARD & & ! is_lag ) ) )
2011-10-12 14:58:29 +02:00
{
for ( int j = first_equation ; j < first_equation + c_Size ; j + + )
{
2018-11-23 17:19:59 +01:00
auto it = dynamic_jacobian . find ( { - 1 , equation_reordered [ eq ] , variable_reordered [ j ] } ) ;
2011-10-12 14:58:29 +02:00
if ( it ! = dynamic_jacobian . end ( ) )
is_lag = true ;
2018-11-23 17:19:59 +01:00
it = dynamic_jacobian . find ( { + 1 , equation_reordered [ eq ] , variable_reordered [ j ] } ) ;
2011-10-12 14:58:29 +02:00
if ( it ! = dynamic_jacobian . end ( ) )
is_lead = true ;
}
}
2019-12-20 16:59:30 +01:00
if ( ( prev_Type = = EVALUATE_FORWARD & & Simulation_Type = = EVALUATE_FORWARD & & ! is_lead )
| | ( prev_Type = = EVALUATE_BACKWARD & & Simulation_Type = = EVALUATE_BACKWARD & & ! is_lag ) )
2009-12-16 14:21:31 +01:00
{
//merge the current block with the previous one
2018-11-23 17:19:59 +01:00
BlockSimulationType c_Type = get < 0 > ( block_type_size_mfs [ block_type_size_mfs . size ( ) - 1 ] ) ;
2010-07-23 11:20:24 +02:00
c_Size + + ;
2018-11-23 17:19:59 +01:00
block_type_size_mfs [ block_type_size_mfs . size ( ) - 1 ] = { c_Type , first_equation , c_Size , c_Size } ;
2009-12-16 18:13:23 +01:00
if ( block_lag_lead [ block_type_size_mfs . size ( ) - 1 ] . first > Lag )
2009-12-16 14:21:31 +01:00
Lag = block_lag_lead [ block_type_size_mfs . size ( ) - 1 ] . first ;
2009-12-16 18:13:23 +01:00
if ( block_lag_lead [ block_type_size_mfs . size ( ) - 1 ] . second > Lead )
2009-12-16 14:21:31 +01:00
Lead = block_lag_lead [ block_type_size_mfs . size ( ) - 1 ] . second ;
2018-06-04 16:36:46 +02:00
block_lag_lead [ block_type_size_mfs . size ( ) - 1 ] = { Lag , Lead } ;
2018-11-23 17:19:59 +01:00
auto tmp = block_col_type [ block_col_type . size ( ) - 1 ] ;
block_col_type [ block_col_type . size ( ) - 1 ] = { get < 0 > ( tmp ) + l_n_static , get < 1 > ( tmp ) + l_n_forward , get < 2 > ( tmp ) + l_n_backward , get < 3 > ( tmp ) + l_n_mixed } ;
2009-12-16 14:21:31 +01:00
}
else
{
2018-11-23 17:19:59 +01:00
block_type_size_mfs . emplace_back ( Simulation_Type , eq , Blck_Size , MFS_Size ) ;
2018-06-04 16:36:46 +02:00
block_lag_lead . emplace_back ( Lag , Lead ) ;
2018-11-23 17:19:59 +01:00
block_col_type . emplace_back ( l_n_static , l_n_forward , l_n_backward , l_n_mixed ) ;
2009-12-16 14:21:31 +01:00
}
}
else
{
2018-11-23 17:19:59 +01:00
block_type_size_mfs . emplace_back ( Simulation_Type , eq , Blck_Size , MFS_Size ) ;
2018-06-04 16:36:46 +02:00
block_lag_lead . emplace_back ( Lag , Lead ) ;
2018-11-23 17:19:59 +01:00
block_col_type . emplace_back ( l_n_static , l_n_forward , l_n_backward , l_n_mixed ) ;
2009-12-16 14:21:31 +01:00
}
}
else
{
2018-11-23 17:19:59 +01:00
block_type_size_mfs . emplace_back ( Simulation_Type , eq , Blck_Size , MFS_Size ) ;
2018-06-04 16:36:46 +02:00
block_lag_lead . emplace_back ( Lag , Lead ) ;
2018-11-23 17:19:59 +01:00
block_col_type . emplace_back ( l_n_static , l_n_forward , l_n_backward , l_n_mixed ) ;
2009-12-16 14:21:31 +01:00
}
prev_Type = Simulation_Type ;
eq + = Blck_Size ;
}
2018-11-23 17:19:59 +01:00
return block_type_size_mfs ;
2009-12-16 14:21:31 +01:00
}
2018-09-21 17:13:19 +02:00
vector < bool >
2018-11-23 17:19:59 +01:00
ModelTree : : equationLinear ( map < tuple < int , int , int > , expr_t > first_order_endo_derivatives ) const
2018-09-21 17:13:19 +02:00
{
vector < bool > is_linear ( symbol_table . endo_nbr ( ) , true ) ;
2018-11-23 17:19:59 +01:00
for ( const auto & it : first_order_endo_derivatives )
2018-09-21 17:13:19 +02:00
{
2019-12-20 16:59:30 +01:00
expr_t Id = it . second ;
set < pair < int , int > > endogenous ;
Id - > collectEndogenous ( endogenous ) ;
if ( endogenous . size ( ) > 0 )
{
int eq = get < 0 > ( it . first ) ;
is_linear [ eq ] = false ;
}
2018-09-21 17:13:19 +02:00
}
return is_linear ;
}
2009-12-16 14:21:31 +01:00
vector < bool >
2010-09-16 19:00:48 +02:00
ModelTree : : BlockLinear ( const blocks_derivatives_t & blocks_derivatives , const vector < int > & variable_reordered ) const
2009-12-16 14:21:31 +01:00
{
unsigned int nb_blocks = getNbBlocks ( ) ;
vector < bool > blocks_linear ( nb_blocks , true ) ;
2009-12-16 18:13:23 +01:00
for ( unsigned int block = 0 ; block < nb_blocks ; block + + )
2009-12-16 14:21:31 +01:00
{
BlockSimulationType simulation_type = getBlockSimulationType ( block ) ;
int block_size = getBlockSize ( block ) ;
2010-09-16 19:00:48 +02:00
block_derivatives_equation_variable_laglead_nodeid_t derivatives_block = blocks_derivatives [ block ] ;
2009-12-16 14:21:31 +01:00
int first_variable_position = getBlockFirstEquation ( block ) ;
2009-12-16 18:13:23 +01:00
if ( simulation_type = = SOLVE_BACKWARD_COMPLETE | | simulation_type = = SOLVE_FORWARD_COMPLETE )
2019-12-16 19:42:59 +01:00
for ( const auto & [ ignore , ignore2 , lag , d1 ] : derivatives_block )
{
if ( lag = = 0 )
{
set < pair < int , int > > endogenous ;
d1 - > collectEndogenous ( endogenous ) ;
if ( endogenous . size ( ) > 0 )
2009-12-16 18:13:23 +01:00
for ( int l = 0 ; l < block_size ; l + + )
2019-12-16 19:42:59 +01:00
if ( endogenous . find ( { variable_reordered [ first_variable_position + l ] , 0 } ) ! = endogenous . end ( ) )
{
blocks_linear [ block ] = false ;
goto the_end ;
}
}
}
else if ( simulation_type = = SOLVE_TWO_BOUNDARIES_COMPLETE | | simulation_type = = SOLVE_TWO_BOUNDARIES_SIMPLE )
for ( const auto & [ ignore , ignore2 , lag , d1 ] : derivatives_block )
{
set < pair < int , int > > endogenous ;
d1 - > collectEndogenous ( endogenous ) ;
if ( endogenous . size ( ) > 0 )
for ( int l = 0 ; l < block_size ; l + + )
if ( endogenous . find ( { variable_reordered [ first_variable_position + l ] , lag } ) ! = endogenous . end ( ) )
{
blocks_linear [ block ] = false ;
goto the_end ;
}
}
2009-12-16 18:13:23 +01:00
the_end :
;
2009-12-16 14:21:31 +01:00
}
2018-11-23 17:19:59 +01:00
return blocks_linear ;
2009-12-16 14:21:31 +01:00
}
2008-02-03 11:28:36 +01:00
int
ModelTree : : equation_number ( ) const
2009-01-23 11:59:37 +01:00
{
2009-12-16 18:13:23 +01:00
return ( equations . size ( ) ) ;
2009-01-23 11:59:37 +01:00
}
2008-02-03 11:28:36 +01:00
void
ModelTree : : writeDerivative ( ostream & output , int eq , int symb_id , int lag ,
ExprNodeOutputType output_type ,
2010-09-16 19:00:48 +02:00
const temporary_terms_t & temporary_terms ) const
2009-01-23 11:59:37 +01:00
{
2019-12-16 19:42:59 +01:00
if ( auto it = derivatives [ 1 ] . find ( { eq , getDerivID ( symb_id , lag ) } ) ;
it ! = derivatives [ 1 ] . end ( ) )
it - > second - > writeOutput ( output , output_type , temporary_terms , { } ) ;
2009-01-23 11:59:37 +01:00
else
output < < 0 ;
}
2008-02-03 11:28:36 +01:00
void
2018-11-22 14:32:40 +01:00
ModelTree : : computeDerivatives ( int order , const set < int > & vars )
2008-02-03 11:28:36 +01:00
{
2019-12-20 16:59:30 +01:00
assert ( order > = 1 ) ;
2008-02-03 11:28:36 +01:00
2018-11-22 14:32:40 +01:00
// Do not shrink the vectors, since they have a minimal size of 4 (see constructor)
2019-04-12 15:41:52 +02:00
derivatives . resize ( max ( static_cast < size_t > ( order + 1 ) , derivatives . size ( ) ) ) ;
NNZDerivatives . resize ( max ( static_cast < size_t > ( order + 1 ) , NNZDerivatives . size ( ) ) , 0 ) ;
2009-04-20 12:48:54 +02:00
2018-11-22 14:32:40 +01:00
// First-order derivatives
for ( int var : vars )
2019-04-12 15:41:52 +02:00
for ( int eq = 0 ; eq < static_cast < int > ( equations . size ( ) ) ; eq + + )
2018-11-22 14:32:40 +01:00
{
expr_t d1 = equations [ eq ] - > getDerivative ( var ) ;
if ( d1 = = Zero )
continue ;
derivatives [ 1 ] [ { eq , var } ] = d1 ;
+ + NNZDerivatives [ 1 ] ;
}
2009-04-20 12:48:54 +02:00
2018-11-22 14:32:40 +01:00
// Higher-order derivatives
for ( int o = 2 ; o < = order ; o + + )
for ( const auto & it : derivatives [ o - 1 ] )
for ( int var : vars )
2008-02-03 11:28:36 +01:00
{
2018-11-22 14:32:40 +01:00
if ( it . first . back ( ) > var )
2009-04-20 12:48:54 +02:00
continue ;
2018-11-22 14:32:40 +01:00
expr_t d = it . second - > getDerivative ( var ) ;
if ( d = = Zero )
2009-04-20 12:48:54 +02:00
continue ;
2018-11-22 14:32:40 +01:00
vector < int > indices { it . first } ;
indices . push_back ( var ) ;
// At this point, indices of endogenous variables are sorted in non-decreasing order
derivatives [ o ] [ indices ] = d ;
2019-06-17 15:28:33 +02:00
// We output symmetric elements at order = 2
if ( o = = 2 & & indices [ 1 ] ! = indices [ 2 ] )
NNZDerivatives [ o ] + = 2 ;
2019-04-12 15:41:52 +02:00
else
2018-11-22 14:32:40 +01:00
NNZDerivatives [ o ] + + ;
2008-02-03 11:28:36 +01:00
}
}
void
2018-09-25 15:56:52 +02:00
ModelTree : : computeTemporaryTerms ( bool is_matlab , bool no_tmp_terms )
2008-02-03 11:28:36 +01:00
{
2018-11-16 16:53:07 +01:00
/* Collect all model local variables appearing in equations (and only those,
because printing unused model local variables can lead to a crash ,
see Dynare / dynare # 101 ) .
Then store them in a dedicated structure ( temporary_terms_mlv ) , that will
be treated as the rest of temporary terms . */
2018-11-30 12:22:13 +01:00
temporary_terms_mlv . clear ( ) ;
2018-05-23 19:19:11 +02:00
set < int > used_local_vars ;
2019-12-20 16:59:30 +01:00
for ( auto & equation : equations )
2018-07-17 18:34:07 +02:00
equation - > collectVariables ( SymbolType : : modelLocalVariable , used_local_vars ) ;
2018-06-04 12:26:16 +02:00
for ( int used_local_var : used_local_vars )
2018-03-27 17:14:30 +02:00
{
2018-06-04 12:26:16 +02:00
VariableNode * v = AddVariable ( used_local_var ) ;
temporary_terms_mlv [ v ] = local_variables_table . find ( used_local_var ) - > second ;
2018-03-27 17:14:30 +02:00
}
2018-11-30 12:22:13 +01:00
// Compute the temporary terms in equations and derivatives
map < pair < int , int > , temporary_terms_t > temp_terms_map ;
2019-08-14 15:28:41 +02:00
map < expr_t , pair < int , pair < int , int > > > reference_count ;
2019-12-20 16:59:30 +01:00
for ( auto & equation : equations )
2019-08-14 15:28:41 +02:00
equation - > computeTemporaryTerms ( { 0 , 0 } ,
temp_terms_map ,
reference_count ,
is_matlab ) ;
for ( int order = 1 ; order < static_cast < int > ( derivatives . size ( ) ) ; order + + )
for ( const auto & it : derivatives [ order ] )
it . second - > computeTemporaryTerms ( { 0 , order } ,
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 )
// The following loop can be simplified with std::erase_if() in C++20
2019-12-20 16:59:30 +01:00
for ( auto it2 = it . second . begin ( ) ; it2 ! = it . second . end ( ) ; )
2019-08-14 15:28:41 +02:00
if ( ! dynamic_cast < AbstractExternalFunctionNode * > ( * it2 ) )
it2 = it . second . erase ( it2 ) ;
else
+ + it2 ;
2015-09-03 13:50:02 +02:00
2018-11-30 12:22:13 +01:00
// Fill the (now obsolete) temporary_terms structure
temporary_terms . clear ( ) ;
for ( const auto & it : temp_terms_map )
temporary_terms . insert ( it . second . begin ( ) , it . second . end ( ) ) ;
2015-09-03 13:50:02 +02:00
2018-11-30 12:22:13 +01:00
// Fill the new structure
temporary_terms_derivatives . clear ( ) ;
temporary_terms_derivatives . resize ( derivatives . size ( ) ) ;
2019-04-23 11:07:32 +02:00
for ( int order = 0 ; order < static_cast < int > ( derivatives . size ( ) ) ; order + + )
2018-11-30 12:22:13 +01:00
temporary_terms_derivatives [ order ] = move ( temp_terms_map [ { 0 , order } ] ) ;
2018-03-27 17:14:30 +02:00
2018-11-30 12:22:13 +01:00
// Compute indices in MATLAB/Julia vector
2018-05-25 15:19:33 +02:00
int idx = 0 ;
2018-11-16 16:53:07 +01:00
for ( auto & it : temporary_terms_mlv )
temporary_terms_idxs [ it . first ] = idx + + ;
2019-04-23 11:07:32 +02:00
for ( int order = 0 ; order < static_cast < int > ( derivatives . size ( ) ) ; order + + )
2018-11-30 12:22:13 +01:00
for ( const auto & it : temporary_terms_derivatives [ order ] )
temporary_terms_idxs [ it ] = idx + + ;
2018-03-27 17:14:30 +02:00
}
void
2018-11-16 16:53:07 +01:00
ModelTree : : writeModelLocalVariableTemporaryTerms ( temporary_terms_t & temp_term_union ,
2018-11-16 18:13:34 +01:00
const temporary_terms_idxs_t & tt_idxs ,
2018-03-27 17:14:30 +02:00
ostream & output , ExprNodeOutputType output_type ,
deriv_node_temp_terms_t & tef_terms ) const
{
2018-11-16 16:53:07 +01:00
temporary_terms_t tto ;
for ( auto it : temporary_terms_mlv )
tto . insert ( it . first ) ;
for ( auto & it : temporary_terms_mlv )
2018-03-27 17:14:30 +02:00
{
2018-09-25 19:15:22 +02:00
if ( isJuliaOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " @inbounds const " ;
2018-11-16 18:13:34 +01:00
it . first - > writeOutput ( output , output_type , tto , tt_idxs , tef_terms ) ;
2018-03-27 17:14:30 +02:00
output < < " = " ;
2018-11-16 18:13:34 +01:00
it . second - > writeOutput ( output , output_type , temp_term_union , tt_idxs , tef_terms ) ;
2018-03-27 17:14:30 +02:00
2018-09-05 18:27:13 +02:00
if ( isCOutput ( output_type ) | | isMatlabOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " ; " ;
output < < endl ;
2018-11-16 16:53:07 +01:00
/* We put in temp_term_union the VariableNode corresponding to the MLV,
not its definition , so that when equations use the MLV ,
T ( XXX ) is printed instead of the MLV name */
temp_term_union . insert ( it . first ) ;
2018-03-27 17:14:30 +02:00
}
}
void
2018-05-24 19:29:53 +02:00
ModelTree : : writeTemporaryTerms ( const temporary_terms_t & tt ,
2018-11-16 16:53:07 +01:00
temporary_terms_t & temp_term_union ,
2018-05-28 15:23:15 +02:00
const temporary_terms_idxs_t & tt_idxs ,
2018-03-27 17:14:30 +02:00
ostream & output , ExprNodeOutputType output_type , deriv_node_temp_terms_t & tef_terms ) const
{
2018-11-16 16:53:07 +01:00
for ( auto it : tt )
2018-03-27 17:14:30 +02:00
{
2019-12-16 19:42:59 +01:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( it ) )
2018-11-16 16:53:07 +01:00
it - > writeExternalFunctionOutput ( output , output_type , temp_term_union , tt_idxs , tef_terms ) ;
2018-03-27 17:14:30 +02:00
2018-09-25 19:15:22 +02:00
if ( isJuliaOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " @inbounds " ;
2018-11-16 16:53:07 +01:00
it - > writeOutput ( output , output_type , tt , tt_idxs , tef_terms ) ;
2018-03-27 17:14:30 +02:00
output < < " = " ;
2018-11-16 16:53:07 +01:00
it - > writeOutput ( output , output_type , temp_term_union , tt_idxs , tef_terms ) ;
2018-03-27 17:14:30 +02:00
2018-09-05 18:27:13 +02:00
if ( isCOutput ( output_type ) | | isMatlabOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " ; " ;
output < < endl ;
2018-11-16 16:53:07 +01:00
temp_term_union . insert ( it ) ;
2018-03-27 17:14:30 +02:00
}
2009-01-23 11:59:37 +01:00
}
2008-02-03 11:28:36 +01:00
2017-02-20 12:18:11 +01:00
void
2019-04-18 14:34:48 +02:00
ModelTree : : writeJsonTemporaryTerms ( const temporary_terms_t & tt ,
temporary_terms_t & temp_term_union ,
ostream & output ,
deriv_node_temp_terms_t & tef_terms , const string & concat ) const
2017-02-20 12:18:11 +01:00
{
// Local var used to keep track of temp nodes already written
bool wrote_term = false ;
2019-04-18 14:34:48 +02:00
temporary_terms_t tt2 = temp_term_union ;
2017-02-20 12:18:11 +01:00
2019-04-03 16:32:52 +02:00
output < < R " ( " external_functions_temporary_terms_ ) " << concat << R " ( " : [) " ;
2018-06-04 12:26:16 +02:00
for ( auto it : tt )
2019-04-18 14:34:48 +02:00
{
2019-12-16 19:42:59 +01:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( it ) )
2019-04-18 14:34:48 +02:00
{
if ( wrote_term )
output < < " , " ;
vector < string > efout ;
it - > writeJsonExternalFunctionOutput ( efout , tt2 , tef_terms ) ;
for ( auto it1 = efout . begin ( ) ; it1 ! = efout . end ( ) ; + + it1 )
{
if ( it1 ! = efout . begin ( ) )
output < < " , " ;
output < < * it1 ;
}
wrote_term = true ;
}
tt2 . insert ( it ) ;
}
2017-02-20 12:18:11 +01:00
wrote_term = false ;
output < < " ] "
2019-04-03 16:32:52 +02:00
< < R " (, " temporary_terms_ ) " << concat << R " ( " : [) " ;
2019-04-18 14:34:48 +02:00
for ( const auto & it : tt )
{
if ( wrote_term )
output < < " , " ;
output < < R " ({ " temporary_term " : " ) " ;
it - > writeJsonOutput ( output , tt , tef_terms ) ;
output < < R " ( " ) "
< < R " (, " value " : " ) " ;
it - > writeJsonOutput ( output , temp_term_union , tef_terms ) ;
output < < R " ( " } ) " << endl;
wrote_term = true ;
temp_term_union . insert ( it ) ;
}
2017-02-20 12:18:11 +01:00
output < < " ] " ;
}
2016-12-28 14:02:50 +01:00
void
2017-01-05 15:19:13 +01:00
ModelTree : : fixNestedParenthesis ( ostringstream & output , map < string , string > & tmp_paren_vars , bool & message_printed ) const
2016-12-28 14:02:50 +01:00
{
string str = output . str ( ) ;
2016-12-30 18:32:20 +01:00
if ( ! testNestedParenthesis ( str ) )
return ;
2016-12-28 14:02:50 +01:00
int open = 0 ;
2016-12-30 11:26:56 +01:00
int first_open_paren = 0 ;
int matching_paren = 0 ;
bool hit_limit = false ;
int i1 = 0 ;
2017-01-02 11:37:15 +01:00
for ( size_t i = 0 ; i < str . length ( ) ; i + + )
2016-12-28 14:02:50 +01:00
{
2016-12-30 11:26:56 +01:00
if ( str . at ( i ) = = ' ( ' )
{
if ( open = = 0 )
first_open_paren = i ;
open + + ;
}
else if ( str . at ( i ) = = ' ) ' )
{
open - - ;
if ( open = = 0 )
matching_paren = i ;
}
2016-12-28 14:02:50 +01:00
if ( open > 32 )
2016-12-30 11:26:56 +01:00
hit_limit = true ;
if ( hit_limit & & open = = 0 )
2016-12-28 14:02:50 +01:00
{
2017-01-05 15:19:13 +01:00
if ( ! message_printed )
{
2019-10-02 11:24:37 +02:00
cerr < < " Warning: A .m file created by Dynare will have more than 32 nested parenthesis. MATLAB cannot support this. " < < endl
2017-01-05 15:19:13 +01:00
< < " We are going to modify, albeit inefficiently, this output to have fewer than 32 nested parenthesis. " < < endl
< < " It would hence behoove you to use the use_dll option of the model block to circumnavigate this problem. " < < endl
2019-10-02 11:24:37 +02:00
< < " If you have not yet set up a compiler on your system, see the MATLAB documentation for doing so. " < < endl
2017-01-05 15:19:13 +01:00
< < " For Windows, see: https://www.mathworks.com/help/matlab/matlab_external/install-mingw-support-package.html " < < endl < < endl ;
message_printed = true ;
}
2016-12-30 11:26:56 +01:00
string str1 = str . substr ( first_open_paren , matching_paren - first_open_paren + 1 ) ;
2019-12-16 19:42:59 +01:00
string repstr , varname ;
2016-12-30 11:26:56 +01:00
while ( testNestedParenthesis ( str1 ) )
{
2019-12-20 16:59:30 +01:00
size_t open_paren_idx = string : : npos ;
2017-01-02 11:37:15 +01:00
size_t match_paren_idx = string : : npos ;
size_t last_open_paren = string : : npos ;
for ( size_t j = 0 ; j < str1 . length ( ) ; j + + )
2016-12-30 11:26:56 +01:00
{
if ( str1 . at ( j ) = = ' ( ' )
{
// don't match, e.g. y(1)
2019-12-16 19:42:59 +01:00
if ( size_t idx = str1 . find_last_of ( " */-+ " , j - 1 ) ;
j = = 0 | | ( idx ! = string : : npos & & idx = = j - 1 ) )
2016-12-30 11:26:56 +01:00
open_paren_idx = j ;
last_open_paren = j ;
}
else if ( str1 . at ( j ) = = ' ) ' )
{
// don't match, e.g. y(1)
2019-12-16 19:42:59 +01:00
if ( size_t idx = str1 . find_last_not_of ( " 0123456789 " , j - 1 ) ;
idx ! = string : : npos & & idx ! = last_open_paren )
2016-12-30 11:26:56 +01:00
match_paren_idx = j ;
}
2017-01-02 11:37:15 +01:00
if ( open_paren_idx ! = string : : npos & & match_paren_idx ! = string : : npos )
2016-12-30 11:26:56 +01:00
{
string val = str1 . substr ( open_paren_idx , match_paren_idx - open_paren_idx + 1 ) ;
2019-12-16 19:42:59 +01:00
if ( auto it = tmp_paren_vars . find ( val ) ;
it = = tmp_paren_vars . end ( ) )
2016-12-30 11:26:56 +01:00
{
2017-01-02 11:21:28 +01:00
ostringstream ptvstr ;
ptvstr < < i1 + + ;
varname = " paren32_tmp_var_ " + ptvstr . str ( ) ;
2016-12-30 11:26:56 +01:00
repstr = repstr + varname + " = " + val + " ; \n " ;
tmp_paren_vars [ val ] = varname ;
}
else
varname = it - > second ;
str1 . replace ( open_paren_idx , match_paren_idx - open_paren_idx + 1 , varname ) ;
break ;
}
}
}
2019-12-16 19:42:59 +01:00
if ( auto it = tmp_paren_vars . find ( str1 ) ;
it = = tmp_paren_vars . end ( ) )
2016-12-30 11:26:56 +01:00
{
2017-01-02 11:21:28 +01:00
ostringstream ptvstr ;
ptvstr < < i1 + + ;
varname = " paren32_tmp_var_ " + ptvstr . str ( ) ;
2016-12-30 11:26:56 +01:00
repstr = repstr + varname + " = " + str1 + " ; \n " ;
}
else
varname = it - > second ;
str . replace ( first_open_paren , matching_paren - first_open_paren + 1 , varname ) ;
size_t insertLoc = str . find_last_of ( " \n " , first_open_paren ) ;
str . insert ( insertLoc + 1 , repstr ) ;
hit_limit = false ;
i = - 1 ;
first_open_paren = matching_paren = open = 0 ;
2016-12-28 14:02:50 +01:00
}
}
2016-12-30 11:26:56 +01:00
output . str ( str ) ;
}
bool
ModelTree : : testNestedParenthesis ( const string & str ) const
{
int open = 0 ;
2018-06-04 12:26:16 +02:00
for ( char i : str )
2016-12-30 11:26:56 +01:00
{
2018-06-04 12:26:16 +02:00
if ( i = = ' ( ' )
2016-12-30 11:26:56 +01:00
open + + ;
2018-06-04 12:26:16 +02:00
else if ( i = = ' ) ' )
2016-12-30 11:26:56 +01:00
open - - ;
if ( open > 32 )
return true ;
}
return false ;
2016-12-28 14:02:50 +01:00
}
2010-01-22 11:03:29 +01:00
void
2010-07-23 11:20:24 +02:00
ModelTree : : compileTemporaryTerms ( ostream & code_file , unsigned int & instruction_number , const temporary_terms_t & tt , map_idx_t map_idx , bool dynamic , bool steady_dynamic ) const
2010-01-22 11:03:29 +01:00
{
// Local var used to keep track of temp nodes already written
2010-09-16 19:00:48 +02:00
temporary_terms_t tt2 ;
2010-12-10 11:50:27 +01:00
// To store the functions that have already been written in the form TEF* = ext_fun();
deriv_node_temp_terms_t tef_terms ;
2018-06-04 12:26:16 +02:00
for ( auto it : tt )
2010-01-22 11:03:29 +01:00
{
2019-12-16 19:42:59 +01:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( it ) )
2010-12-10 11:50:27 +01:00
{
2018-06-04 12:26:16 +02:00
it - > compileExternalFunctionOutput ( code_file , instruction_number , false , tt2 , map_idx , dynamic , steady_dynamic , tef_terms ) ;
2010-12-10 11:50:27 +01:00
}
2019-04-23 11:07:32 +02:00
FNUMEXPR_ fnumexpr ( TemporaryTerm , static_cast < int > ( map_idx . find ( it - > idx ) - > second ) ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2018-06-04 12:26:16 +02:00
it - > compile ( code_file , instruction_number , false , tt2 , map_idx , dynamic , steady_dynamic , tef_terms ) ;
2010-01-22 11:03:29 +01:00
if ( dynamic )
{
2019-04-23 11:07:32 +02:00
FSTPT_ fstpt ( static_cast < int > ( map_idx . find ( it - > idx ) - > second ) ) ;
2010-07-23 11:20:24 +02:00
fstpt . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
}
else
{
2019-04-23 11:07:32 +02:00
FSTPST_ fstpst ( static_cast < int > ( map_idx . find ( it - > idx ) - > second ) ) ;
2010-07-23 11:20:24 +02:00
fstpst . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
}
// Insert current node into tt2
2018-06-04 12:26:16 +02:00
tt2 . insert ( it ) ;
2010-01-22 11:03:29 +01:00
}
}
2017-02-20 12:18:11 +01:00
void
ModelTree : : writeJsonModelLocalVariables ( ostream & output , deriv_node_temp_terms_t & tef_terms ) const
{
/* Collect all model local variables appearing in equations, and print only
them . Printing unused model local variables can lead to a crash ( see
ticket # 101 ) . */
set < int > used_local_vars ;
// Use an empty set for the temporary terms
const temporary_terms_t tt ;
2018-06-04 12:26:16 +02:00
for ( auto equation : equations )
2018-07-17 18:34:07 +02:00
equation - > collectVariables ( SymbolType : : modelLocalVariable , used_local_vars ) ;
2017-02-20 12:18:11 +01:00
2019-04-03 16:32:52 +02:00
output < < R " ( " model_local_variables " : [) " ;
2017-08-28 15:14:11 +02:00
bool printed = false ;
2018-06-04 12:26:16 +02:00
for ( int it : local_variables_vector )
if ( used_local_vars . find ( it ) ! = used_local_vars . end ( ) )
2017-08-28 15:14:11 +02:00
{
if ( printed )
output < < " , " ;
else
printed = true ;
2018-06-04 12:26:16 +02:00
int id = it ;
2017-11-06 15:21:11 +01:00
vector < string > efout ;
expr_t value = local_variables_table . find ( id ) - > second ;
value - > writeJsonExternalFunctionOutput ( efout , tt , tef_terms ) ;
2019-12-16 19:42:59 +01:00
for ( auto it1 = efout . begin ( ) ; it1 ! = efout . end ( ) ; + + it1 )
2017-11-06 15:21:11 +01:00
{
if ( it1 ! = efout . begin ( ) )
output < < " , " ;
output < < * it1 ;
}
if ( ! efout . empty ( ) )
output < < " , " ;
2017-08-28 15:14:11 +02:00
/* We append underscores to avoid name clashes with "g1" or "oo_" (see
also VariableNode : : writeOutput ) */
2019-04-03 16:32:52 +02:00
output < < R " ({ " variable " : " ) " << symbol_table.getName(id) << R " ( __ " ) "
< < R " (, " value " : " ) " ;
2017-08-28 15:14:11 +02:00
value - > writeJsonOutput ( output , tt , tef_terms ) ;
2019-04-03 16:32:52 +02:00
output < < R " ( " } ) " << endl;
2017-08-28 15:14:11 +02:00
}
2017-02-20 12:18:11 +01:00
output < < " ] " ;
}
2008-02-03 11:28:36 +01:00
void
ModelTree : : writeModelEquations ( ostream & output , ExprNodeOutputType output_type ) const
2009-01-23 11:59:37 +01:00
{
2018-03-27 17:14:30 +02:00
temporary_terms_t tt ;
temporary_terms_idxs_t ttidxs ;
2018-05-24 19:29:53 +02:00
writeModelEquations ( output , output_type , tt ) ;
2018-03-27 17:14:30 +02:00
}
2015-09-01 17:39:49 +02:00
2018-03-27 17:14:30 +02:00
void
ModelTree : : writeModelEquations ( ostream & output , ExprNodeOutputType output_type ,
2018-05-24 19:29:53 +02:00
const temporary_terms_t & temporary_terms ) const
2018-03-27 17:14:30 +02:00
{
2019-04-23 11:07:32 +02:00
for ( int eq = 0 ; eq < static_cast < int > ( equations . size ( ) ) ; eq + + )
2009-01-23 11:59:37 +01:00
{
BinaryOpNode * eq_node = equations [ eq ] ;
2018-11-28 14:32:26 +01:00
expr_t lhs = eq_node - > arg1 ;
expr_t rhs = eq_node - > arg2 ;
2008-02-03 11:28:36 +01:00
2009-10-29 18:16:10 +01:00
// Test if the right hand side of the equation is empty.
double vrhs = 1.0 ;
try
{
2010-09-16 19:00:48 +02:00
vrhs = rhs - > eval ( eval_context_t ( ) ) ;
2009-10-29 18:16:10 +01:00
}
2009-12-16 18:13:23 +01:00
catch ( ExprNode : : EvalException & e )
2009-10-29 18:16:10 +01:00
{
}
2009-12-16 18:13:23 +01:00
if ( vrhs ! = 0 ) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
2018-09-05 18:27:13 +02:00
if ( isJuliaOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
{
output < < " @inbounds residual " < < LEFT_ARRAY_SUBSCRIPT ( output_type )
< < eq + ARRAY_SUBSCRIPT_OFFSET ( output_type )
< < RIGHT_ARRAY_SUBSCRIPT ( output_type )
< < " = ( " ;
lhs - > writeOutput ( output , output_type , temporary_terms , temporary_terms_idxs ) ;
output < < " ) - ( " ;
rhs - > writeOutput ( output , output_type , temporary_terms , temporary_terms_idxs ) ;
output < < " ) " < < endl ;
}
else
{
output < < " lhs = " ;
lhs - > writeOutput ( output , output_type , temporary_terms , temporary_terms_idxs ) ;
output < < " ; " < < endl
< < " rhs = " ;
rhs - > writeOutput ( output , output_type , temporary_terms , temporary_terms_idxs ) ;
output < < " ; " < < endl
< < " residual " < < LEFT_ARRAY_SUBSCRIPT ( output_type )
< < eq + ARRAY_SUBSCRIPT_OFFSET ( output_type )
< < RIGHT_ARRAY_SUBSCRIPT ( output_type )
< < " = lhs - rhs; " < < endl ;
}
2009-12-16 18:13:23 +01:00
else // The right hand side of the equation is empty ==> residual=lhs;
2009-10-29 18:16:10 +01:00
{
2018-09-05 18:27:13 +02:00
if ( isJuliaOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " @inbounds " ;
2009-12-16 18:13:23 +01:00
output < < " residual " < < LEFT_ARRAY_SUBSCRIPT ( output_type )
< < eq + ARRAY_SUBSCRIPT_OFFSET ( output_type )
< < RIGHT_ARRAY_SUBSCRIPT ( output_type )
< < " = " ;
2018-03-27 17:14:30 +02:00
lhs - > writeOutput ( output , output_type , temporary_terms , temporary_terms_idxs ) ;
2009-12-16 18:13:23 +01:00
output < < " ; " < < endl ;
2009-10-29 18:16:10 +01:00
}
2009-01-23 11:59:37 +01:00
}
}
2008-02-03 11:28:36 +01:00
2010-01-22 11:03:29 +01:00
void
2010-07-23 11:20:24 +02:00
ModelTree : : compileModelEquations ( ostream & code_file , unsigned int & instruction_number , const temporary_terms_t & tt , const map_idx_t & map_idx , bool dynamic , bool steady_dynamic ) const
2010-01-22 11:03:29 +01:00
{
2019-04-23 11:07:32 +02:00
for ( int eq = 0 ; eq < static_cast < int > ( equations . size ( ) ) ; eq + + )
2010-01-22 11:03:29 +01:00
{
BinaryOpNode * eq_node = equations [ eq ] ;
2018-11-28 14:32:26 +01:00
expr_t lhs = eq_node - > arg1 ;
expr_t rhs = eq_node - > arg2 ;
2010-01-22 17:42:08 +01:00
FNUMEXPR_ fnumexpr ( ModelEquation , eq ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
// Test if the right hand side of the equation is empty.
double vrhs = 1.0 ;
try
{
2010-09-16 19:00:48 +02:00
vrhs = rhs - > eval ( eval_context_t ( ) ) ;
2010-01-22 11:03:29 +01:00
}
catch ( ExprNode : : EvalException & e )
{
}
if ( vrhs ! = 0 ) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
{
2010-07-23 11:20:24 +02:00
lhs - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , dynamic , steady_dynamic ) ;
rhs - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , dynamic , steady_dynamic ) ;
2010-01-22 11:03:29 +01:00
2018-07-18 16:18:26 +02:00
FBINARY_ fbinary { static_cast < int > ( BinaryOpcode : : minus ) } ;
2010-07-23 11:20:24 +02:00
fbinary . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
FSTPR_ fstpr ( eq ) ;
2010-07-23 11:20:24 +02:00
fstpr . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
}
else // The right hand side of the equation is empty ==> residual=lhs;
{
2010-07-23 11:20:24 +02:00
lhs - > compile ( code_file , instruction_number , false , temporary_terms , map_idx , dynamic , steady_dynamic ) ;
2010-01-22 11:03:29 +01:00
FSTPR_ fstpr ( eq ) ;
2010-07-23 11:20:24 +02:00
fstpr . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
}
}
}
void
2018-06-27 15:01:31 +02:00
ModelTree : : Write_Inf_To_Bin_File ( const string & filename ,
2011-02-04 16:25:38 +01:00
int & u_count_int , bool & file_open , bool is_two_boundaries , int block_mfs ) const
2010-01-22 11:03:29 +01:00
{
int j ;
std : : ofstream SaveCode ;
if ( file_open )
2018-06-27 15:01:31 +02:00
SaveCode . open ( filename , ios : : out | ios : : in | ios : : binary | ios : : ate ) ;
2010-01-22 11:03:29 +01:00
else
2018-06-27 15:01:31 +02:00
SaveCode . open ( filename , ios : : out | ios : : binary ) ;
2010-01-22 11:03:29 +01:00
if ( ! SaveCode . is_open ( ) )
{
2019-04-03 16:32:52 +02:00
cerr < < R " (Error : Can't open file " ) " << filename << R " ( " for writing) " < < endl ;
2010-01-22 11:03:29 +01:00
exit ( EXIT_FAILURE ) ;
}
u_count_int = 0 ;
2019-12-16 19:42:59 +01:00
for ( const auto & [ indices , d1 ] : derivatives [ 1 ] )
2010-01-22 11:03:29 +01:00
{
2019-12-16 19:42:59 +01:00
int deriv_id = indices [ 1 ] ;
2018-07-17 18:34:07 +02:00
if ( getTypeByDerivID ( deriv_id ) = = SymbolType : : endogenous )
2010-01-22 11:03:29 +01:00
{
2019-12-16 19:42:59 +01:00
int eq = indices [ 0 ] ;
2010-01-22 11:03:29 +01:00
int symb = getSymbIDByDerivID ( deriv_id ) ;
int var = symbol_table . getTypeSpecificID ( symb ) ;
int lag = getLagByDerivID ( deriv_id ) ;
SaveCode . write ( reinterpret_cast < char * > ( & eq ) , sizeof ( eq ) ) ;
int varr = var + lag * block_mfs ;
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 + + ;
}
}
if ( is_two_boundaries )
2019-12-20 16:59:30 +01:00
u_count_int + = symbol_table . endo_nbr ( ) ;
2019-04-23 11:07:32 +02:00
for ( j = 0 ; j < static_cast < int > ( symbol_table . endo_nbr ( ) ) ; j + + )
2010-01-22 11:03:29 +01:00
SaveCode . write ( reinterpret_cast < char * > ( & j ) , sizeof ( j ) ) ;
2019-04-23 11:07:32 +02:00
for ( j = 0 ; j < static_cast < int > ( symbol_table . endo_nbr ( ) ) ; j + + )
2010-01-22 11:03:29 +01:00
SaveCode . write ( reinterpret_cast < char * > ( & j ) , sizeof ( j ) ) ;
SaveCode . close ( ) ;
}
2009-04-30 15:14:33 +02:00
void
2019-12-16 19:42:59 +01:00
ModelTree : : writeLatexModelFile ( const string & mod_basename , const string & latex_basename , ExprNodeOutputType output_type , bool write_equation_tags ) const
2009-04-30 15:14:33 +02:00
{
2019-04-04 17:01:37 +02:00
filesystem : : create_directories ( mod_basename + " /latex " ) ;
2019-07-11 17:33:53 +02:00
2015-07-15 08:58:15 +02:00
ofstream output , content_output ;
2019-07-11 17:33:53 +02:00
string filename = mod_basename + " /latex/ " + latex_basename + " .tex " ;
string content_filename = mod_basename + " /latex/ " + latex_basename + " _content " + " .tex " ;
2018-06-27 15:12:12 +02:00
output . open ( filename , ios : : out | ios : : binary ) ;
2009-04-30 15:14:33 +02:00
if ( ! output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2018-06-27 15:12:12 +02:00
content_output . open ( content_filename , ios : : out | ios : : binary ) ;
2015-07-15 08:58:15 +02:00
if ( ! content_output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < content_filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2019-04-03 16:32:52 +02:00
output < < R " ( \ documentclass[10pt,a4paper]{article}) " < < endl
< < R " ( \ usepackage[landscape]{geometry}) " < < endl
< < R " ( \ usepackage{fullpage}) " < < endl
< < R " ( \ usepackage{amsfonts}) " < < endl
< < R " ( \ usepackage{breqn}) " < < endl
< < R " ( \b egin{document}) " < < endl
< < R " ( \f ootnotesize) " < < endl ;
2009-04-30 15:14:33 +02:00
// Write model local variables
2018-06-04 12:26:16 +02:00
for ( int id : local_variables_vector )
2009-04-30 15:14:33 +02:00
{
2017-08-28 15:14:11 +02:00
expr_t value = local_variables_table . find ( id ) - > second ;
2009-04-30 15:14:33 +02:00
2019-04-03 16:32:52 +02:00
content_output < < R " ( \b egin{dmath*}) " < < endl
2015-09-16 12:20:20 +02:00
< < symbol_table . getTeXName ( id ) < < " = " ;
2009-04-30 15:14:33 +02:00
// Use an empty set for the temporary terms
2015-07-15 08:58:15 +02:00
value - > writeOutput ( content_output , output_type ) ;
2019-04-03 16:32:52 +02:00
content_output < < endl < < R " ( \ end{dmath*}) " < < endl ;
2009-04-30 15:14:33 +02:00
}
2019-04-23 11:07:32 +02:00
for ( int eq = 0 ; eq < static_cast < int > ( equations . size ( ) ) ; eq + + )
2009-04-30 15:14:33 +02:00
{
2017-04-05 10:28:37 +02:00
content_output < < " % Equation " < < eq + 1 < < endl ;
2017-04-04 15:28:27 +02:00
if ( write_equation_tags )
2017-04-05 10:28:37 +02:00
{
2020-01-07 12:41:27 +01:00
auto escape_special_latex_symbols
= [ ] ( string str )
{
const regex special_latex_chars ( R " ([&%$#_{}]) " ) ;
const regex backslash ( R " ( \\ ) " ) ;
const regex tilde ( R " (~) " ) ;
const regex carrot ( R " ( \ ^) " ) ;
const regex textbackslash ( R " ( \\ textbackslash) " ) ;
str = regex_replace ( str , backslash , R " ( \t extbackslash) " ) ;
str = regex_replace ( str , special_latex_chars , R " ( \ $&) " ) ;
str = regex_replace ( str , carrot , R " ( \ ^{}) " ) ;
str = regex_replace ( str , tilde , R " ( \t extasciitilde{}) " ) ;
return regex_replace ( str , textbackslash , R " ( \t extbackslash{}) " ) ;
} ;
2018-06-05 14:47:36 +02:00
bool wrote_eq_tag = false ;
2019-12-16 19:42:59 +01:00
for ( const auto & [ tagged_eq , tag_pair ] : equation_tags )
if ( tagged_eq = = eq )
2017-04-05 10:28:37 +02:00
{
if ( ! wrote_eq_tag )
2019-04-03 16:32:52 +02:00
content_output < < R " ( \n oindent[) " ;
2017-04-05 10:28:37 +02:00
else
content_output < < " , " ;
2020-01-07 12:41:27 +01:00
content_output < < escape_special_latex_symbols ( tag_pair . first ) ;
2017-04-05 10:28:37 +02:00
2019-12-16 19:42:59 +01:00
if ( ! ( tag_pair . second . empty ( ) ) )
2020-01-07 12:41:27 +01:00
content_output < < " = ` " < < escape_special_latex_symbols ( tag_pair . second ) < < " ' " ;
2017-04-05 10:28:37 +02:00
wrote_eq_tag = true ;
}
2018-06-05 14:47:36 +02:00
if ( wrote_eq_tag )
2020-01-07 12:26:40 +01:00
content_output < < " ] " < < endl ;
2017-04-05 10:28:37 +02:00
}
2017-04-04 15:28:27 +02:00
2019-04-03 16:32:52 +02:00
content_output < < R " ( \b egin{dmath}) " < < endl ;
2010-03-09 12:16:32 +01:00
// Here it is necessary to cast to superclass ExprNode, otherwise the overloaded writeOutput() method is not found
2015-07-15 08:58:15 +02:00
dynamic_cast < ExprNode * > ( equations [ eq ] ) - > writeOutput ( content_output , output_type ) ;
2019-04-03 16:32:52 +02:00
content_output < < endl < < R " ( \ end{dmath}) " < < endl ;
2009-04-30 15:14:33 +02:00
}
2019-07-11 17:33:53 +02:00
output < < R " ( \ include{) " < < latex_basename + " _content " < < " } " < < endl
2019-04-03 16:32:52 +02:00
< < R " ( \ end{document}) " < < endl ;
2009-04-30 15:14:33 +02:00
output . close ( ) ;
2015-07-15 08:58:15 +02:00
content_output . close ( ) ;
2009-04-30 15:14:33 +02:00
}
2008-02-03 11:28:36 +01:00
void
2014-01-27 16:41:43 +01:00
ModelTree : : addEquation ( expr_t eq , int lineno )
2009-01-23 11:59:37 +01:00
{
2019-12-16 19:42:59 +01:00
auto beq = dynamic_cast < BinaryOpNode * > ( eq ) ;
assert ( beq & & beq - > op_code = = BinaryOpcode : : equal ) ;
2008-02-03 11:28:36 +01:00
equations . push_back ( beq ) ;
2014-01-27 16:41:43 +01:00
equations_lineno . push_back ( lineno ) ;
2008-02-03 11:28:36 +01:00
}
2009-09-02 16:37:59 +02:00
2019-11-18 17:13:49 +01:00
vector < int >
ModelTree : : includeExcludeEquations ( set < pair < string , string > > & eqs , bool exclude_eqs ,
vector < BinaryOpNode * > & equations , vector < int > & equations_lineno ,
vector < pair < int , pair < string , string > > > & equation_tags ,
multimap < pair < string , string > , int > & equation_tags_xref , bool static_equations ) const
{
vector < int > excluded_vars ;
if ( equations . empty ( ) )
return excluded_vars ;
// Get equation numbers of tags
set < int > tag_eqns ;
2019-12-20 16:59:30 +01:00
for ( auto & it : eqs )
2019-11-18 17:13:49 +01:00
if ( equation_tags_xref . find ( it ) ! = equation_tags_xref . end ( ) )
{
auto range = equation_tags_xref . equal_range ( it ) ;
2019-12-20 16:59:30 +01:00
for_each ( range . first , range . second , [ & tag_eqns ] ( auto & x ) { tag_eqns . insert ( x . second ) ; } ) ;
2019-11-18 17:13:49 +01:00
eqs . erase ( it ) ;
}
if ( tag_eqns . empty ( ) )
return excluded_vars ;
set < int > eqns ;
if ( exclude_eqs )
eqns = tag_eqns ;
else
for ( size_t i = 0 ; i < equations . size ( ) ; i + + )
if ( tag_eqns . find ( i ) = = tag_eqns . end ( ) )
eqns . insert ( i ) ;
// remove from equations, equations_lineno, equation_tags, equation_tags_xref
vector < BinaryOpNode * > new_eqns ;
vector < int > new_equations_lineno ;
map < int , int > old_eqn_num_2_new ;
for ( size_t i = 0 ; i < equations . size ( ) ; i + + )
if ( eqns . find ( i ) ! = eqns . end ( ) )
{
bool found = false ;
2019-12-16 19:42:59 +01:00
for ( const auto & [ tagged_eq , tag_pair ] : equation_tags )
if ( tagged_eq = = static_cast < int > ( i ) & & tag_pair . first = = " endogenous " )
2019-11-18 17:13:49 +01:00
{
found = true ;
2019-12-16 19:42:59 +01:00
excluded_vars . push_back ( symbol_table . getID ( tag_pair . second ) ) ;
2019-11-18 17:13:49 +01:00
break ;
}
if ( ! found )
{
set < pair < int , int > > result ;
equations [ i ] - > arg1 - > collectDynamicVariables ( SymbolType : : endogenous , result ) ;
if ( result . size ( ) = = 1 )
excluded_vars . push_back ( result . begin ( ) - > first ) ;
else
{
cerr < < " ERROR: Equation " < < i
< < " has been excluded but does not have a single variable on LHS or `endogenous` tag " < < endl ;
exit ( EXIT_FAILURE ) ;
}
}
}
else
{
new_eqns . emplace_back ( equations [ i ] ) ;
old_eqn_num_2_new [ i ] = new_eqns . size ( ) - 1 ;
new_equations_lineno . emplace_back ( equations_lineno [ i ] ) ;
}
int n_excl = equations . size ( ) - new_eqns . size ( ) ;
equations = new_eqns ;
equations_lineno = new_equations_lineno ;
equation_tags . erase ( remove_if ( equation_tags . begin ( ) , equation_tags . end ( ) ,
2019-12-20 16:59:30 +01:00
[ & ] ( const auto & it ) { return eqns . find ( it . first ) ! = eqns . end ( ) ; } ) ,
2019-11-18 17:13:49 +01:00
equation_tags . end ( ) ) ;
2019-12-20 16:59:30 +01:00
for ( auto & it : old_eqn_num_2_new )
for ( auto & it1 : equation_tags )
2019-11-18 17:13:49 +01:00
if ( it1 . first = = it . first )
it1 . first = it . second ;
equation_tags_xref . clear ( ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & it : equation_tags )
2019-11-18 17:13:49 +01:00
equation_tags_xref . emplace ( it . second , it . first ) ;
if ( ! static_equations )
for ( size_t i = 0 ; i < excluded_vars . size ( ) ; i + + )
for ( size_t j = i + 1 ; j < excluded_vars . size ( ) ; j + + )
if ( excluded_vars [ i ] = = excluded_vars [ j ] )
{
cerr < < " Error: Variable " < < symbol_table . getName ( i ) < < " was excluded twice "
< < " via in/exclude_eqs option " < < endl ;
exit ( EXIT_FAILURE ) ;
}
cout < < " Excluded " < < n_excl < < ( static_equations ? " static " : " dynamic " )
2019-12-20 16:59:30 +01:00
< < " equation " < < ( n_excl > 1 ? " s " : " " ) < < " via in/exclude_eqs option " < < endl ;
2019-11-18 17:13:49 +01:00
return excluded_vars ;
}
2019-01-28 14:57:30 +01:00
void
ModelTree : : simplifyEquations ( )
{
size_t last_subst_table_size = 0 ;
2019-01-29 17:29:24 +01:00
map < VariableNode * , NumConstNode * > subst_table ;
2019-01-28 14:57:30 +01:00
findConstantEquations ( subst_table ) ;
while ( subst_table . size ( ) ! = last_subst_table_size )
{
last_subst_table_size = subst_table . size ( ) ;
2019-12-20 16:59:30 +01:00
for ( auto & equation : equations )
2019-01-28 14:57:30 +01:00
equation = dynamic_cast < BinaryOpNode * > ( equation - > replaceVarsInEquation ( subst_table ) ) ;
subst_table . clear ( ) ;
findConstantEquations ( subst_table ) ;
}
}
void
2019-01-29 17:29:24 +01:00
ModelTree : : findConstantEquations ( map < VariableNode * , NumConstNode * > & subst_table ) const
2019-01-28 14:57:30 +01:00
{
2019-12-20 16:59:30 +01:00
for ( auto & equation : equations )
2019-01-28 14:57:30 +01:00
equation - > findConstantEquations ( subst_table ) ;
}
2009-09-02 16:37:59 +02:00
void
2018-06-04 14:17:36 +02:00
ModelTree : : addEquation ( expr_t eq , int lineno , const vector < pair < string , string > > & eq_tags )
2009-09-02 16:37:59 +02:00
{
2017-02-24 11:20:54 +01:00
int n = equations . size ( ) ;
2019-12-20 16:59:30 +01:00
for ( const auto & eq_tag : eq_tags )
2019-10-29 12:56:53 +01:00
{
equation_tags . emplace_back ( n , eq_tag ) ;
equation_tags_xref . emplace ( eq_tag , n ) ;
}
2014-01-27 16:41:43 +01:00
addEquation ( eq , lineno ) ;
2009-09-02 16:37:59 +02:00
}
2009-09-30 17:10:31 +02:00
void
2010-09-16 19:18:45 +02:00
ModelTree : : addAuxEquation ( expr_t eq )
2009-09-30 17:10:31 +02:00
{
2019-12-16 19:42:59 +01:00
auto beq = dynamic_cast < BinaryOpNode * > ( eq ) ;
assert ( beq & & beq - > op_code = = BinaryOpcode : : equal ) ;
2009-09-30 17:10:31 +02:00
aux_equations . push_back ( beq ) ;
}
2010-10-15 19:05:16 +02:00
void
2019-07-05 18:36:10 +02:00
ModelTree : : addTrendVariables ( const vector < int > & trend_vars , expr_t growth_factor ) noexcept ( false )
2010-10-15 19:05:16 +02:00
{
2019-07-05 18:36:10 +02:00
for ( int id : trend_vars )
if ( trend_symbols_map . find ( id ) ! = trend_symbols_map . end ( ) )
throw TrendException ( symbol_table . getName ( id ) ) ;
2010-10-15 19:05:16 +02:00
else
2019-07-05 18:36:10 +02:00
trend_symbols_map [ id ] = growth_factor ;
2010-10-15 19:05:16 +02:00
}
void
2019-07-05 18:36:10 +02:00
ModelTree : : addNonstationaryVariables ( const vector < int > & nonstationary_vars , bool log_deflator , expr_t deflator ) noexcept ( false )
2010-10-15 19:05:16 +02:00
{
2019-07-05 18:36:10 +02:00
for ( int id : nonstationary_vars )
if ( nonstationary_symbols_map . find ( id ) ! = nonstationary_symbols_map . end ( ) )
throw TrendException ( symbol_table . getName ( id ) ) ;
2010-10-15 19:05:16 +02:00
else
2019-07-05 18:36:10 +02:00
nonstationary_symbols_map [ id ] = { log_deflator , deflator } ;
2010-10-15 19:05:16 +02:00
}
2011-06-22 11:56:07 +02:00
void
ModelTree : : initializeVariablesAndEquations ( )
{
2017-06-14 23:49:10 +02:00
for ( size_t j = 0 ; j < equations . size ( ) ; j + + )
2018-10-10 17:06:19 +02:00
equation_reordered . push_back ( j ) ;
for ( int j = 0 ; j < symbol_table . endo_nbr ( ) ; j + + )
variable_reordered . push_back ( j ) ;
2011-06-22 11:56:07 +02:00
}
void
ModelTree : : set_cutoff_to_zero ( )
{
cutoff = 0 ;
}
2011-08-19 15:05:57 +02:00
void
ModelTree : : jacobianHelper ( ostream & output , int eq_nb , int col_nb , ExprNodeOutputType output_type ) const
{
2018-09-05 18:27:13 +02:00
if ( isJuliaOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " @inbounds " ;
2015-08-20 12:12:17 +02:00
output < < " g1 " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) ;
2018-09-05 18:27:13 +02:00
if ( isMatlabOutput ( output_type ) | | isJuliaOutput ( output_type ) )
2011-08-19 15:05:57 +02:00
output < < eq_nb + 1 < < " , " < < col_nb + 1 ;
else
output < < eq_nb + col_nb * equations . size ( ) ;
output < < RIGHT_ARRAY_SUBSCRIPT ( output_type ) ;
}
void
ModelTree : : sparseHelper ( int order , ostream & output , int row_nb , int col_nb , ExprNodeOutputType output_type ) const
{
2018-03-27 17:14:30 +02:00
output < < " v " < < order < < LEFT_ARRAY_SUBSCRIPT ( output_type ) ;
2018-09-05 18:27:13 +02:00
if ( isMatlabOutput ( output_type ) | | isJuliaOutput ( output_type ) )
2011-08-19 15:05:57 +02:00
output < < row_nb + 1 < < " , " < < col_nb + 1 ;
else
2018-11-15 16:39:53 +01:00
output < < row_nb + col_nb * NNZDerivatives [ order ] ;
2011-08-19 15:05:57 +02:00
output < < RIGHT_ARRAY_SUBSCRIPT ( output_type ) ;
}
2012-11-29 18:07:48 +01:00
void
2016-05-18 12:26:19 +02:00
ModelTree : : computeParamsDerivatives ( int paramsDerivsOrder )
2012-11-29 18:07:48 +01:00
{
2018-11-22 15:41:11 +01:00
assert ( paramsDerivsOrder > = 1 ) ;
2012-11-29 18:07:48 +01:00
set < int > deriv_id_set ;
addAllParamDerivId ( deriv_id_set ) ;
2016-03-25 15:38:49 +01:00
2018-11-22 15:41:11 +01:00
// First-order derivatives w.r.t. params
2018-06-04 12:26:16 +02:00
for ( int param : deriv_id_set )
2012-11-29 18:07:48 +01:00
{
2019-04-23 11:07:32 +02:00
for ( int eq = 0 ; eq < static_cast < int > ( equations . size ( ) ) ; eq + + )
2012-11-29 18:07:48 +01:00
{
2018-11-22 15:41:11 +01:00
expr_t d = equations [ eq ] - > getDerivative ( param ) ;
if ( d = = Zero )
2012-11-29 18:07:48 +01:00
continue ;
2018-11-22 15:41:11 +01:00
params_derivatives [ { 0 , 1 } ] [ { eq , param } ] = d ;
2012-11-29 18:07:48 +01:00
}
2019-04-23 11:07:32 +02:00
for ( int endoOrd = 1 ; endoOrd < static_cast < int > ( derivatives . size ( ) ) ; endoOrd + + )
2019-12-16 19:42:59 +01:00
for ( const auto & [ indices , dprev ] : derivatives [ endoOrd ] )
2016-05-12 12:02:34 +02:00
{
2019-12-16 19:42:59 +01:00
expr_t d = dprev - > getDerivative ( param ) ;
2018-11-22 15:41:11 +01:00
if ( d = = Zero )
2016-05-12 12:02:34 +02:00
continue ;
2019-12-16 19:42:59 +01:00
vector < int > new_indices = indices ;
new_indices . push_back ( param ) ;
params_derivatives [ { endoOrd , 1 } ] [ new_indices ] = d ;
2016-05-12 12:02:34 +02:00
}
2018-11-22 15:41:11 +01:00
}
2012-11-29 18:07:48 +01:00
2018-11-22 15:41:11 +01:00
// Higher-order derivatives w.r.t. parameters
2019-04-23 11:07:32 +02:00
for ( int endoOrd = 0 ; endoOrd < static_cast < int > ( derivatives . size ( ) ) ; endoOrd + + )
2018-11-22 15:41:11 +01:00
for ( int paramOrd = 2 ; paramOrd < = paramsDerivsOrder ; paramOrd + + )
2019-12-16 19:42:59 +01:00
for ( const auto & [ indices , dprev ] : params_derivatives [ { endoOrd , paramOrd - 1 } ] )
2018-11-22 15:41:11 +01:00
for ( int param : deriv_id_set )
{
2019-12-16 19:42:59 +01:00
if ( indices . back ( ) > param )
2018-11-22 15:41:11 +01:00
continue ;
2016-05-18 10:33:45 +02:00
2019-12-16 19:42:59 +01:00
expr_t d = dprev - > getDerivative ( param ) ;
2018-11-22 15:41:11 +01:00
if ( d = = Zero )
continue ;
2019-12-16 19:42:59 +01:00
vector < int > new_indices = indices ;
new_indices . push_back ( param ) ;
2018-11-22 15:41:11 +01:00
// At this point, indices of both endogenous and parameters are sorted in non-decreasing order
2019-12-16 19:42:59 +01:00
params_derivatives [ { endoOrd , paramOrd } ] [ new_indices ] = d ;
2018-11-22 15:41:11 +01:00
}
2012-11-29 18:07:48 +01:00
}
void
ModelTree : : computeParamsDerivativesTemporaryTerms ( )
{
2018-11-30 12:22:13 +01:00
map < expr_t , pair < int , pair < int , int > > > reference_count ;
2018-11-16 18:24:06 +01:00
/* The temp terms should be constructed in the same order as the for loops in
{ Static , Dynamic } Model : : write { Json , } ParamsDerivativesFile ( ) */
2018-11-30 12:22:13 +01:00
params_derivs_temporary_terms . clear ( ) ;
2019-12-16 19:42:59 +01:00
for ( const auto & [ order , derivs ] : params_derivatives )
for ( const auto & [ indices , d ] : derivs )
d - > computeTemporaryTerms ( order , params_derivs_temporary_terms ,
reference_count , true ) ;
2018-05-28 15:23:15 +02:00
int idx = 0 ;
2019-12-16 19:42:59 +01:00
for ( auto & [ mlv , value ] : temporary_terms_mlv )
params_derivs_temporary_terms_idxs [ mlv ] = idx + + ;
for ( const auto & [ order , tts ] : params_derivs_temporary_terms )
for ( const auto & tt : tts )
2018-11-30 12:22:13 +01:00
params_derivs_temporary_terms_idxs [ tt ] = idx + + ;
2012-11-29 18:07:48 +01:00
}
2013-10-29 11:47:59 +01:00
2017-06-14 07:01:31 +02:00
bool
ModelTree : : isNonstationary ( int symb_id ) const
2013-10-29 11:47:59 +01:00
{
2019-12-16 19:42:59 +01:00
return nonstationary_symbols_map . find ( symb_id ) ! = nonstationary_symbols_map . end ( ) ;
2013-10-29 11:47:59 +01:00
}
2017-02-02 15:09:43 +01:00
void
2017-02-20 12:18:11 +01:00
ModelTree : : writeJsonModelEquations ( ostream & output , bool residuals ) const
2017-02-02 15:09:43 +01:00
{
2017-02-20 12:18:11 +01:00
if ( residuals )
2019-04-03 16:32:52 +02:00
output < < endl < < R " ( " residuals " :[) " < < endl ;
2017-02-20 12:18:11 +01:00
else
2019-04-03 16:32:52 +02:00
output < < endl < < R " ( " model " :[) " < < endl ;
2019-04-23 11:07:32 +02:00
for ( int eq = 0 ; eq < static_cast < int > ( equations . size ( ) ) ; eq + + )
2017-02-02 15:09:43 +01:00
{
2017-02-20 12:18:11 +01:00
if ( eq > 0 )
output < < " , " ;
2017-03-15 12:52:55 +01:00
BinaryOpNode * eq_node = equations [ eq ] ;
2018-11-28 14:32:26 +01:00
expr_t lhs = eq_node - > arg1 ;
expr_t rhs = eq_node - > arg2 ;
2017-03-15 12:52:55 +01:00
2017-02-20 12:18:11 +01:00
if ( residuals )
2017-02-02 15:09:43 +01:00
{
2019-04-03 16:32:52 +02:00
output < < R " ({ " residual " : {) "
< < R " ( " lhs " : " ) " ;
2018-05-29 11:59:42 +02:00
lhs - > writeJsonOutput ( output , temporary_terms , { } ) ;
2019-04-03 16:32:52 +02:00
output < < R " ( " ) " ;
2017-02-20 12:18:11 +01:00
2019-04-03 16:32:52 +02:00
output < < R " (, " rhs " : " ) " ;
2018-05-29 11:59:42 +02:00
rhs - > writeJsonOutput ( output , temporary_terms , { } ) ;
2019-04-03 16:32:52 +02:00
output < < R " ( " ) " ;
2017-02-20 12:18:11 +01:00
try
{
// Test if the right hand side of the equation is empty.
if ( rhs - > eval ( eval_context_t ( ) ) ! = 0 )
{
2019-04-03 16:32:52 +02:00
output < < R " (, " rhs " : " ) " ;
2018-05-29 11:59:42 +02:00
rhs - > writeJsonOutput ( output , temporary_terms , { } ) ;
2019-04-03 16:32:52 +02:00
output < < R " ( " ) " ;
2017-02-20 12:18:11 +01:00
}
}
catch ( ExprNode : : EvalException & e )
2017-02-02 15:09:43 +01:00
{
}
output < < " } " ;
}
2017-02-20 12:18:11 +01:00
else
{
2019-04-03 16:32:52 +02:00
output < < R " ({ " lhs " : " ) " ;
2019-12-16 19:42:59 +01:00
lhs - > writeJsonOutput ( output , { } , { } ) ;
2019-04-03 16:32:52 +02:00
output < < R " ( " , " rhs " : " ) " ;
2019-12-16 19:42:59 +01:00
rhs - > writeJsonOutput ( output , { } , { } ) ;
2019-04-03 16:32:52 +02:00
output < < R " ( " ) "
< < R " (, " line " : ) " < < equations_lineno [ eq ] ;
2017-02-20 12:18:11 +01:00
2019-12-16 19:42:59 +01:00
if ( auto eqtags = getEquationTags ( eq ) ;
! eqtags . empty ( ) )
2017-02-20 12:18:11 +01:00
{
2019-04-03 16:32:52 +02:00
output < < R " (, " tags " : {) " ;
2017-02-20 12:18:11 +01:00
int i = 0 ;
2019-12-16 19:42:59 +01:00
for ( const auto & [ name , value ] : eqtags )
2017-02-20 12:18:11 +01:00
{
if ( i ! = 0 )
output < < " , " ;
2019-12-16 19:42:59 +01:00
output < < R " ( " ) " << name << R " ( " : " ) " << value << R " ( " ) " ;
i + + ;
2017-02-20 12:18:11 +01:00
}
output < < " } " ;
eqtags . clear ( ) ;
}
}
output < < " } " < < endl ;
2017-02-02 15:09:43 +01:00
}
output < < endl < < " ] " < < endl ;
}
2018-10-26 11:44:26 +02:00
string
ModelTree : : matlab_arch ( const string & mexext )
{
if ( mexext = = " mexglx " )
return " glnx86 " ;
else if ( mexext = = " mexa64 " )
return " glnxa64 " ;
if ( mexext = = " mexw32 " )
return " win32 " ;
else if ( mexext = = " mexw64 " )
return " win64 " ;
else if ( mexext = = " mexmaci " )
2019-11-04 15:50:26 +01:00
{
cerr < < " 32-bit MATLAB not supported on macOS " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2018-10-26 11:44:26 +02:00
else if ( mexext = = " mexmaci64 " )
return " maci64 " ;
else
{
cerr < < " ERROR: 'mexext' option to preprocessor incorrectly set, needed with 'use_dll' " < < endl ;
exit ( EXIT_FAILURE ) ;
}
}
void
2019-12-02 19:21:14 +01:00
ModelTree : : compileDll ( const string & basename , const string & static_or_dynamic , const string & mexext , const filesystem : : path & matlabroot , const filesystem : : path & dynareroot ) const
2018-10-26 11:44:26 +02:00
{
const string opt_flags = " -O3 -g0 --param ira-max-conflict-table-size=1 -fno-forward-propagate -fno-gcse -fno-dce -fno-dse -fno-tree-fre -fno-tree-pre -fno-tree-cselim -fno-tree-dse -fno-tree-dce -fno-tree-pta -fno-gcse-after-reload " ;
2019-04-04 17:01:37 +02:00
filesystem : : path compiler ;
2018-10-26 11:44:26 +02:00
ostringstream flags ;
string libs ;
if ( mexext = = " mex " )
{
// Octave
compiler = matlabroot / " bin " / " mkoctfile " ;
flags < < " --mex " ;
}
else
{
// MATLAB
compiler = " gcc " ;
string arch = matlab_arch ( mexext ) ;
auto include_dir = matlabroot / " extern " / " include " ;
flags < < " -I " < < include_dir ;
auto bin_dir = matlabroot / " bin " / arch ;
flags < < " -L " < < bin_dir ;
flags < < " -fexceptions -DNDEBUG " ;
2019-07-11 16:11:24 +02:00
libs = " -lmex -lmx " ;
2018-10-26 11:44:26 +02:00
if ( mexext = = " mexglx " | | mexext = = " mexa64 " )
{
// GNU/Linux
flags < < " -D_GNU_SOURCE -fPIC -pthread "
< < " -shared -Wl,--no-undefined -Wl,-rpath-link, " < < bin_dir ;
libs + = " -lm -lstdc++ " ;
if ( mexext = = " mexglx " )
flags < < " -D_FILE_OFFSET_BITS=64 -m32 " ;
else
flags < < " -fno-omit-frame-pointer " ;
}
else if ( mexext = = " mexw32 " | | mexext = = " mexw64 " )
{
// Windows
flags < < " -static-libgcc -static-libstdc++ -shared " ;
// Put the MinGW environment shipped with Dynare in the path
2019-12-16 19:42:59 +01:00
auto mingwpath = dynareroot / ( string { " mingw " } + ( mexext = = " mexw32 " ? " 32 " : " 64 " ) ) / " bin " ;
2018-10-26 11:44:26 +02:00
string newpath = " PATH= " + mingwpath . string ( ) + ' ; ' + string { getenv ( " PATH " ) } ;
if ( putenv ( const_cast < char * > ( newpath . c_str ( ) ) ) ! = 0 )
{
cerr < < " Can't set PATH " < < endl ;
exit ( EXIT_FAILURE ) ;
}
2019-12-20 16:59:30 +01:00
}
2018-10-26 11:44:26 +02:00
else
{
// macOS
2019-11-04 15:50:26 +01:00
# ifdef __APPLE__
char dynare_m_path [ PATH_MAX ] ;
uint32_t size = PATH_MAX ;
string gcc_relative_path = " " ;
if ( _NSGetExecutablePath ( dynare_m_path , & size ) = = 0 )
{
string str = dynare_m_path ;
gcc_relative_path = str . substr ( 0 , str . find_last_of ( " / " ) ) + " /../../.brew/bin/gcc-9 " ;
}
if ( filesystem : : exists ( gcc_relative_path ) )
compiler = gcc_relative_path ;
else if ( filesystem : : exists ( " /usr/local/bin/gcc-9 " ) )
compiler = " /usr/local/bin/gcc-9 " ;
else
{
cerr < < " ERROR: You must install gcc-9 on your system before using the `use_dll` option of Dynare. "
< < " You can do this via the Dynare installation package. " < < endl ;
exit ( EXIT_FAILURE ) ;
}
# endif
flags < < " -fno-common -arch x86_64 -mmacosx-version-min=10.9 -Wl,-twolevel_namespace -undefined error -bundle " ;
2018-10-26 11:44:26 +02:00
libs + = " -lm -lstdc++ " ;
}
}
2019-04-04 17:01:37 +02:00
auto model_dir = filesystem : : path { basename } / " model " / " src " ;
filesystem : : path main_src { model_dir / ( static_or_dynamic + " .c " ) } ,
2018-10-26 11:44:26 +02:00
mex_src { model_dir / ( static_or_dynamic + " _mex.c " ) } ;
2019-04-04 17:01:37 +02:00
filesystem : : path mex_dir { " + " + basename } ;
filesystem : : path binary { mex_dir / ( static_or_dynamic + " . " + mexext ) } ;
2018-10-26 11:44:26 +02:00
ostringstream cmd ;
# ifdef _WIN32
/* On Windows, system() hands the command over to "cmd.exe /C". We need to
enclose the whole command line within double quotes if we want the inner
quotes to be correctly handled . See " cmd /? " for more details . */
cmd < < ' " ' ;
# endif
2019-12-02 19:21:14 +01:00
if ( user_set_compiler . empty ( ) )
cmd < < compiler < < " " ;
else
if ( ! filesystem : : exists ( user_set_compiler ) )
{
cerr < < " Error: The specified compiler ' " < < user_set_compiler < < " ' cannot be found on your system " < < endl ;
exit ( EXIT_FAILURE ) ;
}
else
cmd < < user_set_compiler < < " " ;
if ( user_set_subst_flags . empty ( ) )
cmd < < opt_flags < < " " < < flags . str ( ) < < " " ;
else
cmd < < user_set_subst_flags < < " " ;
if ( ! user_set_add_flags . empty ( ) )
cmd < < user_set_add_flags < < " " ;
cmd < < main_src < < " " < < mex_src < < " -o " < < binary < < " " ;
if ( user_set_subst_libs . empty ( ) )
cmd < < libs ;
else
cmd < < user_set_subst_libs ;
if ( ! user_set_add_libs . empty ( ) )
cmd < < " " < < user_set_add_libs ;
2018-10-26 11:44:26 +02:00
# ifdef _WIN32
cmd < < ' " ' ;
# endif
cout < < " Compiling " < < static_or_dynamic < < " MEX... " < < endl < < cmd . str ( ) < < endl ;
if ( system ( cmd . str ( ) . c_str ( ) ) )
{
cerr < < " Compilation failed " < < endl ;
exit ( EXIT_FAILURE ) ;
}
}
2019-12-03 14:19:32 +01:00
void
ModelTree : : reorderAuxiliaryEquations ( )
{
using namespace boost ;
// Create the mapping between auxiliary variables and auxiliary equations
int n = static_cast < int > ( aux_equations . size ( ) ) ;
map < int , int > auxEndoToEq ;
for ( int i = 0 ; i < n ; i + + )
{
auto varexpr = dynamic_cast < VariableNode * > ( aux_equations [ i ] - > arg1 ) ;
assert ( varexpr & & symbol_table . getType ( varexpr - > symb_id ) = = SymbolType : : endogenous ) ;
auxEndoToEq [ varexpr - > symb_id ] = i ;
}
assert ( static_cast < int > ( auxEndoToEq . size ( ) ) = = n ) ;
/* Construct the directed acyclic graph where auxiliary equations are
vertices and edges represent dependency relationships . */
using Graph = adjacency_list < vecS , vecS , directedS > ;
Graph g ( n ) ;
for ( int i = 0 ; i < n ; i + + )
{
set < int > endos ;
aux_equations [ i ] - > collectVariables ( SymbolType : : endogenous , endos ) ;
for ( int endo : endos )
2019-12-16 19:42:59 +01:00
if ( auto it = auxEndoToEq . find ( endo ) ;
it ! = auxEndoToEq . end ( ) & & it - > second ! = i )
add_edge ( i , it - > second , g ) ;
2019-12-03 14:19:32 +01:00
}
// Topological sort of the graph
using Vertex = graph_traits < Graph > : : vertex_descriptor ;
vector < Vertex > ordered ;
topological_sort ( g , back_inserter ( ordered ) ) ;
// Reorder auxiliary equations accordingly
auto aux_equations_old = aux_equations ;
auto index = get ( vertex_index , g ) ; // Maps vertex descriptors to their index
for ( int i = 0 ; i < n ; i + + )
aux_equations [ i ] = aux_equations_old [ index [ ordered [ i ] ] ] ;
}