29 #ifndef __ADTypeOrder2_H
30 #define __ADTypeOrder2_H 1
32 #include <boost/numeric/ublas/symmetric.hpp>
33 #include <boost/numeric/ublas/vector.hpp>
35 #include <boost/numeric/ublas/io.hpp>
69 template<
typename T,
int Nvar,
int Var>
70 class ADType<T, Nvar, 2, Var>
76 enum { depvar = Var };
78 typedef ADVariable<Var> variable_type;
82 typedef boost::numeric::ublas::vector<value_type> gradient_type;
85 typedef std::vector<bool> deps_type;
88 typedef boost::numeric::ublas::symmetric_matrix<value_type, boost::numeric::ublas::upper> hessian_type;
90 typedef ADType<T,Nvar,2, Var> This;
92 typedef std::set<int> dependency_type;
94 template<
typename NumT,
int NumVar,
int Order,
int VarNum>
friend class ADType;
100 ADType( value_type
val = value_type( 0 ) )
104 M_hessian( nvar, nvar ),
108 M_grad = boost::numeric::ublas::zero_vector<double>( nvar );
109 M_hessian = boost::numeric::ublas::zero_matrix<double>( nvar, nvar );
113 M_grad[ depvar ] = 1.;
114 __dep[ depvar ] =
true;
119 ADType( ADType<T,Nvar,2,VarNum>
const& sad )
122 M_grad( sad.M_grad ),
123 M_hessian( sad.M_hessian ),
129 template<
typename ExprT>
130 ADType (
const ADExpr<ExprT>& expr )
134 M_hessian( nvar, nvar ),
147 value_type value()
const
153 bool deps(
int __i )
const
160 value_type grad(
int __i )
const
163 return M_grad( __i );
166 gradient_type grad()
const
177 value_type
hessian(
int __i,
int __j )
const
180 return M_hessian( __i, __j );
192 bool& deps(
int __i )
198 value_type& grad(
int __i )
201 return M_grad( __i );
204 gradient_type& grad()
215 value_type&
hessian(
int __i,
int __j )
218 return M_hessian( __i, __j );
228 template <
class ExprT> This&
operator=(
const ADExpr<ExprT>& expr );
231 ADExpr< ADUnaryPlus< This > > operator+ ()
const
233 typedef ADUnaryPlus<This> expr_t;
234 return ADExpr<expr_t> ( expr_t ( *
this ) );
238 ADExpr< ADUnaryMinus< This > > operator- ()
const
240 typedef ADUnaryMinus<This> expr_t;
241 return ADExpr<expr_t> ( expr_t ( *
this ) );
244 #define AD_UNARY_OP( op ) \
245 This& operator op ( value_type val ) \
255 This& operator += ( This
const& sad )
258 M_grad += sad.M_grad;
259 M_hessian += sad.M_hessian;
262 This& operator -= ( This
const& sad )
265 M_grad -= sad.M_grad;
266 M_hessian -= sad.M_hessian;
271 This& operator *= ( This
const& sad )
274 M_grad = M_grad * sad.M_val + M_val * sad.M_grad;
278 This& operator /= ( This
const& sad )
281 M_grad = ( M_grad * sad.M_val - M_val * sad.M_grad ) / ( sad.M_val * sad.M_val );
285 template<
typename Expr>
286 This& operator += ( ADExpr<Expr>
const& sad )
292 template<
typename Expr>
293 This& operator -= ( ADExpr<Expr>
const& sad )
299 template<
typename Expr>
300 This& operator *= ( ADExpr<Expr>
const& sad )
306 template<
typename Expr>
307 This& operator /= ( ADExpr<Expr>
const& sad )
319 gradient_type M_grad;
320 hessian_type M_hessian;
323 dependency_type M_deps;
326 template<
typename T,
int Nvar,
int Var>
327 ADType<T, Nvar, 2, Var>&
332 M_grad = boost::numeric::ublas::zero_vector<double>( nvar );
333 M_hessian = boost::numeric::ublas::zero_matrix<double>( nvar, nvar );;
337 template<
typename T,
int Nvar,
int Var>
338 ADType<T, Nvar, 2, Var>&
344 M_hessian = sad.M_hessian;
347 template<
typename T,
int Nvar,
int Var>
348 template <
class ExprT>
349 ADType<T,Nvar, 2, Var> &
354 M_val = expr.value();
356 typedef typename Feel::STL::SNoDuplicates<typename Feel::STL::SFlatten<typename ADExpr<ExprT>::Var>::Result>::Result Expr;
357 typedef typename Feel::Meta::IF<Feel::STL::SIndexOf<Expr,Variable<-1> >::value == -1,
360 >::Result ForEachExpr;
363 M_grad = boost::numeric::ublas::zero_vector<double>( nvar );
364 ForEachExpr( *
this, expr );
369 M_val = expr.value();
373 for (
int __i = 0; __i < nvar; ++__i )
375 __dep( __i ) = expr.deps( __i );
380 M_deps.insert( __i );
384 if ( ! M_deps.empty() && M_deps.size() < nvar )
389 dependency_type::iterator __begin = M_deps.begin();
390 dependency_type::iterator __end = M_deps.end();
392 dependency_type::iterator __it = __begin;
393 dependency_type::iterator __it2;
395 while ( __it != __end )
397 M_grad( *__it )= expr.grad( *__it );
400 while ( __it2 != __it )
402 M_hessian( *__it, *__it2 ) = expr.hessian( *__it, *__it2 );
408 M_hessian( *__it, *__it ) = expr.hessian( *__it, *__it );
417 for (
int __i = 0; __i < nvar; ++__i )
419 M_grad( __i )= expr.grad( __i );
420 M_hessian( __i, __i ) = expr.hessian( __i, __i );
422 for (
int __j = 0; __j < __i; ++__j )
424 M_hessian( __i, __j ) = expr.hessian( __i, __j );
441 template <
class T,
int Nvar,
int Var>
443 operator << ( std::ostream& os, const Feel::ADType<T, Nvar, 2, Var>& a )
445 os.setf( std::ios::fixed,std::ios::floatfield );
447 os <<
"value = " << a.value() <<
" \n"
448 <<
"gradient = " << a.grad() <<
"\n"
449 <<
"hessian = " << a.hessian() <<
"\n";
452 for (
int __i = 0 ; __i < Nvar; ++__i )
454 for (
int __j = 0 ; __j < Nvar; ++__j )
455 os << a.hessian( __i, __j ) <<
" ";