source: trunk/cppad/local/checkpoint.hpp @ 3301

Last change on this file since 3301 was 3301, checked in by bradbell, 6 years ago

merge in multiple forward direcitons from branches/forward_dir

  • Property svn:keywords set to Id
File size: 14.5 KB
Line 
1/* $Id: checkpoint.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_CHECKPOINT_INCLUDED
3# define CPPAD_CHECKPOINT_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16namespace CppAD { // BEGIN_CPPAD_NAMESPACE
17/*!
18\file checkpoint.hpp
19defining checkpoint functions.
20*/
21
22/*
23$begin checkpoint$$
24$spell
25        cppad.hpp
26        CppAD
27        checkpoint
28        checkpointing
29        algo
30        afun
31        const
32$$
33
34$section Checkpointing Functions$$
35$index function, checkpoint$$
36$index checkpoint, function$$
37
38$head Syntax$$
39$codei%checkpoint<%Base%> %afun%(%name%, %algo%, %ax%, %ay%)
40%afun%.option(%option_value%)
41%algo%(%ax%, %ay%)
42%afun%(%ax%, %ay%)
43checkpoint<%Base%>::clear()%$$
44
45$head Purpose$$
46You can reduce the size of the tape and memory required for AD by
47checkpointing functions of the form $latex y = f(x)$$ where
48$latex f : B^n \rightarrow B^m$$.
49
50$head Method$$
51The $code checkpoint$$ class is derived from $code atomic_base$$
52and makes this easy.
53It implements all the $code atomic_base$$
54$cref/virtual functions/atomic_base/Virtual Functions/$$
55and hence its source code $code cppad/local/checkpoint.hpp$$
56provides an example implementation of $cref atomic_base$$.
57The difference is that $code checkpoint.hpp$$ uses AD
58instead of user provided derivatives.
59
60$head constructor$$
61The constructor
62$codei%
63        checkpoint<%Base%> %afun%(%name%, %algo%, %ax%, %ay%)
64%$$
65cannot be called in $cref/parallel/ta_in_parallel/$$ mode.
66In addition, you cannot currently be recording
67$codei%AD<%Base%>%$$ operations when the constructor is called.
68This class is implemented as a derived class of
69$cref/atomic_base/atomic_ctor/atomic_base/$$ and hence
70some of its error message will refer to $code atomic_base$$.
71
72$head Base$$
73The type $icode Base$$ specifies the base type for AD operations.
74
75$head ADVector$$
76The type $icode ADVector$$ must be a
77$cref/simple vector class/SimpleVector/$$ with elements of type
78$codei%AD<%Base%>%$$.
79
80$head name$$
81This $icode checkpoint$$ constructor argument has prototype
82$codei%
83        const char* %name%
84%$$
85It is the name used for error reporting.
86The suggested value for $icode name$$ is $icode afun$$; i.e.,
87the same name as used for the function.
88
89$head ax$$
90This argument has prototype
91$codei%
92        const %ADVector%& %ax%
93%$$
94and size must be equal to $icode n$$.
95It specifies vector $latex x \in B^n$$
96at which an $codei%AD<%Base%>%$$ version of
97$latex y = f(x)$$ is to be evaluated.
98
99$head ay$$
100This argument has prototype
101$codei%
102        %ADVector%& %ay%
103%$$
104Its input size must be equal to $icode m$$ and does not change.
105The input values of its elements do not matter.
106Upon return, it is an $codei%AD<%Base%>%$$ version of
107$latex y = f(x)$$.
108
109$head option$$
110The $code option$$ syntax can be used to set the type of sparsity
111pattern used by $icode afun$$.
112This is an $codei%atomic_base<%Base%>%$$ function and its documentation
113can be found at $cref atomic_option$$.
114
115$head algo$$
116The type of $icode algo$$ is arbitrary, except for the fact that
117the syntax
118$codei%
119        %algo%(%ax%, %ay%)
120%$$
121must evaluate the function $latex y = f(x)$$ using
122$codei%AD<%Base%>%$$ operations.
123In addition, we assume that the
124$cref/operation sequence/glossary/Operation/Sequence/$$
125does not depend on the value of $icode ax$$.
126
127$head afun$$
128Given $icode ax$$ it computes the corresponding value of $icode ay$$
129using the operation sequence corresponding to $icode algo$$.
130If $codei%AD<%Base%>%$$ operations are being recorded,
131it enters the computation as single operation in the recording
132see $cref/start recording/Independent/Start Recording/$$.
133(Currently each use of $icode afun$$ actually corresponds to
134$icode%m%+%n%+2%$$ operations and creates $icode m$$ new variables,
135but this is not part of the CppAD specifications and my change.)
136
137$head clear$$
138The $code atomic_base$$ class holds onto static work space in order to
139increase speed by avoiding system memory allocation calls.
140This call makes to work space $cref/available/ta_available/$$ to
141for other uses by the same thread.
142This should be called when you are done using the
143user atomic functions for a specific value of $icode Base$$.
144
145$subhead Restriction$$
146The $code clear$$ routine cannot be called
147while in $cref/parallel/ta_in_parallel/$$ execution mode.
148
149$children%
150        example/atomic/checkpoint.cpp
151%$$
152$head Example$$
153The file $cref checkpoint.cpp$$ contains an example and test
154of these operations.
155It returns true if it succeeds and false if it fails.
156
157$end
158*/
159template <class Base>
160class checkpoint : public atomic_base<Base> {
161private:
162        ADFun<Base> f_;
163public:
164        /*!
165        Constructor of a checkpoint object
166
167        \param name [in]
168        is the user's name for the AD version of this atomic operation.
169
170        \param algo [in/out]
171        user routine that compute AD function values
172        (not const because state may change during evaluation).
173
174        \param ax [in]
175        argument value where algo operation sequence is taped.
176
177        \param ay [out]
178        function value at specified argument value.
179        */
180        template <class Algo, class ADVector>
181        checkpoint(const char* name, 
182                Algo& algo, const ADVector& ax, ADVector& ay)
183        : atomic_base<Base>(name)
184        {       CheckSimpleVector< CppAD::AD<Base> , ADVector>();
185
186                // make a copy of ax because Independent modifies AD information
187                ADVector x_tmp(ax);
188                // delcare x_tmp as the independent variables
189                Independent(x_tmp);
190                // record mapping from x_tmp to ay
191                algo(x_tmp, ay); 
192                // create function f_ : x -> y
193                f_.Dependent(ay);
194                // suppress checking for nan in f_ results
195                // (see optimize documentation for atomic functions)
196                f_.check_for_nan(false);
197                // now optimize (we expect to use this function many times).
198                f_.optimize();
199        }
200        /*!
201        Implement the user call to <tt>afun(ax, ay)</tt>.
202       
203        \tparam ADVector
204        A simple vector class with elements of type <code>AD<Base></code>.
205       
206        \param id
207        optional parameter which must be zero if present.
208       
209        \param ax
210        is the argument vector for this call,
211        <tt>ax.size()</tt> determines the number of arguments.
212       
213        \param ay
214        is the result vector for this call,
215        <tt>ay.size()</tt> determines the number of results.
216        */
217        template <class ADVector>
218        void operator()(const ADVector& ax, ADVector& ay, size_t id = 0)
219        {       CPPAD_ASSERT_KNOWN(
220                        id == 0,
221                        "checkpoint: id is non-zero in afun(ax, ay, id)"
222                );
223                this->atomic_base<Base>::operator()(ax, ay, id);
224        }
225        /*!
226        Link from user_atomic to forward mode
227
228        \copydetails atomic_base::forward
229        */
230        virtual bool forward(
231                size_t                    p ,
232                size_t                    q ,
233                const vector<bool>&      vx , 
234                      vector<bool>&      vy , 
235                const vector<Base>&      tx ,
236                      vector<Base>&      ty )
237        {
238                CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
239                CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
240                CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
241                size_t n = tx.size() / (q+1);
242                size_t m = ty.size() / (q+1);
243                bool ok  = true;       
244                size_t i, j;
245
246                // 2DO: test both forward and reverse vy information
247                if( vx.size() > 0 )
248                {       //Compute Jacobian sparsity pattern.
249                        vector< std::set<size_t> > s(m);
250                        if( n <= m )
251                        {       vector< std::set<size_t> > r(n);
252                                for(j = 0; j < n; j++)
253                                        r[j].insert(j);
254                                s = f_.ForSparseJac(n, r);
255                        }
256                        else
257                        {       vector< std::set<size_t> > r(m);
258                                for(i = 0; i < m; i++)
259                                        r[i].insert(i);
260                                s = f_.RevSparseJac(m, r);
261                        }
262                        std::set<size_t>::const_iterator itr;
263                        for(i = 0; i < m; i++)
264                        {       vy[i] = false;
265                                for(itr = s[i].begin(); itr != s[i].end(); itr++)
266                                {       j = *itr;
267                                        assert( j < n );
268                                        // y[i] depends on the value of x[j]
269                                        vy[i] |= vx[j];
270                                }
271                        }
272                }
273                ty = f_.Forward(q, tx);
274
275                // no longer need the Taylor coefficients in f_
276                // (have to reconstruct them every time)
277                size_t c = 0;
278                size_t r = 0;
279                f_.capacity_order(c, r);
280                return ok;
281        }
282        /*!
283        Link from user_atomic to reverse mode
284
285        \copydetails atomic_base::reverse
286        */
287        virtual bool reverse(
288                size_t                    q  ,
289                const vector<Base>&       tx ,
290                const vector<Base>&       ty ,
291                      vector<Base>&       px ,
292                const vector<Base>&       py )
293        {
294                CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
295                CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
296                CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
297                bool ok  = true;       
298
299                // put proper forward mode coefficients in f_
300# ifdef NDEBUG
301                f_.Forward(q, tx);
302# else
303                size_t n = tx.size() / (q+1);
304                size_t m = ty.size() / (q+1);
305                CPPAD_ASSERT_UNKNOWN( px.size() == n * (q+1) );
306                CPPAD_ASSERT_UNKNOWN( py.size() == m * (q+1) );
307                size_t i, j, k;
308                //
309                vector<Base> check_ty = f_.Forward(q, tx);
310                for(i = 0; i < m; i++)
311                {       for(k = 0; k <= q; k++)
312                        {       j = i * (q+1) + k;
313                                CPPAD_ASSERT_UNKNOWN( check_ty[j] == ty[j] );
314                        }
315                }
316# endif
317                // now can run reverse mode
318                px = f_.Reverse(q+1, py);
319
320                // no longer need the Taylor coefficients in f_
321                // (have to reconstruct them every time)
322                size_t c = 0;
323                size_t r = 0;
324                f_.capacity_order(c, r);
325                return ok;
326        }
327        /*!
328        Link from user_atomic to forward sparse Jacobian
329
330        \copydetails atomic_base::for_sparse_jac
331        */
332        virtual bool for_sparse_jac(
333                size_t                                  q  ,
334                const vector< std::set<size_t> >&       r  ,
335                      vector< std::set<size_t> >&       s  )
336        {
337                bool ok = true;
338                s = f_.ForSparseJac(q, r);
339
340                // no longer need the forward mode sparsity pattern
341                // (have to reconstruct them every time)
342                f_.size_forward_set(0);
343               
344                return ok; 
345        }
346        /*!
347        Link from user_atomic to forward sparse Jacobian
348
349        \copydetails atomic_base::for_sparse_jac
350        */
351        virtual bool for_sparse_jac(
352                size_t                                  q  ,
353                const vector<bool>&                     r  ,
354                      vector<bool>&                     s  )
355        {
356                bool ok = true;
357                s = f_.ForSparseJac(q, r);
358
359                // no longer need the forward mode sparsity pattern
360                // (have to reconstruct them every time)
361                f_.size_forward_bool(0);
362               
363                return ok; 
364        }
365        /*!
366        Link from user_atomic to forward sparse Jacobian
367
368        \copydetails atomic_base::rev_sparse_jac
369        */
370        virtual bool rev_sparse_jac(
371                size_t                                  q  ,
372                const vector< std::set<size_t> >&       rt ,
373                      vector< std::set<size_t> >&       st )
374        {
375                bool ok  = true;
376
377                // compute rt
378                // 2DO: remove need for nz_compare all the time. It is only really
379                // necessary when optimizer calls this member function.
380                bool transpose = true;
381                bool nz_compare = true;
382                st = f_.RevSparseJac(q, rt, transpose, nz_compare);
383
384                return ok; 
385        }
386        /*!
387        Link from user_atomic to forward sparse Jacobian
388
389        \copydetails atomic_base::rev_sparse_jac
390        */
391        virtual bool rev_sparse_jac(
392                size_t                                  q  ,
393                const vector<bool>&                     rt ,
394                      vector<bool>&                     st )
395        {
396                bool ok  = true;
397
398                // compute rt
399                bool transpose  = true;
400                bool nz_compare = true;
401                // 2DO: remove need for nz_compare all the time. It is only really
402                // necessary when optimizer calls this member function.
403                st = f_.RevSparseJac(q, rt, transpose, nz_compare);
404
405                return ok; 
406        }
407        /*!
408        Link from user_atomic to forward sparse Jacobian
409
410        \copydetails atomic_base::rev_sparse_hes
411        */
412        virtual bool rev_sparse_hes(
413                const vector<bool>&                     vx ,
414                const vector<bool>&                     s  ,
415                      vector<bool>&                     t  ,
416                size_t                                  q  ,
417                const vector< std::set<size_t> >&       r  ,
418                const vector< std::set<size_t> >&       u  ,
419                      vector< std::set<size_t> >&       v  )
420        {       size_t n       = v.size();
421                size_t m       = u.size();
422                CPPAD_ASSERT_UNKNOWN( r.size() == v.size() );
423                CPPAD_ASSERT_UNKNOWN( s.size() == m );
424                CPPAD_ASSERT_UNKNOWN( t.size() == n );
425                bool ok        = true;
426                bool transpose = true;
427                std::set<size_t>::const_iterator itr;
428
429                // compute sparsity pattern for T(x) = S(x) * f'(x)
430                t = f_.RevSparseJac(1, s);
431# ifndef NDEBUG
432                for(size_t j = 0; j < n; j++)
433                        CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
434# endif
435
436                // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
437                // U(x) = g''(y) * f'(x) * R
438                // S(x) = g'(y)
439               
440                // compute sparsity pattern for A(x) = f'(x)^T * U(x)
441                vector< std::set<size_t> > a(n);
442                a = f_.RevSparseJac(q, u, transpose);
443
444                // set version of s
445                vector< std::set<size_t> > set_s(1);
446                CPPAD_ASSERT_UNKNOWN( set_s[0].empty() );
447                size_t i;
448                for(i = 0; i < m; i++)
449                        if( s[i] )
450                                set_s[0].insert(i);
451
452                // compute sparsity pattern for H(x) = (S(x) * F)''(x) * R
453                // (store it in v)
454                f_.ForSparseJac(q, r);
455                v = f_.RevSparseHes(q, set_s, transpose);
456
457                // compute sparsity pattern for V(x) = A(x) + H(x)
458                for(i = 0; i < n; i++)
459                {       for(itr = a[i].begin(); itr != a[i].end(); itr++)
460                        {       size_t j = *itr;
461                                CPPAD_ASSERT_UNKNOWN( j < q );
462                                v[i].insert(j);
463                        }
464                }
465
466                // no longer need the forward mode sparsity pattern
467                // (have to reconstruct them every time)
468                f_.size_forward_set(0);
469
470                return ok;
471        }
472        /*!
473        Link from user_atomic to forward sparse Jacobian
474
475        \copydetails atomic_base::rev_sparse_hes
476        */
477        virtual bool rev_sparse_hes(
478                const vector<bool>&                     vx ,
479                const vector<bool>&                     s  ,
480                      vector<bool>&                     t  ,
481                size_t                                  q  ,
482                const vector<bool>&                     r  ,
483                const vector<bool>&                     u  ,
484                      vector<bool>&                     v  )
485        {
486                CPPAD_ASSERT_UNKNOWN( r.size() == v.size() );
487                CPPAD_ASSERT_UNKNOWN( s.size() == u.size() / q );
488                CPPAD_ASSERT_UNKNOWN( t.size() == v.size() / q );
489                size_t n       = t.size();
490                bool ok        = true;
491                bool transpose = true;
492                std::set<size_t>::const_iterator itr;
493                size_t i, j;
494
495                // compute sparsity pattern for T(x) = S(x) * f'(x)
496                t = f_.RevSparseJac(1, s);
497# ifndef NDEBUG
498                for(j = 0; j < n; j++)
499                        CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
500# endif
501
502                // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
503                // U(x) = g''(y) * f'(x) * R
504                // S(x) = g'(y)
505
506                // compute sparsity pattern for A(x) = f'(x)^T * U(x)
507                vector<bool> a(n * q);
508                a = f_.RevSparseJac(q, u, transpose);
509
510                // compute sparsity pattern for H(x) =(S(x) * F)''(x) * R
511                // (store it in v)
512                f_.ForSparseJac(q, r);
513                v = f_.RevSparseHes(q, s, transpose);
514
515                // compute sparsity pattern for V(x) = A(x) + H(x)
516                for(i = 0; i < n; i++)
517                {       for(j = 0; j < q; j++)
518                                v[ i * q + j ] |= a[ i * q + j];
519                }
520
521                // no longer need the forward mode sparsity pattern
522                // (have to reconstruct them every time)
523                f_.size_forward_set(0);
524
525                return ok;
526        }
527};
528
529} // END_CPPAD_NAMESPACE
530# endif
Note: See TracBrowser for help on using the repository browser.