29 #ifndef __DirScalingMatrix_H
30 #define __DirScalingMatrix_H 1
35 #include <boost/numeric/ublas/vector.hpp>
36 #include <boost/numeric/ublas/banded.hpp>
37 #include <boost/numeric/ublas/io.hpp>
41 using namespace boost::numeric::ublas;
50 template<
typename NumType>
60 typedef enum mode_type { NO_JACOBIAN, WITH_JACOBIAN };
62 typedef NumType value_type;
64 typedef vector<value_type> vector_type;
65 typedef banded_matrix<value_type> matrix_type;
80 M_trust_region_active(
false )
87 M_lb_ub( M_ub - M_lb ),
88 M_value( __lb.size(), __lb.size(), 0, 0 ),
89 M_jacobian( __lb.size(), __lb.size(), 0, 0 ),
90 M_trust_region_active(
false )
92 GST_SMART_ASSERT( __lb.size() == __ub.size() )( __lb )( __ub )(
"inconsistent bounds definition" );
93 GST_SMART_ASSERT( *std::min_element( M_lb_ub.begin(), M_lb_ub.end() ) >= 0 )
94 ( M_lb )( M_ub )(
"lower and upper bounds are not properly defined" );
114 value_type zeta( vector_type
const& __x )
const;
115 value_type zeta()
const
120 matrix_type
const& operator()()
const
125 matrix_type
const& jacobian()
const
130 bool isTrustRegionActive()
const
132 return M_trust_region_active;
140 void update( value_type
const&, vector_type
const&, vector_type
const&, mode_type = WITH_JACOBIAN );
142 void setBounds( vector_type
const& __lb, vector_type
const& __up )
144 GST_SMART_ASSERT( __lb.size() == __up.size() )( __lb )( __up )(
"inconsistent bounds definition" );
147 M_lb_ub = __up - __lb;
149 GST_SMART_ASSERT( *std::min_element( M_lb_ub.begin(), M_lb_ub.end() ) >= 0 )
150 ( M_lb )( M_ub )(
"lower and upper bounds are not properly defined" );
164 vector_type distanceToLB( vector_type
const& __x )
const
166 GST_SMART_ASSERT( __x.size() == M_lb.size() )( __x )( M_lb )(
"inconsistent bounds definition" );
170 vector_type distanceToUB( vector_type
const& __x )
const
172 GST_SMART_ASSERT( __x.size() == M_ub.size() )( __x )( M_ub )(
"inconsistent bounds definition" );
183 matrix_type M_jacobian;
185 mutable value_type M_zeta;
187 bool M_trust_region_active;
190 template<
typename NumType>
191 typename DirScalingMatrix<NumType>::value_type
196 M_zeta = std::max( M_zeta, std::max( ublas::norm_inf( distanceToLB( __x ) ),
197 ublas::norm_inf( distanceToUB( __x ) ) ) );
200 template<
typename NumType>
202 DirScalingMatrix<NumType>::update( value_type
const& __Delta,
203 vector_type
const& __x,
204 vector_type
const& __s,
207 GST_SMART_ASSERT( __x.size() == M_lb.size() )( __x )( M_lb )(
"inconsistent bounds definition" );
208 GST_SMART_ASSERT( __x.size() == M_ub.size() )( __x )( M_ub )(
"inconsistent bounds definition" );
210 M_value.resize( __x.size(), __x.size(), 0, 0 );
211 M_jacobian.resize( __x.size(), __x.size(), 0, 0 );
214 M_zeta = zeta( __x );
216 vector_type __dl = distanceToLB( __x );
217 vector_type __du = distanceToUB( __x );
219 if ( M_zeta * std::min( norm_inf( __dl ), norm_inf( __du ) ) > __Delta )
222 M_value = identity_matrix<value_type>( M_value.size1(), M_value.size2() );
227 M_trust_region_active =
true;
229 for (
size_t __i = 0; __i < __x.size(); ++__i )
231 if ( __s ( __i ) < 0 )
232 M_value ( __i, __i ) = M_zeta * std::min( 1. , __dl( __i ) ) / __Delta;
235 M_value ( __i, __i ) = M_zeta * std::min( 1. , __du( __i ) ) / __Delta;
239 if ( __mode == WITH_JACOBIAN )
241 for (
size_t i = 0; i < __x.size(); i++ )
243 if ( ( __s( i ) < 0 ) && ( __x( i ) < M_lb( i ) + __Delta ) )
245 M_jacobian ( i, i ) = M_zeta / __Delta;
248 else if ( ( __s ( i ) > 0 ) && ( __x ( i ) > M_ub( i ) - __Delta ) )
250 M_jacobian ( i, i ) = -M_zeta / __Delta;
255 M_jacobian ( i, i ) = 0;
Definition: solverlinear.hpp:33
implements the directional Scaling Matrix for directionally-scaled trust region
Definition: dirscalingmatrix.hpp:51
ElementType element_div(ElementType const &v1, ElementType const &v2)
Definition: functionspace.hpp:6268