Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
adtypeorder2.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2008-02-14
7 
8  Copyright (C) 2008 Université Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __ADTypeOrder2_H
30 #define __ADTypeOrder2_H 1
31 
32 #include <boost/numeric/ublas/symmetric.hpp>
33 #include <boost/numeric/ublas/vector.hpp>
34 
35 #include <boost/numeric/ublas/io.hpp>
36 
37 
38 namespace Feel
39 {
69 template<typename T, int Nvar, int Var>
70 class ADType<T, Nvar, 2, Var>
71 {
72 public:
73 
74  enum { nvar = Nvar };
75  enum { order = 2 };
76  enum { depvar = Var };
77 
78  typedef ADVariable<Var> variable_type;
79  //typedef typename Feel::STL::SFlatten<typename SListGenerator<Nvar,ADVariable>::list_type>::Result variable_list_type;
80 
81  typedef T value_type;
82  typedef boost::numeric::ublas::vector<value_type> gradient_type;
83  //typedef typename STinyVector<value_type, Nvar>::type gradient_type;
84  //typedef typename STinyVector<bool, Nvar>::type deps_type;
85  typedef std::vector<bool> deps_type;
86 
87  //typedef typename STinyVector<typename STinyVector<value_type,Nvar>::type, Nvar>::type hessian_type;
88  typedef boost::numeric::ublas::symmetric_matrix<value_type, boost::numeric::ublas::upper> hessian_type;
89 
90  typedef ADType<T,Nvar,2, Var> This;
91 
92  typedef std::set<int> dependency_type;
93 
94  template<typename NumT, int NumVar, int Order, int VarNum> friend class ADType;
95 
98  //@ {
99 
100  ADType( value_type val = value_type( 0 ) )
101  :
102  M_val( val ),
103  M_grad( nvar ),
104  M_hessian( nvar, nvar ),
105  __dep( false ),
106  M_deps()
107  {
108  M_grad = boost::numeric::ublas::zero_vector<double>( nvar );
109  M_hessian = boost::numeric::ublas::zero_matrix<double>( nvar, nvar );
110 
111  if ( depvar != -1 )
112  {
113  M_grad[ depvar ] = 1.;
114  __dep[ depvar ] = true;
115  }
116  }
117 
118  template<int VarNum>
119  ADType( ADType<T,Nvar,2,VarNum> const& sad )
120  :
121  M_val( sad.M_val ),
122  M_grad( sad.M_grad ),
123  M_hessian( sad.M_hessian ),
124  __dep( sad.__dep ),
125  M_deps()
126  {
127 
128  }
129  template<typename ExprT>
130  ADType ( const ADExpr<ExprT>& expr )
131  :
132  M_val( 0 ),
133  M_grad( nvar ),
134  M_hessian( nvar, nvar ),
135  __dep( false ),
136  M_deps()
137  {
138  *this = expr;
139  }
141 
145 
147  value_type value() const
148  {
149  return M_val;
150  }
151 
153  bool deps( int __i ) const
154  {
155  //ENSURE( __i >= 0 && __i < nvar );
156  return __dep[ __i ];
157  }
158 
160  value_type grad( int __i ) const
161  {
162  //ENSURE( __i >= 0 && __i < nvar );
163  return M_grad( __i );
164  }
165 
166  gradient_type grad() const
167  {
168  return M_grad;
169  }
170 
171  hessian_type hessian() const
172  {
173  return M_hessian;
174  }
175 
177  value_type hessian( int __i, int __j ) const
178  {
179  //ENSURE( __i >= 0 && __i < nvar && __j >= 0 && __j < nvar );
180  return M_hessian( __i, __j );
181  }
183 
187  value_type& value()
188  {
189  return M_val;
190  }
192  bool& deps( int __i )
193  {
194  //ENSURE( __i >= 0 && __i < nvar );
195  return __dep[ __i ];
196  }
198  value_type& grad( int __i )
199  {
200  //ENSURE( __i >= 0 && __i < nvar );
201  return M_grad( __i );
202  }
203 
204  gradient_type& grad()
205  {
206  return M_grad;
207  }
208 
209  hessian_type& hessian()
210  {
211  return M_hessian;
212  }
213 
215  value_type& hessian( int __i, int __j )
216  {
217  //ENSURE( __i >= 0 && __i < nvar && __j >= 0 && __j < nvar );
218  return M_hessian( __i, __j );
219  }
221 
225 
226  This& operator=( T const& );
227  This& operator=( This const& );
228  template <class ExprT> This& operator=( const ADExpr<ExprT>& expr );
229 
231  ADExpr< ADUnaryPlus< This > > operator+ () const
232  {
233  typedef ADUnaryPlus<This> expr_t;
234  return ADExpr<expr_t> ( expr_t ( *this ) );
235  }
236 
238  ADExpr< ADUnaryMinus< This > > operator- () const
239  {
240  typedef ADUnaryMinus<This> expr_t;
241  return ADExpr<expr_t> ( expr_t ( *this ) );
242  }
243 
244 #define AD_UNARY_OP( op ) \
245  This& operator op ( value_type val ) \
246  { \
247  M_val op val; \
248  return *this; \
249  }
250  AD_UNARY_OP( += );
251  AD_UNARY_OP( -= );
252  AD_UNARY_OP( *= );
253  AD_UNARY_OP( /= );
254 
255  This& operator += ( This const& sad )
256  {
257  M_val += sad.M_val;
258  M_grad += sad.M_grad;
259  M_hessian += sad.M_hessian;
260  return *this;
261  }
262  This& operator -= ( This const& sad )
263  {
264  M_val -= sad.M_val;
265  M_grad -= sad.M_grad;
266  M_hessian -= sad.M_hessian;
267  return *this;
268  }
269 
270 
271  This& operator *= ( This const& sad )
272  {
273  M_val *= sad.M_val;
274  M_grad = M_grad * sad.M_val + M_val * sad.M_grad;
275  return *this;
276  }
277 
278  This& operator /= ( This const& sad )
279  {
280  M_val /= sad.M_val;
281  M_grad = ( M_grad * sad.M_val - M_val * sad.M_grad ) / ( sad.M_val * sad.M_val );
282  return *this;
283  }
284 
285  template<typename Expr>
286  This& operator += ( ADExpr<Expr> const& sad )
287  {
288  *this = *this + sad;
289  return *this;
290  }
291 
292  template<typename Expr>
293  This& operator -= ( ADExpr<Expr> const& sad )
294  {
295  *this = *this - sad;
296  return *this;
297  }
298 
299  template<typename Expr>
300  This& operator *= ( ADExpr<Expr> const& sad )
301  {
302  *this = *this * sad;
303  return *this;
304  }
305 
306  template<typename Expr>
307  This& operator /= ( ADExpr<Expr> const& sad )
308  {
309  *this = *this / sad;
310  return *this;
311  }
313 
314 
315 protected:
316 
317 
318  value_type M_val;
319  gradient_type M_grad;
320  hessian_type M_hessian;
321 
322  deps_type __dep;
323  dependency_type M_deps;
324 };
325 
326 template<typename T,int Nvar, int Var>
327 ADType<T, Nvar, 2, Var>&
328 ADType<T, Nvar, 2, Var>::operator=( value_type const& val )
329 {
330  __dep = false;
331  M_val = val;
332  M_grad = boost::numeric::ublas::zero_vector<double>( nvar );
333  M_hessian = boost::numeric::ublas::zero_matrix<double>( nvar, nvar );;
334  return *this;
335 }
336 
337 template<typename T,int Nvar, int Var>
338 ADType<T, Nvar, 2, Var>&
339 ADType<T, Nvar, 2, Var>::operator=( This const& sad )
340 {
341  __dep = sad.__dep;
342  M_val = sad.M_val;
343  M_grad = sad.M_grad;
344  M_hessian = sad.M_hessian;
345  return *this;
346 }
347 template<typename T,int Nvar, int Var>
348 template <class ExprT>
349 ADType<T,Nvar, 2, Var> &
350 ADType<T,Nvar, 2, Var>::operator=( const ADExpr<ExprT>& expr )
351 {
352 
353 #if 0
354  M_val = expr.value();
355 
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,
358  FOR_EACH<Expr>,
359  depManager
360  >::Result ForEachExpr;
361 
362  // should generate the complement to VarList of Expr
363  M_grad = boost::numeric::ublas::zero_vector<double>( nvar );
364  ForEachExpr( *this, expr );
365 
366  //FOR_EACH<FOR_EACH<Expr> >( M_grad, expr );
367  return *this;
368 #endif
369  M_val = expr.value();
370 #if 0
371  M_deps.clear();
372 
373  for ( int __i = 0; __i < nvar; ++__i )
374  {
375  __dep( __i ) = expr.deps( __i );
376 
377  if ( __dep( __i ) )
378  {
379  //std::cerr << "depends on " << __i << "\n";
380  M_deps.insert( __i );
381  }
382  }
383 
384  if ( ! M_deps.empty() && M_deps.size() < nvar )
385  {
386  //M_grad = 0;
387  //M_hessian = M_grad;
388 
389  dependency_type::iterator __begin = M_deps.begin();
390  dependency_type::iterator __end = M_deps.end();
391 
392  dependency_type::iterator __it = __begin;
393  dependency_type::iterator __it2;
394 
395  while ( __it != __end )
396  {
397  M_grad( *__it )= expr.grad( *__it );
398  __it2 = __begin;
399 
400  while ( __it2 != __it )
401  {
402  M_hessian( *__it, *__it2 ) = expr.hessian( *__it, *__it2 );
403  //M_hessian( *__it2, *__it ) = M_hessian( *__it][ *__it2 ];
404 
405  ++__it2;
406  }
407 
408  M_hessian( *__it, *__it ) = expr.hessian( *__it, *__it );
409  ++__it;
410  }
411  }
412 
413  else
414  {
415 #endif
416 
417  for ( int __i = 0; __i < nvar; ++__i )
418  {
419  M_grad( __i )= expr.grad( __i );
420  M_hessian( __i, __i ) = expr.hessian( __i, __i );
421 
422  for ( int __j = 0; __j < __i; ++__j )
423  {
424  M_hessian( __i, __j ) = expr.hessian( __i, __j );
425  //M_hessian( __j, __i ) = M_hessian[ __i ][ __j ];
426  }
427  }
428 
429 #if 0
430  }
431 
432 #endif
433  return *this;
434 
435 }
436 
437 } // Feel
438 
439 //------------------------------- AD ostream operator ------------------------------------------
440 
441 template <class T, int Nvar, int Var>
442 std::ostream&
443 operator << ( std::ostream& os, const Feel::ADType<T, Nvar, 2, Var>& a )
444 {
445  os.setf( std::ios::fixed,std::ios::floatfield );
446  //os.width(12);
447  os << "value = " << a.value() << " \n"
448  << "gradient = " << a.grad() << "\n"
449  << "hessian = " << a.hessian() << "\n";
450 #if 0
451 
452  for ( int __i = 0 ; __i < Nvar; ++__i )
453  {
454  for ( int __j = 0 ; __j < Nvar; ++__j )
455  os << a.hessian( __i, __j ) << " ";
456 
457  os << "\n";
458  }
459 
460  os << "]\n";
461 #endif
462  return os;
463 }
464 
465 #endif

Generated on Sun Oct 20 2013 08:24:54 for Feel++ by doxygen 1.8.4