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

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

optimize.hpp: change order to make it easier to handle cexp_connected.

  • Property svn:keywords set to Id
File size: 61.5 KB
Line 
1/* $Id: optimize.hpp 2956 2013-10-16 15:46:09Z 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 i_var is sum_connected to and
1277                                // it is the result of a summation.
1278                                if( tape[i_var].connect_type == sum_connected )
1279                                        tape[i_var].connect_type = csum_connected;
1280
1281                                // check for case where arg[0] is sum_connected
1282                                if( tape[arg[0]].connect_type == not_connected )
1283                                        tape[arg[0]].connect_type = sum_connected;
1284                                else    tape[arg[0]].connect_type = yes_connected;
1285
1286                        }
1287                        break; // --------------------------------------------
1288               
1289                        // Special case for AddpvOp and SubpvOp
1290                        case AddpvOp:
1291                        case SubpvOp:
1292                        if( tape[i_var].connect_type != not_connected )
1293                        {
1294                                // check for case where i_var is sum_connected to and
1295                                // it is the result of a summation.
1296                                if( tape[i_var].connect_type == sum_connected )
1297                                        tape[i_var].connect_type = csum_connected;
1298
1299                                // check for case where arg[1] is sum_connected
1300                                if( tape[arg[1]].connect_type == not_connected )
1301                                        tape[arg[1]].connect_type = sum_connected;
1302                                else    tape[arg[1]].connect_type = yes_connected;
1303
1304                        }
1305                        break; // --------------------------------------------
1306
1307               
1308                        // Special case for AddvvOp and SubvvOp
1309                        case AddvvOp:
1310                        case SubvvOp:
1311                        if( tape[i_var].connect_type != not_connected )
1312                        {
1313                                // check for case where i_var is sum_connected to and
1314                                // it is the result of a summation.
1315                                if( tape[i_var].connect_type == sum_connected )
1316                                        tape[i_var].connect_type = csum_connected;
1317
1318                                // check for case where arg[0] is sum_connected
1319                                if( tape[arg[0]].connect_type == not_connected )
1320                                        tape[arg[0]].connect_type = sum_connected;
1321                                else    tape[arg[0]].connect_type = yes_connected;
1322
1323                                // check for case where arg[1] is sum_connected
1324                                if( tape[arg[1]].connect_type == not_connected )
1325                                        tape[arg[1]].connect_type = sum_connected;
1326                                else    tape[arg[1]].connect_type = yes_connected;
1327                        }
1328                        break; // --------------------------------------------
1329
1330                        // Other binary operators
1331                        // where operands are arg[0], arg[1]
1332                        case DivvvOp:
1333                        case MulvvOp:
1334                        case PowvvOp:
1335                        if( tape[i_var].connect_type != not_connected )
1336                        {
1337                                tape[arg[0]].connect_type = yes_connected;
1338                                tape[arg[1]].connect_type = yes_connected;
1339                        }
1340                        break; // --------------------------------------------
1341
1342                        // Conditional expression operators
1343                        case CExpOp:
1344                        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
1345                        if( tape[i_var].connect_type != not_connected )
1346                        {
1347                                mask = 1;
1348                                for(i = 2; i < 6; i++)
1349                                {       if( arg[1] & mask )
1350                                        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[i]) < i_var );
1351                                                tape[arg[i]].connect_type = yes_connected;
1352                                        }
1353                                        mask = mask << 1;
1354                                }
1355                        }
1356                        break;  // --------------------------------------------
1357
1358                        // Operations where there is noting to do
1359                        case BeginOp:
1360                        case ComOp:
1361                        case EndOp:
1362                        case InvOp:
1363                        case ParOp:
1364                        case PriOp:
1365                        break;  // --------------------------------------------
1366
1367                        // Load using a parameter index
1368                        case LdpOp:
1369                        if( tape[i_var].connect_type != not_connected )
1370                        {
1371                                i                = vecad[ arg[0] - 1 ];
1372                                vecad_connect[i] = yes_connected;
1373                        }
1374                        break; // --------------------------------------------
1375
1376                        // Load using a variable index
1377                        case LdvOp:
1378                        if( tape[i_var].connect_type != not_connected )
1379                        {
1380                                i                    = vecad[ arg[0] - 1 ];
1381                                vecad_connect[i]     = yes_connected;
1382                                tape[arg[1]].connect_type = yes_connected;
1383                        }
1384                        break; // --------------------------------------------
1385
1386                        // Store a variable using a parameter index
1387                        case StpvOp:
1388                        i = vecad[ arg[0] - 1 ];
1389                        if( vecad_connect[i] != not_connected )
1390                                tape[arg[2]].connect_type = yes_connected;
1391                        break; // --------------------------------------------
1392
1393                        // Store a variable using a variable index
1394                        case StvvOp:
1395                        i = vecad[ arg[0] - 1 ];
1396                        if( vecad_connect[i] )
1397                        {       tape[arg[1]].connect_type = yes_connected;
1398                                tape[arg[2]].connect_type = yes_connected;
1399                        }
1400                        break; 
1401                        // ============================================================
1402                        case UserOp:
1403                        // start or end atomic operation sequence
1404                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
1405                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
1406                        if( user_state == user_end )
1407                        {       user_index = arg[0];
1408                                user_id    = arg[1];
1409                                user_n     = arg[2];
1410                                user_m     = arg[3];
1411                                user_q     = 1;
1412                                if(user_s.size() != user_n )
1413                                        user_s.resize(user_n);
1414                                if(user_r.size() != user_m )
1415                                        user_r.resize(user_m);
1416                                user_j     = user_n;
1417                                user_i     = user_m;
1418                                user_state = user_ret;
1419                                user_keep.push(false);
1420                        }
1421                        else
1422                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
1423                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
1424                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
1425                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
1426                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
1427                                user_state = user_end;
1428               }
1429                        break;
1430
1431                        case UsrapOp:
1432                        // parameter argument in an atomic operation sequence
1433                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
1434                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
1435                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1436                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
1437                        --user_j;
1438                        if( user_j == 0 )
1439                                user_state = user_start;
1440                        break;
1441
1442                        case UsravOp:
1443                        // variable argument in an atomic operation sequence
1444                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
1445                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
1446                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1447                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
1448                        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
1449                        --user_j;
1450                        if( ! user_s[user_j].empty() )
1451                        {       tape[arg[0]].connect_type = yes_connected;
1452                                user_keep.top() = true;
1453                        }
1454                        if( user_j == 0 )
1455                                user_state = user_start;
1456                        break;
1457
1458                        case UsrrpOp:
1459                        // parameter result in an atomic operation sequence
1460                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
1461                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
1462                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1463                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
1464                        --user_i;
1465                        user_r[user_i].clear();
1466                        if( user_i == 0 )
1467                        {       // call users function for this operation
1468                                atomic_base<Base>* atom = 
1469                                        atomic_base<Base>::class_object(user_index);
1470                                atom->set_id(user_id);
1471                                atom->rev_sparse_jac(
1472                                        user_q, user_r, user_s
1473                                );
1474                                user_state = user_arg;
1475                        }
1476                        break;
1477
1478                        case UsrrvOp:
1479                        // variable result in an atomic operation sequence
1480                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
1481                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
1482                        --user_i;
1483                        user_r[user_i].clear();
1484                        if( tape[i_var].connect_type != not_connected )
1485                                user_r[user_i].insert(0);
1486                        if( user_i == 0 )
1487                        {       // call users function for this operation
1488                                atomic_base<Base>* atom = 
1489                                        atomic_base<Base>::class_object(user_index);
1490                                atom->set_id(user_id);
1491                                atom->rev_sparse_jac(
1492                                        user_q, user_r, user_s
1493                                );
1494                                user_state = user_arg;
1495                        }
1496                        break;
1497                        // ============================================================
1498
1499                        // all cases should be handled above
1500                        default:
1501                        CPPAD_ASSERT_UNKNOWN(0);
1502                }
1503        }
1504        // values corresponding to BeginOp
1505        CPPAD_ASSERT_UNKNOWN( i_op == 0 && i_var == 0 && op == BeginOp );
1506        tape[i_var].op = op;
1507        // -------------------------------------------------------------
1508
1509        // Erase all information in the recording
1510        rec->free();
1511
1512        // Initilaize table mapping hash code to variable index in tape
1513        // as pointing to the BeginOp at the beginning of the tape
1514        CppAD::vector<size_t>  hash_table_var(CPPAD_HASH_TABLE_SIZE);
1515        for(i = 0; i < CPPAD_HASH_TABLE_SIZE; i++)
1516                hash_table_var[i] = 0;
1517        CPPAD_ASSERT_UNKNOWN( tape[0].op == BeginOp );
1518
1519        // initialize mapping from old variable index to new variable index
1520        for(i = 0; i < num_var; i++)
1521                tape[i].new_var = num_var; // invalid index
1522       
1523
1524        // initialize mapping from old VecAD index to new VecAD index
1525        CppAD::vector<size_t> new_vecad_ind(num_vecad_ind);
1526        for(i = 0; i < num_vecad_ind; i++)
1527                new_vecad_ind[i] = num_vecad_ind; // invalid index
1528
1529        j = 0;     // index into the old set of indices
1530        for(i = 0; i < num_vecad_vec; i++)
1531        {       // length of this VecAD
1532                size_t length = play->GetVecInd(j);
1533                if( vecad_connect[i] != not_connected )
1534                {       // Put this VecAD vector in new recording
1535                        CPPAD_ASSERT_UNKNOWN(length < num_vecad_ind);
1536                        new_vecad_ind[j] = rec->PutVecInd(length);
1537                        for(k = 1; k <= length; k++) new_vecad_ind[j+k] =
1538                                rec->PutVecInd(
1539                                        rec->PutPar(
1540                                                play->GetPar( 
1541                                                        play->GetVecInd(j+k)
1542                        ) ) );
1543                }
1544                // start of next VecAD
1545                j       += length + 1;
1546        }
1547        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
1548
1549        // start playing the operations in the forward direction
1550        play->start_forward(op, arg, i_op, i_var);
1551
1552        // playing forward skips BeginOp at the beginning, but not EndOp at
1553        // the end.  Put BeginOp at beginning of recording
1554        CPPAD_ASSERT_UNKNOWN( op == BeginOp );
1555        CPPAD_ASSERT_NARG_NRES(BeginOp, 0, 1);
1556        tape[i_var].new_var = rec->PutOp(BeginOp);
1557
1558        // temporary buffer for new argument values
1559        addr_t new_arg[6];
1560
1561        // temporary work space used by optimize_record_csum
1562        // (decalared here to avoid realloaction of memory)
1563        optimize_csum_stacks csum_work;
1564
1565        user_state = user_start;
1566        while(op != EndOp)
1567        {       // next op
1568                play->next_forward(op, arg, i_op, i_var);
1569                CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
1570                CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
1571
1572                // determine if we should keep this operation in the new
1573                // operation sequence
1574                bool keep;
1575                switch( op )
1576                {       case ComOp:
1577                        case PriOp:
1578                        keep = false;
1579                        break;
1580
1581                        case InvOp:
1582                        case EndOp:
1583                        keep = true;
1584                        break;
1585
1586                        case StppOp:
1587                        case StvpOp:
1588                        case StpvOp:
1589                        case StvvOp:
1590                        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
1591                        i = vecad[ arg[0] - 1 ];
1592                        keep = vecad_connect[i] != not_connected;
1593                        break;
1594
1595                        case AddpvOp:
1596                        case AddvvOp:
1597                        case SubpvOp:
1598                        case SubvpOp:
1599                        case SubvvOp:
1600                        keep  = tape[i_var].connect_type != not_connected;
1601                        keep &= tape[i_var].connect_type != csum_connected;
1602                        break; 
1603
1604                        case UserOp:
1605                        case UsrapOp:
1606                        case UsravOp:
1607                        case UsrrpOp:
1608                        case UsrrvOp:
1609                        keep = user_keep.top();
1610                        break;
1611
1612                        default:
1613                        keep = tape[i_var].connect_type != not_connected;
1614                        break;
1615                }
1616
1617                size_t         match_var    = 0;
1618                unsigned short code         = 0;
1619                bool           replace_hash = false;
1620                if( keep ) switch( op )
1621                {
1622                        // Unary operator where operand is arg[0]
1623                        case AbsOp:
1624                        case AcosOp:
1625                        case AsinOp:
1626                        case AtanOp:
1627                        case CosOp:
1628                        case CoshOp:
1629                        case ExpOp:
1630                        case LogOp:
1631                        case SinOp:
1632                        case SinhOp:
1633                        case SqrtOp:
1634                        case TanOp:
1635                        case TanhOp:
1636                        match_var = optimize_unary_match(
1637                                tape                ,  // inputs
1638                                i_var               ,
1639                                play->num_rec_par() ,
1640                                play->GetPar()      ,
1641                                hash_table_var      ,
1642                                code                  // outputs
1643                        );
1644                        if( match_var > 0 )
1645                                tape[i_var].new_var = match_var;
1646                        else
1647                        {
1648                                replace_hash = true;
1649                                new_arg[0]   = tape[ arg[0] ].new_var;
1650                                rec->PutArg( new_arg[0] );
1651                                i                   = rec->PutOp(op);
1652                                tape[i_var].new_var = i;
1653                                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
1654                        }
1655                        break;
1656                        // ---------------------------------------------------
1657                        // Binary operators where
1658                        // left is a variable and right is a parameter
1659                        case SubvpOp:
1660                        if( tape[arg[0]].connect_type == csum_connected )
1661                        {
1662                                // convert to a sequence of summation operators
1663                                tape[i_var].new_var = optimize_record_csum(
1664                                        tape                , // inputs
1665                                        i_var               ,
1666                                        play->num_rec_par() ,
1667                                        play->GetPar()      ,
1668                                        rec                 ,
1669                                        csum_work
1670                                );
1671                                // abort rest of this case
1672                                break;
1673                        }
1674                        case DivvpOp:
1675                        case PowvpOp:
1676                        match_var = optimize_binary_match(
1677                                tape                ,  // inputs
1678                                i_var               ,
1679                                play->num_rec_par() ,
1680                                play->GetPar()      ,
1681                                hash_table_var      ,
1682                                code                  // outputs
1683                        );
1684                        if( match_var > 0 )
1685                                tape[i_var].new_var = match_var;
1686                        else
1687                        {       tape[i_var].new_var = optimize_record_vp(
1688                                        tape                , // inputs
1689                                        i_var               ,
1690                                        play->num_rec_par() ,
1691                                        play->GetPar()      ,
1692                                        rec                 ,
1693                                        op                  ,
1694                                        arg
1695                                );
1696                                replace_hash = true;
1697                        }
1698                        break;
1699                        // ---------------------------------------------------
1700                        // Binary operators where
1701                        // left is a parameter and right is a variable
1702                        case SubpvOp:
1703                        case AddpvOp:
1704                        if( tape[arg[1]].connect_type == csum_connected )
1705                        {
1706                                // convert to a sequence of summation operators
1707                                tape[i_var].new_var = optimize_record_csum(
1708                                        tape                , // inputs
1709                                        i_var               ,
1710                                        play->num_rec_par() ,
1711                                        play->GetPar()      ,
1712                                        rec                 ,
1713                                        csum_work
1714                                );
1715                                // abort rest of this case
1716                                break;
1717                        }
1718                        case DivpvOp:
1719                        case MulpvOp:
1720                        case PowpvOp:
1721                        match_var = optimize_binary_match(
1722                                tape                ,  // inputs
1723                                i_var               ,
1724                                play->num_rec_par() ,
1725                                play->GetPar()      ,
1726                                hash_table_var      ,
1727                                code                  // outputs
1728                        );
1729                        if( match_var > 0 )
1730                                tape[i_var].new_var = match_var;
1731                        else
1732                        {       tape[i_var].new_var = optimize_record_pv(
1733                                        tape                , // inputs
1734                                        i_var               ,
1735                                        play->num_rec_par() ,
1736                                        play->GetPar()      ,
1737                                        rec                 ,
1738                                        op                  ,
1739                                        arg
1740                                );
1741                                replace_hash = true;
1742                        }
1743                        break;
1744                        // ---------------------------------------------------
1745                        // Binary operator where
1746                        // both operators are variables
1747                        case AddvvOp:
1748                        case SubvvOp:
1749                        if( (tape[arg[0]].connect_type == csum_connected) |
1750                            (tape[arg[1]].connect_type == csum_connected)
1751                        )
1752                        {
1753                                // convert to a sequence of summation operators
1754                                tape[i_var].new_var = optimize_record_csum(
1755                                        tape                , // inputs
1756                                        i_var               ,
1757                                        play->num_rec_par() ,
1758                                        play->GetPar()      ,
1759                                        rec                 ,
1760                                        csum_work
1761                                );
1762                                // abort rest of this case
1763                                break;
1764                        }
1765                        case DivvvOp:
1766                        case MulvvOp:
1767                        case PowvvOp:
1768                        match_var = optimize_binary_match(
1769                                tape                ,  // inputs
1770                                i_var               ,
1771                                play->num_rec_par() ,
1772                                play->GetPar()      ,
1773                                hash_table_var      ,
1774                                code                  // outputs
1775                        );
1776                        if( match_var > 0 )
1777                                tape[i_var].new_var = match_var;
1778                        else
1779                        {       tape[i_var].new_var = optimize_record_vv(
1780                                        tape                , // inputs
1781                                        i_var               ,
1782                                        play->num_rec_par() ,
1783                                        play->GetPar()      ,
1784                                        rec                 ,
1785                                        op                  ,
1786                                        arg
1787                                );
1788                                replace_hash = true;
1789                        }
1790                        break;
1791                        // ---------------------------------------------------
1792                        // Conditional expression operators
1793                        case CExpOp:
1794                        CPPAD_ASSERT_NARG_NRES(op, 6, 1);
1795                        new_arg[0] = arg[0];
1796                        new_arg[1] = arg[1];
1797                        mask = 1;
1798                        for(i = 2; i < 6; i++)
1799                        {       if( arg[1] & mask )
1800                                {       new_arg[i] = tape[arg[i]].new_var;
1801                                        CPPAD_ASSERT_UNKNOWN( 
1802                                                size_t(new_arg[i]) < num_var
1803                                        );
1804                                }
1805                                else    new_arg[i] = rec->PutPar( 
1806                                                play->GetPar( arg[i] )
1807                                );
1808                                mask = mask << 1;
1809                        }
1810                        rec->PutArg(
1811                                new_arg[0] ,
1812                                new_arg[1] ,
1813                                new_arg[2] ,
1814                                new_arg[3] ,
1815                                new_arg[4] ,
1816                                new_arg[5] 
1817                        );
1818                        tape[i_var].new_var = rec->PutOp(op);
1819                        break;
1820                        // ---------------------------------------------------
1821                        // Operations with no arguments and no results
1822                        case EndOp:
1823                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
1824                        rec->PutOp(op);
1825                        break;
1826                        // ---------------------------------------------------
1827                        // Operations with no arguments and one result
1828                        case InvOp:
1829                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
1830                        tape[i_var].new_var = rec->PutOp(op);
1831                        break;
1832                        // ---------------------------------------------------
1833                        // Operations with one argument that is a parameter
1834                        case ParOp:
1835                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
1836                        new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
1837
1838                        rec->PutArg( new_arg[0] );
1839                        tape[i_var].new_var = rec->PutOp(op);
1840                        break;
1841                        // ---------------------------------------------------
1842                        // Load using a parameter index
1843                        case LdpOp:
1844                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1845                        new_arg[0] = new_vecad_ind[ arg[0] ];
1846                        new_arg[1] = arg[1];
1847                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1848                        rec->PutArg( 
1849                                new_arg[0], 
1850                                new_arg[1], 
1851                                0
1852                        );
1853                        tape[i_var].new_var = rec->PutOp(op);
1854                        break;
1855                        // ---------------------------------------------------
1856                        // Load using a variable index
1857                        case LdvOp:
1858                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1859                        new_arg[0] = new_vecad_ind[ arg[0] ];
1860                        new_arg[1] = tape[arg[1]].new_var;
1861                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1862                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1863                        rec->PutArg( 
1864                                new_arg[0], 
1865                                new_arg[1], 
1866                                0
1867                        );
1868                        tape[i_var].new_var = rec->PutOp(op);
1869                        break;
1870                        // ---------------------------------------------------
1871                        // Store a parameter using a parameter index
1872                        case StppOp:
1873                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1874                        new_arg[0] = new_vecad_ind[ arg[0] ];
1875                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1876                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1877                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1878                        rec->PutArg(
1879                                new_arg[0], 
1880                                new_arg[1], 
1881                                new_arg[2]
1882                        );
1883                        rec->PutOp(op);
1884                        break;
1885                        // ---------------------------------------------------
1886                        // Store a parameter using a variable index
1887                        case StvpOp:
1888                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1889                        new_arg[0] = new_vecad_ind[ arg[0] ];
1890                        new_arg[1] = tape[arg[1]].new_var;
1891                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1892                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1893                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1894                        rec->PutArg(
1895                                new_arg[0], 
1896                                new_arg[1], 
1897                                new_arg[2]
1898                        );
1899                        rec->PutOp(op);
1900                        break;
1901                        // ---------------------------------------------------
1902                        // Store a variable using a parameter index
1903                        case StpvOp:
1904                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1905                        new_arg[0] = new_vecad_ind[ arg[0] ];
1906                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1907                        new_arg[2] = tape[arg[2]].new_var;
1908                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1909                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1910                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
1911                        rec->PutArg(
1912                                new_arg[0], 
1913                                new_arg[1], 
1914                                new_arg[2]
1915                        );
1916                        rec->PutOp(op);
1917                        break;
1918                        // ---------------------------------------------------
1919                        // Store a variable using a variable index
1920                        case StvvOp:
1921                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1922                        new_arg[0] = new_vecad_ind[ arg[0] ];
1923                        new_arg[1] = tape[arg[1]].new_var;
1924                        new_arg[2] = tape[arg[2]].new_var;
1925                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1926                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1927                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
1928                        rec->PutArg(
1929                                new_arg[0], 
1930                                new_arg[1], 
1931                                new_arg[2]
1932                        );
1933                        rec->PutOp(op);
1934                        break;
1935
1936                        // -----------------------------------------------------------
1937                        case UserOp:
1938                        CPPAD_ASSERT_NARG_NRES(op, 4, 0);
1939                        if( user_state == user_start )
1940                                user_state = user_arg;
1941                        else
1942                        {       user_state = user_start;
1943                                user_keep.pop();       
1944                        }
1945                        // user_index, user_id, user_n, user_m
1946                        rec->PutArg(arg[0], arg[1], arg[2], arg[3]);
1947                        rec->PutOp(UserOp);
1948                        break;
1949
1950                        case UsrapOp:
1951                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1952                        new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
1953                        rec->PutArg(new_arg[0]);
1954                        rec->PutOp(UsrapOp);
1955                        break;
1956
1957                        case UsravOp:
1958                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1959                        new_arg[0] = tape[arg[0]].new_var;
1960                        if( size_t(new_arg[0]) < num_var )
1961                        {       rec->PutArg(new_arg[0]);
1962                                rec->PutOp(UsravOp);
1963                        }
1964                        else
1965                        {       // This argument does not affect the result and
1966                                // has been optimized out so use nan in its place.
1967                                new_arg[0] = rec->PutPar( nan(Base(0)) );
1968                                rec->PutArg(new_arg[0]);
1969                                rec->PutOp(UsrapOp);
1970                        }
1971                        break;
1972
1973                        case UsrrpOp:
1974                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1975                        new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
1976                        rec->PutArg(new_arg[0]);
1977                        rec->PutOp(UsrrpOp);
1978                        break;
1979                       
1980                        case UsrrvOp:
1981                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
1982                        tape[i_var].new_var = rec->PutOp(UsrrvOp);
1983                        break;
1984                        // ---------------------------------------------------
1985
1986                        // all cases should be handled above
1987                        default:
1988                        CPPAD_ASSERT_UNKNOWN(false);
1989
1990                }
1991                if( replace_hash )
1992                {       // The old variable index i_var corresponds to the
1993                        // new variable index tape[i_var].new_var. In addition
1994                        // this is the most recent variable that has this code.
1995                        hash_table_var[code] = i_var;
1996                }
1997
1998        }
1999        // modify the dependent variable vector to new indices
2000        for(i = 0; i < dep_taddr.size(); i++ )
2001        {       CPPAD_ASSERT_UNKNOWN( size_t(tape[ dep_taddr[i] ].new_var) < num_var );
2002                dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
2003        }
2004}
2005
2006/*!
2007Optimize a player object operation sequence
2008
2009The operation sequence for this object is replaced by one with fewer operations
2010but the same funcition and derivative values.
2011
2012\tparam Base
2013base type for the operator; i.e., this operation was recorded
2014using AD< \a Base > and computations by this routine are done using type
2015\a Base.
2016
2017*/
2018template <class Base>
2019void ADFun<Base>::optimize(void)
2020{       // place to store the optimized version of the recording
2021        recorder<Base> rec;
2022
2023        // number of independent variables
2024        size_t n = ind_taddr_.size();
2025
2026# ifndef NDEBUG
2027        size_t i, j, m = dep_taddr_.size();
2028        CppAD::vector<Base> x(n), y(m), check(m);
2029        bool check_zero_order = taylor_per_var_ > 0;
2030        Base max_taylor(0);
2031        if( check_zero_order )
2032        {       // zero order coefficients for independent vars
2033                for(j = 0; j < n; j++)
2034                {       CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
2035                        CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
2036                        x[j] = taylor_[ ind_taddr_[j] * taylor_col_dim_ + 0];
2037                }
2038                // zero order coefficients for dependent vars
2039                for(i = 0; i < m; i++)
2040                {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
2041                        y[i] = taylor_[ dep_taddr_[i] * taylor_col_dim_ + 0];
2042                }
2043                // maximum zero order coefficient not counting BeginOp at beginning
2044                // (which is correpsonds to uninitialized memory).
2045                for(i = 1; i < total_num_var_; i++)
2046                {       if(  abs_geq(taylor_[i*taylor_col_dim_+0] , max_taylor) )
2047                                max_taylor = taylor_[i*taylor_col_dim_+0];
2048                }
2049        }
2050# endif
2051
2052        // create the optimized recording
2053        CppAD::optimize<Base>(n, dep_taddr_, &play_, &rec);
2054
2055        // number of variables in the recording
2056        total_num_var_ = rec.num_rec_var();
2057
2058        // now replace the recording
2059        play_.get(rec);
2060
2061        // free memory allocated for sparse Jacobian calculation
2062        // (the results are no longer valid)
2063        for_jac_sparse_pack_.resize(0, 0);
2064        for_jac_sparse_set_.resize(0,0);
2065
2066        // free old Taylor coefficient memory
2067        taylor_.free();
2068        taylor_per_var_ = 0;
2069        taylor_col_dim_ = 0;
2070
2071# ifndef NDEBUG
2072        if( check_zero_order )
2073        {
2074                // zero order forward calculation using new operation sequence
2075                check = Forward(0, x);
2076
2077                // check results
2078                Base eps = 10. * epsilon<Base>();
2079                for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
2080                        abs_geq( eps * max_taylor , check[i] - y[i] ) ,
2081                        "Error during check of f.optimize()."
2082                );
2083
2084                // Erase memory that this calculation was done so NDEBUG gives
2085                // same final state for this object (from users perspective)
2086                taylor_per_var_ = 0;
2087        }
2088# endif
2089}
2090
2091/*! \} */
2092} // END_CPPAD_NAMESPACE
2093# endif
Note: See TracBrowser for help on using the repository browser.