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

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

change plan: use a bool for each operation that specifies skip or not

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