2008-02-03 11:28:36 +01:00
/*
2014-01-27 16:41:43 +01:00
* Copyright ( C ) 2003 - 2014 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 boost ;
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
{
const int n = equation_number ( ) ;
assert ( n = = symbol_table . endo_nbr ( ) ) ;
typedef adjacency_list < vecS , vecS , undirectedS > BipartiteGraph ;
/*
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
set < pair < int , int > > endo ;
2010-09-16 19:00:48 +02:00
for ( jacob_map_t : : const_iterator it = contemporaneous_jacobian . begin ( ) ; it ! = contemporaneous_jacobian . end ( ) ; it + + )
2009-12-21 11:29:21 +01:00
add_edge ( it - > first . first + n , it - > first . second , g ) ;
2009-12-16 14:21:31 +01:00
// Compute maximum cardinality matching
vector < int > mate_map ( 2 * n ) ;
# if 1
bool check = checked_edmonds_maximum_cardinality_matching ( g , & mate_map [ 0 ] ) ;
# else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
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
endo2eq . resize ( equation_number ( ) ) ;
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
vector < int > : : const_iterator it = find ( mate_map . begin ( ) , mate_map . begin ( ) + n , graph_traits < BipartiteGraph > : : null_vertex ( ) ) ;
if ( it ! = mate_map . begin ( ) + n )
2009-12-21 11:29:21 +01:00
{
if ( verbose )
cerr < < " ERROR: Could not normalize the model. Variable "
< < symbol_table . getName ( symbol_table . getID ( eEndogenous , it - mate_map . begin ( ) ) )
< < " 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 ;
2009-12-21 11:29:21 +01:00
int n = equation_number ( ) ;
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
2010-09-16 19:00:48 +02:00
for ( jacob_map_t : : iterator iter = normalized_contemporaneous_jacobian . begin ( ) ; iter ! = normalized_contemporaneous_jacobian . end ( ) ; iter + + )
2009-12-21 11:29:21 +01:00
iter - > second / = max_val [ iter - > first . first ] ;
//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 ;
2010-09-16 19:00:48 +02:00
for ( jacob_map_t : : iterator iter = normalized_contemporaneous_jacobian . begin ( ) ; iter ! = normalized_contemporaneous_jacobian . end ( ) ; iter + + )
2009-12-21 11:29:21 +01:00
if ( fabs ( iter - > second ) > max ( current_cutoff , cutoff ) )
tmp_normalized_contemporaneous_jacobian [ make_pair ( iter - > first . first , iter - > first . second ) ] = iter - > second ;
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 ;
2009-12-16 14:21:31 +01: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 ) ;
2009-12-16 18:13:23 +01:00
for ( set < pair < int , int > > : : const_iterator it = endo . begin ( ) ; it ! = endo . end ( ) ; it + + )
2009-12-21 11:29:21 +01:00
tmp_normalized_contemporaneous_jacobian [ make_pair ( 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
{
if ( static_jacobian . find ( make_pair ( it - > first . first , it - > first . second ) ) = = static_jacobian . end ( ) )
static_jacobian [ make_pair ( it - > first . first , it - > first . second ) ] = 0 ;
if ( dynamic_jacobian . find ( make_pair ( 0 , make_pair ( it - > first . first , it - > first . second ) ) ) = = dynamic_jacobian . end ( ) )
dynamic_jacobian [ make_pair ( 0 , make_pair ( it - > first . first , it - > first . second ) ) ] = 0 ;
if ( contemporaneous_jacobian . find ( make_pair ( it - > first . first , it - > first . second ) ) = = contemporaneous_jacobian . end ( ) )
contemporaneous_jacobian [ make_pair ( it - > first . first , it - > first . second ) ] = 0 ;
2010-07-23 11:20:24 +02:00
if ( first_derivatives . find ( make_pair ( it - > first . first , getDerivID ( symbol_table . getID ( eEndogenous , it - > first . second ) , 0 ) ) ) = = first_derivatives . end ( ) )
first_derivatives [ make_pair ( it - > first . first , getDerivID ( symbol_table . getID ( eEndogenous , it - > first . second ) , 0 ) ) ] = Zero ;
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
{
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < equation_number ( ) ; i + + )
2009-12-16 14:21:31 +01:00
{
VariableNode * lhs = dynamic_cast < VariableNode * > ( equations [ i ] - > get_arg1 ( ) ) ;
if ( lhs = = NULL )
continue ;
int symb_id = lhs - > get_symb_id ( ) ;
if ( symbol_table . getType ( symb_id ) ! = eEndogenous )
continue ;
set < pair < int , int > > endo ;
equations [ i ] - > get_arg2 ( ) - > collectEndogenous ( endo ) ;
if ( endo . find ( make_pair ( symbol_table . getTypeSpecificID ( symb_id ) , 0 ) ) ! = endo . end ( ) )
continue ;
endo2eqs . insert ( make_pair ( symbol_table . getTypeSpecificID ( symb_id ) , i ) ) ;
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 ;
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 ;
if ( getTypeByDerivID ( deriv_id ) = = eEndogenous )
{
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 ;
2009-12-16 14:21:31 +01:00
Id - > writeOutput ( cerr , oMatlabDynamicModelSparse , temporary_terms ) ;
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 ;
jacobian_elements_to_delete . insert ( make_pair ( eq , deriv_id ) ) ;
}
else
{
if ( lag = = 0 )
{
nb_elements_contemparenous_Jacobian + + ;
2009-12-16 18:13:23 +01:00
contemporaneous_jacobian [ make_pair ( eq , var ) ] = val ;
2009-12-16 14:21:31 +01:00
}
if ( static_jacobian . find ( make_pair ( eq , var ) ) ! = static_jacobian . end ( ) )
static_jacobian [ make_pair ( eq , var ) ] + = val ;
else
static_jacobian [ make_pair ( eq , var ) ] = val ;
dynamic_jacobian [ make_pair ( lag , make_pair ( eq , var ) ) ] = Id ;
}
}
}
// Get rid of the elements of the Jacobian matrix below the cutoff
2009-12-16 18:13:23 +01:00
for ( set < pair < int , int > > : : const_iterator it = jacobian_elements_to_delete . begin ( ) ; it ! = jacobian_elements_to_delete . end ( ) ; it + + )
2009-12-16 14:21:31 +01:00
first_derivatives . erase ( * it ) ;
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 ;
}
}
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
{
2009-12-16 18:13:23 +01:00
vector < int > eq2endo ( equation_number ( ) , 0 ) ;
2009-12-16 14:21:31 +01:00
equation_reordered . resize ( equation_number ( ) ) ;
variable_reordered . resize ( equation_number ( ) ) ;
bool * IM ;
int n = equation_number ( ) ;
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 ;
}
2010-09-16 19:00:48 +02:00
for ( jacob_map_t : : const_iterator it = static_jacobian_arg . begin ( ) ; it ! = static_jacobian_arg . end ( ) ; it + + )
2009-12-16 14:21:31 +01:00
IM [ it - > first . first * n + endo2eq [ it - > first . second ] ] = true ;
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
2010-09-16 19:18:45 +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 ;
2010-09-16 19:18:45 +02:00
map < pair < int , pair < int , int > > , expr_t > : : const_iterator derivative = first_order_endo_derivatives . find ( make_pair ( eq , make_pair ( var , 0 ) ) ) ;
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
{
set < pair < int , int > > result ;
derivative - > second - > collectEndogenous ( result ) ;
set < pair < int , int > > : : const_iterator d_endo_variable = result . find ( make_pair ( var , 0 ) ) ;
//Determine whether the equation could be evaluated rather than to be solved
2010-04-16 16:54:55 +02:00
if ( lhs - > isVariableNodeEqualTo ( eEndogenous , Index_Var_IM [ i ] , 0 ) & & derivative - > second - > isNumConstNodeEqualTo ( 1 ) )
2009-12-16 14:21:31 +01:00
{
Equation_Simulation_Type = E_EVALUATE ;
}
else
{
2010-09-16 19:18:45 +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 ;
}
}
}
V_Equation_Simulation_Type [ eq ] = make_pair ( Equation_Simulation_Type , dynamic_cast < BinaryOpNode * > ( res . second ) ) ;
}
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 ( ) ;
2010-09-16 19:00:48 +02:00
variable_lead_lag = lag_lead_vector_t ( nb_endo , make_pair ( 0 , 0 ) ) ;
equation_lead_lag = lag_lead_vector_t ( nb_endo , make_pair ( 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 ) ;
}
}
2010-09-16 19:00:48 +02:00
for ( dynamic_jacob_map_t : : const_iterator it = dynamic_jacobian . begin ( ) ; it ! = dynamic_jacobian . end ( ) ; it + + )
2009-12-16 14:21:31 +01:00
{
int lag = it - > first . first ;
2010-07-23 11:20:24 +02:00
int j_1 = it - > first . second . first ;
2009-12-16 14:21:31 +01:00
int i_1 = it - > first . second . second ;
if ( variable_blck [ i_1 ] = = equation_blck [ j_1 ] )
{
if ( lag > variable_lead_lag [ i_1 ] . second )
variable_lead_lag [ i_1 ] = make_pair ( variable_lead_lag [ i_1 ] . first , lag ) ;
if ( lag < - variable_lead_lag [ i_1 ] . first )
variable_lead_lag [ i_1 ] = make_pair ( - lag , variable_lead_lag [ i_1 ] . second ) ;
if ( lag > equation_lead_lag [ j_1 ] . second )
equation_lead_lag [ j_1 ] = make_pair ( equation_lead_lag [ j_1 ] . first , lag ) ;
if ( lag < - equation_lead_lag [ j_1 ] . first )
equation_lead_lag [ j_1 ] = make_pair ( - lag , equation_lead_lag [ j_1 ] . second ) ;
}
}
}
void
2010-07-23 11:20:24 +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
property_map < AdjacencyList_t , vertex_index_t > : : type v_index = get ( vertex_index , G2 ) ;
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 ;
}
2010-09-16 19:00:48 +02:00
for ( jacob_map_t : : const_iterator it = static_jacobian . begin ( ) ; it ! = static_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 ) ) ;
2011-02-04 16:25:38 +01:00
iterator_property_map < int * , property_map < AdjacencyList_t , vertex_index_t > : : type , int , int & > endo2block_map ( & endo2block [ 0 ] , get ( 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
blocks = vector < pair < int , int > > ( num , make_pair ( 0 , 0 ) ) ;
// Create directed acyclic graph associated to the strongly connected components
typedef adjacency_list < vecS , vecS , directedS > DirectedGraph ;
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.
vector < pair < set < int > , pair < set < int > , vector < int > > > > components_set ( num ) ;
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 ) ;
property_map < AdjacencyList_t , vertex_index_t > : : type v_index = get ( 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
{
2010-07-23 11:20:24 +02:00
for ( vector < int > : : iterator its = Reordered_Vertice . begin ( ) ; its ! = Reordered_Vertice . end ( ) ; its + + )
{
bool something_done = false ;
2010-11-20 15:33:56 +01: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 ;
}
2010-11-20 15:33:56 +01: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 ;
}
2010-11-20 15:33:56 +01: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 ;
}
2010-11-20 15:33:56 +01: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 )
{
equation_reordered [ order ] = tmp_equation_reordered [ * its + prologue ] ;
variable_reordered [ order ] = tmp_variable_reordered [ * its + prologue ] ;
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
{
2010-07-23 11:20:24 +02:00
for ( set < int > : : iterator its = feed_back_vertices . begin ( ) ; its ! = feed_back_vertices . end ( ) ; its + + )
{
bool something_done = false ;
2010-11-20 15:33:56 +01:00
if ( j = = 2 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ] . second ! = 0 )
2010-07-23 11:20:24 +02:00
{
n_mixed [ prologue + i ] + + ;
something_done = true ;
}
2010-11-20 15:33:56 +01:00
else if ( j = = 3 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ] . second ! = 0 )
2010-07-23 11:20:24 +02:00
{
n_forward [ prologue + i ] + + ;
something_done = true ;
}
2010-11-20 15:33:56 +01:00
else if ( j = = 1 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ] . first ! = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ] . second = = 0 )
2010-07-23 11:20:24 +02:00
{
n_backward [ prologue + i ] + + ;
something_done = true ;
}
2010-11-20 15:33:56 +01:00
else if ( j = = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ] . first = = 0 & & variable_lag_lead [ tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ] . second = = 0 )
2010-07-23 11:20:24 +02:00
{
n_static [ prologue + i ] + + ;
something_done = true ;
}
if ( something_done )
{
equation_reordered [ order ] = tmp_equation_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ;
variable_reordered [ order ] = tmp_variable_reordered [ v_index [ vertex ( * its , G ) ] + prologue ] ;
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
2010-09-16 17:51:50 +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
2011-02-04 16:25:38 +01:00
ModelTree : : reduceBlocksAndTypeDetermination ( const dynamic_jacob_map_t & dynamic_jacobian , vector < pair < int , int > > & blocks , const equation_type_and_normalized_equation_t & Equation_Type , const vector < int > & variable_reordered , const vector < int > & equation_reordered , vector < unsigned int > & n_static , vector < unsigned int > & n_forward , vector < unsigned int > & n_backward , vector < unsigned int > & n_mixed , vector < pair < pair < int , int > , pair < int , int > > > & block_col_type )
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 ;
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 ) ;
for ( set < pair < int , int > > : : const_iterator it = endo . begin ( ) ; it ! = endo . end ( ) ; it + + )
{
int curr_variable = it - > first ;
int curr_lag = it - > second ;
2012-02-17 10:50:54 +01:00
vector < int > : : const_iterator it1 = find ( variable_reordered . begin ( ) + first_count_equ , variable_reordered . begin ( ) + ( first_count_equ + Blck_Size ) , curr_variable ) ;
if ( it1 ! = variable_reordered . begin ( ) + ( first_count_equ + Blck_Size ) )
2009-12-16 14:21:31 +01:00
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 ;
}
}
}
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 )
| | ( prev_Type = = EVALUATE_BACKWARD & & Simulation_Type = = EVALUATE_BACKWARD & & ! is_lag ) ) )
{
for ( int j = first_equation ; j < first_equation + c_Size ; j + + )
{
dynamic_jacob_map_t : : const_iterator it = dynamic_jacobian . find ( make_pair ( - 1 , make_pair ( equation_reordered [ eq ] , variable_reordered [ j ] ) ) ) ;
if ( it ! = dynamic_jacobian . end ( ) )
is_lag = true ;
it = dynamic_jacobian . find ( make_pair ( + 1 , make_pair ( equation_reordered [ eq ] , variable_reordered [ j ] ) ) ) ;
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 + + ;
block_type_size_mfs [ block_type_size_mfs . size ( ) - 1 ] = make_pair ( make_pair ( c_Type , first_equation ) , make_pair ( 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 ;
block_lag_lead [ block_type_size_mfs . size ( ) - 1 ] = make_pair ( Lag , Lead ) ;
2010-07-23 11:20:24 +02:00
pair < pair < unsigned int , unsigned int > , pair < unsigned int , unsigned int > > tmp = block_col_type [ block_col_type . size ( ) - 1 ] ;
2011-02-04 16:25:38 +01:00
block_col_type [ block_col_type . size ( ) - 1 ] = make_pair ( make_pair ( tmp . first . first + l_n_static , tmp . first . second + l_n_forward ) , make_pair ( tmp . second . first + l_n_backward , tmp . second . second + l_n_mixed ) ) ;
2009-12-16 14:21:31 +01:00
}
else
{
block_type_size_mfs . push_back ( make_pair ( make_pair ( Simulation_Type , eq ) , make_pair ( Blck_Size , MFS_Size ) ) ) ;
block_lag_lead . push_back ( make_pair ( Lag , Lead ) ) ;
2011-02-04 16:25:38 +01:00
block_col_type . push_back ( make_pair ( 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
{
block_type_size_mfs . push_back ( make_pair ( make_pair ( Simulation_Type , eq ) , make_pair ( Blck_Size , MFS_Size ) ) ) ;
block_lag_lead . push_back ( make_pair ( Lag , Lead ) ) ;
2011-02-04 16:25:38 +01:00
block_col_type . push_back ( make_pair ( 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
{
block_type_size_mfs . push_back ( make_pair ( make_pair ( Simulation_Type , eq ) , make_pair ( Blck_Size , MFS_Size ) ) ) ;
block_lag_lead . push_back ( make_pair ( Lag , Lead ) ) ;
2011-02-04 16:25:38 +01:00
block_col_type . push_back ( make_pair ( 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 ) ;
}
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 ;
2009-12-16 14:21:31 +01:00
set < pair < int , int > > endogenous ;
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
{
if ( endogenous . find ( make_pair ( variable_reordered [ first_variable_position + l ] , 0 ) ) ! = endogenous . end ( ) )
{
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 ; //
2009-12-16 14:21:31 +01:00
set < pair < int , int > > endogenous ;
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
{
if ( endogenous . find ( make_pair ( variable_reordered [ first_variable_position + l ] , lag ) ) ! = endogenous . end ( ) )
{
blocks_linear [ block ] = false ;
goto the_end ;
}
}
}
}
}
2009-12-16 18:13:23 +01:00
the_end :
;
2009-12-16 14:21:31 +01:00
}
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 ,
ExternalFunctionsTable & external_functions_table_arg ) :
2011-06-22 11:56:07 +02:00
DataTree ( symbol_table_arg , num_constants_arg , external_functions_table_arg ) ,
cutoff ( 1e-15 ) ,
mfs ( 0 )
2008-02-03 11:28:36 +01:00
{
2009-12-16 18:13:23 +01:00
for ( int i = 0 ; i < 3 ; i + + )
2009-05-12 22:32:27 +02:00
NNZDerivatives [ i ] = 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
{
2010-09-16 19:00:48 +02:00
first_derivatives_t : : const_iterator it = first_derivatives . find ( make_pair ( eq , getDerivID ( symb_id , lag ) ) ) ;
2009-01-23 11:59:37 +01:00
if ( it ! = first_derivatives . end ( ) )
( it - > second ) - > writeOutput ( output , output_type , temporary_terms ) ;
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
{
2009-12-16 18:13:23 +01:00
for ( set < int > : : const_iterator it = vars . begin ( ) ;
it ! = vars . end ( ) ; it + + )
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 + + )
{
expr_t d1 = equations [ eq ] - > getDerivative ( * it ) ;
if ( d1 = = Zero )
continue ;
first_derivatives [ make_pair ( eq , * it ) ] = d1 ;
+ + NNZDerivatives [ 0 ] ;
}
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 )
{
2010-09-16 19:00:48 +02:00
for ( first_derivatives_t : : const_iterator it = first_derivatives . begin ( ) ;
2009-04-20 12:48:54 +02:00
it ! = first_derivatives . end ( ) ; it + + )
2008-02-03 11:28:36 +01:00
{
2009-04-20 12:48:54 +02:00
int eq = it - > first . first ;
int var1 = it - > first . second ;
2010-09-16 19:18:45 +02:00
expr_t d1 = it - > second ;
2009-04-20 12:48:54 +02:00
// Store only second derivatives with var2 <= var1
2009-12-16 18:13:23 +01:00
for ( set < int > : : const_iterator it2 = vars . begin ( ) ;
it2 ! = vars . end ( ) ; it2 + + )
2008-02-03 11:28:36 +01:00
{
2009-04-20 12:48:54 +02:00
int var2 = * it2 ;
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 ;
second_derivatives [ make_pair ( eq , make_pair ( 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 )
{
2010-09-16 19:00:48 +02:00
for ( second_derivatives_t : : const_iterator it = second_derivatives . begin ( ) ;
2009-04-20 12:48:54 +02:00
it ! = second_derivatives . end ( ) ; it + + )
2008-02-03 11:28:36 +01:00
{
2009-04-20 12:48:54 +02:00
int eq = it - > first . first ;
int var1 = it - > first . second . first ;
int var2 = it - > first . second . second ;
// By construction, var2 <= var1
2010-09-16 19:18:45 +02:00
expr_t d2 = it - > second ;
2009-04-20 12:48:54 +02:00
// Store only third derivatives such that var3 <= var2 <= var1
2009-12-16 18:13:23 +01:00
for ( set < int > : : const_iterator it2 = vars . begin ( ) ;
it2 ! = vars . end ( ) ; it2 + + )
2008-02-03 11:28:36 +01:00
{
2009-04-20 12:48:54 +02:00
int var3 = * it2 ;
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 ;
third_derivatives [ make_pair ( eq , make_pair ( var1 , make_pair ( 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
2009-07-07 16:20:48 +02:00
ModelTree : : computeTemporaryTerms ( bool is_matlab )
2008-02-03 11:28:36 +01:00
{
2010-09-16 19:18:45 +02:00
map < expr_t , int > reference_count ;
2008-02-03 11:28:36 +01:00
temporary_terms . clear ( ) ;
2008-12-24 16:49:01 +01:00
for ( vector < BinaryOpNode * > : : iterator it = equations . begin ( ) ;
it ! = equations . end ( ) ; it + + )
2008-02-03 11:28:36 +01:00
( * it ) - > computeTemporaryTerms ( reference_count , temporary_terms , is_matlab ) ;
2010-09-16 19:00:48 +02:00
for ( first_derivatives_t : : iterator it = first_derivatives . begin ( ) ;
2008-12-24 16:49:01 +01:00
it ! = first_derivatives . end ( ) ; it + + )
2008-02-03 11:28:36 +01:00
it - > second - > computeTemporaryTerms ( reference_count , temporary_terms , is_matlab ) ;
2010-09-16 19:00:48 +02:00
for ( second_derivatives_t : : iterator it = second_derivatives . begin ( ) ;
2009-04-20 12:48:54 +02:00
it ! = second_derivatives . end ( ) ; it + + )
it - > second - > computeTemporaryTerms ( reference_count , temporary_terms , is_matlab ) ;
2008-02-03 11:28:36 +01:00
2010-09-16 19:00:48 +02:00
for ( third_derivatives_t : : iterator it = third_derivatives . begin ( ) ;
2009-04-20 12:48:54 +02:00
it ! = third_derivatives . end ( ) ; it + + )
it - > second - > computeTemporaryTerms ( reference_count , temporary_terms , is_matlab ) ;
2008-02-03 11:28:36 +01:00
}
void
2010-09-16 19:00:48 +02:00
ModelTree : : writeTemporaryTerms ( const temporary_terms_t & tt , ostream & output ,
2011-04-12 16:41:29 +02:00
ExprNodeOutputType output_type , deriv_node_temp_terms_t & tef_terms ) const
2009-01-23 11:59:37 +01:00
{
2009-04-20 15:58:15 +02:00
// Local var used to keep track of temp nodes already written
2010-09-16 19:00:48 +02:00
temporary_terms_t tt2 ;
2008-02-03 11:28:36 +01:00
2010-09-16 19:00:48 +02:00
for ( temporary_terms_t : : const_iterator it = tt . begin ( ) ;
2009-04-20 15:58:15 +02:00
it ! = tt . end ( ) ; it + + )
2009-01-23 11:59:37 +01:00
{
2010-03-04 16:40:07 +01:00
if ( dynamic_cast < ExternalFunctionNode * > ( * it ) ! = NULL )
( * it ) - > writeExternalFunctionOutput ( output , output_type , tt2 , tef_terms ) ;
2010-10-05 18:48:15 +02:00
if ( IS_C ( output_type ) )
output < < " double " ;
2010-03-04 16:40:07 +01:00
( * it ) - > writeOutput ( output , output_type , tt , tef_terms ) ;
2009-01-23 11:59:37 +01:00
output < < " = " ;
2010-03-04 16:40:07 +01:00
( * it ) - > writeOutput ( output , output_type , tt2 , tef_terms ) ;
2008-02-03 11:28:36 +01:00
2010-10-05 18:48:15 +02:00
if ( IS_C ( output_type ) )
output < < " ; " < < endl ;
2009-01-23 11:59:37 +01:00
// Insert current node into tt2
tt2 . insert ( * it ) ;
2008-02-03 11:28:36 +01:00
2009-04-30 15:14:33 +02:00
if ( IS_MATLAB ( output_type ) )
2009-01-23 11:59:37 +01:00
output < < " ; " < < endl ;
}
}
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 : : 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 ;
2010-09-16 19:00:48 +02:00
for ( temporary_terms_t : : const_iterator it = tt . begin ( ) ;
2010-01-22 11:03:29 +01:00
it ! = tt . end ( ) ; it + + )
{
2010-12-10 11:50:27 +01:00
if ( dynamic_cast < ExternalFunctionNode * > ( * it ) ! = NULL )
{
( * it ) - > compileExternalFunctionOutput ( code_file , instruction_number , false , tt2 , map_idx , dynamic , steady_dynamic , tef_terms ) ;
}
2011-02-04 16:25:38 +01:00
FNUMEXPR_ fnumexpr ( TemporaryTerm , ( int ) ( map_idx . find ( ( * it ) - > idx ) - > second ) ) ;
2010-07-23 11:20:24 +02:00
fnumexpr . write ( code_file , instruction_number ) ;
2010-12-10 11:50:27 +01:00
( * it ) - > compile ( code_file , instruction_number , false , tt2 , map_idx , dynamic , steady_dynamic , tef_terms ) ;
2010-01-22 11:03:29 +01:00
if ( dynamic )
{
2011-02-04 16:25:38 +01: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
{
2011-02-04 16:25:38 +01:00
FSTPST_ fstpst ( ( int ) ( map_idx . find ( ( * it ) - > idx ) - > second ) ) ;
2010-07-23 11:20:24 +02:00
fstpst . write ( code_file , instruction_number ) ;
2010-01-22 11:03:29 +01:00
}
// Insert current node into tt2
tt2 . insert ( * it ) ;
}
}
2008-02-03 11:28:36 +01:00
void
2011-04-12 16:41:29 +02:00
ModelTree : : writeModelLocalVariables ( ostream & output , ExprNodeOutputType output_type , deriv_node_temp_terms_t & tef_terms ) const
2009-01-23 11:59:37 +01:00
{
2010-09-23 12:13:30 +02:00
/* 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 ;
2011-04-12 16:41:29 +02:00
// Use an empty set for the temporary terms
const temporary_terms_t tt ;
2010-09-23 12:13:30 +02:00
for ( size_t i = 0 ; i < equations . size ( ) ; i + + )
2013-11-29 15:32:49 +01:00
equations [ i ] - > collectVariables ( eModelLocalVariable , used_local_vars ) ;
2010-09-23 12:13:30 +02:00
for ( set < int > : : const_iterator it = used_local_vars . begin ( ) ;
it ! = used_local_vars . end ( ) ; + + it )
2009-01-23 11:59:37 +01:00
{
2010-09-23 12:13:30 +02:00
int id = * it ;
expr_t value = local_variables_table . find ( id ) - > second ;
2011-04-12 16:41:29 +02:00
value - > writeExternalFunctionOutput ( output , output_type , tt , tef_terms ) ;
2008-02-03 11:28:36 +01:00
2009-04-30 15:14:33 +02:00
if ( IS_C ( output_type ) )
2009-01-23 11:59:37 +01:00
output < < " double " ;
2008-02-03 11:28:36 +01:00
2011-02-21 16:32:21 +01:00
/* We append underscores to avoid name clashes with "g1" or "oo_" (see
also VariableNode : : writeOutput ) */
output < < symbol_table . getName ( id ) < < " __ = " ;
2011-04-12 16:41:29 +02:00
value - > writeOutput ( output , output_type , tt , tef_terms ) ;
2009-01-23 11:59:37 +01:00
output < < " ; " < < endl ;
}
}
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
{
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;
2009-10-29 18:16:10 +01:00
{
output < < " lhs = " ;
lhs - > writeOutput ( output , output_type , temporary_terms ) ;
output < < " ; " < < endl ;
output < < " rhs = " ;
rhs - > writeOutput ( output , output_type , temporary_terms ) ;
output < < " ; " < < endl ;
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 )
2009-10-29 18:16:10 +01:00
< < " = 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
{
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 )
< < " = " ;
2009-10-29 18:16:10 +01:00
lhs - > writeOutput ( output , output_type , temporary_terms ) ;
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
FBINARY_ fbinary ( oMinus ) ;
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
ModelTree : : Write_Inf_To_Bin_File ( const string & basename ,
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 ;
const string bin_basename = basename + " .bin " ;
if ( file_open )
SaveCode . open ( bin_basename . c_str ( ) , ios : : out | ios : : in | ios : : binary | ios : : ate ) ;
else
SaveCode . open ( bin_basename . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! SaveCode . is_open ( ) )
{
cout < < " Error : Can't open file \" " < < bin_basename < < " \" for writing \n " ;
exit ( EXIT_FAILURE ) ;
}
u_count_int = 0 ;
2011-02-04 16:25:38 +01:00
for ( first_derivatives_t : : const_iterator it = first_derivatives . begin ( ) ; it ! = first_derivatives . end ( ) ; it + + )
2010-01-22 11:03:29 +01:00
{
int deriv_id = it - > first . second ;
if ( getTypeByDerivID ( deriv_id ) = = eEndogenous )
{
int eq = it - > first . first ;
int symb = getSymbIDByDerivID ( deriv_id ) ;
int var = symbol_table . getTypeSpecificID ( symb ) ;
int lag = getLagByDerivID ( deriv_id ) ;
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
ModelTree : : writeLatexModelFile ( const string & filename , ExprNodeOutputType output_type ) const
{
ofstream output ;
output . open ( filename . c_str ( ) , ios : : out | ios : : binary ) ;
if ( ! output . is_open ( ) )
{
cerr < < " ERROR: Can't open file " < < filename < < " for writing " < < endl ;
exit ( EXIT_FAILURE ) ;
}
output < < " \\ documentclass[10pt,a4paper]{article} " < < endl
< < " \\ usepackage[landscape]{geometry} " < < endl
< < " \\ usepackage{fullpage} " < < 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
2010-09-16 19:18:45 +02:00
for ( map < int , expr_t > : : const_iterator it = local_variables_table . begin ( ) ;
2009-04-30 15:14:33 +02:00
it ! = local_variables_table . end ( ) ; it + + )
{
int id = it - > first ;
2010-09-16 19:18:45 +02:00
expr_t value = it - > second ;
2009-04-30 15:14:33 +02:00
2011-10-09 18:24:39 +02:00
output < < " \\ begin{dmath*} " < < endl
2009-04-30 15:14:33 +02:00
< < symbol_table . getName ( id ) < < " = " ;
// Use an empty set for the temporary terms
2010-03-09 12:16:32 +01:00
value - > writeOutput ( output , output_type ) ;
2011-10-09 18:24:39 +02:00
output < < endl < < " \\ end{dmath*} " < < endl ;
2009-04-30 15:14:33 +02:00
}
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
{
2011-10-09 18:24:39 +02:00
output < < " \\ begin{dmath} " < < endl
2009-04-30 15:14:33 +02:00
< < " % Equation " < < eq + 1 < < 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
dynamic_cast < ExprNode * > ( equations [ eq ] ) - > writeOutput ( output , output_type ) ;
2011-10-09 18:24:39 +02:00
output < < endl < < " \\ end{dmath} " < < endl ;
2009-04-30 15:14:33 +02:00
}
output < < " \\ end{document} " < < endl ;
output . close ( ) ;
}
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
{
2009-04-14 16:39:53 +02:00
BinaryOpNode * beq = dynamic_cast < BinaryOpNode * > ( eq ) ;
2009-04-28 18:21:39 +02:00
assert ( beq ! = NULL & & beq - > get_op_code ( ) = = oEqual ) ;
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
2014-01-27 16:41:43 +01:00
ModelTree : : addEquation ( expr_t eq , int lineno , vector < pair < string , string > > & eq_tags )
2009-09-02 16:37:59 +02:00
{
2013-04-11 17:07:39 +02:00
int n = equation_number ( ) ;
for ( size_t i = 0 ; i < eq_tags . size ( ) ; i + + )
equation_tags . push_back ( make_pair ( n , eq_tags [ i ] ) ) ;
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
{
BinaryOpNode * beq = dynamic_cast < BinaryOpNode * > ( eq ) ;
assert ( beq ! = NULL & & beq - > get_op_code ( ) = = oEqual ) ;
aux_equations . push_back ( beq ) ;
}
2010-10-15 19:05:16 +02:00
void
ModelTree : : addTrendVariables ( vector < int > trend_vars , expr_t growth_factor ) throw ( TrendException )
{
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
2013-03-26 16:46:18 +01:00
ModelTree : : addNonstationaryVariables ( vector < int > nonstationary_vars , bool log_deflator , expr_t deflator ) throw ( TrendException )
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
{
2013-03-26 16:46:18 +01:00
nonstationary_symbols_map [ nonstationary_vars . back ( ) ] = make_pair ( 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 ( )
{
for ( int j = 0 ; j < equation_number ( ) ; j + + )
{
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
{
output < < " g1 " < < LEFT_ARRAY_SUBSCRIPT ( output_type ) ;
if ( IS_MATLAB ( output_type ) )
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
{
output < < " v " < < order < < LEFT_ARRAY_SUBSCRIPT ( output_type ) ;
if ( IS_MATLAB ( output_type ) )
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
ModelTree : : computeParamsDerivatives ( )
{
set < int > deriv_id_set ;
addAllParamDerivId ( deriv_id_set ) ;
for ( set < int > : : const_iterator it = deriv_id_set . begin ( ) ;
it ! = deriv_id_set . end ( ) ; it + + )
{
const int param = * it ;
for ( int eq = 0 ; eq < ( int ) equations . size ( ) ; eq + + )
{
expr_t d1 = equations [ eq ] - > getDerivative ( param ) ;
if ( d1 = = Zero )
continue ;
residuals_params_derivatives [ make_pair ( eq , param ) ] = d1 ;
}
for ( first_derivatives_t : : const_iterator it2 = residuals_params_derivatives . begin ( ) ;
it2 ! = residuals_params_derivatives . end ( ) ; it2 + + )
{
int eq = it2 - > first . first ;
int param1 = it2 - > first . second ;
expr_t d1 = it2 - > second ;
expr_t d2 = d1 - > getDerivative ( param ) ;
if ( d2 = = Zero )
continue ;
residuals_params_second_derivatives [ make_pair ( eq , make_pair ( param1 , param ) ) ] = d2 ;
}
for ( first_derivatives_t : : const_iterator it2 = first_derivatives . begin ( ) ;
it2 ! = first_derivatives . end ( ) ; it2 + + )
{
int eq = it2 - > first . first ;
int var = it2 - > first . second ;
expr_t d1 = it2 - > second ;
expr_t d2 = d1 - > getDerivative ( param ) ;
if ( d2 = = Zero )
continue ;
jacobian_params_derivatives [ make_pair ( eq , make_pair ( var , param ) ) ] = d2 ;
}
for ( second_derivatives_t : : const_iterator it2 = jacobian_params_derivatives . begin ( ) ;
it2 ! = jacobian_params_derivatives . end ( ) ; it2 + + )
{
int eq = it2 - > first . first ;
int var = it2 - > first . second . first ;
int param1 = it2 - > first . second . second ;
expr_t d1 = it2 - > second ;
expr_t d2 = d1 - > getDerivative ( param ) ;
if ( d2 = = Zero )
continue ;
jacobian_params_second_derivatives [ make_pair ( eq , make_pair ( var , make_pair ( param1 , param ) ) ) ] = d2 ;
}
for ( second_derivatives_t : : const_iterator it2 = second_derivatives . begin ( ) ;
it2 ! = second_derivatives . end ( ) ; it2 + + )
{
int eq = it2 - > first . first ;
int var1 = it2 - > first . second . first ;
int var2 = it2 - > first . second . second ;
expr_t d1 = it2 - > second ;
expr_t d2 = d1 - > getDerivative ( param ) ;
if ( d2 = = Zero )
continue ;
hessian_params_derivatives [ make_pair ( eq , make_pair ( var1 , make_pair ( var2 , param ) ) ) ] = d2 ;
}
}
}
void
ModelTree : : computeParamsDerivativesTemporaryTerms ( )
{
map < expr_t , int > reference_count ;
params_derivs_temporary_terms . clear ( ) ;
for ( first_derivatives_t : : iterator it = residuals_params_derivatives . begin ( ) ;
it ! = residuals_params_derivatives . end ( ) ; it + + )
it - > second - > computeTemporaryTerms ( reference_count , params_derivs_temporary_terms , true ) ;
for ( second_derivatives_t : : iterator it = jacobian_params_derivatives . begin ( ) ;
it ! = jacobian_params_derivatives . end ( ) ; it + + )
it - > second - > computeTemporaryTerms ( reference_count , params_derivs_temporary_terms , true ) ;
for ( second_derivatives_t : : const_iterator it = residuals_params_second_derivatives . begin ( ) ;
it ! = residuals_params_second_derivatives . end ( ) ; + + it )
it - > second - > computeTemporaryTerms ( reference_count , params_derivs_temporary_terms , true ) ;
for ( third_derivatives_t : : const_iterator it = jacobian_params_second_derivatives . begin ( ) ;
it ! = jacobian_params_second_derivatives . end ( ) ; + + it )
it - > second - > computeTemporaryTerms ( reference_count , params_derivs_temporary_terms , true ) ;
for ( third_derivatives_t : : const_iterator it = hessian_params_derivatives . begin ( ) ;
it ! = hessian_params_derivatives . end ( ) ; + + it )
it - > second - > computeTemporaryTerms ( reference_count , params_derivs_temporary_terms , true ) ;
}
2013-10-29 11:47:59 +01:00
bool ModelTree : : isNonstationary ( int symb_id ) const
{
return ( nonstationary_symbols_map . find ( symb_id )
! = nonstationary_symbols_map . end ( ) ) ;
}