2008-02-03 11:28:36 +01:00
/*
2018-03-27 17:14:30 +02:00
* Copyright ( C ) 2003 - 2018 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/>.
*/
2009-04-30 15:14:33 +02:00
# include <cstdlib>
2009-04-28 18:21:39 +02:00
# include <cassert>
2010-01-28 15:18:40 +01:00
# include <cmath>
2008-02-03 11:28:36 +01:00
# include <iostream>
2009-07-06 12:36:36 +02:00
# include <fstream>
2008-02-03 11:28:36 +01:00
# include "ModelTree.hh"
2009-12-16 14:21:31 +01:00
# include "MinimumFeedbackSet.hh"
# 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>
using namespace MFS ;
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
2018-06-04 12:26:16 +02:00
for ( const auto & it : contemporaneous_jacobian )
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
fill ( mate_map . begin ( ) , mate_map . end ( ) , graph_traits < BipartiteGraph > : : null_vertex ( ) ) ;
multimap < int , int > natural_endo2eqs ;
computeNormalizedEquations ( natural_endo2eqs ) ;
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 ) ;
}
edmonds_augmenting_path_finder < BipartiteGraph , size_t * , property_map < BipartiteGraph , vertex_index_t > : : type > augmentor ( g , & mate_map [ 0 ] , get ( vertex_index , g ) ) ;
bool not_maximum_yet = true ;
2009-12-16 18:13:23 +01:00
while ( not_maximum_yet )
2009-12-16 14:21:31 +01:00
{
not_maximum_yet = augmentor . augment_matching ( ) ;
}
augmentor . get_current_matching ( & mate_map [ 0 ] ) ;
bool check = maximum_cardinality_matching_verifier < BipartiteGraph , size_t * , property_map < BipartiteGraph , vertex_index_t > : : type > : : verify_matching ( g , & mate_map [ 0 ] , get ( vertex_index , g ) ) ;
# 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 ( ) ) ;
2009-12-16 14:21:31 +01:00
transform ( mate_map . begin ( ) , mate_map . begin ( ) + n , endo2eq . begin ( ) , bind2nd ( minus < int > ( ) , n ) ) ;
# ifdef DEBUG
multimap < int , int > natural_endo2eqs ;
computeNormalizedEquations ( natural_endo2eqs ) ;
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 + + ;
pair < multimap < int , int > : : const_iterator , multimap < int , int > : : const_iterator > x = natural_endo2eqs . equal_range ( i ) ;
if ( find_if ( x . first , x . second , compose1 ( bind2nd ( equal_to < int > ( ) , endo2eq [ i ] ) , select2nd < multimap < int , int > : : value_type > ( ) ) ) = = x . second )
cout < < " Natural normalization of variable " < < symbol_table . getName ( symbol_table . getID ( eEndogenous , i ) )
< < " 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
2018-06-05 14:56:34 +02:00
auto it = find ( mate_map . begin ( ) , mate_map . begin ( ) + n , boost : : graph_traits < BipartiteGraph > : : null_vertex ( ) ) ;
2009-12-16 14:21:31 +01:00
if ( 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 ) ;
2010-09-16 19:00:48 +02:00
for ( jacob_map_t : : const_iterator iter = contemporaneous_jacobian . begin ( ) ; iter ! = contemporaneous_jacobian . end ( ) ; iter + + )
2009-12-21 11:29:21 +01:00
if ( fabs ( iter - > second ) > max_val [ iter - > first . first ] )
max_val [ iter - > first . first ] = fabs ( iter - > second ) ;
2009-12-16 14:21:31 +01:00
2018-06-04 12:26:16 +02:00
for ( auto & iter : normalized_contemporaneous_jacobian )
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 ;
2018-06-04 12:26:16 +02:00
for ( auto & iter : normalized_contemporaneous_jacobian )
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 ) ;
2018-06-04 12:26:16 +02: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
2010-09-16 19:00:48 +02:00
for ( jacob_map_t : : const_iterator it = tmp_normalized_contemporaneous_jacobian . begin ( ) ; it ! = tmp_normalized_contemporaneous_jacobian . end ( ) ; it + + )
2009-12-21 11:29:21 +01:00
{
2018-06-04 16:36:46 +02:00
if ( static_jacobian . find ( { it - > first . first , it - > first . second } ) = = static_jacobian . end ( ) )
static_jacobian [ { it - > first . first , it - > first . second } ] = 0 ;
if ( dynamic_jacobian . find ( { 0 , { it - > first . first , it - > first . second } } ) = = dynamic_jacobian . end ( ) )
dynamic_jacobian [ { 0 , { it - > first . first , it - > first . second } } ] = nullptr ;
if ( contemporaneous_jacobian . find ( { it - > first . first , it - > first . second } ) = = contemporaneous_jacobian . end ( ) )
contemporaneous_jacobian [ { it - > first . first , it - > first . second } ] = 0 ;
2014-07-15 17:32:12 +02:00
try
{
2018-07-17 18:34:07 +02:00
if ( first_derivatives . find ( { it - > first . first , getDerivID ( symbol_table . getID ( SymbolType : : endogenous , it - > first . second ) , 0 ) } ) = = first_derivatives . end ( ) )
first_derivatives [ { it - > first . first , getDerivID ( symbol_table . getID ( SymbolType : : endogenous , it - > first . 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
{
2018-07-17 18:34:07 +02:00
cerr < < " The variable " < < symbol_table . getName ( symbol_table . getID ( SymbolType : : endogenous , it - > first . 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
}
void
ModelTree : : computeNormalizedEquations ( multimap < int , int > & endo2eqs ) const
{
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
{
2018-06-04 15:03:26 +02:00
auto * lhs = dynamic_cast < VariableNode * > ( equations [ i ] - > get_arg1 ( ) ) ;
2018-06-04 12:52:14 +02:00
if ( lhs = = nullptr )
2009-12-16 14:21:31 +01:00
continue ;
int symb_id = lhs - > get_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 ;
2009-12-16 14:21:31 +01:00
equations [ i ] - > get_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 ;
2018-06-04 16:36:46 +02:00
endo2eqs . emplace ( symbol_table . getTypeSpecificID ( symb_id ) , ( int ) i ) ;
2009-12-16 14:21:31 +01:00
cout < < " Endogenous " < < symbol_table . getName ( symb_id ) < < " normalized in equation " < < ( i + 1 ) < < endl ;
}
}
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-06-04 14:17:36 +02:00
set < pair < int , int > > jacobian_elements_to_delete ;
2010-09-16 19:00:48 +02:00
for ( first_derivatives_t : : const_iterator it = first_derivatives . begin ( ) ;
2009-12-16 14:21:31 +01:00
it ! = first_derivatives . end ( ) ; it + + )
{
int deriv_id = it - > first . second ;
2018-07-17 18:34:07 +02:00
if ( getTypeByDerivID ( deriv_id ) = = SymbolType : : endogenous )
2009-12-16 14:21:31 +01:00
{
2010-09-16 19:18:45 +02:00
expr_t Id = it - > second ;
2009-12-16 14:21:31 +01:00
int eq = it - > first . first ;
int symb = getSymbIDByDerivID ( deriv_id ) ;
int var = symbol_table . getTypeSpecificID ( symb ) ;
int lag = getLagByDerivID ( deriv_id ) ;
double val = 0 ;
try
{
val = Id - > eval ( eval_context ) ;
}
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 ;
2018-09-05 18:27:13 +02:00
Id - > 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 ;
2018-06-04 16:36:46 +02:00
jacobian_elements_to_delete . emplace ( 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 ;
dynamic_jacobian [ { lag , { eq , var } } ] = Id ;
2009-12-16 14:21:31 +01:00
}
}
}
// Get rid of the elements of the Jacobian matrix below the cutoff
2018-06-04 12:26:16 +02:00
for ( const auto & it : jacobian_elements_to_delete )
first_derivatives . 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
{
cout < < jacobian_elements_to_delete . size ( ) < < " elements among " < < first_derivatives . size ( ) < < " in the incidence matrices are below the cutoff ( " < < cutoff < < " ) and are discarded " < < endl
< < " The contemporaneous incidence matrix has " < < nb_elements_contemparenous_Jacobian < < " elements " < < endl ;
}
}
2018-09-21 17:13:19 +02:00
vector < pair < int , int > >
ModelTree : : select_non_linear_equations_and_variables ( vector < bool > is_equation_linear , const dynamic_jacob_map_t & dynamic_jacobian , vector < int > & equation_reordered , vector < int > & variable_reordered ,
vector < int > & inv_equation_reordered , vector < int > & inv_variable_reordered ,
lag_lead_vector_t & equation_lag_lead , lag_lead_vector_t & variable_lag_lead ,
vector < unsigned int > & n_static , vector < unsigned int > & n_forward , vector < unsigned int > & n_backward , vector < unsigned int > & n_mixed )
{
vector < int > eq2endo ( equations . size ( ) , 0 ) ;
/*equation_reordered.resize(equations.size());
variable_reordered . resize ( equations . size ( ) ) ; */
unsigned int num = 0 ;
for ( vector < int > : : const_iterator it = endo2eq . begin ( ) ; it ! = endo2eq . end ( ) ; it + + )
if ( ! is_equation_linear [ * it ] )
num + + ;
vector < int > endo2block = vector < int > ( endo2eq . size ( ) , 1 ) ;
vector < pair < set < int > , pair < set < int > , vector < int > > > > components_set ( num ) ;
int i = 0 , j = 0 ;
for ( vector < int > : : const_iterator it = endo2eq . begin ( ) ; it ! = endo2eq . end ( ) ; it + + , j + + )
{
if ( ! is_equation_linear [ * it ] )
{
equation_reordered [ i ] = * it ;
variable_reordered [ i ] = j ;
endo2block [ j ] = 0 ;
components_set [ endo2block [ j ] ] . first . insert ( i ) ;
/*cout << " -----------------------------------------" << endl;
cout . flush ( ) ;
cout < < " equation_reordered[ " < < * it < < " ] = " < < i < < endl ;
cout . flush ( ) ;
cout < < " variable_reordered[ " < < j < < " ] = " < < i < < endl ;
cout . flush ( ) ;
cout < < " components_set[ " < < endo2block [ j ] < < " ].first[ " < < i < < " ] = " < < i < < endl ;
cout < < " endo2block[ " < < j < < " ] = " < < 0 < < endl ;
cout . flush ( ) ;
*/
i + + ;
}
}
/* for (unsigned int j = 0; j < is_equation_linear.size() ; j++)
cout < < " endo2block[ " < < j < < " ] = " < < endo2block [ j ] < < endl ; */
/*cout << "before getVariableLeadLAgByBlock\n";
cout . flush ( ) ; */
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 + + )
{
if ( variable_lag_lead [ variable_reordered [ i ] ] . first ! = 0 & & variable_lag_lead [ variable_reordered [ i ] ] . second ! = 0 )
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 ( ) ;
vector < pair < int , int > > blocks = vector < pair < int , int > > ( 1 , make_pair ( i , i ) ) ;
inv_equation_reordered = vector < int > ( nb_endo ) ;
inv_variable_reordered = vector < int > ( nb_endo ) ;
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 ;
set < pair < int , int > > result ;
endo2eq . resize ( equations . size ( ) ) ;
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
if ( ! is_equation_linear [ eq ] )
{
BinaryOpNode * eq_node = equations [ eq ] ;
expr_t lhs = eq_node - > get_arg1 ( ) ;
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 !
vector < int > : : iterator it = find ( endo2eq . begin ( ) , endo2eq . end ( ) , result . begin ( ) - > first ) ;
if ( it = = endo2eq . end ( ) )
endo2eq [ result . begin ( ) - > first ] = eq ;
else
{
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 ( ) ) ;
2009-12-16 14:21:31 +01:00
bool * IM ;
2017-02-24 11:20:54 +01:00
int n = equations . size ( ) ;
2009-12-16 18:13:23 +01:00
IM = ( bool * ) calloc ( n * n , sizeof ( bool ) ) ;
2009-12-16 14:21:31 +01:00
int i = 0 ;
2009-12-16 18:13:23 +01:00
for ( vector < int > : : const_iterator it = endo2eq . begin ( ) ; it ! = endo2eq . end ( ) ; it + + , i + + )
2009-12-16 14:21:31 +01:00
{
eq2endo [ * it ] = i ;
equation_reordered [ i ] = i ;
variable_reordered [ * it ] = 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 ) ;
2018-06-04 12:26:16 +02:00
for ( const auto & it : endo )
IM [ i * n + endo2eq [ it . first ] ] = true ;
2017-06-14 07:01:31 +02:00
}
}
else
2018-06-04 12:26:16 +02:00
for ( const auto & it : static_jacobian_arg )
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 ;
2009-12-16 18:13:23 +01:00
for ( int i = prologue ; i < n - ( 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 ;
}
free ( IM ) ;
}
2010-09-16 19:00:48 +02:00
equation_type_and_normalized_equation_t
2018-06-04 14:17:36 +02:00
ModelTree : : equationTypeDetermination ( const map < pair < int , pair < 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 ] ;
lhs = eq_node - > get_arg1 ( ) ;
Equation_Simulation_Type = E_SOLVE ;
2018-06-04 16:36:46 +02: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 ) )
2009-12-16 14:21:31 +01:00
{
Equation_Simulation_Type = E_EVALUATE ;
}
else
{
2018-06-04 14:17:36 +02:00
vector < pair < int , pair < expr_t , expr_t > > > List_of_Op_RHS ;
2009-12-16 14:21:31 +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
}
return ( V_Equation_Simulation_Type ) ;
}
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 + + )
{
2011-01-13 19:03:55 +01:00
if ( i < ( int ) prologue )
2009-12-16 14:21:31 +01:00
{
variable_blck [ variable_reordered [ i ] ] = i ;
equation_blck [ equation_reordered [ i ] ] = i ;
}
2011-01-13 19:03:55 +01:00
else if ( i < ( 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 ) ;
}
}
2018-06-04 12:26:16 +02:00
for ( const auto & it : dynamic_jacobian )
2009-12-16 14:21:31 +01:00
{
2018-06-04 12:26:16 +02:00
int lag = it . first . first ;
int j_1 = it . first . second . first ;
int i_1 = it . first . second . second ;
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 ) ;
2018-06-04 12:26:16 +02: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 ;
2016-03-25 15:38:49 +01:00
for ( jacob_map_t : : const_iterator it = tmp_normalized_contemporaneous_jacobian . begin ( ) ; it ! = tmp_normalized_contemporaneous_jacobian . end ( ) ; it + + )
2011-01-13 19:03:55 +01:00
if ( reverse_equation_reordered [ it - > first . first ] > = ( int ) prologue & & reverse_equation_reordered [ it - > first . first ] < ( int ) ( nb_var - epilogue )
& & reverse_variable_reordered [ it - > first . second ] > = ( int ) prologue & & reverse_variable_reordered [ it - > first . second ] < ( int ) ( nb_var - epilogue )
2009-12-16 18:13:23 +01:00
& & it - > first . first ! = endo2eq [ it - > first . second ] )
2010-12-13 14:23:04 +01:00
add_edge ( vertex ( reverse_equation_reordered [ endo2eq [ it - > first . second ] ] - prologue , G2 ) ,
vertex ( reverse_equation_reordered [ it - > first . first ] - prologue , G2 ) ,
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-06-04 14:17:36 +02:00
vector < pair < set < int > , pair < 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 + + ;
components_set [ endo2block [ i ] ] . first . insert ( i ) ;
}
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
{
for ( int i = 0 ; i < n ; i + + )
if ( Equation_Type [ equation_reordered [ i + prologue ] ] . first = = E_SOLVE | | 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
}
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 ) ;
2011-01-13 19:03:55 +01:00
for ( int i = 0 ; i < ( int ) prologue ; i + + )
2010-07-23 11:20:24 +02: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 + + )
{
2010-12-13 14:23:04 +01:00
AdjacencyList_t G = extract_subgraph ( G2 , components_set [ i ] . first ) ;
2009-12-16 14:21:31 +01:00
set < int > feed_back_vertices ;
//Print(G);
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 ) ;
2009-12-16 14:21:31 +01:00
components_set [ i ] . second . first = feed_back_vertices ;
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 + + )
2009-12-16 14:21:31 +01:00
{
2018-06-04 12:26:16 +02:00
for ( int & its : Reordered_Vertice )
2010-07-23 11:20:24 +02:00
{
bool something_done = false ;
2018-06-04 12:26:16 +02:00
if ( j = = 2 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . second ! = 0 )
2010-07-23 11:20:24 +02:00
{
n_mixed [ prologue + i ] + + ;
something_done = true ;
}
2018-06-04 12:26:16 +02:00
else if ( j = = 3 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . second ! = 0 )
2010-07-23 11:20:24 +02:00
{
n_forward [ prologue + i ] + + ;
something_done = true ;
}
2018-06-04 12:26:16 +02:00
else if ( j = = 1 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . second = = 0 )
2010-07-23 11:20:24 +02:00
{
n_backward [ prologue + i ] + + ;
something_done = true ;
}
2018-06-04 12:26:16 +02:00
else if ( j = = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ its + prologue ] ] . second = = 0 )
2010-07-23 11:20:24 +02:00
{
n_static [ prologue + i ] + + ;
something_done = true ;
}
if ( something_done )
{
2018-06-04 12:26:16 +02:00
equation_reordered [ order ] = tmp_equation_reordered [ its + prologue ] ;
variable_reordered [ order ] = tmp_variable_reordered [ its + prologue ] ;
2010-07-23 11:20:24 +02:00
order + + ;
}
}
2009-12-16 14:21:31 +01:00
}
components_set [ i ] . second . second = Reordered_Vertice ;
//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 + + )
2009-12-16 14:21:31 +01:00
{
2018-06-04 12:26:16 +02:00
for ( int feed_back_vertice : feed_back_vertices )
2010-07-23 11:20:24 +02:00
{
bool something_done = false ;
2018-06-04 12:26:16 +02:00
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 )
2010-07-23 11:20:24 +02:00
{
n_mixed [ prologue + i ] + + ;
something_done = true ;
}
2018-06-04 12:26:16 +02:00
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 )
2010-07-23 11:20:24 +02:00
{
n_forward [ prologue + i ] + + ;
something_done = true ;
}
2018-06-04 12:26:16 +02:00
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 )
2010-07-23 11:20:24 +02:00
{
n_backward [ prologue + i ] + + ;
something_done = true ;
}
2018-06-04 12:26:16 +02:00
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 )
2010-07-23 11:20:24 +02:00
{
n_static [ prologue + i ] + + ;
something_done = true ;
}
if ( something_done )
{
2018-06-04 12:26:16 +02:00
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 ] ;
2010-07-23 11:20:24 +02:00
order + + ;
}
}
2009-12-16 14:21:31 +01:00
}
}
2010-07-23 11:20:24 +02:00
2011-01-13 19:03:55 +01:00
for ( int i = 0 ; i < ( int ) epilogue ; i + + )
2010-07-23 11:20:24 +02:00
{
2010-10-27 15:34:48 +02:00
if ( variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . second ! = 0 )
2010-07-23 11:20:24 +02:00
n_mixed [ prologue + num + i ] + + ;
2010-10-27 15:34:48 +02:00
else if ( variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . second ! = 0 )
2010-07-23 11:20:24 +02:00
n_forward [ prologue + num + i ] + + ;
2010-10-27 15:34:48 +02:00
else if ( variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . second = = 0 )
2010-07-23 11:20:24 +02:00
n_backward [ prologue + num + i ] + + ;
2010-10-27 15:34:48 +02:00
else if ( variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ prologue + n + i ] ] . second = = 0 )
2010-07-23 11:20:24 +02:00
n_static [ prologue + num + i ] + + ;
}
2009-12-16 14:21:31 +01:00
inv_equation_reordered = vector < int > ( nb_var ) ;
inv_variable_reordered = vector < int > ( 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
{
int largest_block = 0 ;
int Nb_SimulBlocks = 0 ;
int Nb_feedback_variable = 0 ;
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-09-21 17:13:19 +02: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 < pair < pair < int , int > , pair < 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 ;
2010-07-23 11:20:24 +02:00
unsigned int l_n_static = 0 ;
unsigned int l_n_forward = 0 ;
unsigned int l_n_backward = 0 ;
unsigned int l_n_mixed = 0 ;
2011-01-13 19:03:55 +01:00
for ( i = 0 ; i < ( int ) ( prologue + blocks . size ( ) + epilogue ) ; i + + )
2009-12-16 14:21:31 +01:00
{
int first_count_equ = count_equ ;
2011-01-13 19:03:55 +01:00
if ( i < ( int ) prologue )
2009-12-16 14:21:31 +01:00
{
Blck_Size = 1 ;
MFS_Size = 1 ;
}
2011-01-13 19:03:55 +01:00
else if ( i < ( 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 + + ;
}
2011-01-13 19:03:55 +01:00
else if ( i < ( 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 ;
2009-12-16 18:13:23 +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 ) ;
2018-06-04 12:26:16 +02: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-06-04 15:03:26 +02:00
auto it1 = find ( variable_reordered . begin ( ) + first_count_equ , variable_reordered . begin ( ) + ( first_count_equ + Blck_Size ) , curr_variable ) ;
2018-09-21 17:13:19 +02:00
if ( linear_decomposition )
{
if ( dynamic_jacobian . find ( make_pair ( curr_lag , make_pair ( equation_reordered [ count_equ ] , curr_variable ) ) ) ! = dynamic_jacobian . end ( ) )
{
if ( curr_lag > Lead )
Lead = curr_lag ;
else if ( - curr_lag > Lag )
Lag = - curr_lag ;
}
}
else
{
if ( it1 ! = variable_reordered . begin ( ) + ( first_count_equ + Blck_Size ) )
if ( 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 ;
}
}
2009-12-16 14:21:31 +01:00
}
}
if ( ( Lag > 0 ) & & ( Lead > 0 ) )
{
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 ;
int c_Size = ( block_type_size_mfs [ block_type_size_mfs . size ( ) - 1 ] ) . second . first ;
int first_equation = ( block_type_size_mfs [ block_type_size_mfs . size ( ) - 1 ] ) . first . second ;
if ( c_Size > 0 & & ( ( prev_Type = = EVALUATE_FORWARD & & Simulation_Type = = EVALUATE_FORWARD & & ! is_lead )
2017-06-14 07:01:31 +02:00
| | ( 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-06-04 16:36:46 +02: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-06-04 16:36:46 +02: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 ;
}
}
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
BlockSimulationType c_Type = ( block_type_size_mfs [ block_type_size_mfs . size ( ) - 1 ] ) . first . first ;
2010-07-23 11:20:24 +02:00
c_Size + + ;
2018-06-04 16:36:46 +02: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-06-04 14:17:36 +02:00
pair < pair < unsigned int , unsigned int > , pair < unsigned int , unsigned int > > tmp = block_col_type [ block_col_type . size ( ) - 1 ] ;
2018-06-04 16:36:46 +02:00
block_col_type [ block_col_type . size ( ) - 1 ] = { { tmp . first . first + l_n_static , tmp . first . second + l_n_forward } , { tmp . second . first + l_n_backward , tmp . second . second + l_n_mixed } } ;
2009-12-16 14:21:31 +01:00
}
else
{
2018-06-04 16:36:46 +02:00
block_type_size_mfs . emplace_back ( make_pair ( Simulation_Type , eq ) , make_pair ( Blck_Size , MFS_Size ) ) ;
block_lag_lead . emplace_back ( Lag , Lead ) ;
2018-06-04 12:48:09 +02:00
block_col_type . emplace_back ( make_pair ( l_n_static , l_n_forward ) , make_pair ( l_n_backward , l_n_mixed ) ) ;
2009-12-16 14:21:31 +01:00
}
}
else
{
2018-06-04 16:36:46 +02:00
block_type_size_mfs . emplace_back ( make_pair ( Simulation_Type , eq ) , make_pair ( Blck_Size , MFS_Size ) ) ;
block_lag_lead . emplace_back ( Lag , Lead ) ;
2018-06-04 12:48:09 +02:00
block_col_type . emplace_back ( make_pair ( l_n_static , l_n_forward ) , make_pair ( l_n_backward , l_n_mixed ) ) ;
2009-12-16 14:21:31 +01:00
}
}
else
{
2018-06-04 16:36:46 +02:00
block_type_size_mfs . emplace_back ( make_pair ( Simulation_Type , eq ) , make_pair ( Blck_Size , MFS_Size ) ) ;
block_lag_lead . emplace_back ( Lag , Lead ) ;
2018-06-04 12:48:09 +02:00
block_col_type . emplace_back ( make_pair ( l_n_static , l_n_forward ) , make_pair ( l_n_backward , l_n_mixed ) ) ;
2009-12-16 14:21:31 +01:00
}
prev_Type = Simulation_Type ;
eq + = Blck_Size ;
}
return ( block_type_size_mfs ) ;
}
2018-09-21 17:13:19 +02:00
vector < bool >
ModelTree : : equationLinear ( map < pair < int , pair < int , int > > , expr_t > first_order_endo_derivatives ) const
{
vector < bool > is_linear ( symbol_table . endo_nbr ( ) , true ) ;
for ( map < pair < int , pair < int , int > > , expr_t > : : const_iterator it = first_order_endo_derivatives . begin ( ) ; it ! = first_order_endo_derivatives . end ( ) ; it + + )
{
expr_t Id = it - > second ;
set < pair < int , int > > endogenous ;
Id - > collectEndogenous ( endogenous ) ;
if ( endogenous . size ( ) > 0 )
{
int eq = it - > first . first ;
is_linear [ eq ] = false ;
}
}
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 )
2009-12-16 14:21:31 +01:00
{
2010-09-16 19:00:48 +02:00
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = derivatives_block . begin ( ) ; it ! = derivatives_block . end ( ) ; it + + )
2009-12-16 14:21:31 +01:00
{
int lag = it - > second . first ;
if ( lag = = 0 )
{
2010-09-16 19:18:45 +02:00
expr_t Id = it - > second . second ;
2018-06-04 14:17:36 +02:00
set < pair < int , int > > endogenous ;
2009-12-16 14:21:31 +01:00
Id - > collectEndogenous ( endogenous ) ;
if ( endogenous . size ( ) > 0 )
{
2009-12-16 18:13:23 +01:00
for ( int l = 0 ; l < block_size ; l + + )
2009-12-16 14:21:31 +01:00
{
2018-06-04 16:36:46 +02:00
if ( endogenous . find ( { variable_reordered [ first_variable_position + l ] , 0 } ) ! = endogenous . end ( ) )
2009-12-16 14:21:31 +01:00
{
blocks_linear [ block ] = false ;
goto the_end ;
}
}
}
}
}
}
2009-12-16 18:13:23 +01:00
else if ( simulation_type = = SOLVE_TWO_BOUNDARIES_COMPLETE | | simulation_type = = SOLVE_TWO_BOUNDARIES_SIMPLE )
2009-12-16 14:21:31 +01:00
{
2010-09-16 19:00:48 +02:00
for ( block_derivatives_equation_variable_laglead_nodeid_t : : const_iterator it = derivatives_block . begin ( ) ; it ! = derivatives_block . end ( ) ; it + + )
2009-12-16 14:21:31 +01:00
{
int lag = it - > second . first ;
2010-09-16 19:18:45 +02:00
expr_t Id = it - > second . second ; //
2018-06-04 14:17:36 +02:00
set < pair < int , int > > endogenous ;
2009-12-16 14:21:31 +01:00
Id - > collectEndogenous ( endogenous ) ;
if ( endogenous . size ( ) > 0 )
{
2009-12-16 18:13:23 +01:00
for ( int l = 0 ; l < block_size ; l + + )
2009-12-16 14:21:31 +01:00
{
2018-06-04 16:36:46 +02:00
if ( endogenous . find ( { variable_reordered [ first_variable_position + l ] , lag } ) ! = endogenous . end ( ) )
2009-12-16 14:21:31 +01:00
{
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
}
2009-12-16 18:13:23 +01:00
return ( blocks_linear ) ;
2009-12-16 14:21:31 +01:00
}
2008-02-03 11:28:36 +01:00
ModelTree : : ModelTree ( SymbolTable & symbol_table_arg ,
2010-02-22 17:33:38 +01:00
NumericalConstants & num_constants_arg ,
2018-09-14 17:04:06 +02:00
ExternalFunctionsTable & external_functions_table_arg ) :
DataTree ( symbol_table_arg , num_constants_arg , external_functions_table_arg ) ,
2011-06-22 11:56:07 +02:00
cutoff ( 1e-15 ) ,
mfs ( 0 )
2008-02-03 11:28:36 +01:00
{
2018-06-04 12:26:16 +02:00
for ( int & NNZDerivative : NNZDerivatives )
NNZDerivative = 0 ;
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
{
2018-06-04 16:36:46 +02:00
auto it = first_derivatives . find ( { eq , getDerivID ( symb_id , lag ) } ) ;
2009-01-23 11:59:37 +01:00
if ( it ! = first_derivatives . end ( ) )
2018-05-29 11:12:22 +02:00
( 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
2009-04-20 12:48:54 +02:00
ModelTree : : computeJacobian ( const set < int > & vars )
2008-02-03 11:28:36 +01:00
{
2018-06-04 12:26:16 +02:00
for ( int var : vars )
2013-03-22 16:34:50 +01:00
{
2013-03-26 15:33:28 +01:00
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
{
2018-06-04 12:26:16 +02:00
expr_t d1 = equations [ eq ] - > getDerivative ( var ) ;
2013-03-26 15:33:28 +01:00
if ( d1 = = Zero )
continue ;
2018-06-04 16:36:46 +02:00
first_derivatives [ { eq , var } ] = d1 ;
2013-03-26 15:33:28 +01:00
+ + NNZDerivatives [ 0 ] ;
2016-03-25 15:38:49 +01:00
}
2013-03-22 16:34:50 +01:00
}
2009-04-20 12:48:54 +02:00
}
2008-02-03 11:28:36 +01:00
2009-04-20 12:48:54 +02:00
void
ModelTree : : computeHessian ( const set < int > & vars )
{
2018-06-05 15:34:34 +02:00
for ( const auto & it : first_derivatives )
2008-02-03 11:28:36 +01:00
{
2018-06-05 15:34:34 +02:00
int eq , var1 ;
tie ( eq , var1 ) = it . first ;
expr_t d1 = it . second ;
2009-04-20 12:48:54 +02:00
// Store only second derivatives with var2 <= var1
2018-06-04 12:26:16 +02:00
for ( int var2 : vars )
2008-02-03 11:28:36 +01:00
{
2009-04-20 12:48:54 +02:00
if ( var2 > var1 )
continue ;
2010-09-16 19:18:45 +02:00
expr_t d2 = d1 - > getDerivative ( var2 ) ;
2009-04-20 12:48:54 +02:00
if ( d2 = = Zero )
continue ;
2018-06-05 15:34:34 +02:00
second_derivatives [ { eq , var1 , var2 } ] = d2 ;
2009-12-16 18:13:23 +01:00
if ( var2 = = var1 )
+ + NNZDerivatives [ 1 ] ;
else
NNZDerivatives [ 1 ] + = 2 ;
2008-02-03 11:28:36 +01:00
}
}
2009-04-20 12:48:54 +02:00
}
2008-02-03 11:28:36 +01:00
2009-04-20 12:48:54 +02:00
void
ModelTree : : computeThirdDerivatives ( const set < int > & vars )
{
2018-06-05 15:34:34 +02:00
for ( const auto & it : second_derivatives )
2008-02-03 11:28:36 +01:00
{
2018-06-05 15:34:34 +02:00
int eq , var1 , var2 ;
tie ( eq , var1 , var2 ) = it . first ;
2009-04-20 12:48:54 +02:00
// By construction, var2 <= var1
2018-06-05 15:34:34 +02:00
expr_t d2 = it . second ;
2009-04-20 12:48:54 +02:00
// Store only third derivatives such that var3 <= var2 <= var1
2018-06-04 12:26:16 +02:00
for ( int var3 : vars )
2008-02-03 11:28:36 +01:00
{
2009-04-20 12:48:54 +02:00
if ( var3 > var2 )
continue ;
2010-09-16 19:18:45 +02:00
expr_t d3 = d2 - > getDerivative ( var3 ) ;
2009-04-20 12:48:54 +02:00
if ( d3 = = Zero )
continue ;
2018-06-05 15:34:34 +02:00
third_derivatives [ { eq , var1 , var2 , var3 } ] = d3 ;
2009-12-16 18:13:23 +01:00
if ( var3 = = var2 & & var2 = = var1 )
+ + NNZDerivatives [ 2 ] ;
else if ( var3 = = var2 | | var2 = = var1 )
NNZDerivatives [ 2 ] + = 3 ;
else
NNZDerivatives [ 2 ] + = 6 ;
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-06-04 14:17:36 +02:00
map < expr_t , pair < int , NodeTreeReference > > reference_count ;
2008-02-03 11:28:36 +01:00
temporary_terms . clear ( ) ;
2018-03-27 17:14:30 +02:00
temporary_terms_mlv . clear ( ) ;
2015-09-01 17:39:49 +02:00
temporary_terms_res . clear ( ) ;
2015-09-02 18:50:09 +02:00
temporary_terms_g1 . clear ( ) ;
temporary_terms_g2 . clear ( ) ;
temporary_terms_g3 . clear ( ) ;
2018-03-27 17:14:30 +02:00
// Collect all model local variables appearing in equations. See #101
// All used model local variables are automatically set as temporary variables
2018-05-23 19:19:11 +02:00
set < int > used_local_vars ;
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 ) ;
2018-03-27 17:14:30 +02:00
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-07-18 16:35:19 +02:00
reference_count [ v ] = { ExprNode : : min_cost ( is_matlab ) + 1 , NodeTreeReference : : residuals } ;
2018-03-27 17:14:30 +02:00
}
2018-09-25 15:56:52 +02:00
/* When option notmpterms is set, we only need to process model local
variables ( and turn them into temporary terms ) ; no need to go further */
if ( no_tmp_terms )
return ;
2015-09-03 13:50:02 +02:00
map < NodeTreeReference , temporary_terms_t > temp_terms_map ;
2018-07-18 16:35:19 +02:00
temp_terms_map [ NodeTreeReference : : residuals ] = temporary_terms_res ;
temp_terms_map [ NodeTreeReference : : firstDeriv ] = temporary_terms_g1 ;
temp_terms_map [ NodeTreeReference : : secondDeriv ] = temporary_terms_g2 ;
temp_terms_map [ NodeTreeReference : : thirdDeriv ] = temporary_terms_g3 ;
2008-02-03 11:28:36 +01:00
2018-06-04 12:26:16 +02:00
for ( auto & equation : equations )
equation - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
is_matlab , NodeTreeReference : : residuals ) ;
2008-02-03 11:28:36 +01:00
2018-06-04 12:26:16 +02:00
for ( auto & first_derivative : first_derivatives )
first_derivative . second - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
is_matlab , NodeTreeReference : : firstDeriv ) ;
2008-02-03 11:28:36 +01:00
2018-06-04 12:26:16 +02:00
for ( auto & second_derivative : second_derivatives )
second_derivative . second - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
is_matlab , NodeTreeReference : : secondDeriv ) ;
2008-02-03 11:28:36 +01:00
2018-06-04 12:26:16 +02:00
for ( auto & third_derivative : third_derivatives )
third_derivative . second - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
is_matlab , NodeTreeReference : : thirdDeriv ) ;
2015-09-03 13:50:02 +02:00
for ( map < NodeTreeReference , temporary_terms_t > : : const_iterator it = temp_terms_map . begin ( ) ;
it ! = temp_terms_map . end ( ) ; it + + )
2018-03-27 17:14:30 +02:00
temporary_terms . insert ( it - > second . begin ( ) , it - > second . end ( ) ) ;
2015-09-03 13:50:02 +02:00
2018-07-18 16:35:19 +02:00
temporary_terms_res = temp_terms_map [ NodeTreeReference : : residuals ] ;
temporary_terms_g1 = temp_terms_map [ NodeTreeReference : : firstDeriv ] ;
temporary_terms_g2 = temp_terms_map [ NodeTreeReference : : secondDeriv ] ;
temporary_terms_g3 = temp_terms_map [ NodeTreeReference : : thirdDeriv ] ;
2018-03-27 17:14:30 +02:00
2018-05-25 15:19:33 +02:00
int idx = 0 ;
2018-05-28 15:50:29 +02:00
for ( map < expr_t , expr_t , ExprNodeLess > : : const_iterator it = temporary_terms_mlv . begin ( ) ;
2018-03-27 17:14:30 +02:00
it ! = temporary_terms_mlv . end ( ) ; it + + )
2018-05-24 19:29:53 +02:00
temporary_terms_idxs [ it - > first ] = idx + + ;
2018-03-27 17:14:30 +02:00
2018-06-04 12:26:16 +02:00
for ( auto temporary_terms_re : temporary_terms_res )
temporary_terms_idxs [ temporary_terms_re ] = idx + + ;
2018-03-27 17:14:30 +02:00
2018-06-04 12:26:16 +02:00
for ( auto it : temporary_terms_g1 )
temporary_terms_idxs [ it ] = idx + + ;
2018-03-27 17:14:30 +02:00
2018-06-04 12:26:16 +02:00
for ( auto it : temporary_terms_g2 )
temporary_terms_idxs [ it ] = idx + + ;
2018-03-27 17:14:30 +02:00
2018-06-04 12:26:16 +02:00
for ( auto it : temporary_terms_g3 )
temporary_terms_idxs [ it ] = idx + + ;
2018-03-27 17:14:30 +02:00
}
void
2018-05-28 15:50:29 +02:00
ModelTree : : writeModelLocalVariableTemporaryTerms ( const temporary_terms_t & tto , const map < expr_t , expr_t , ExprNodeLess > & tt ,
2018-03-27 17:14:30 +02:00
ostream & output , ExprNodeOutputType output_type ,
deriv_node_temp_terms_t & tef_terms ) const
{
temporary_terms_t tt2 ;
2018-06-04 12:26:16 +02:00
for ( auto it : tt )
2018-03-27 17:14:30 +02:00
{
2018-09-05 18:27:13 +02:00
if ( isCOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " double " ;
2018-09-05 18:27:13 +02:00
else if ( isJuliaOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " @inbounds const " ;
2018-06-04 12:26:16 +02:00
it . first - > writeOutput ( output , output_type , tto , temporary_terms_idxs , tef_terms ) ;
2018-03-27 17:14:30 +02:00
output < < " = " ;
2018-06-04 12:26:16 +02:00
it . second - > writeOutput ( output , output_type , tt2 , temporary_terms_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 ;
// Insert current node into tt2
2018-06-04 12:26:16 +02:00
tt2 . 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 ,
const temporary_terms_t & ttm1 ,
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
{
// Local var used to keep track of temp nodes already written
temporary_terms_t tt2 = ttm1 ;
2018-06-04 15:03:26 +02:00
for ( auto it = tt . begin ( ) ;
2018-05-24 19:29:53 +02:00
it ! = tt . end ( ) ; it + + )
2018-03-27 17:14:30 +02:00
{
2018-06-04 12:52:14 +02:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( * it ) ! = nullptr )
2018-05-28 15:23:15 +02:00
( * it ) - > writeExternalFunctionOutput ( output , output_type , tt2 , tt_idxs , tef_terms ) ;
2018-03-27 17:14:30 +02:00
2018-09-05 18:27:13 +02:00
if ( isCOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " double " ;
2018-09-05 18:27:13 +02:00
else if ( isJuliaOutput ( output_type ) )
2018-03-27 17:14:30 +02:00
output < < " @inbounds " ;
2018-05-28 15:23:15 +02:00
( * it ) - > writeOutput ( output , output_type , tt , tt_idxs , tef_terms ) ;
2018-03-27 17:14:30 +02:00
output < < " = " ;
2018-05-28 15:23:15 +02:00
( * it ) - > writeOutput ( output , output_type , tt2 , 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 ;
// Insert current node into tt2
tt2 . insert ( * it ) ;
}
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
ModelTree : : writeJsonTemporaryTerms ( const temporary_terms_t & tt , const temporary_terms_t & ttm1 , ostream & output ,
deriv_node_temp_terms_t & tef_terms , string & concat ) const
{
// Local var used to keep track of temp nodes already written
bool wrote_term = false ;
temporary_terms_t tt2 = ttm1 ;
output < < " \" external_functions_temporary_terms_ " < < concat < < " \" : [ " ;
2018-06-04 12:26:16 +02:00
for ( auto it : tt )
if ( ttm1 . find ( it ) = = ttm1 . end ( ) )
2017-02-20 12:18:11 +01:00
{
2018-06-04 12:52:14 +02:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( it ) ! = nullptr )
2017-02-20 12:18:11 +01:00
{
if ( wrote_term )
output < < " , " ;
vector < string > efout ;
2018-06-04 12:26:16 +02:00
it - > writeJsonExternalFunctionOutput ( efout , tt2 , tef_terms ) ;
2017-02-20 12:18:11 +01:00
for ( vector < string > : : const_iterator it1 = efout . begin ( ) ; it1 ! = efout . end ( ) ; it1 + + )
{
if ( it1 ! = efout . begin ( ) )
output < < " , " ;
output < < * it1 ;
}
wrote_term = true ;
}
2018-06-04 12:26:16 +02:00
tt2 . insert ( it ) ;
2017-02-20 12:18:11 +01:00
}
tt2 = ttm1 ;
wrote_term = false ;
output < < " ] "
< < " , \" temporary_terms_ " < < concat < < " \" : [ " ;
2018-06-04 15:03:26 +02:00
for ( auto it = tt . begin ( ) ;
2017-02-20 12:18:11 +01:00
it ! = tt . end ( ) ; it + + )
if ( ttm1 . find ( * it ) = = ttm1 . end ( ) )
{
if ( wrote_term )
output < < " , " ;
output < < " { \" temporary_term \" : \" " ;
( * it ) - > writeJsonOutput ( output , tt , tef_terms ) ;
2017-02-27 12:23:39 +01:00
output < < " \" "
< < " , \" value \" : \" " ;
2017-02-20 12:18:11 +01:00
( * it ) - > writeJsonOutput ( output , tt2 , tef_terms ) ;
output < < " \" } " < < endl ;
wrote_term = true ;
// Insert current node into tt2
tt2 . insert ( * it ) ;
}
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 ;
map < string , string > : : iterator it ;
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 )
{
cerr < < " Warning: A .m file created by Dynare will have more than 32 nested parenthesis. Matlab cannot support this. " < < endl
< < " 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
< < " If you have not yet set up a compiler on your system, see the Matlab documentation for doing so. " < < endl
< < " 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 ) ;
string repstr = " " ;
string varname ;
while ( testNestedParenthesis ( str1 ) )
{
2017-01-02 11:37:15 +01:00
size_t open_paren_idx = string : : npos ;
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)
size_t idx = str1 . find_last_of ( " */-+ " , j - 1 ) ;
if ( j = = 0 | | ( idx ! = string : : npos & & idx = = j - 1 ) )
open_paren_idx = j ;
last_open_paren = j ;
}
else if ( str1 . at ( j ) = = ' ) ' )
{
// don't match, e.g. y(1)
size_t idx = str1 . find_last_not_of ( " 0123456789 " , j - 1 ) ;
if ( idx ! = string : : npos & & idx ! = last_open_paren )
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 ) ;
it = tmp_paren_vars . find ( val ) ;
if ( it = = tmp_paren_vars . end ( ) )
{
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 ;
}
}
}
it = tmp_paren_vars . find ( str1 ) ;
if ( it = = tmp_paren_vars . end ( ) )
{
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
{
2018-06-04 12:52:14 +02:00
if ( dynamic_cast < AbstractExternalFunctionNode * > ( it ) ! = nullptr )
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
}
2018-06-04 12:26:16 +02:00
FNUMEXPR_ fnumexpr ( TemporaryTerm , ( int ) ( map_idx . find ( it - > idx ) - > second ) ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
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 )
{
2018-06-04 12:26:16 +02:00
FSTPT_ fstpt ( ( 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
{
2018-06-04 12:26:16 +02:00
FSTPST_ fstpst ( ( 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
2017-11-06 15:21:11 +01:00
output < < " \" 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 ) ;
for ( vector < string > : : const_iterator it1 = efout . begin ( ) ; it1 ! = efout . end ( ) ; it1 + + )
{
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 ) */
output < < " { \" variable \" : \" " < < symbol_table . getName ( id ) < < " __ \" "
< < " , \" value \" : \" " ;
value - > writeJsonOutput ( output , tt , tef_terms ) ;
output < < " \" } " < < endl ;
}
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
{
2009-01-23 11:59:37 +01:00
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
{
BinaryOpNode * eq_node = equations [ eq ] ;
2010-09-16 19:18:45 +02:00
expr_t lhs = eq_node - > get_arg1 ( ) ;
expr_t rhs = eq_node - > get_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
{
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
{
BinaryOpNode * eq_node = equations [ eq ] ;
2010-09-16 19:18:45 +02:00
expr_t lhs = eq_node - > get_arg1 ( ) ;
expr_t rhs = eq_node - > get_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 ( ) )
{
2018-06-27 15:01:31 +02:00
cerr < < " Error : Can't open file \" " < < filename < < " \" for writing " < < endl ;
2010-01-22 11:03:29 +01:00
exit ( EXIT_FAILURE ) ;
}
u_count_int = 0 ;
2018-06-04 12:26:16 +02:00
for ( const auto & first_derivative : first_derivatives )
2010-01-22 11:03:29 +01:00
{
2018-06-04 12:26:16 +02:00
int deriv_id = first_derivative . first . second ;
2018-07-17 18:34:07 +02:00
if ( getTypeByDerivID ( deriv_id ) = = SymbolType : : endogenous )
2010-01-22 11:03:29 +01:00
{
2018-06-04 12:26:16 +02:00
int eq = first_derivative . first . first ;
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 )
u_count_int + = symbol_table . endo_nbr ( ) ;
for ( j = 0 ; j < ( int ) symbol_table . endo_nbr ( ) ; j + + )
SaveCode . write ( reinterpret_cast < char * > ( & j ) , sizeof ( j ) ) ;
for ( j = 0 ; j < ( int ) symbol_table . endo_nbr ( ) ; j + + )
SaveCode . write ( reinterpret_cast < char * > ( & j ) , sizeof ( j ) ) ;
SaveCode . close ( ) ;
}
2009-04-30 15:14:33 +02:00
void
2017-04-04 15:28:27 +02:00
ModelTree : : writeLatexModelFile ( const string & basename , ExprNodeOutputType output_type , const bool write_equation_tags ) const
2009-04-30 15:14:33 +02:00
{
2015-07-15 08:58:15 +02:00
ofstream output , content_output ;
string filename = basename + " .tex " ;
string content_basename = basename + " _content " ;
string content_filename = content_basename + " .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 ) ;
}
2009-04-30 15:14:33 +02:00
output < < " \\ documentclass[10pt,a4paper]{article} " < < endl
< < " \\ usepackage[landscape]{geometry} " < < endl
< < " \\ usepackage{fullpage} " < < endl
2015-02-16 08:31:30 +01:00
< < " \\ usepackage{amsfonts} " < < endl
2011-10-09 18:24:39 +02:00
< < " \\ usepackage{breqn} " < < endl
2009-04-30 15:14:33 +02:00
< < " \\ begin{document} " < < endl
< < " \\ footnotesize " < < endl ;
// 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
2015-07-15 08:58:15 +02:00
content_output < < " \\ begin{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 ) ;
content_output < < endl < < " \\ end{dmath*} " < < endl ;
2009-04-30 15:14:33 +02:00
}
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
{
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
{
2018-06-05 14:47:36 +02:00
bool wrote_eq_tag = false ;
2018-06-04 12:26:16 +02:00
for ( const auto & equation_tag : equation_tags )
if ( equation_tag . first = = eq )
2017-04-05 10:28:37 +02:00
{
if ( ! wrote_eq_tag )
content_output < < " \\ noindent[ " ;
else
content_output < < " , " ;
2018-06-04 12:26:16 +02:00
content_output < < equation_tag . second . first ;
2017-04-05 10:28:37 +02:00
2018-06-04 12:26:16 +02:00
if ( ! ( equation_tag . second . second . empty ( ) ) )
content_output < < " = ` " < < equation_tag . second . 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 )
content_output < < " ] " ;
2017-04-05 10:28:37 +02:00
}
2017-04-04 15:28:27 +02:00
2017-04-05 10:28:37 +02:00
content_output < < " \\ begin{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 ) ;
content_output < < endl < < " \\ end{dmath} " < < endl ;
2009-04-30 15:14:33 +02:00
}
2015-07-15 08:58:15 +02:00
output < < " \\ include{ " < < content_basename < < " } " < < endl
< < " \\ 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
{
2018-06-04 15:03:26 +02:00
auto * beq = dynamic_cast < BinaryOpNode * > ( eq ) ;
2018-07-18 16:18:26 +02:00
assert ( beq ! = nullptr & & beq - > get_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
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 ( ) ;
2018-06-04 12:26:16 +02:00
for ( const auto & eq_tag : eq_tags )
2018-06-04 12:48:09 +02:00
equation_tags . emplace_back ( n , eq_tag ) ;
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
{
2018-06-04 15:03:26 +02:00
auto * beq = dynamic_cast < BinaryOpNode * > ( eq ) ;
2018-07-18 16:18:26 +02:00
assert ( beq ! = nullptr & & beq - > get_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
2018-06-04 12:50:53 +02:00
ModelTree : : addTrendVariables ( vector < int > trend_vars , expr_t growth_factor ) noexcept ( false )
2010-10-15 19:05:16 +02:00
{
while ( ! trend_vars . empty ( ) )
if ( trend_symbols_map . find ( trend_vars . back ( ) ) ! = trend_symbols_map . end ( ) )
throw TrendException ( symbol_table . getName ( trend_vars . back ( ) ) ) ;
else
{
trend_symbols_map [ trend_vars . back ( ) ] = growth_factor ;
trend_vars . pop_back ( ) ;
}
}
void
2018-06-04 12:50:53 +02:00
ModelTree : : addNonstationaryVariables ( vector < int > nonstationary_vars , bool log_deflator , expr_t deflator ) noexcept ( false )
2010-10-15 19:05:16 +02:00
{
while ( ! nonstationary_vars . empty ( ) )
if ( nonstationary_symbols_map . find ( nonstationary_vars . back ( ) ) ! = nonstationary_symbols_map . end ( ) )
throw TrendException ( symbol_table . getName ( nonstationary_vars . back ( ) ) ) ;
else
{
2018-06-04 16:36:46 +02:00
nonstationary_symbols_map [ nonstationary_vars . back ( ) ] = { log_deflator , deflator } ;
2010-10-15 19:05:16 +02:00
nonstationary_vars . pop_back ( ) ;
}
}
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 + + )
2011-06-22 11:56:07 +02:00
{
equation_reordered . push_back ( j ) ;
variable_reordered . push_back ( j ) ;
}
}
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
output < < row_nb + col_nb * NNZDerivatives [ order - 1 ] ;
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
{
2016-05-18 12:26:19 +02:00
if ( ! ( paramsDerivsOrder = = 1 | | paramsDerivsOrder = = 2 ) )
2016-05-12 12:02:34 +02:00
return ;
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-06-04 12:26:16 +02:00
for ( int param : deriv_id_set )
2012-11-29 18:07:48 +01:00
{
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
{
expr_t d1 = equations [ eq ] - > getDerivative ( param ) ;
if ( d1 = = Zero )
continue ;
2018-06-04 16:36:46 +02:00
residuals_params_derivatives [ { eq , param } ] = d1 ;
2012-11-29 18:07:48 +01:00
}
2016-05-18 12:26:19 +02:00
if ( paramsDerivsOrder = = 2 )
2018-06-05 15:34:34 +02:00
for ( const auto & it : residuals_params_derivatives )
2016-05-12 12:02:34 +02:00
{
2018-06-05 15:34:34 +02:00
int eq , param1 ;
tie ( eq , param1 ) = it . first ;
expr_t d1 = it . second ;
2012-11-29 18:07:48 +01:00
2016-05-12 12:02:34 +02:00
expr_t d2 = d1 - > getDerivative ( param ) ;
if ( d2 = = Zero )
continue ;
2018-06-05 15:34:34 +02:00
residuals_params_second_derivatives [ { eq , param1 , param } ] = d2 ;
2016-05-12 12:02:34 +02:00
}
2012-11-29 18:07:48 +01:00
2018-06-05 15:34:34 +02:00
for ( const auto & it : first_derivatives )
2012-11-29 18:07:48 +01:00
{
2018-06-05 15:34:34 +02:00
int eq , var ;
tie ( eq , var ) = it . first ;
expr_t d1 = it . second ;
2012-11-29 18:07:48 +01:00
expr_t d2 = d1 - > getDerivative ( param ) ;
if ( d2 = = Zero )
continue ;
2018-06-05 15:34:34 +02:00
jacobian_params_derivatives [ { eq , var , param } ] = d2 ;
2012-11-29 18:07:48 +01:00
}
2016-05-18 12:26:19 +02:00
if ( paramsDerivsOrder = = 2 )
{
2018-06-05 15:34:34 +02:00
for ( const auto & it : jacobian_params_derivatives )
2016-05-18 12:26:19 +02:00
{
2018-06-05 15:34:34 +02:00
int eq , var , param1 ;
tie ( eq , var , param1 ) = it . first ;
expr_t d1 = it . second ;
2016-05-18 12:26:19 +02:00
expr_t d2 = d1 - > getDerivative ( param ) ;
if ( d2 = = Zero )
continue ;
2018-06-05 15:34:34 +02:00
jacobian_params_second_derivatives [ { eq , var , param1 , param } ] = d2 ;
2016-05-18 12:26:19 +02:00
}
2016-05-18 10:33:45 +02:00
2018-06-05 15:34:34 +02:00
for ( const auto & it : second_derivatives )
2016-05-18 12:26:19 +02:00
{
2018-06-05 15:34:34 +02:00
int eq , var1 , var2 ;
tie ( eq , var1 , var2 ) = it . first ;
expr_t d1 = it . second ;
2016-05-18 12:26:19 +02:00
expr_t d2 = d1 - > getDerivative ( param ) ;
if ( d2 = = Zero )
continue ;
2018-06-05 15:34:34 +02:00
hessian_params_derivatives [ { eq , var1 , var2 , param } ] = d2 ;
2016-05-18 12:26:19 +02:00
}
}
2012-11-29 18:07:48 +01:00
}
}
void
ModelTree : : computeParamsDerivativesTemporaryTerms ( )
{
2018-06-04 14:17:36 +02:00
map < expr_t , pair < int , NodeTreeReference > > reference_count ;
2012-11-29 18:07:48 +01:00
params_derivs_temporary_terms . clear ( ) ;
2015-09-03 13:50:02 +02:00
map < NodeTreeReference , temporary_terms_t > temp_terms_map ;
2018-07-18 16:35:19 +02:00
temp_terms_map [ NodeTreeReference : : residualsParamsDeriv ] = params_derivs_temporary_terms_res ;
temp_terms_map [ NodeTreeReference : : jacobianParamsDeriv ] = params_derivs_temporary_terms_g1 ;
temp_terms_map [ NodeTreeReference : : residualsParamsSecondDeriv ] = params_derivs_temporary_terms_res2 ;
temp_terms_map [ NodeTreeReference : : jacobianParamsSecondDeriv ] = params_derivs_temporary_terms_g12 ;
temp_terms_map [ NodeTreeReference : : hessianParamsDeriv ] = params_derivs_temporary_terms_g2 ;
2012-11-29 18:07:48 +01:00
2018-06-04 12:26:16 +02:00
for ( auto & residuals_params_derivative : residuals_params_derivatives )
residuals_params_derivative . second - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
true , NodeTreeReference : : residualsParamsDeriv ) ;
2012-11-29 18:07:48 +01:00
2018-06-04 12:26:16 +02:00
for ( auto & jacobian_params_derivative : jacobian_params_derivatives )
jacobian_params_derivative . second - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
true , NodeTreeReference : : jacobianParamsDeriv ) ;
2012-11-29 18:07:48 +01:00
for ( second_derivatives_t : : const_iterator it = residuals_params_second_derivatives . begin ( ) ;
it ! = residuals_params_second_derivatives . end ( ) ; + + it )
2015-09-02 18:50:09 +02:00
it - > second - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
true , NodeTreeReference : : residualsParamsSecondDeriv ) ;
2012-11-29 18:07:48 +01:00
for ( third_derivatives_t : : const_iterator it = jacobian_params_second_derivatives . begin ( ) ;
it ! = jacobian_params_second_derivatives . end ( ) ; + + it )
2015-09-02 18:50:09 +02:00
it - > second - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
true , NodeTreeReference : : jacobianParamsSecondDeriv ) ;
2012-11-29 18:07:48 +01:00
for ( third_derivatives_t : : const_iterator it = hessian_params_derivatives . begin ( ) ;
it ! = hessian_params_derivatives . end ( ) ; + + it )
2015-09-02 18:50:09 +02:00
it - > second - > computeTemporaryTerms ( reference_count ,
2015-09-03 13:50:02 +02:00
temp_terms_map ,
2018-07-18 16:35:19 +02:00
true , NodeTreeReference : : hessianParamsDeriv ) ;
2015-09-03 13:50:02 +02:00
for ( map < NodeTreeReference , temporary_terms_t > : : const_iterator it = temp_terms_map . begin ( ) ;
it ! = temp_terms_map . end ( ) ; it + + )
2015-09-03 15:47:53 +02:00
params_derivs_temporary_terms . insert ( it - > second . begin ( ) , it - > second . end ( ) ) ;
2015-09-03 13:50:02 +02:00
2018-07-18 16:35:19 +02:00
params_derivs_temporary_terms_res = temp_terms_map [ NodeTreeReference : : residualsParamsDeriv ] ;
params_derivs_temporary_terms_g1 = temp_terms_map [ NodeTreeReference : : jacobianParamsDeriv ] ;
params_derivs_temporary_terms_res2 = temp_terms_map [ NodeTreeReference : : residualsParamsSecondDeriv ] ;
params_derivs_temporary_terms_g12 = temp_terms_map [ NodeTreeReference : : jacobianParamsSecondDeriv ] ;
params_derivs_temporary_terms_g2 = temp_terms_map [ NodeTreeReference : : hessianParamsDeriv ] ;
2018-05-28 15:23:15 +02:00
int idx = 0 ;
for ( auto tt : params_derivs_temporary_terms_res )
params_derivs_temporary_terms_idxs [ tt ] = idx + + ;
for ( auto tt : params_derivs_temporary_terms_g1 )
params_derivs_temporary_terms_idxs [ tt ] = idx + + ;
for ( auto tt : params_derivs_temporary_terms_res2 )
params_derivs_temporary_terms_idxs [ tt ] = idx + + ;
for ( auto tt : params_derivs_temporary_terms_g12 )
params_derivs_temporary_terms_idxs [ tt ] = idx + + ;
for ( auto tt : params_derivs_temporary_terms_g2 )
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
{
return ( nonstationary_symbols_map . find ( symb_id )
! = nonstationary_symbols_map . end ( ) ) ;
}
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
{
2018-06-04 14:17:36 +02:00
vector < pair < string , string > > eqtags ;
2017-02-20 12:18:11 +01:00
temporary_terms_t tt_empty ;
if ( residuals )
2017-02-27 15:40:34 +01:00
output < < endl < < " \" residuals \" :[ " < < endl ;
2017-02-20 12:18:11 +01:00
else
2017-02-27 15:40:34 +01:00
output < < endl < < " \" model \" :[ " < < endl ;
2017-02-02 15:09:43 +01:00
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
{
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 ] ;
expr_t lhs = eq_node - > get_arg1 ( ) ;
expr_t rhs = eq_node - > get_arg2 ( ) ;
2017-02-20 12:18:11 +01:00
if ( residuals )
2017-02-02 15:09:43 +01:00
{
2017-02-20 12:18:11 +01:00
output < < " { \" residual \" : { "
< < " \" lhs \" : \" " ;
2018-05-29 11:59:42 +02:00
lhs - > writeJsonOutput ( output , temporary_terms , { } ) ;
2017-02-20 12:18:11 +01:00
output < < " \" " ;
output < < " , \" rhs \" : \" " ;
2018-05-29 11:59:42 +02:00
rhs - > writeJsonOutput ( output , temporary_terms , { } ) ;
2017-02-20 12:18:11 +01:00
output < < " \" " ;
try
{
// Test if the right hand side of the equation is empty.
if ( rhs - > eval ( eval_context_t ( ) ) ! = 0 )
{
output < < " , \" rhs \" : \" " ;
2018-05-29 11:59:42 +02:00
rhs - > writeJsonOutput ( output , temporary_terms , { } ) ;
2017-02-20 12:18:11 +01:00
output < < " \" " ;
}
}
catch ( ExprNode : : EvalException & e )
2017-02-02 15:09:43 +01:00
{
}
output < < " } " ;
}
2017-02-20 12:18:11 +01:00
else
{
2017-03-15 12:52:55 +01:00
output < < " { \" lhs \" : \" " ;
2018-05-29 11:59:42 +02:00
lhs - > writeJsonOutput ( output , tt_empty , { } ) ;
2017-03-15 12:52:55 +01:00
output < < " \" , \" rhs \" : \" " ;
2018-05-29 11:59:42 +02:00
rhs - > writeJsonOutput ( output , tt_empty , { } ) ;
2017-02-20 12:18:11 +01:00
output < < " \" "
< < " , \" line \" : " < < equations_lineno [ eq ] ;
2018-06-04 12:26:16 +02:00
for ( const auto & equation_tag : equation_tags )
if ( equation_tag . first = = eq )
eqtags . push_back ( equation_tag . second ) ;
2017-02-20 12:18:11 +01:00
if ( ! eqtags . empty ( ) )
{
output < < " , \" tags \" : { " ;
int i = 0 ;
2018-06-04 14:17:36 +02:00
for ( vector < pair < string , string > > : : const_iterator it = eqtags . begin ( ) ; it ! = eqtags . end ( ) ; it + + , i + + )
2017-02-20 12:18:11 +01:00
{
if ( i ! = 0 )
output < < " , " ;
output < < " \" " < < it - > first < < " \" : \" " < < it - > second < < " \" " ;
}
output < < " } " ;
eqtags . clear ( ) ;
}
}
output < < " } " < < endl ;
2017-02-02 15:09:43 +01:00
}
output < < endl < < " ] " < < endl ;
}