source: branches/opt_cond_exp/cppad/local/optimize.hpp @ 2955

Last change on this file since 2955 was 2955, checked in by bradbell, 7 years ago

op_code.hpp: fix bug in printing atomic op codes.
optimize.hpp: improve comments about cummulative summation.

  • Property svn:keywords set to Id
File size: 61.5 KB
Line 
1/* $Id: optimize.hpp 2955 2013-10-16 14:56:16Z bradbell $ */
2# ifndef CPPAD_OPTIMIZE_INCLUDED
3# define CPPAD_OPTIMIZE_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 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
16/*
17$begin optimize$$
18$spell
19        Taylor
20        var
21        CppAD
22        cppad
23        std
24        CondExpEq
25$$
26
27$section Optimize an ADFun Object Tape$$
28
29$index optimize$$
30$index tape, optimize$$
31$index sequence, optimize operations$$
32$index operations, optimize sequence$$
33$index speed, optimize$$
34$index memory, optimize$$
35
36$head Syntax$$
37$icode%f%.optimize()%$$
38
39
40$head Purpose$$
41The operation sequence corresponding to an $cref ADFun$$ object can
42be very large and involve many operations; see the
43size functions in $cref seq_property$$.
44The $icode%f%.optimize%$$ procedure reduces the number of operations,
45and thereby the time and the memory, required to
46compute function and derivative values.
47
48$head f$$
49The object $icode f$$ has prototype
50$codei%
51        ADFun<%Base%> %f%
52%$$
53
54$head Improvements$$
55You can see the reduction in number of variables in the operation sequence
56by calling the function $cref/f.size_var()/seq_property/size_var/$$
57before and after the optimization procedure.
58Given that the optimization procedure takes time,
59it may be faster to skip this optimize procedure and just compute
60derivatives using the original operation sequence.
61
62$subhead Testing$$
63You can run the CppAD $cref/speed/speed_main/$$ tests and see
64the corresponding changes in number of variables and execution time;
65see $cref cppad_test$$.
66
67$head Efficiency$$
68The $code optimize$$ member function
69may greatly reduce the number of variables
70in the operation sequence; see $cref/size_var/seq_property/size_var/$$.
71If a $cref/zero order forward/ForwardZero/$$ calculation is done during
72the construction of $icode f$$, it will require more memory
73and time than required after the optimization procedure.
74In addition, it will need to be redone.
75For this reason, it is more efficient to use
76$codei%
77        ADFun<%Base%> %f%;
78        %f%.Dependent(%x%, %y%);
79        %f%.optimize();
80%$$
81instead of
82$codei%
83        ADFun<%Base%> %f%(%x%, %y%)
84        %f%.optimize();
85%$$
86See the discussion about
87$cref/sequence constructors/FunConstruct/Sequence Constructor/$$.
88
89$head Comparison Operators$$
90Any comparison operators that are in the tape are removed by this operation.
91Hence the return value of $cref CompareChange$$ will always be zero
92for an optimized tape (even if $code NDEBUG$$ is not defined).
93
94$head Atomic Functions$$
95There are some subtitle issue with optimized $cref atomic$$ functions
96$latex v = g(u)$$.
97$list number$$
98The $cref atomic_rev_sparse_jac$$ function is be used to determine
99which components of $icode u$$ affect the dependent variables of $icode f$$.
100Currently this always uses $code std::set<size_t>$$ for the sparsity patterns.
101It should use the current setting of
102$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ for the
103atomic function $latex g(u)$$.
104$lnext
105If $icode%u%[%i%]%$$ does not affect the value of
106the dependent variables for $icode f$$,
107the value of $icode%u%[%i%]%$$ is set to $cref nan$$.
108$lend
109
110
111$head Checking Optimization$$
112$index NDEBUG$$
113If $cref/NDEBUG/Faq/Speed/NDEBUG/$$ is not defined,
114and $cref/f.size_taylor()/size_taylor/$$ is greater than zero,
115a $cref ForwardZero$$ calculation is done using the optimized version
116of $icode f$$ and the results are checked to see that they are
117the same as before.
118If they are not the same, the
119$cref ErrorHandler$$ is called with a known error message
120related to $icode%f%.optimize()%$$.
121
122$head Example$$
123$children%
124        example/optimize.cpp
125%$$
126The file
127$cref optimize.cpp$$
128contains an example and test of this operation.
129It returns true if it succeeds and false otherwise.
130
131$end
132-----------------------------------------------------------------------------
133*/
134# include <stack>
135
136namespace CppAD { // BEGIN_CPPAD_NAMESPACE
137/*!
138\defgroup optimize_hpp optimize.hpp
139\{
140\file optimize.hpp
141Routines for optimizing a tape
142*/
143
144
145/*!
146State for this variable set during reverse sweep.
147*/
148enum optimize_connection_type {
149        /// There is no operation that connects this variable to the
150        /// independent variables.
151        not_connected     ,
152
153        /// There is one or more operations that connects this variable to the
154        /// independent variables.
155        yes_connected      ,
156
157        /// There is only one parrent that connects this variable to the
158        /// independent variables and the parent is a summation operation; i.e.,
159        /// AddvvOp, AddpvOp, SubpvOp, SubvpOp, or SubvvOp.
160        sum_connected      ,
161
162        /// Satisfies the sum_connected assumptions above and in addition
163        /// this variable is the result of summation operator.
164        csum_connected
165};
166
167
168/*!
169Structure used by \c optimize to hold information about one variable.
170in the old operation seqeunce.
171*/
172struct optimize_old_variable {
173        /// Operator for which this variable is the result, \c NumRes(op) > 0.
174        /// Set by the reverse sweep at beginning of optimization.
175        OpCode              op;       
176
177        /// Pointer to first argument (child) for this operator.
178        /// Set by the reverse sweep at beginning of optimization.
179        const addr_t*       arg;
180
181        /// How is this variable connected to the independent variables
182        optimize_connection_type connect_type; 
183
184        /// Set during forward sweep to the index in the
185        /// new operation sequence corresponding to this old varable.
186        addr_t new_var;
187};
188
189/*!
190Structures used by \c optimize_record_csum
191to hold information about one variable.
192*/
193struct optimize_csum_variable {
194        /// Operator for which this variable is the result, \c NumRes(op) > 0.
195        OpCode              op;       
196
197        /// Pointer to first argument (child) for this operator.
198        /// Set by the reverse sweep at beginning of optimization.
199        const addr_t*       arg;
200
201        /// Is this variable added to the summation
202        /// (if not it is subtracted)
203        bool                add;
204};
205
206/*!
207Structure used to pass work space from \c optimize to \c optimize_record_csum
208(so that stacks do not start from zero size every time).
209*/
210struct optimize_csum_stacks {
211        /// stack of operations in the cummulative summation
212        std::stack<struct optimize_csum_variable>   op_stack;
213        /// stack of variables to be added
214        std::stack<size_t >                         add_stack;
215        /// stack of variables to be subtracted
216        std::stack<size_t >                         sub_stack;
217};
218
219
220/*!
221Shared documentation for optimization helper functions (not called).
222
223<!-- define optimize_prototype -->
224\param tape
225is a vector that maps a variable index, in the old operation sequence,
226to an <tt>optimize_old_variable</tt> information record.
227Note that the index for this vector must be greater than or equal zero and
228less than <tt>tape.size()</tt>.
229
230\li <tt>tape[i].op</tt>
231is the operator in the old operation sequence
232corresponding to the old variable index \c i.
233Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
234
235\li <tt>tape[i].arg</tt>
236for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
237is the j-th the argument, in the old operation sequence,
238corresponding to the old variable index \c i.
239Assertion: <tt>tape[i].arg[j] < i</tt>.
240
241\li <tt>tape[i].new_var</tt>
242Suppose
243<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
244and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
245It follows that <tt>tape[k].new_var</tt>
246has alread been set to the variable in the new operation sequence
247corresponding to the old variable index \c k.
248This means that the \c new_var value has been set
249for all the possible arguments that come before \a current.
250
251\param current
252is the index in the old operation sequence for
253the variable corresponding to the result for the current operator.
254Assertions:
255<tt>
256current < tape.size(),
257NumRes( tape[current].op ) > 0.
258</tt>
259
260\param npar
261is the number of parameters corresponding to this operation sequence.
262
263\param par
264is a vector of length \a npar containing the parameters
265for this operation sequence; i.e.,
266given a parameter index \c i, the corresponding parameter value is
267<tt>par[i]</tt>.
268<!-- end optimize_prototype -->
269*/
270template <class Base>
271void optimize_prototype(
272        const CppAD::vector<struct optimize_old_variable>& tape           ,
273        size_t                                             current        ,
274        size_t                                             npar           ,
275        const Base*                                        par            )
276{       CPPAD_ASSERT_UNKNOWN(false); }
277
278/*!
279Check a unary operator for a complete match with a previous operator.
280
281A complete match means that the result of the previous operator
282can be used inplace of the result for current operator.
283
284<!-- replace optimize_prototype -->
285\param tape
286is a vector that maps a variable index, in the old operation sequence,
287to an <tt>optimize_old_variable</tt> information record.
288Note that the index for this vector must be greater than or equal zero and
289less than <tt>tape.size()</tt>.
290
291\li <tt>tape[i].op</tt>
292is the operator in the old operation sequence
293corresponding to the old variable index \c i.
294Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
295
296\li <tt>tape[i].arg</tt>
297for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
298is the j-th the argument, in the old operation sequence,
299corresponding to the old variable index \c i.
300Assertion: <tt>tape[i].arg[j] < i</tt>.
301
302\li <tt>tape[i].new_var</tt>
303Suppose
304<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
305and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
306It follows that <tt>tape[k].new_var</tt>
307has alread been set to the variable in the new operation sequence
308corresponding to the old variable index \c k.
309This means that the \c new_var value has been set
310for all the possible arguments that come before \a current.
311
312\param current
313is the index in the old operation sequence for
314the variable corresponding to the result for the current operator.
315Assertions:
316<tt>
317current < tape.size(),
318NumRes( tape[current].op ) > 0.
319</tt>
320
321\param npar
322is the number of parameters corresponding to this operation sequence.
323
324\param par
325is a vector of length \a npar containing the parameters
326for this operation sequence; i.e.,
327given a parameter index \c i, the corresponding parameter value is
328<tt>par[i]</tt>.
329<!-- end optimize_prototype -->
330
331\param hash_table_var
332is a vector with size CPPAD_HASH_TABLE_SIZE
333that maps a hash code to the corresponding
334variable index in the old operation sequence.
335All the values in this table must be less than \a current.
336
337\param code
338The input value of code does not matter.
339The output value of code is the hash code corresponding to
340this operation in the new operation sequence.
341
342\return
343If the return value is zero,
344no match was found.
345If the return value is greater than zero,
346it is the index of a new variable that can be used to replace the
347old variable.
348
349\par Restrictions:
350NumArg( tape[current].op ) == 1
351*/
352template <class Base>
353size_t optimize_unary_match(
354        const CppAD::vector<struct optimize_old_variable>& tape           ,
355        size_t                                             current        ,
356        size_t                                             npar           ,
357        const Base*                                        par            ,
358        const CppAD::vector<size_t>&                       hash_table_var ,
359        unsigned short&                                    code           )
360{       const addr_t* arg = tape[current].arg;
361        OpCode        op  = tape[current].op;
362        addr_t new_arg[1];
363       
364        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
365        CPPAD_ASSERT_UNKNOWN( NumRes(op) > 0  );
366        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
367        new_arg[0] = tape[arg[0]].new_var;
368        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < current );
369        code = hash_code(
370                op                  , 
371                new_arg             ,
372                npar                ,
373                par
374        );
375        size_t  i               = hash_table_var[code];
376        CPPAD_ASSERT_UNKNOWN( i < current );
377        if( op == tape[i].op )
378        {       size_t k = tape[i].arg[0];
379                CPPAD_ASSERT_UNKNOWN( k < i );
380                if (new_arg[0] == tape[k].new_var )
381                        return tape[i].new_var;
382        }
383        return 0;
384} 
385
386/*!
387Check a binary operator for a complete match with a previous operator,
388
389<!-- replace optimize_prototype -->
390\param tape
391is a vector that maps a variable index, in the old operation sequence,
392to an <tt>optimize_old_variable</tt> information record.
393Note that the index for this vector must be greater than or equal zero and
394less than <tt>tape.size()</tt>.
395
396\li <tt>tape[i].op</tt>
397is the operator in the old operation sequence
398corresponding to the old variable index \c i.
399Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
400
401\li <tt>tape[i].arg</tt>
402for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
403is the j-th the argument, in the old operation sequence,
404corresponding to the old variable index \c i.
405Assertion: <tt>tape[i].arg[j] < i</tt>.
406
407\li <tt>tape[i].new_var</tt>
408Suppose
409<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
410and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
411It follows that <tt>tape[k].new_var</tt>
412has alread been set to the variable in the new operation sequence
413corresponding to the old variable index \c k.
414This means that the \c new_var value has been set
415for all the possible arguments that come before \a current.
416
417\param current
418is the index in the old operation sequence for
419the variable corresponding to the result for the current operator.
420Assertions:
421<tt>
422current < tape.size(),
423NumRes( tape[current].op ) > 0.
424</tt>
425
426\param npar
427is the number of parameters corresponding to this operation sequence.
428
429\param par
430is a vector of length \a npar containing the parameters
431for this operation sequence; i.e.,
432given a parameter index \c i, the corresponding parameter value is
433<tt>par[i]</tt>.
434<!-- end optimize_prototype -->
435
436\param hash_table_var
437is a vector with size CPPAD_HASH_TABLE_SIZE
438that maps a hash code to the corresponding
439variable index in the old operation sequence.
440All the values in this table must be less than \a current.
441
442\param code
443The input value of code does not matter.
444The output value of code is the hash code corresponding to
445this operation in the new operation sequence.
446
447\return
448If the return value is zero,
449no match was found.
450If the return value is greater than zero,
451it is the index of a new variable that can be used to replace the
452old variable.
453
454
455\par Restrictions:
456The binary operator must be an addition, subtraction, multiplication, division
457or power operator.  NumArg( tape[current].op ) == 1.
458*/
459template <class Base>
460inline size_t optimize_binary_match(
461        const CppAD::vector<struct optimize_old_variable>& tape           ,
462        size_t                                             current        ,
463        size_t                                             npar           ,
464        const Base*                                        par            ,
465        const CppAD::vector<size_t>&                       hash_table_var ,
466        unsigned short&                                    code           )
467{       OpCode        op         = tape[current].op;
468        const addr_t* arg        = tape[current].arg;
469        addr_t        new_arg[2];
470        bool          parameter[2];
471
472        // initialize return value
473        size_t  match_var = 0;
474
475        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
476        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
477        switch(op)
478        {       // parameter op variable ----------------------------------
479                case AddpvOp:
480                case MulpvOp:
481                case DivpvOp:
482                case PowpvOp:
483                case SubpvOp:
484                // arg[0]
485                parameter[0] = true;
486                new_arg[0]   = arg[0];
487                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar );
488                // arg[1]
489                parameter[1] = false;
490                new_arg[1]   = tape[arg[1]].new_var;
491                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
492                break;
493
494                // variable op parameter -----------------------------------
495                case DivvpOp:
496                case PowvpOp:
497                case SubvpOp:
498                // arg[0]
499                parameter[0] = false;
500                new_arg[0]   = tape[arg[0]].new_var;
501                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
502                // arg[1]
503                parameter[1] = true;
504                new_arg[1]   = arg[1];
505                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar );
506                break;
507
508                // variable op variable -----------------------------------
509                case AddvvOp:
510                case MulvvOp:
511                case DivvvOp:
512                case PowvvOp:
513                case SubvvOp:
514                // arg[0]
515                parameter[0] = false;
516                new_arg[0]   = tape[arg[0]].new_var;
517                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
518                // arg[1]
519                parameter[1] = false;
520                new_arg[1]   = tape[arg[1]].new_var;
521                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
522                break;
523
524                // must be one of the cases above
525                default:
526                CPPAD_ASSERT_UNKNOWN(false);
527        }
528        code = hash_code(
529                op                  , 
530                new_arg             ,
531                npar                ,
532                par
533        );
534        size_t  i  = hash_table_var[code];
535        CPPAD_ASSERT_UNKNOWN( i < current );
536        if( op == tape[i].op )
537        {       bool match = true;
538                size_t j;
539                for(j = 0; j < 2; j++)
540                {       size_t k = tape[i].arg[j];
541                        if( parameter[j] )
542                        {       CPPAD_ASSERT_UNKNOWN( k < npar );
543                                match &= IdenticalEqualPar(
544                                        par[ arg[j] ], par[k]
545                                );
546                        }
547                        else
548                        {       CPPAD_ASSERT_UNKNOWN( k < i );
549                                match &= (new_arg[j] == tape[k].new_var);
550                        }
551                }
552                if( match )
553                        match_var = tape[i].new_var;
554        }
555        if( (match_var > 0) | ( (op != AddvvOp) & (op != MulvvOp ) ) )
556                return match_var;
557
558        // check for match with argument order switched ----------------------
559        CPPAD_ASSERT_UNKNOWN( op == AddvvOp || op == MulvvOp );
560        i          = new_arg[0];
561        new_arg[0] = new_arg[1];
562        new_arg[1] = i;
563        unsigned short code_switch = hash_code(
564                op                  , 
565                new_arg             ,
566                npar                ,
567                par
568        );
569        i  = hash_table_var[code_switch];
570        CPPAD_ASSERT_UNKNOWN( i < current );
571        if( op == tape[i].op )
572        {       bool match = true;
573                size_t j;
574                for(j = 0; j < 2; j++)
575                {       size_t k = tape[i].arg[j];
576                        CPPAD_ASSERT_UNKNOWN( k < i );
577                        match &= (new_arg[j] == tape[k].new_var);
578                }
579                if( match )
580                        match_var = tape[i].new_var;
581        }
582        return match_var;
583} 
584
585/*!
586Record an operation of the form (parameter op variable).
587
588<!-- replace optimize_prototype -->
589\param tape
590is a vector that maps a variable index, in the old operation sequence,
591to an <tt>optimize_old_variable</tt> information record.
592Note that the index for this vector must be greater than or equal zero and
593less than <tt>tape.size()</tt>.
594
595\li <tt>tape[i].op</tt>
596is the operator in the old operation sequence
597corresponding to the old variable index \c i.
598Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
599
600\li <tt>tape[i].arg</tt>
601for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
602is the j-th the argument, in the old operation sequence,
603corresponding to the old variable index \c i.
604Assertion: <tt>tape[i].arg[j] < i</tt>.
605
606\li <tt>tape[i].new_var</tt>
607Suppose
608<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
609and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
610It follows that <tt>tape[k].new_var</tt>
611has alread been set to the variable in the new operation sequence
612corresponding to the old variable index \c k.
613This means that the \c new_var value has been set
614for all the possible arguments that come before \a current.
615
616\param current
617is the index in the old operation sequence for
618the variable corresponding to the result for the current operator.
619Assertions:
620<tt>
621current < tape.size(),
622NumRes( tape[current].op ) > 0.
623</tt>
624
625\param npar
626is the number of parameters corresponding to this operation sequence.
627
628\param par
629is a vector of length \a npar containing the parameters
630for this operation sequence; i.e.,
631given a parameter index \c i, the corresponding parameter value is
632<tt>par[i]</tt>.
633<!-- end optimize_prototype -->
634
635\param rec
636is the object that will record the operations.
637
638\param op
639is the operator that we are recording which must be one of the following:
640AddpvOp, DivpvOp, MulpvOp, PowvpOp, SubpvOp.
641 
642\param arg
643is the vector of arguments for this operator.
644
645\return
646the result value is the index corresponding to the current
647operation in the new operation sequence.
648*/
649template <class Base>
650size_t optimize_record_pv(
651        const CppAD::vector<struct optimize_old_variable>& tape           ,
652        size_t                                             current        ,
653        size_t                                             npar           ,
654        const Base*                                        par            ,
655        recorder<Base>*                                    rec            ,
656        OpCode                                             op             ,
657        const addr_t*                                      arg            )
658{
659# ifndef NDEBUG
660        switch(op)
661        {       case AddpvOp:
662                case DivpvOp:
663                case MulpvOp:
664                case PowpvOp:
665                case SubpvOp:
666                break;
667
668                default:
669                CPPAD_ASSERT_UNKNOWN(false);
670        }
671# endif
672        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar    );
673        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
674        addr_t new_arg[2];
675        new_arg[0]   = rec->PutPar( par[arg[0]] );
676        new_arg[1]   = tape[ arg[1] ].new_var;
677        rec->PutArg( new_arg[0], new_arg[1] );
678        size_t i     = rec->PutOp(op);
679        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
680        return i;
681}
682
683
684/*!
685Record an operation of the form (variable op parameter).
686
687<!-- replace optimize_prototype -->
688\param tape
689is a vector that maps a variable index, in the old operation sequence,
690to an <tt>optimize_old_variable</tt> information record.
691Note that the index for this vector must be greater than or equal zero and
692less than <tt>tape.size()</tt>.
693
694\li <tt>tape[i].op</tt>
695is the operator in the old operation sequence
696corresponding to the old variable index \c i.
697Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
698
699\li <tt>tape[i].arg</tt>
700for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
701is the j-th the argument, in the old operation sequence,
702corresponding to the old variable index \c i.
703Assertion: <tt>tape[i].arg[j] < i</tt>.
704
705\li <tt>tape[i].new_var</tt>
706Suppose
707<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
708and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
709It follows that <tt>tape[k].new_var</tt>
710has alread been set to the variable in the new operation sequence
711corresponding to the old variable index \c k.
712This means that the \c new_var value has been set
713for all the possible arguments that come before \a current.
714
715\param current
716is the index in the old operation sequence for
717the variable corresponding to the result for the current operator.
718Assertions:
719<tt>
720current < tape.size(),
721NumRes( tape[current].op ) > 0.
722</tt>
723
724\param npar
725is the number of parameters corresponding to this operation sequence.
726
727\param par
728is a vector of length \a npar containing the parameters
729for this operation sequence; i.e.,
730given a parameter index \c i, the corresponding parameter value is
731<tt>par[i]</tt>.
732<!-- end optimize_prototype -->
733
734\param rec
735is the object that will record the operations.
736
737\param op
738is the operator that we are recording which must be one of the following:
739DivvpOp, PowvpOp, SubvpOp.
740 
741\param arg
742is the vector of arguments for this operator.
743
744\return
745the result value is the index corresponding to the current
746operation in the new operation sequence.
747*/
748template <class Base>
749size_t optimize_record_vp(
750        const CppAD::vector<struct optimize_old_variable>& tape           ,
751        size_t                                             current        ,
752        size_t                                             npar           ,
753        const Base*                                        par            ,
754        recorder<Base>*                                    rec            ,
755        OpCode                                             op             ,
756        const addr_t*                                      arg            )
757{
758# ifndef NDEBUG
759        switch(op)
760        {       case DivvpOp:
761                case PowvpOp:
762                case SubvpOp:
763                break;
764
765                default:
766                CPPAD_ASSERT_UNKNOWN(false);
767        }
768# endif
769        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
770        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar    );
771        addr_t new_arg[2];
772        new_arg[0]   = tape[ arg[0] ].new_var;
773        new_arg[1]   = rec->PutPar( par[arg[1]] );
774        rec->PutArg( new_arg[0], new_arg[1] );
775        size_t i     = rec->PutOp(op);
776        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
777        return i;
778}
779
780/*!
781Record an operation of the form (variable op variable).
782
783<!-- replace optimize_prototype -->
784\param tape
785is a vector that maps a variable index, in the old operation sequence,
786to an <tt>optimize_old_variable</tt> information record.
787Note that the index for this vector must be greater than or equal zero and
788less than <tt>tape.size()</tt>.
789
790\li <tt>tape[i].op</tt>
791is the operator in the old operation sequence
792corresponding to the old variable index \c i.
793Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
794
795\li <tt>tape[i].arg</tt>
796for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
797is the j-th the argument, in the old operation sequence,
798corresponding to the old variable index \c i.
799Assertion: <tt>tape[i].arg[j] < i</tt>.
800
801\li <tt>tape[i].new_var</tt>
802Suppose
803<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
804and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
805It follows that <tt>tape[k].new_var</tt>
806has alread been set to the variable in the new operation sequence
807corresponding to the old variable index \c k.
808This means that the \c new_var value has been set
809for all the possible arguments that come before \a current.
810
811\param current
812is the index in the old operation sequence for
813the variable corresponding to the result for the current operator.
814Assertions:
815<tt>
816current < tape.size(),
817NumRes( tape[current].op ) > 0.
818</tt>
819
820\param npar
821is the number of parameters corresponding to this operation sequence.
822
823\param par
824is a vector of length \a npar containing the parameters
825for this operation sequence; i.e.,
826given a parameter index \c i, the corresponding parameter value is
827<tt>par[i]</tt>.
828<!-- end optimize_prototype -->
829
830\param rec
831is the object that will record the operations.
832
833\param op
834is the operator that we are recording which must be one of the following:
835AddvvOp, DivvvOp, MulvvOp, PowvpOp, SubvvOp.
836 
837\param arg
838is the vector of arguments for this operator.
839
840\return
841the result value is the index corresponding to the current
842operation in the new operation sequence.
843*/
844template <class Base>
845size_t optimize_record_vv(
846        const CppAD::vector<struct optimize_old_variable>& tape           ,
847        size_t                                             current        ,
848        size_t                                             npar           ,
849        const Base*                                        par            ,
850        recorder<Base>*                                    rec            ,
851        OpCode                                             op             ,
852        const addr_t*                                      arg            )
853{
854# ifndef NDEBUG
855        switch(op)
856        {       case AddvvOp:
857                case DivvvOp:
858                case MulvvOp:
859                case PowvvOp:
860                case SubvvOp:
861                break;
862
863                default:
864                CPPAD_ASSERT_UNKNOWN(false);
865        }
866# endif
867        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
868        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
869        addr_t new_arg[2];
870        new_arg[0]   = tape[ arg[0] ].new_var;
871        new_arg[1]   = tape[ arg[1] ].new_var;
872        rec->PutArg( new_arg[0], new_arg[1] );
873        size_t i     = rec->PutOp(op);
874        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
875        return i;
876}
877
878// ==========================================================================
879
880/*!
881Recording a cummulative cummulative summation starting at its highest parrent.
882
883<!-- replace optimize_prototype -->
884\param tape
885is a vector that maps a variable index, in the old operation sequence,
886to an <tt>optimize_old_variable</tt> information record.
887Note that the index for this vector must be greater than or equal zero and
888less than <tt>tape.size()</tt>.
889
890\li <tt>tape[i].op</tt>
891is the operator in the old operation sequence
892corresponding to the old variable index \c i.
893Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
894
895\li <tt>tape[i].arg</tt>
896for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
897is the j-th the argument, in the old operation sequence,
898corresponding to the old variable index \c i.
899Assertion: <tt>tape[i].arg[j] < i</tt>.
900
901\li <tt>tape[i].new_var</tt>
902Suppose
903<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
904and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
905It follows that <tt>tape[k].new_var</tt>
906has alread been set to the variable in the new operation sequence
907corresponding to the old variable index \c k.
908This means that the \c new_var value has been set
909for all the possible arguments that come before \a current.
910
911\param current
912is the index in the old operation sequence for
913the variable corresponding to the result for the current operator.
914Assertions:
915<tt>
916current < tape.size(),
917NumRes( tape[current].op ) > 0.
918</tt>
919
920\param npar
921is the number of parameters corresponding to this operation sequence.
922
923\param par
924is a vector of length \a npar containing the parameters
925for this operation sequence; i.e.,
926given a parameter index \c i, the corresponding parameter value is
927<tt>par[i]</tt>.
928<!-- end optimize_prototype -->
929
930\param rec
931is the object that will record the operations.
932
933\param work
934Is used for computation. On input and output,
935<tt>work.op_stack.empty()</tt>,
936<tt>work.add_stack.empty()</tt>, and
937<tt>work.sub_stack.empty()</tt>,
938are all true true.
939These stacks are passed in so that elements can be allocated once
940and then the elements can be reused with calls to \c optimize_record_csum.
941
942\par Exception
943<tt>tape[i].new_var</tt>
944is not yet defined for any node \c i that is \c csum_connected
945to the \a current node
946(or that is \c sum_connected to a node that is \c csum_connected).
947For example; suppose that index \c j corresponds to a variable
948in the current operator,
949<tt>i = tape[current].arg[j]</tt>,
950and
951<tt>tape[arg[j]].connect_type == csum_connected</tt>.
952It then follows that
953<tt>tape[i].new_var == tape.size()</tt>.
954
955\par Restriction:
956\li <tt>tape[current].op</tt>
957must be one of <tt>AddpvOp, AddvvOp, SubpvOp, SubvpOp, SubvvOp</tt>.
958
959\li <tt>tape[current].connect_type</tt> must be \c yes_connected.
960
961\li <tt>tape[j].connect_type == csum_connected</tt> for some index
962j that is a variable operand for the current operation.
963*/
964
965
966template <class Base>
967size_t optimize_record_csum(
968        const CppAD::vector<struct optimize_old_variable>& tape           ,
969        size_t                                             current        ,
970        size_t                                             npar           ,
971        const Base*                                        par            ,
972        recorder<Base>*                                    rec            ,
973        optimize_csum_stacks&                              work           )
974{
975       
976        CPPAD_ASSERT_UNKNOWN( work.op_stack.empty() );
977        CPPAD_ASSERT_UNKNOWN( work.add_stack.empty() );
978        CPPAD_ASSERT_UNKNOWN( work.sub_stack.empty() );
979        CPPAD_ASSERT_UNKNOWN( tape[current].connect_type == yes_connected );
980
981        size_t                        i;
982        OpCode                        op;
983        const addr_t*                 arg;
984        bool                          add;
985        struct optimize_csum_variable var;
986
987        var.op  = tape[current].op;
988        var.arg = tape[current].arg;
989        var.add = true; 
990        work.op_stack.push( var );
991        Base sum_par(0);
992
993# ifndef NDEBUG
994        bool ok = false;
995        if( var.op == SubvpOp )
996                ok = tape[ tape[current].arg[0] ].connect_type == csum_connected;
997        if( var.op == AddpvOp || var.op == SubpvOp )
998                ok = tape[ tape[current].arg[1] ].connect_type == csum_connected;
999        if( var.op == AddvvOp || var.op == SubvvOp )
1000        {       ok  = tape[ tape[current].arg[0] ].connect_type == csum_connected;
1001                ok |= tape[ tape[current].arg[1] ].connect_type == csum_connected;
1002        }
1003        CPPAD_ASSERT_UNKNOWN( ok );
1004# endif
1005        while( ! work.op_stack.empty() )
1006        {       var     = work.op_stack.top();
1007                work.op_stack.pop();
1008                op      = var.op;
1009                arg     = var.arg;
1010                add     = var.add;
1011                // process first argument to this operator
1012                switch(op)
1013                {       case AddpvOp:
1014                        case SubpvOp:
1015                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar );
1016                        if( add )
1017                                sum_par += par[arg[0]];
1018                        else    sum_par -= par[arg[0]];
1019                        break;
1020
1021                        case AddvvOp:
1022                        case SubvpOp:
1023                        case SubvvOp:
1024                        if( tape[arg[0]].connect_type == csum_connected )
1025                        {       CPPAD_ASSERT_UNKNOWN(
1026                                        size_t(tape[arg[0]].new_var) == tape.size()
1027                                );
1028                                var.op  = tape[arg[0]].op;
1029                                var.arg = tape[arg[0]].arg;
1030                                var.add = add; 
1031                                work.op_stack.push( var );
1032                        }
1033                        else if( add )
1034                                work.add_stack.push(arg[0]);
1035                        else    work.sub_stack.push(arg[0]);
1036                        break;
1037
1038                        default:
1039                        CPPAD_ASSERT_UNKNOWN(false);
1040                }
1041                // process second argument to this operator
1042                switch(op)
1043                {
1044                        case SubvpOp:
1045                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar );
1046                        if( add )
1047                                sum_par -= par[arg[1]];
1048                        else    sum_par += par[arg[1]];
1049                        break;
1050
1051                        case SubvvOp:
1052                        case SubpvOp:
1053                        add = ! add;
1054
1055                        case AddvvOp:
1056                        case AddpvOp:
1057                        if( tape[arg[1]].connect_type == csum_connected )
1058                        {       CPPAD_ASSERT_UNKNOWN(
1059                                        size_t(tape[arg[1]].new_var) == tape.size()
1060                                );
1061                                var.op   = tape[arg[1]].op;
1062                                var.arg  = tape[arg[1]].arg;
1063                                var.add  = add;
1064                                work.op_stack.push( var );
1065                        }
1066                        else if( add )
1067                                work.add_stack.push(arg[1]);
1068                        else    work.sub_stack.push(arg[1]);
1069                        break;
1070
1071                        default:
1072                        CPPAD_ASSERT_UNKNOWN(false);
1073                }
1074        }
1075        // number of variables in this cummulative sum operator
1076        size_t n_add = work.add_stack.size();
1077        size_t n_sub = work.sub_stack.size();
1078        size_t old_arg, new_arg;
1079        rec->PutArg(n_add);                // arg[0]
1080        rec->PutArg(n_sub);                // arg[1]
1081        new_arg = rec->PutPar( sum_par );
1082        rec->PutArg(new_arg);              // arg[2]
1083        for(i = 0; i < n_add; i++)
1084        {       CPPAD_ASSERT_UNKNOWN( ! work.add_stack.empty() );
1085                old_arg = work.add_stack.top();
1086                new_arg = tape[old_arg].new_var;
1087                CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
1088                rec->PutArg(new_arg);      // arg[3+i]
1089                work.add_stack.pop();
1090        }
1091        for(i = 0; i < n_sub; i++)
1092        {       CPPAD_ASSERT_UNKNOWN( ! work.sub_stack.empty() );
1093                old_arg = work.sub_stack.top();
1094                new_arg = tape[old_arg].new_var;
1095                CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
1096                rec->PutArg(new_arg);      // arg[3 + arg[0] + i]
1097                work.sub_stack.pop();
1098        }
1099        rec->PutArg(n_add + n_sub);        // arg[3 + arg[0] + arg[1]]
1100        i = rec->PutOp(CSumOp);
1101        CPPAD_ASSERT_UNKNOWN(new_arg < tape.size());
1102
1103        return i;
1104}
1105// ==========================================================================
1106/*!
1107Convert a player object to an optimized recorder object
1108
1109\tparam Base
1110base type for the operator; i.e., this operation was recorded
1111using AD< \a Base > and computations by this routine are done using type
1112\a Base.
1113
1114\param n
1115is the number of independent variables on the tape.
1116
1117\param dep_taddr
1118On input this vector contains the indices for each of the dependent
1119variable values in the operation sequence corresponding to \a play.
1120Upon return it contains the indices for the same variables but in
1121the operation sequence corresponding to \a rec.
1122
1123\param play
1124This is the operation sequence that we are optimizing.
1125It is essentially const, except for play back state which
1126changes while it plays back the operation seqeunce.
1127
1128\param rec
1129The input contents of this recording does not matter.
1130Upon return, it contains an optimized verison of the
1131operation sequence corresponding to \a play.
1132*/
1133
1134template <class Base>
1135void optimize(
1136        size_t                       n         ,
1137        CppAD::vector<size_t>&       dep_taddr ,
1138        player<Base>*                play      ,
1139        recorder<Base>*              rec       ) 
1140{
1141        // temporary indices
1142        size_t i, j, k;
1143
1144        // temporary variables
1145        OpCode        op;   // current operator
1146        const addr_t* arg;  // operator arguments
1147        size_t        i_var;  // index of first result for current operator
1148
1149        // range and domain dimensions for F
1150        size_t m = dep_taddr.size();
1151
1152        // number of variables in the player
1153        const size_t num_var = play->num_rec_var(); 
1154
1155# ifndef NDEBUG
1156        // number of parameters in the player
1157        const size_t num_par = play->num_rec_par();
1158# endif
1159
1160        // number of  VecAD indices
1161        size_t num_vecad_ind   = play->num_rec_vecad_ind();
1162
1163        // number of VecAD vectors
1164        size_t num_vecad_vec   = play->num_rec_vecad_vec();
1165
1166        // -------------------------------------------------------------
1167        // data structure that maps variable index in original operation
1168        // sequence to corresponding operator information
1169        CppAD::vector<struct optimize_old_variable> tape(num_var);
1170        // -------------------------------------------------------------
1171        // Determine how each variable is connected to the dependent variables
1172
1173        // initialize all variables has having no connections
1174        for(i = 0; i < num_var; i++)
1175                tape[i].connect_type = not_connected;
1176
1177        for(j = 0; j < m; j++)
1178        {       // mark dependent variables as having one or more connections
1179                tape[ dep_taddr[j] ].connect_type = yes_connected;
1180        }
1181
1182        // vecad_connect contains a value for each VecAD object.
1183        // vecad maps a VecAD index (which corresponds to the beginning of the
1184        // VecAD object) to the vecad_connect falg for the VecAD object.
1185        CppAD::vector<optimize_connection_type>   vecad_connect(num_vecad_vec);
1186        CppAD::vector<size_t> vecad(num_vecad_ind);
1187        j = 0;
1188        for(i = 0; i < num_vecad_vec; i++)
1189        {       vecad_connect[i] = not_connected;
1190                // length of this VecAD
1191                size_t length = play->GetVecInd(j);
1192                // set to proper index for this VecAD
1193                vecad[j] = i; 
1194                for(k = 1; k <= length; k++)
1195                        vecad[j+k] = num_vecad_vec; // invalid index
1196                // start of next VecAD
1197                j       += length + 1;
1198        }
1199        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
1200
1201        // work space used by UserOp.
1202        typedef std::set<size_t> size_set;
1203        size_t user_q     = 0;       // maximum set element plus one
1204        vector< size_set > user_r;   // sparsity pattern for the argument x
1205        vector< size_set > user_s;   // sparisty pattern for the result y
1206        size_t user_index = 0;       // indentifier for this user_atomic operation
1207        size_t user_id    = 0;       // user identifier for this call to operator
1208        size_t user_i     = 0;       // index in result vector
1209        size_t user_j     = 0;       // index in argument vector
1210        size_t user_m     = 0;       // size of result vector
1211        size_t user_n     = 0;       // size of arugment vector
1212        // next expected operator in a UserOp sequence
1213        enum { user_start, user_arg, user_ret, user_end } user_state;
1214
1215        // During reverse mode, push true if user operation is connected,
1216        // push false otherwise. During forward mode, use to determine if
1217        // we are keeping this operation and then pop.
1218        std::stack<bool> user_keep;
1219
1220        // Initialize a reverse mode sweep through the operation sequence
1221        size_t i_op;
1222        play->start_reverse(op, arg, i_op, i_var);
1223        CPPAD_ASSERT_UNKNOWN( op == EndOp );
1224        size_t mask;
1225        user_state = user_end;
1226        while(op != BeginOp)
1227        {       // next op
1228                play->next_reverse(op, arg, i_op, i_var);
1229                // This if is not necessary becasue last assignment
1230                // with this value of i_var will have NumRes(op) > 0
1231                if( NumRes(op) > 0 )
1232                {       tape[i_var].op = op;
1233                        tape[i_var].arg = arg;
1234                }
1235# ifndef NDEBUG
1236                if( i_op <= n )
1237                {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
1238                }
1239                else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
1240# endif
1241                switch( op )
1242                {
1243                        // Unary operator where operand is arg[0]
1244                        case AbsOp:
1245                        case AcosOp:
1246                        case AsinOp:
1247                        case AtanOp:
1248                        case CosOp:
1249                        case CoshOp:
1250                        case DivvpOp:
1251                        case ExpOp:
1252                        case LogOp:
1253                        case PowvpOp:
1254                        case SinOp:
1255                        case SinhOp:
1256                        case SqrtOp:
1257                        case TanOp:
1258                        case TanhOp:
1259                        if( tape[i_var].connect_type != not_connected )
1260                                tape[arg[0]].connect_type = yes_connected;
1261                        break; // --------------------------------------------
1262
1263                        // Unary operator where operand is arg[1]
1264                        case DisOp:
1265                        case DivpvOp:
1266                        case MulpvOp:
1267                        case PowpvOp:
1268                        if( tape[i_var].connect_type != not_connected )
1269                                tape[arg[1]].connect_type = yes_connected;
1270                        break; // --------------------------------------------
1271               
1272                        // Special case for SubvpOp
1273                        case SubvpOp:
1274                        if( tape[i_var].connect_type != not_connected )
1275                        {
1276                                // check for case where arg[0] is sum_connected
1277                                if( tape[arg[0]].connect_type == not_connected )
1278                                        tape[arg[0]].connect_type = sum_connected;
1279                                else    tape[arg[0]].connect_type = yes_connected;
1280
1281                                // check for case where i_var is sum_connected to and
1282                                // it is the result of a summation.
1283                                if( tape[i_var].connect_type == sum_connected )
1284                                        tape[i_var].connect_type = csum_connected;
1285                        }
1286                        break; // --------------------------------------------
1287               
1288                        // Special case for AddpvOp and SubpvOp
1289                        case AddpvOp:
1290                        case SubpvOp:
1291                        if( tape[i_var].connect_type != not_connected )
1292                        {
1293                                // check for case where arg[1] is sum_connected
1294                                if( tape[arg[1]].connect_type == not_connected )
1295                                        tape[arg[1]].connect_type = sum_connected;
1296                                else    tape[arg[1]].connect_type = yes_connected;
1297
1298                                // check for case where i_var is sum_connected to and
1299                                // it is the result of a summation.
1300                                if( tape[i_var].connect_type == sum_connected )
1301                                        tape[i_var].connect_type = csum_connected;
1302                        }
1303                        break; // --------------------------------------------
1304
1305               
1306                        // Special case for AddvvOp and SubvvOp
1307                        case AddvvOp:
1308                        case SubvvOp:
1309                        if( tape[i_var].connect_type != not_connected )
1310                        {
1311                                // check for case where arg[0] is sum_connected
1312                                if( tape[arg[0]].connect_type == not_connected )
1313                                        tape[arg[0]].connect_type = sum_connected;
1314                                else    tape[arg[0]].connect_type = yes_connected;
1315
1316                                // check for case where arg[1] is sum_connected
1317                                if( tape[arg[1]].connect_type == not_connected )
1318                                        tape[arg[1]].connect_type = sum_connected;
1319                                else    tape[arg[1]].connect_type = yes_connected;
1320
1321                                // check for case where i_var is sum_connected to and
1322                                // it is the result of a summation.
1323                                if( tape[i_var].connect_type == sum_connected )
1324                                        tape[i_var].connect_type = csum_connected;
1325                        }
1326                        break; // --------------------------------------------
1327
1328                        // Other binary operators
1329                        // where operands are arg[0], arg[1]
1330                        case DivvvOp:
1331                        case MulvvOp:
1332                        case PowvvOp:
1333                        if( tape[i_var].connect_type != not_connected )
1334                        {
1335                                tape[arg[0]].connect_type = yes_connected;
1336                                tape[arg[1]].connect_type = yes_connected;
1337                        }
1338                        break; // --------------------------------------------
1339
1340                        // Conditional expression operators
1341                        case CExpOp:
1342                        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
1343                        if( tape[i_var].connect_type != not_connected )
1344                        {
1345                                mask = 1;
1346                                for(i = 2; i < 6; i++)
1347                                {       if( arg[1] & mask )
1348                                        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[i]) < i_var );
1349                                                tape[arg[i]].connect_type = yes_connected;
1350                                        }
1351                                        mask = mask << 1;
1352                                }
1353                        }
1354                        break;  // --------------------------------------------
1355
1356                        // Operations where there is noting to do
1357                        case BeginOp:
1358                        case ComOp:
1359                        case EndOp:
1360                        case InvOp:
1361                        case ParOp:
1362                        case PriOp:
1363                        break;  // --------------------------------------------
1364
1365                        // Load using a parameter index
1366                        case LdpOp:
1367                        if( tape[i_var].connect_type != not_connected )
1368                        {
1369                                i                = vecad[ arg[0] - 1 ];
1370                                vecad_connect[i] = yes_connected;
1371                        }
1372                        break; // --------------------------------------------
1373
1374                        // Load using a variable index
1375                        case LdvOp:
1376                        if( tape[i_var].connect_type != not_connected )
1377                        {
1378                                i                    = vecad[ arg[0] - 1 ];
1379                                vecad_connect[i]     = yes_connected;
1380                                tape[arg[1]].connect_type = yes_connected;
1381                        }
1382                        break; // --------------------------------------------
1383
1384                        // Store a variable using a parameter index
1385                        case StpvOp:
1386                        i = vecad[ arg[0] - 1 ];
1387                        if( vecad_connect[i] != not_connected )
1388                                tape[arg[2]].connect_type = yes_connected;
1389                        break; // --------------------------------------------
1390
1391                        // Store a variable using a variable index
1392                        case StvvOp:
1393                        i = vecad[ arg[0] - 1 ];
1394                        if( vecad_connect[i] )
1395                        {       tape[arg[1]].connect_type = yes_connected;
1396                                tape[arg[2]].connect_type = yes_connected;
1397                        }
1398                        break; 
1399                        // ============================================================
1400                        case UserOp:
1401                        // start or end atomic operation sequence
1402                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
1403                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
1404                        if( user_state == user_end )
1405                        {       user_index = arg[0];
1406                                user_id    = arg[1];
1407                                user_n     = arg[2];
1408                                user_m     = arg[3];
1409                                user_q     = 1;
1410                                if(user_s.size() != user_n )
1411                                        user_s.resize(user_n);
1412                                if(user_r.size() != user_m )
1413                                        user_r.resize(user_m);
1414                                user_j     = user_n;
1415                                user_i     = user_m;
1416                                user_state = user_ret;
1417                                user_keep.push(false);
1418                        }
1419                        else
1420                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
1421                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
1422                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
1423                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
1424                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
1425                                user_state = user_end;
1426               }
1427                        break;
1428
1429                        case UsrapOp:
1430                        // parameter argument in an atomic operation sequence
1431                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
1432                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
1433                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1434                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
1435                        --user_j;
1436                        if( user_j == 0 )
1437                                user_state = user_start;
1438                        break;
1439
1440                        case UsravOp:
1441                        // variable argument in an atomic operation sequence
1442                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
1443                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
1444                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1445                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
1446                        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
1447                        --user_j;
1448                        if( ! user_s[user_j].empty() )
1449                        {       tape[arg[0]].connect_type = yes_connected;
1450                                user_keep.top() = true;
1451                        }
1452                        if( user_j == 0 )
1453                                user_state = user_start;
1454                        break;
1455
1456                        case UsrrpOp:
1457                        // parameter result in an atomic operation sequence
1458                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
1459                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
1460                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1461                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
1462                        --user_i;
1463                        user_r[user_i].clear();
1464                        if( user_i == 0 )
1465                        {       // call users function for this operation
1466                                atomic_base<Base>* atom = 
1467                                        atomic_base<Base>::class_object(user_index);
1468                                atom->set_id(user_id);
1469                                atom->rev_sparse_jac(
1470                                        user_q, user_r, user_s
1471                                );
1472                                user_state = user_arg;
1473                        }
1474                        break;
1475
1476                        case UsrrvOp:
1477                        // variable result in an atomic operation sequence
1478                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
1479                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
1480                        --user_i;
1481                        user_r[user_i].clear();
1482                        if( tape[i_var].connect_type != not_connected )
1483                                user_r[user_i].insert(0);
1484                        if( user_i == 0 )
1485                        {       // call users function for this operation
1486                                atomic_base<Base>* atom = 
1487                                        atomic_base<Base>::class_object(user_index);
1488                                atom->set_id(user_id);
1489                                atom->rev_sparse_jac(
1490                                        user_q, user_r, user_s
1491                                );
1492                                user_state = user_arg;
1493                        }
1494                        break;
1495                        // ============================================================
1496
1497                        // all cases should be handled above
1498                        default:
1499                        CPPAD_ASSERT_UNKNOWN(0);
1500                }
1501        }
1502        // values corresponding to BeginOp
1503        CPPAD_ASSERT_UNKNOWN( i_op == 0 && i_var == 0 && op == BeginOp );
1504        tape[i_var].op = op;
1505        // -------------------------------------------------------------
1506
1507        // Erase all information in the recording
1508        rec->free();
1509
1510        // Initilaize table mapping hash code to variable index in tape
1511        // as pointing to the BeginOp at the beginning of the tape
1512        CppAD::vector<size_t>  hash_table_var(CPPAD_HASH_TABLE_SIZE);
1513        for(i = 0; i < CPPAD_HASH_TABLE_SIZE; i++)
1514                hash_table_var[i] = 0;
1515        CPPAD_ASSERT_UNKNOWN( tape[0].op == BeginOp );
1516
1517        // initialize mapping from old variable index to new variable index
1518        for(i = 0; i < num_var; i++)
1519                tape[i].new_var = num_var; // invalid index
1520       
1521
1522        // initialize mapping from old VecAD index to new VecAD index
1523        CppAD::vector<size_t> new_vecad_ind(num_vecad_ind);
1524        for(i = 0; i < num_vecad_ind; i++)
1525                new_vecad_ind[i] = num_vecad_ind; // invalid index
1526
1527        j = 0;     // index into the old set of indices
1528        for(i = 0; i < num_vecad_vec; i++)
1529        {       // length of this VecAD
1530                size_t length = play->GetVecInd(j);
1531                if( vecad_connect[i] != not_connected )
1532                {       // Put this VecAD vector in new recording
1533                        CPPAD_ASSERT_UNKNOWN(length < num_vecad_ind);
1534                        new_vecad_ind[j] = rec->PutVecInd(length);
1535                        for(k = 1; k <= length; k++) new_vecad_ind[j+k] =
1536                                rec->PutVecInd(
1537                                        rec->PutPar(
1538                                                play->GetPar( 
1539                                                        play->GetVecInd(j+k)
1540                        ) ) );
1541                }
1542                // start of next VecAD
1543                j       += length + 1;
1544        }
1545        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
1546
1547        // start playing the operations in the forward direction
1548        play->start_forward(op, arg, i_op, i_var);
1549
1550        // playing forward skips BeginOp at the beginning, but not EndOp at
1551        // the end.  Put BeginOp at beginning of recording
1552        CPPAD_ASSERT_UNKNOWN( op == BeginOp );
1553        CPPAD_ASSERT_NARG_NRES(BeginOp, 0, 1);
1554        tape[i_var].new_var = rec->PutOp(BeginOp);
1555
1556        // temporary buffer for new argument values
1557        addr_t new_arg[6];
1558
1559        // temporary work space used by optimize_record_csum
1560        // (decalared here to avoid realloaction of memory)
1561        optimize_csum_stacks csum_work;
1562
1563        user_state = user_start;
1564        while(op != EndOp)
1565        {       // next op
1566                play->next_forward(op, arg, i_op, i_var);
1567                CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
1568                CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
1569
1570                // determine if we should keep this operation in the new
1571                // operation sequence
1572                bool keep;
1573                switch( op )
1574                {       case ComOp:
1575                        case PriOp:
1576                        keep = false;
1577                        break;
1578
1579                        case InvOp:
1580                        case EndOp:
1581                        keep = true;
1582                        break;
1583
1584                        case StppOp:
1585                        case StvpOp:
1586                        case StpvOp:
1587                        case StvvOp:
1588                        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
1589                        i = vecad[ arg[0] - 1 ];
1590                        keep = vecad_connect[i] != not_connected;
1591                        break;
1592
1593                        case AddpvOp:
1594                        case AddvvOp:
1595                        case SubpvOp:
1596                        case SubvpOp:
1597                        case SubvvOp:
1598                        keep  = tape[i_var].connect_type != not_connected;
1599                        keep &= tape[i_var].connect_type != csum_connected;
1600                        break; 
1601
1602                        case UserOp:
1603                        case UsrapOp:
1604                        case UsravOp:
1605                        case UsrrpOp:
1606                        case UsrrvOp:
1607                        keep = user_keep.top();
1608                        break;
1609
1610                        default:
1611                        keep = tape[i_var].connect_type != not_connected;
1612                        break;
1613                }
1614
1615                size_t         match_var    = 0;
1616                unsigned short code         = 0;
1617                bool           replace_hash = false;
1618                if( keep ) switch( op )
1619                {
1620                        // Unary operator where operand is arg[0]
1621                        case AbsOp:
1622                        case AcosOp:
1623                        case AsinOp:
1624                        case AtanOp:
1625                        case CosOp:
1626                        case CoshOp:
1627                        case ExpOp:
1628                        case LogOp:
1629                        case SinOp:
1630                        case SinhOp:
1631                        case SqrtOp:
1632                        case TanOp:
1633                        case TanhOp:
1634                        match_var = optimize_unary_match(
1635                                tape                ,  // inputs
1636                                i_var               ,
1637                                play->num_rec_par() ,
1638                                play->GetPar()      ,
1639                                hash_table_var      ,
1640                                code                  // outputs
1641                        );
1642                        if( match_var > 0 )
1643                                tape[i_var].new_var = match_var;
1644                        else
1645                        {
1646                                replace_hash = true;
1647                                new_arg[0]   = tape[ arg[0] ].new_var;
1648                                rec->PutArg( new_arg[0] );
1649                                i                   = rec->PutOp(op);
1650                                tape[i_var].new_var = i;
1651                                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
1652                        }
1653                        break;
1654                        // ---------------------------------------------------
1655                        // Binary operators where
1656                        // left is a variable and right is a parameter
1657                        case SubvpOp:
1658                        if( tape[arg[0]].connect_type == csum_connected )
1659                        {
1660                                // convert to a sequence of summation operators
1661                                tape[i_var].new_var = optimize_record_csum(
1662                                        tape                , // inputs
1663                                        i_var               ,
1664                                        play->num_rec_par() ,
1665                                        play->GetPar()      ,
1666                                        rec                 ,
1667                                        csum_work
1668                                );
1669                                // abort rest of this case
1670                                break;
1671                        }
1672                        case DivvpOp:
1673                        case PowvpOp:
1674                        match_var = optimize_binary_match(
1675                                tape                ,  // inputs
1676                                i_var               ,
1677                                play->num_rec_par() ,
1678                                play->GetPar()      ,
1679                                hash_table_var      ,
1680                                code                  // outputs
1681                        );
1682                        if( match_var > 0 )
1683                                tape[i_var].new_var = match_var;
1684                        else
1685                        {       tape[i_var].new_var = optimize_record_vp(
1686                                        tape                , // inputs
1687                                        i_var               ,
1688                                        play->num_rec_par() ,
1689                                        play->GetPar()      ,
1690                                        rec                 ,
1691                                        op                  ,
1692                                        arg
1693                                );
1694                                replace_hash = true;
1695                        }
1696                        break;
1697                        // ---------------------------------------------------
1698                        // Binary operators where
1699                        // left is a parameter and right is a variable
1700                        case SubpvOp:
1701                        case AddpvOp:
1702                        if( tape[arg[1]].connect_type == csum_connected )
1703                        {
1704                                // convert to a sequence of summation operators
1705                                tape[i_var].new_var = optimize_record_csum(
1706                                        tape                , // inputs
1707                                        i_var               ,
1708                                        play->num_rec_par() ,
1709                                        play->GetPar()      ,
1710                                        rec                 ,
1711                                        csum_work
1712                                );
1713                                // abort rest of this case
1714                                break;
1715                        }
1716                        case DivpvOp:
1717                        case MulpvOp:
1718                        case PowpvOp:
1719                        match_var = optimize_binary_match(
1720                                tape                ,  // inputs
1721                                i_var               ,
1722                                play->num_rec_par() ,
1723                                play->GetPar()      ,
1724                                hash_table_var      ,
1725                                code                  // outputs
1726                        );
1727                        if( match_var > 0 )
1728                                tape[i_var].new_var = match_var;
1729                        else
1730                        {       tape[i_var].new_var = optimize_record_pv(
1731                                        tape                , // inputs
1732                                        i_var               ,
1733                                        play->num_rec_par() ,
1734                                        play->GetPar()      ,
1735                                        rec                 ,
1736                                        op                  ,
1737                                        arg
1738                                );
1739                                replace_hash = true;
1740                        }
1741                        break;
1742                        // ---------------------------------------------------
1743                        // Binary operator where
1744                        // both operators are variables
1745                        case AddvvOp:
1746                        case SubvvOp:
1747                        if( (tape[arg[0]].connect_type == csum_connected) |
1748                            (tape[arg[1]].connect_type == csum_connected)
1749                        )
1750                        {
1751                                // convert to a sequence of summation operators
1752                                tape[i_var].new_var = optimize_record_csum(
1753                                        tape                , // inputs
1754                                        i_var               ,
1755                                        play->num_rec_par() ,
1756                                        play->GetPar()      ,
1757                                        rec                 ,
1758                                        csum_work
1759                                );
1760                                // abort rest of this case
1761                                break;
1762                        }
1763                        case DivvvOp:
1764                        case MulvvOp:
1765                        case PowvvOp:
1766                        match_var = optimize_binary_match(
1767                                tape                ,  // inputs
1768                                i_var               ,
1769                                play->num_rec_par() ,
1770                                play->GetPar()      ,
1771                                hash_table_var      ,
1772                                code                  // outputs
1773                        );
1774                        if( match_var > 0 )
1775                                tape[i_var].new_var = match_var;
1776                        else
1777                        {       tape[i_var].new_var = optimize_record_vv(
1778                                        tape                , // inputs
1779                                        i_var               ,
1780                                        play->num_rec_par() ,
1781                                        play->GetPar()      ,
1782                                        rec                 ,
1783                                        op                  ,
1784                                        arg
1785                                );
1786                                replace_hash = true;
1787                        }
1788                        break;
1789                        // ---------------------------------------------------
1790                        // Conditional expression operators
1791                        case CExpOp:
1792                        CPPAD_ASSERT_NARG_NRES(op, 6, 1);
1793                        new_arg[0] = arg[0];
1794                        new_arg[1] = arg[1];
1795                        mask = 1;
1796                        for(i = 2; i < 6; i++)
1797                        {       if( arg[1] & mask )
1798                                {       new_arg[i] = tape[arg[i]].new_var;
1799                                        CPPAD_ASSERT_UNKNOWN( 
1800                                                size_t(new_arg[i]) < num_var
1801                                        );
1802                                }
1803                                else    new_arg[i] = rec->PutPar( 
1804                                                play->GetPar( arg[i] )
1805                                );
1806                                mask = mask << 1;
1807                        }
1808                        rec->PutArg(
1809                                new_arg[0] ,
1810                                new_arg[1] ,
1811                                new_arg[2] ,
1812                                new_arg[3] ,
1813                                new_arg[4] ,
1814                                new_arg[5] 
1815                        );
1816                        tape[i_var].new_var = rec->PutOp(op);
1817                        break;
1818                        // ---------------------------------------------------
1819                        // Operations with no arguments and no results
1820                        case EndOp:
1821                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
1822                        rec->PutOp(op);
1823                        break;
1824                        // ---------------------------------------------------
1825                        // Operations with no arguments and one result
1826                        case InvOp:
1827                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
1828                        tape[i_var].new_var = rec->PutOp(op);
1829                        break;
1830                        // ---------------------------------------------------
1831                        // Operations with one argument that is a parameter
1832                        case ParOp:
1833                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
1834                        new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
1835
1836                        rec->PutArg( new_arg[0] );
1837                        tape[i_var].new_var = rec->PutOp(op);
1838                        break;
1839                        // ---------------------------------------------------
1840                        // Load using a parameter index
1841                        case LdpOp:
1842                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1843                        new_arg[0] = new_vecad_ind[ arg[0] ];
1844                        new_arg[1] = arg[1];
1845                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1846                        rec->PutArg( 
1847                                new_arg[0], 
1848                                new_arg[1], 
1849                                0
1850                        );
1851                        tape[i_var].new_var = rec->PutOp(op);
1852                        break;
1853                        // ---------------------------------------------------
1854                        // Load using a variable index
1855                        case LdvOp:
1856                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1857                        new_arg[0] = new_vecad_ind[ arg[0] ];
1858                        new_arg[1] = tape[arg[1]].new_var;
1859                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1860                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1861                        rec->PutArg( 
1862                                new_arg[0], 
1863                                new_arg[1], 
1864                                0
1865                        );
1866                        tape[i_var].new_var = rec->PutOp(op);
1867                        break;
1868                        // ---------------------------------------------------
1869                        // Store a parameter using a parameter index
1870                        case StppOp:
1871                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1872                        new_arg[0] = new_vecad_ind[ arg[0] ];
1873                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1874                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1875                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1876                        rec->PutArg(
1877                                new_arg[0], 
1878                                new_arg[1], 
1879                                new_arg[2]
1880                        );
1881                        rec->PutOp(op);
1882                        break;
1883                        // ---------------------------------------------------
1884                        // Store a parameter using a variable index
1885                        case StvpOp:
1886                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1887                        new_arg[0] = new_vecad_ind[ arg[0] ];
1888                        new_arg[1] = tape[arg[1]].new_var;
1889                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1890                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1891                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1892                        rec->PutArg(
1893                                new_arg[0], 
1894                                new_arg[1], 
1895                                new_arg[2]
1896                        );
1897                        rec->PutOp(op);
1898                        break;
1899                        // ---------------------------------------------------
1900                        // Store a variable using a parameter index
1901                        case StpvOp:
1902                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1903                        new_arg[0] = new_vecad_ind[ arg[0] ];
1904                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1905                        new_arg[2] = tape[arg[2]].new_var;
1906                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1907                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1908                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
1909                        rec->PutArg(
1910                                new_arg[0], 
1911                                new_arg[1], 
1912                                new_arg[2]
1913                        );
1914                        rec->PutOp(op);
1915                        break;
1916                        // ---------------------------------------------------
1917                        // Store a variable using a variable index
1918                        case StvvOp:
1919                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1920                        new_arg[0] = new_vecad_ind[ arg[0] ];
1921                        new_arg[1] = tape[arg[1]].new_var;
1922                        new_arg[2] = tape[arg[2]].new_var;
1923                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1924                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1925                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
1926                        rec->PutArg(
1927                                new_arg[0], 
1928                                new_arg[1], 
1929                                new_arg[2]
1930                        );
1931                        rec->PutOp(op);
1932                        break;
1933
1934                        // -----------------------------------------------------------
1935                        case UserOp:
1936                        CPPAD_ASSERT_NARG_NRES(op, 4, 0);
1937                        if( user_state == user_start )
1938                                user_state = user_arg;
1939                        else
1940                        {       user_state = user_start;
1941                                user_keep.pop();       
1942                        }
1943                        // user_index, user_id, user_n, user_m
1944                        rec->PutArg(arg[0], arg[1], arg[2], arg[3]);
1945                        rec->PutOp(UserOp);
1946                        break;
1947
1948                        case UsrapOp:
1949                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1950                        new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
1951                        rec->PutArg(new_arg[0]);
1952                        rec->PutOp(UsrapOp);
1953                        break;
1954
1955                        case UsravOp:
1956                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1957                        new_arg[0] = tape[arg[0]].new_var;
1958                        if( size_t(new_arg[0]) < num_var )
1959                        {       rec->PutArg(new_arg[0]);
1960                                rec->PutOp(UsravOp);
1961                        }
1962                        else
1963                        {       // This argument does not affect the result and
1964                                // has been optimized out so use nan in its place.
1965                                new_arg[0] = rec->PutPar( nan(Base(0)) );
1966                                rec->PutArg(new_arg[0]);
1967                                rec->PutOp(UsrapOp);
1968                        }
1969                        break;
1970
1971                        case UsrrpOp:
1972                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1973                        new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
1974                        rec->PutArg(new_arg[0]);
1975                        rec->PutOp(UsrrpOp);
1976                        break;
1977                       
1978                        case UsrrvOp:
1979                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
1980                        tape[i_var].new_var = rec->PutOp(UsrrvOp);
1981                        break;
1982                        // ---------------------------------------------------
1983
1984                        // all cases should be handled above
1985                        default:
1986                        CPPAD_ASSERT_UNKNOWN(false);
1987
1988                }
1989                if( replace_hash )
1990                {       // The old variable index i_var corresponds to the
1991                        // new variable index tape[i_var].new_var. In addition
1992                        // this is the most recent variable that has this code.
1993                        hash_table_var[code] = i_var;
1994                }
1995
1996        }
1997        // modify the dependent variable vector to new indices
1998        for(i = 0; i < dep_taddr.size(); i++ )
1999        {       CPPAD_ASSERT_UNKNOWN( size_t(tape[ dep_taddr[i] ].new_var) < num_var );
2000                dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
2001        }
2002}
2003
2004/*!
2005Optimize a player object operation sequence
2006
2007The operation sequence for this object is replaced by one with fewer operations
2008but the same funcition and derivative values.
2009
2010\tparam Base
2011base type for the operator; i.e., this operation was recorded
2012using AD< \a Base > and computations by this routine are done using type
2013\a Base.
2014
2015*/
2016template <class Base>
2017void ADFun<Base>::optimize(void)
2018{       // place to store the optimized version of the recording
2019        recorder<Base> rec;
2020
2021        // number of independent variables
2022        size_t n = ind_taddr_.size();
2023
2024# ifndef NDEBUG
2025        size_t i, j, m = dep_taddr_.size();
2026        CppAD::vector<Base> x(n), y(m), check(m);
2027        bool check_zero_order = taylor_per_var_ > 0;
2028        Base max_taylor(0);
2029        if( check_zero_order )
2030        {       // zero order coefficients for independent vars
2031                for(j = 0; j < n; j++)
2032                {       CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
2033                        CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
2034                        x[j] = taylor_[ ind_taddr_[j] * taylor_col_dim_ + 0];
2035                }
2036                // zero order coefficients for dependent vars
2037                for(i = 0; i < m; i++)
2038                {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
2039                        y[i] = taylor_[ dep_taddr_[i] * taylor_col_dim_ + 0];
2040                }
2041                // maximum zero order coefficient not counting BeginOp at beginning
2042                // (which is correpsonds to uninitialized memory).
2043                for(i = 1; i < total_num_var_; i++)
2044                {       if(  abs_geq(taylor_[i*taylor_col_dim_+0] , max_taylor) )
2045                                max_taylor = taylor_[i*taylor_col_dim_+0];
2046                }
2047        }
2048# endif
2049
2050        // create the optimized recording
2051        CppAD::optimize<Base>(n, dep_taddr_, &play_, &rec);
2052
2053        // number of variables in the recording
2054        total_num_var_ = rec.num_rec_var();
2055
2056        // now replace the recording
2057        play_.get(rec);
2058
2059        // free memory allocated for sparse Jacobian calculation
2060        // (the results are no longer valid)
2061        for_jac_sparse_pack_.resize(0, 0);
2062        for_jac_sparse_set_.resize(0,0);
2063
2064        // free old Taylor coefficient memory
2065        taylor_.free();
2066        taylor_per_var_ = 0;
2067        taylor_col_dim_ = 0;
2068
2069# ifndef NDEBUG
2070        if( check_zero_order )
2071        {
2072                // zero order forward calculation using new operation sequence
2073                check = Forward(0, x);
2074
2075                // check results
2076                Base eps = 10. * epsilon<Base>();
2077                for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
2078                        abs_geq( eps * max_taylor , check[i] - y[i] ) ,
2079                        "Error during check of f.optimize()."
2080                );
2081
2082                // Erase memory that this calculation was done so NDEBUG gives
2083                // same final state for this object (from users perspective)
2084                taylor_per_var_ = 0;
2085        }
2086# endif
2087}
2088
2089/*! \} */
2090} // END_CPPAD_NAMESPACE
2091# endif
Note: See TracBrowser for help on using the repository browser.