source: trunk/cppad/local/optimize.hpp @ 3367

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

Fix bug in conditional skipping of calls to atomic functions.

opt_atomic.sh: test fix here.
optimize.hpp: fix bug here.

  • Property svn:keywords set to Id
File size: 79.7 KB
Line 
1/* $Id: optimize.hpp 3367 2014-09-27 11:55:27Z bradbell $ */
2# ifndef CPPAD_OPTIMIZE_INCLUDED
3# define CPPAD_OPTIMIZE_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16/*
17$begin optimize$$
18$spell
19        jac
20        bool
21        Taylor
22        var
23        CppAD
24        cppad
25        std
26        CondExpEq
27$$
28
29$section Optimize an ADFun Object Tape$$
30
31$index optimize$$
32$index tape, optimize$$
33$index sequence, optimize operations$$
34$index operations, optimize sequence$$
35$index speed, optimize$$
36$index memory, optimize$$
37
38$head Syntax$$
39$icode%f%.optimize()%$$
40
41
42$head Purpose$$
43The operation sequence corresponding to an $cref ADFun$$ object can
44be very large and involve many operations; see the
45size functions in $cref seq_property$$.
46The $icode%f%.optimize%$$ procedure reduces the number of operations,
47and thereby the time and the memory, required to
48compute function and derivative values.
49
50$head f$$
51The object $icode f$$ has prototype
52$codei%
53        ADFun<%Base%> %f%
54%$$
55
56$head Improvements$$
57You can see the reduction in number of variables in the operation sequence
58by calling the function $cref/f.size_var()/seq_property/size_var/$$
59before and after the optimization procedure.
60Given that the optimization procedure takes time,
61it may be faster to skip this optimize procedure and just compute
62derivatives using the original operation sequence.
63
64$subhead Testing$$
65You can run the CppAD $cref/speed/speed_main/$$ tests and see
66the corresponding changes in number of variables and execution time;
67see $cref cmake_check$$.
68
69$head Efficiency$$
70The $code optimize$$ member function
71may greatly reduce the number of variables
72in the operation sequence; see $cref/size_var/seq_property/size_var/$$.
73If a $cref/zero order forward/forward_zero/$$ calculation is done during
74the construction of $icode f$$, it will require more memory
75and time than required after the optimization procedure.
76In addition, it will need to be redone.
77For this reason, it is more efficient to use
78$codei%
79        ADFun<%Base%> %f%;
80        %f%.Dependent(%x%, %y%);
81        %f%.optimize();
82%$$
83instead of
84$codei%
85        ADFun<%Base%> %f%(%x%, %y%)
86        %f%.optimize();
87%$$
88See the discussion about
89$cref/sequence constructors/FunConstruct/Sequence Constructor/$$.
90
91$head Comparison Operators$$
92Any comparison operators that are in the tape are removed by this operation.
93Hence the return value of $cref CompareChange$$ will always be zero
94for an optimized tape (even if $code NDEBUG$$ is not defined).
95
96$head Atomic Functions$$
97There are some subtitle issue with optimized $cref atomic$$ functions
98$latex v = g(u)$$:
99
100$subhead rev_sparse_jac$$
101The $cref atomic_rev_sparse_jac$$ function is be used to determine
102which components of $icode u$$ affect the dependent variables of $icode f$$.
103The current setting of the
104$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for each
105atomic function is used to determine if the $code bool$$ or
106$code std::set<size_t>$$ version of $cref atomic_rev_sparse_jac$$ is used.
107
108$subhead nan$$
109If $icode%u%[%i%]%$$ does not affect the value of
110the dependent variables for $icode f$$,
111the value of $icode%u%[%i%]%$$ is set to $cref nan$$.
112
113
114$head Checking Optimization$$
115$index NDEBUG$$
116If $cref/NDEBUG/Faq/Speed/NDEBUG/$$ is not defined,
117and $cref/f.size_order()/size_order/$$ is greater than zero,
118a $cref forward_zero$$ calculation is done using the optimized version
119of $icode f$$ and the results are checked to see that they are
120the same as before.
121If they are not the same, the
122$cref ErrorHandler$$ is called with a known error message
123related to $icode%f%.optimize()%$$.
124
125$head Example$$
126$children%
127        example/optimize.cpp
128%$$
129The file
130$cref optimize.cpp$$
131contains an example and test of this operation.
132It returns true if it succeeds and false otherwise.
133
134$end
135-----------------------------------------------------------------------------
136*/
137# include <stack>
138
139namespace CppAD { // BEGIN_CPPAD_NAMESPACE
140namespace optimize { // BEGIN_CPPAD_OPTIMIZE_NAMESPACE
141/*!
142\file optimize.hpp
143Routines for optimizing a tape
144*/
145
146
147/*!
148State for this variable set during reverse sweep.
149*/
150enum enum_connect_type {
151        /// There is no operation that connects this variable to the
152        /// independent variables.
153        not_connected        ,
154
155        /// There is one or more operations that connects this variable to the
156        /// independent variables.
157        yes_connected        ,
158
159        /// There is only one parrent that connects this variable to the
160        /// independent variables and the parent is a summation operation; i.e.,
161        /// AddvvOp, AddpvOp, SubpvOp, SubvpOp, or SubvvOp.
162        sum_connected        ,
163
164        /// Satisfies the sum_connected assumptions above and in addition
165        /// this variable is the result of summation operator.
166        csum_connected       ,
167
168        /// This node is only connected in the case where the comparision is
169        /// true for the conditional expression with index \c connect_index.
170        cexp_connected
171
172};
173
174/*!
175Class used to hold information about one conditional expression.
176*/
177class class_cexp_pair {
178public:
179        /// If this is true (false) this connection is only for the case where
180        /// the comparision in the conditional expression is true (false)
181        bool compare_; 
182        /// This is the index of the conditional expression (in cksip_info)
183        /// for this connection
184        size_t index_;
185        /// constructor
186        class_cexp_pair(const bool& compare, const size_t& index)
187        : compare_(compare), index_(index)
188        { }
189        /// assignment operator
190        void operator=(const class_cexp_pair& right)
191        {       index_   = right.index_;
192                compare_ = right.compare_;
193        }
194        /// not equal operator
195        bool operator!=(const class_cexp_pair& right)
196        {       return (index_ != right.index_) | (compare_ != right.compare_);
197        }
198        /// Less than operator
199        /// (required for intersection of two sets of class_cexp_pair elements).
200        bool operator<(const class_cexp_pair& right) const
201        {       if( index_ == right.index_ )
202                        return size_t(compare_) < size_t(right.compare_);       
203                return index_ < right.index_;
204        }
205};
206/*!
207Compute intersection of two sets of class_cexp_pair elements.
208
209\param left
210first operand of the intersection
211
212\param right
213second operand of the intersection
214
215\result
216the intersection of left and right
217*/
218inline std::set<class_cexp_pair> intersection(
219        std::set<class_cexp_pair>& left  ,
220        std::set<class_cexp_pair>& right )
221{       std::set<class_cexp_pair> result;
222        std::set_intersection(
223                left.begin()  ,
224                left.end()    ,
225                right.begin() ,
226                right.end()   ,
227                std::inserter(result, result.begin())
228        );
229        return result; 
230}
231
232/*!
233Structure used by \c optimize to hold information about one variable.
234in the old operation seqeunce.
235*/
236struct struct_old_variable {
237        /// Operator for which this variable is the result, \c NumRes(op) > 0.
238        /// Set by the reverse sweep at beginning of optimization.
239        OpCode              op;       
240
241        /// Pointer to first argument (child) for this operator.
242        /// Set by the reverse sweep at beginning of optimization.
243        const addr_t*       arg;
244
245        /// How is this variable connected to the independent variables
246        enum_connect_type connect_type; 
247
248        /*!
249        If \c connect_type is \c cexp_connected,
250        this is the corresponding infromation for the conditional connections.
251        */
252        std::set<class_cexp_pair> cexp_set;
253
254        /// New operation sequence corresponding to this old varable.
255        /// Set during forward sweep to the index in the new tape
256        addr_t new_var;
257
258        /// New operator index for this varable.
259        /// Set during forward sweep to the index in the new tape
260        size_t new_op;
261};
262
263struct struct_size_pair {
264        size_t i_op;  // an operator index
265        size_t i_var; // a variable index
266};
267
268/*!
269Structures used by \c record_csum
270to hold information about one variable.
271*/
272struct struct_csum_variable {
273        /// Operator for which this variable is the result, \c NumRes(op) > 0.
274        OpCode              op;       
275
276        /// Pointer to first argument (child) for this operator.
277        /// Set by the reverse sweep at beginning of optimization.
278        const addr_t*       arg;
279
280        /// Is this variable added to the summation
281        /// (if not it is subtracted)
282        bool                add;
283};
284
285/*!
286Structure used to pass work space from \c optimize to \c record_csum
287(so that stacks do not start from zero size every time).
288*/
289struct struct_csum_stacks {
290        /// stack of operations in the cummulative summation
291        std::stack<struct struct_csum_variable>   op_stack;
292        /// stack of variables to be added
293        std::stack<size_t >                         add_stack;
294        /// stack of variables to be subtracted
295        std::stack<size_t >                         sub_stack;
296};
297
298/*!
299CExpOp information that is copied to corresponding CSkipOp
300*/
301struct struct_cskip_info {
302        /// comparision operator
303        CompareOp cop;
304        /// (flag & 1) is true if and only if left is a variable
305        /// (flag & 2) is true if and only if right is a variable
306        size_t flag;
307        /// index for left comparison operand
308        size_t left; 
309        /// index for right comparison operand
310        size_t right; 
311        /// maximum variable index between left and right
312        size_t max_left_right;
313        /// set of variables to skip on true
314        CppAD::vector<size_t> skip_var_true;
315        /// set of variables to skip on false
316        CppAD::vector<size_t> skip_var_false;
317        /// set of operations to skip on true
318        CppAD::vector<size_t> skip_op_true;
319        /// set of operations to skip on false
320        CppAD::vector<size_t> skip_op_false;
321        /// size of skip_op_true
322        size_t n_op_true;
323        /// size of skip_op_false
324        size_t n_op_false;
325        /// index in the argument recording of first argument for this CSkipOp
326        size_t i_arg;
327};
328/*!
329Connection information for a user atomic function
330*/
331struct struct_user_info {
332        /// type of connection for this atomic function
333        enum_connect_type connect_type;
334        /// If this is an conditional connection, this is the information
335        /// of the correpsonding CondExpOp operators
336        std::set<class_cexp_pair> cexp_set;
337        /// If this is a conditional connection, this is the operator
338        /// index of the beginning of the atomic call sequence; i.e.,
339        /// the first UserOp.
340        size_t op_begin;
341        /// If this is a conditional connection, this is one more than the
342        ///  operator index of the ending of the atomic call sequence; i.e.,
343        /// the second UserOp.
344        size_t op_end;
345};
346
347/*!
348Shared documentation for optimization helper functions (not called).
349
350<!-- define prototype -->
351\param tape
352is a vector that maps a variable index, in the old operation sequence,
353to an <tt>struct_old_variable</tt> information record.
354Note that the index for this vector must be greater than or equal zero and
355less than <tt>tape.size()</tt>.
356
357\li <tt>tape[i].op</tt>
358is the operator in the old operation sequence
359corresponding to the old variable index \c i.
360Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
361
362\li <tt>tape[i].arg</tt>
363for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
364is the j-th the argument, in the old operation sequence,
365corresponding to the old variable index \c i.
366Assertion: <tt>tape[i].arg[j] < i</tt>.
367
368\li <tt>tape[i].new_var</tt>
369Suppose
370<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
371and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
372It follows that <tt>tape[k].new_var</tt>
373has alread been set to the variable in the new operation sequence
374corresponding to the old variable index \c k.
375This means that the \c new_var value has been set
376for all the possible arguments that come before \a current.
377
378\param current
379is the index in the old operation sequence for
380the variable corresponding to the result for the current operator.
381Assertions:
382<tt>
383current < tape.size(),
384NumRes( tape[current].op ) > 0.
385</tt>
386
387\param npar
388is the number of parameters corresponding to this operation sequence.
389
390\param par
391is a vector of length \a npar containing the parameters
392for this operation sequence; i.e.,
393given a parameter index \c i, the corresponding parameter value is
394<tt>par[i]</tt>.
395<!-- end prototype -->
396*/
397template <class Base>
398void prototype(
399        const CppAD::vector<struct struct_old_variable>& tape           ,
400        size_t                                             current        ,
401        size_t                                             npar           ,
402        const Base*                                        par            )
403{       CPPAD_ASSERT_UNKNOWN(false); }
404
405/*!
406Check a unary operator for a complete match with a previous operator.
407
408A complete match means that the result of the previous operator
409can be used inplace of the result for current operator.
410
411<!-- replace prototype -->
412\param tape
413is a vector that maps a variable index, in the old operation sequence,
414to an <tt>struct_old_variable</tt> information record.
415Note that the index for this vector must be greater than or equal zero and
416less than <tt>tape.size()</tt>.
417
418\li <tt>tape[i].op</tt>
419is the operator in the old operation sequence
420corresponding to the old variable index \c i.
421Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
422
423\li <tt>tape[i].arg</tt>
424for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
425is the j-th the argument, in the old operation sequence,
426corresponding to the old variable index \c i.
427Assertion: <tt>tape[i].arg[j] < i</tt>.
428
429\li <tt>tape[i].new_var</tt>
430Suppose
431<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
432and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
433It follows that <tt>tape[k].new_var</tt>
434has alread been set to the variable in the new operation sequence
435corresponding to the old variable index \c k.
436This means that the \c new_var value has been set
437for all the possible arguments that come before \a current.
438
439\param current
440is the index in the old operation sequence for
441the variable corresponding to the result for the current operator.
442Assertions:
443<tt>
444current < tape.size(),
445NumRes( tape[current].op ) > 0.
446</tt>
447
448\param npar
449is the number of parameters corresponding to this operation sequence.
450
451\param par
452is a vector of length \a npar containing the parameters
453for this operation sequence; i.e.,
454given a parameter index \c i, the corresponding parameter value is
455<tt>par[i]</tt>.
456<!-- end prototype -->
457
458\param hash_table_var
459is a vector with size CPPAD_HASH_TABLE_SIZE
460that maps a hash code to the corresponding
461variable index in the old operation sequence.
462All the values in this table must be less than \a current.
463
464\param code
465The input value of code does not matter.
466The output value of code is the hash code corresponding to
467this operation in the new operation sequence.
468
469\return
470If the return value is zero,
471no match was found.
472If the return value is greater than zero,
473it is the index of a new variable that can be used to replace the
474old variable.
475
476\par Restrictions:
477NumArg( tape[current].op ) == 1
478*/
479template <class Base>
480size_t unary_match(
481        const CppAD::vector<struct struct_old_variable>& tape           ,
482        size_t                                             current        ,
483        size_t                                             npar           ,
484        const Base*                                        par            ,
485        const CppAD::vector<size_t>&                       hash_table_var ,
486        unsigned short&                                    code           )
487{       const addr_t* arg = tape[current].arg;
488        OpCode        op  = tape[current].op;
489        addr_t new_arg[1];
490       
491        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
492        CPPAD_ASSERT_UNKNOWN( NumRes(op) > 0  );
493        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
494        new_arg[0] = tape[arg[0]].new_var;
495        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < current );
496        code = hash_code(
497                op                  , 
498                new_arg             ,
499                npar                ,
500                par
501        );
502        size_t  i               = hash_table_var[code];
503        CPPAD_ASSERT_UNKNOWN( i < current );
504        if( op == tape[i].op )
505        {       size_t k = tape[i].arg[0];
506                CPPAD_ASSERT_UNKNOWN( k < i );
507                if (new_arg[0] == tape[k].new_var )
508                        return tape[i].new_var;
509        }
510        return 0;
511} 
512
513/*!
514Check a binary operator for a complete match with a previous operator,
515
516<!-- replace prototype -->
517\param tape
518is a vector that maps a variable index, in the old operation sequence,
519to an <tt>struct_old_variable</tt> information record.
520Note that the index for this vector must be greater than or equal zero and
521less than <tt>tape.size()</tt>.
522
523\li <tt>tape[i].op</tt>
524is the operator in the old operation sequence
525corresponding to the old variable index \c i.
526Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
527
528\li <tt>tape[i].arg</tt>
529for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
530is the j-th the argument, in the old operation sequence,
531corresponding to the old variable index \c i.
532Assertion: <tt>tape[i].arg[j] < i</tt>.
533
534\li <tt>tape[i].new_var</tt>
535Suppose
536<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
537and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
538It follows that <tt>tape[k].new_var</tt>
539has alread been set to the variable in the new operation sequence
540corresponding to the old variable index \c k.
541This means that the \c new_var value has been set
542for all the possible arguments that come before \a current.
543
544\param current
545is the index in the old operation sequence for
546the variable corresponding to the result for the current operator.
547Assertions:
548<tt>
549current < tape.size(),
550NumRes( tape[current].op ) > 0.
551</tt>
552
553\param npar
554is the number of parameters corresponding to this operation sequence.
555
556\param par
557is a vector of length \a npar containing the parameters
558for this operation sequence; i.e.,
559given a parameter index \c i, the corresponding parameter value is
560<tt>par[i]</tt>.
561<!-- end prototype -->
562
563\param hash_table_var
564is a vector with size CPPAD_HASH_TABLE_SIZE
565that maps a hash code to the corresponding
566variable index in the old operation sequence.
567All the values in this table must be less than \a current.
568
569\param code
570The input value of code does not matter.
571The output value of code is the hash code corresponding to
572this operation in the new operation sequence.
573
574\return
575If the return value is zero,
576no match was found.
577If the return value is greater than zero,
578it is the index of a new variable that can be used to replace the
579old variable.
580
581
582\par Restrictions:
583The binary operator must be an addition, subtraction, multiplication, division
584or power operator.  NumArg( tape[current].op ) == 1.
585*/
586template <class Base>
587inline size_t binary_match(
588        const CppAD::vector<struct struct_old_variable>&   tape           ,
589        size_t                                             current        ,
590        size_t                                             npar           ,
591        const Base*                                        par            ,
592        const CppAD::vector<size_t>&                       hash_table_var ,
593        unsigned short&                                    code           )
594{       OpCode        op         = tape[current].op;
595        const addr_t* arg        = tape[current].arg;
596        addr_t        new_arg[2];
597        bool          parameter[2];
598
599        // initialize return value
600        size_t  match_var = 0;
601
602        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
603        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
604        switch(op)
605        {       // index op variable
606                case DisOp:
607                // parameter not defined for this case
608                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
609                new_arg[0]   = arg[0];
610                new_arg[1]   = tape[arg[1]].new_var;
611                break;
612
613                // parameter op variable ----------------------------------
614                case AddpvOp:
615                case MulpvOp:
616                case DivpvOp:
617                case PowpvOp:
618                case SubpvOp:
619                // arg[0]
620                parameter[0] = true;
621                new_arg[0]   = arg[0];
622                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar );
623                // arg[1]
624                parameter[1] = false;
625                new_arg[1]   = tape[arg[1]].new_var;
626                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
627                break;
628
629                // variable op parameter -----------------------------------
630                case DivvpOp:
631                case PowvpOp:
632                case SubvpOp:
633                // arg[0]
634                parameter[0] = false;
635                new_arg[0]   = tape[arg[0]].new_var;
636                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
637                // arg[1]
638                parameter[1] = true;
639                new_arg[1]   = arg[1];
640                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar );
641                break;
642
643                // variable op variable -----------------------------------
644                case AddvvOp:
645                case MulvvOp:
646                case DivvvOp:
647                case PowvvOp:
648                case SubvvOp:
649                // arg[0]
650                parameter[0] = false;
651                new_arg[0]   = tape[arg[0]].new_var;
652                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
653                // arg[1]
654                parameter[1] = false;
655                new_arg[1]   = tape[arg[1]].new_var;
656                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
657                break;
658
659                // must be one of the cases above
660                default:
661                CPPAD_ASSERT_UNKNOWN(false);
662        }
663        code = hash_code(
664                op                  , 
665                new_arg             ,
666                npar                ,
667                par
668        );
669        size_t  i  = hash_table_var[code];
670        CPPAD_ASSERT_UNKNOWN( i < current );
671        if( op == tape[i].op )
672        {       bool match = true;
673                if( op == DisOp )
674                {       match   &= new_arg[0] == tape[i].arg[0];
675                        size_t k = tape[i].arg[1];
676                        match   &= new_arg[1] == tape[k].new_var;
677                }
678                else
679                {       for(size_t j = 0; j < 2; j++)
680                        {       size_t k = tape[i].arg[j];
681                                if( parameter[j] )
682                                {       CPPAD_ASSERT_UNKNOWN( k < npar );
683                                        match &= IdenticalEqualPar(
684                                                par[ arg[j] ], par[k]
685                                        );
686                                }
687                                else
688                                {       CPPAD_ASSERT_UNKNOWN( k < i );
689                                        match &= (new_arg[j] == tape[k].new_var);
690                                }
691                        }
692                }
693                if( match )
694                        match_var = tape[i].new_var;
695        }
696        if( (match_var > 0) | ( (op != AddvvOp) & (op != MulvvOp ) ) )
697                return match_var;
698
699        // check for match with argument order switched ----------------------
700        CPPAD_ASSERT_UNKNOWN( op == AddvvOp || op == MulvvOp );
701        i          = new_arg[0];
702        new_arg[0] = new_arg[1];
703        new_arg[1] = i;
704        unsigned short code_switch = hash_code(
705                op                  , 
706                new_arg             ,
707                npar                ,
708                par
709        );
710        i  = hash_table_var[code_switch];
711        CPPAD_ASSERT_UNKNOWN( i < current );
712        if( op == tape[i].op )
713        {       bool match = true;
714                size_t j;
715                for(j = 0; j < 2; j++)
716                {       size_t k = tape[i].arg[j];
717                        CPPAD_ASSERT_UNKNOWN( k < i );
718                        match &= (new_arg[j] == tape[k].new_var);
719                }
720                if( match )
721                        match_var = tape[i].new_var;
722        }
723        return match_var;
724} 
725
726/*!
727Record an operation of the form (parameter op variable).
728
729<!-- replace prototype -->
730\param tape
731is a vector that maps a variable index, in the old operation sequence,
732to an <tt>struct_old_variable</tt> information record.
733Note that the index for this vector must be greater than or equal zero and
734less than <tt>tape.size()</tt>.
735
736\li <tt>tape[i].op</tt>
737is the operator in the old operation sequence
738corresponding to the old variable index \c i.
739Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
740
741\li <tt>tape[i].arg</tt>
742for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
743is the j-th the argument, in the old operation sequence,
744corresponding to the old variable index \c i.
745Assertion: <tt>tape[i].arg[j] < i</tt>.
746
747\li <tt>tape[i].new_var</tt>
748Suppose
749<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
750and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
751It follows that <tt>tape[k].new_var</tt>
752has alread been set to the variable in the new operation sequence
753corresponding to the old variable index \c k.
754This means that the \c new_var value has been set
755for all the possible arguments that come before \a current.
756
757\param current
758is the index in the old operation sequence for
759the variable corresponding to the result for the current operator.
760Assertions:
761<tt>
762current < tape.size(),
763NumRes( tape[current].op ) > 0.
764</tt>
765
766\param npar
767is the number of parameters corresponding to this operation sequence.
768
769\param par
770is a vector of length \a npar containing the parameters
771for this operation sequence; i.e.,
772given a parameter index \c i, the corresponding parameter value is
773<tt>par[i]</tt>.
774<!-- end prototype -->
775
776\param rec
777is the object that will record the operations.
778
779\param op
780is the operator that we are recording which must be one of the following:
781AddpvOp, DivpvOp, MulpvOp, PowvpOp, SubpvOp.
782 
783\param arg
784is the vector of arguments for this operator.
785
786\return
787the result is the operaiton and variable index corresponding to the current
788operation in the new operation sequence.
789*/
790template <class Base>
791struct_size_pair record_pv(
792        const CppAD::vector<struct struct_old_variable>& tape           ,
793        size_t                                             current        ,
794        size_t                                             npar           ,
795        const Base*                                        par            ,
796        recorder<Base>*                                    rec            ,
797        OpCode                                             op             ,
798        const addr_t*                                      arg            )
799{
800# ifndef NDEBUG
801        switch(op)
802        {       case AddpvOp:
803                case DivpvOp:
804                case MulpvOp:
805                case PowpvOp:
806                case SubpvOp:
807                break;
808
809                default:
810                CPPAD_ASSERT_UNKNOWN(false);
811        }
812# endif
813        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar    );
814        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
815        addr_t new_arg[2];
816        new_arg[0]   = rec->PutPar( par[arg[0]] );
817        new_arg[1]   = tape[ arg[1] ].new_var;
818        rec->PutArg( new_arg[0], new_arg[1] );
819
820        struct_size_pair ret;
821        ret.i_op  = rec->num_op_rec();
822        ret.i_var = rec->PutOp(op);
823        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < ret.i_var );
824        return ret;
825}
826
827
828/*!
829Record an operation of the form (variable op parameter).
830
831<!-- replace prototype -->
832\param tape
833is a vector that maps a variable index, in the old operation sequence,
834to an <tt>struct_old_variable</tt> information record.
835Note that the index for this vector must be greater than or equal zero and
836less than <tt>tape.size()</tt>.
837
838\li <tt>tape[i].op</tt>
839is the operator in the old operation sequence
840corresponding to the old variable index \c i.
841Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
842
843\li <tt>tape[i].arg</tt>
844for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
845is the j-th the argument, in the old operation sequence,
846corresponding to the old variable index \c i.
847Assertion: <tt>tape[i].arg[j] < i</tt>.
848
849\li <tt>tape[i].new_var</tt>
850Suppose
851<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
852and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
853It follows that <tt>tape[k].new_var</tt>
854has alread been set to the variable in the new operation sequence
855corresponding to the old variable index \c k.
856This means that the \c new_var value has been set
857for all the possible arguments that come before \a current.
858
859\param current
860is the index in the old operation sequence for
861the variable corresponding to the result for the current operator.
862Assertions:
863<tt>
864current < tape.size(),
865NumRes( tape[current].op ) > 0.
866</tt>
867
868\param npar
869is the number of parameters corresponding to this operation sequence.
870
871\param par
872is a vector of length \a npar containing the parameters
873for this operation sequence; i.e.,
874given a parameter index \c i, the corresponding parameter value is
875<tt>par[i]</tt>.
876<!-- end prototype -->
877
878\param rec
879is the object that will record the operations.
880
881\param op
882is the operator that we are recording which must be one of the following:
883DivvpOp, PowvpOp, SubvpOp.
884 
885\param arg
886is the vector of arguments for this operator.
887
888\return
889the result operation and variable index corresponding to the current
890operation in the new operation sequence.
891*/
892template <class Base>
893struct_size_pair record_vp(
894        const CppAD::vector<struct struct_old_variable>& tape           ,
895        size_t                                             current        ,
896        size_t                                             npar           ,
897        const Base*                                        par            ,
898        recorder<Base>*                                    rec            ,
899        OpCode                                             op             ,
900        const addr_t*                                      arg            )
901{
902# ifndef NDEBUG
903        switch(op)
904        {       case DivvpOp:
905                case PowvpOp:
906                case SubvpOp:
907                break;
908
909                default:
910                CPPAD_ASSERT_UNKNOWN(false);
911        }
912# endif
913        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
914        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar    );
915        addr_t new_arg[2];
916        new_arg[0]   = tape[ arg[0] ].new_var;
917        new_arg[1]   = rec->PutPar( par[arg[1]] );
918        rec->PutArg( new_arg[0], new_arg[1] );
919
920        struct_size_pair ret;
921        ret.i_op  = rec->num_op_rec();
922        ret.i_var = rec->PutOp(op);
923        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < ret.i_var );
924        return ret;
925}
926
927/*!
928Record an operation of the form (variable op variable).
929
930<!-- replace prototype -->
931\param tape
932is a vector that maps a variable index, in the old operation sequence,
933to an <tt>struct_old_variable</tt> information record.
934Note that the index for this vector must be greater than or equal zero and
935less than <tt>tape.size()</tt>.
936
937\li <tt>tape[i].op</tt>
938is the operator in the old operation sequence
939corresponding to the old variable index \c i.
940Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
941
942\li <tt>tape[i].arg</tt>
943for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
944is the j-th the argument, in the old operation sequence,
945corresponding to the old variable index \c i.
946Assertion: <tt>tape[i].arg[j] < i</tt>.
947
948\li <tt>tape[i].new_var</tt>
949Suppose
950<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
951and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
952It follows that <tt>tape[k].new_var</tt>
953has alread been set to the variable in the new operation sequence
954corresponding to the old variable index \c k.
955This means that the \c new_var value has been set
956for all the possible arguments that come before \a current.
957
958\param current
959is the index in the old operation sequence for
960the variable corresponding to the result for the current operator.
961Assertions:
962<tt>
963current < tape.size(),
964NumRes( tape[current].op ) > 0.
965</tt>
966
967\param npar
968is the number of parameters corresponding to this operation sequence.
969
970\param par
971is a vector of length \a npar containing the parameters
972for this operation sequence; i.e.,
973given a parameter index \c i, the corresponding parameter value is
974<tt>par[i]</tt>.
975<!-- end prototype -->
976
977\param rec
978is the object that will record the operations.
979
980\param op
981is the operator that we are recording which must be one of the following:
982AddvvOp, DivvvOp, MulvvOp, PowvpOp, SubvvOp.
983 
984\param arg
985is the vector of arguments for this operator.
986
987\return
988the result is the operation and variable index corresponding to the current
989operation in the new operation sequence.
990*/
991template <class Base>
992struct_size_pair record_vv(
993        const CppAD::vector<struct struct_old_variable>& tape           ,
994        size_t                                             current        ,
995        size_t                                             npar           ,
996        const Base*                                        par            ,
997        recorder<Base>*                                    rec            ,
998        OpCode                                             op             ,
999        const addr_t*                                      arg            )
1000{
1001# ifndef NDEBUG
1002        switch(op)
1003        {       case AddvvOp:
1004                case DivvvOp:
1005                case MulvvOp:
1006                case PowvvOp:
1007                case SubvvOp:
1008                break;
1009
1010                default:
1011                CPPAD_ASSERT_UNKNOWN(false);
1012        }
1013# endif
1014        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
1015        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
1016        addr_t new_arg[2];
1017        new_arg[0]   = tape[ arg[0] ].new_var;
1018        new_arg[1]   = tape[ arg[1] ].new_var;
1019        rec->PutArg( new_arg[0], new_arg[1] );
1020
1021        struct_size_pair ret;
1022        ret.i_op  = rec->num_op_rec();
1023        ret.i_var = rec->PutOp(op);
1024        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < ret.i_var );
1025        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < ret.i_var );
1026        return ret;
1027}
1028
1029// ==========================================================================
1030
1031/*!
1032Recording a cummulative cummulative summation starting at its highest parrent.
1033
1034<!-- replace prototype -->
1035\param tape
1036is a vector that maps a variable index, in the old operation sequence,
1037to an <tt>struct_old_variable</tt> information record.
1038Note that the index for this vector must be greater than or equal zero and
1039less than <tt>tape.size()</tt>.
1040
1041\li <tt>tape[i].op</tt>
1042is the operator in the old operation sequence
1043corresponding to the old variable index \c i.
1044Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
1045
1046\li <tt>tape[i].arg</tt>
1047for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
1048is the j-th the argument, in the old operation sequence,
1049corresponding to the old variable index \c i.
1050Assertion: <tt>tape[i].arg[j] < i</tt>.
1051
1052\li <tt>tape[i].new_var</tt>
1053Suppose
1054<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
1055and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
1056It follows that <tt>tape[k].new_var</tt>
1057has alread been set to the variable in the new operation sequence
1058corresponding to the old variable index \c k.
1059This means that the \c new_var value has been set
1060for all the possible arguments that come before \a current.
1061
1062\param current
1063is the index in the old operation sequence for
1064the variable corresponding to the result for the current operator.
1065Assertions:
1066<tt>
1067current < tape.size(),
1068NumRes( tape[current].op ) > 0.
1069</tt>
1070
1071\param npar
1072is the number of parameters corresponding to this operation sequence.
1073
1074\param par
1075is a vector of length \a npar containing the parameters
1076for this operation sequence; i.e.,
1077given a parameter index \c i, the corresponding parameter value is
1078<tt>par[i]</tt>.
1079<!-- end prototype -->
1080
1081\param rec
1082is the object that will record the operations.
1083
1084\param work
1085Is used for computation. On input and output,
1086<tt>work.op_stack.empty()</tt>,
1087<tt>work.add_stack.empty()</tt>, and
1088<tt>work.sub_stack.empty()</tt>,
1089are all true true.
1090These stacks are passed in so that elements can be allocated once
1091and then the elements can be reused with calls to \c record_csum.
1092
1093\par Exception
1094<tt>tape[i].new_var</tt>
1095is not yet defined for any node \c i that is \c csum_connected
1096to the \a current node
1097(or that is \c sum_connected to a node that is \c csum_connected).
1098For example; suppose that index \c j corresponds to a variable
1099in the current operator,
1100<tt>i = tape[current].arg[j]</tt>,
1101and
1102<tt>tape[arg[j]].connect_type == csum_connected</tt>.
1103It then follows that
1104<tt>tape[i].new_var == tape.size()</tt>.
1105
1106\par Restriction:
1107\li <tt>tape[current].op</tt>
1108must be one of <tt>AddpvOp, AddvvOp, SubpvOp, SubvpOp, SubvvOp</tt>.
1109
1110\li <tt>tape[current].connect_type</tt> must be \c yes_connected.
1111
1112\li <tt>tape[j].connect_type == csum_connected</tt> for some index
1113j that is a variable operand for the current operation.
1114*/
1115
1116
1117template <class Base>
1118struct_size_pair record_csum(
1119        const CppAD::vector<struct struct_old_variable>& tape           ,
1120        size_t                                             current        ,
1121        size_t                                             npar           ,
1122        const Base*                                        par            ,
1123        recorder<Base>*                                    rec            ,
1124        struct_csum_stacks&                              work           )
1125{
1126       
1127        CPPAD_ASSERT_UNKNOWN( work.op_stack.empty() );
1128        CPPAD_ASSERT_UNKNOWN( work.add_stack.empty() );
1129        CPPAD_ASSERT_UNKNOWN( work.sub_stack.empty() );
1130        CPPAD_ASSERT_UNKNOWN( tape[current].connect_type == yes_connected );
1131
1132        size_t                        i;
1133        OpCode                        op;
1134        const addr_t*                 arg;
1135        bool                          add;
1136        struct struct_csum_variable var;
1137
1138        var.op  = tape[current].op;
1139        var.arg = tape[current].arg;
1140        var.add = true; 
1141        work.op_stack.push( var );
1142        Base sum_par(0);
1143
1144# ifndef NDEBUG
1145        bool ok = false;
1146        if( var.op == SubvpOp )
1147                ok = tape[ tape[current].arg[0] ].connect_type == csum_connected;
1148        if( var.op == AddpvOp || var.op == SubpvOp )
1149                ok = tape[ tape[current].arg[1] ].connect_type == csum_connected;
1150        if( var.op == AddvvOp || var.op == SubvvOp )
1151        {       ok  = tape[ tape[current].arg[0] ].connect_type == csum_connected;
1152                ok |= tape[ tape[current].arg[1] ].connect_type == csum_connected;
1153        }
1154        CPPAD_ASSERT_UNKNOWN( ok );
1155# endif
1156        while( ! work.op_stack.empty() )
1157        {       var     = work.op_stack.top();
1158                work.op_stack.pop();
1159                op      = var.op;
1160                arg     = var.arg;
1161                add     = var.add;
1162                // process first argument to this operator
1163                switch(op)
1164                {       case AddpvOp:
1165                        case SubpvOp:
1166                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar );
1167                        if( add )
1168                                sum_par += par[arg[0]];
1169                        else    sum_par -= par[arg[0]];
1170                        break;
1171
1172                        case AddvvOp:
1173                        case SubvpOp:
1174                        case SubvvOp:
1175                        if( tape[arg[0]].connect_type == csum_connected )
1176                        {       CPPAD_ASSERT_UNKNOWN(
1177                                        size_t(tape[arg[0]].new_var) == tape.size()
1178                                );
1179                                var.op  = tape[arg[0]].op;
1180                                var.arg = tape[arg[0]].arg;
1181                                var.add = add; 
1182                                work.op_stack.push( var );
1183                        }
1184                        else if( add )
1185                                work.add_stack.push(arg[0]);
1186                        else    work.sub_stack.push(arg[0]);
1187                        break;
1188
1189                        default:
1190                        CPPAD_ASSERT_UNKNOWN(false);
1191                }
1192                // process second argument to this operator
1193                switch(op)
1194                {
1195                        case SubvpOp:
1196                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar );
1197                        if( add )
1198                                sum_par -= par[arg[1]];
1199                        else    sum_par += par[arg[1]];
1200                        break;
1201
1202                        case SubvvOp:
1203                        case SubpvOp:
1204                        add = ! add;
1205
1206                        case AddvvOp:
1207                        case AddpvOp:
1208                        if( tape[arg[1]].connect_type == csum_connected )
1209                        {       CPPAD_ASSERT_UNKNOWN(
1210                                        size_t(tape[arg[1]].new_var) == tape.size()
1211                                );
1212                                var.op   = tape[arg[1]].op;
1213                                var.arg  = tape[arg[1]].arg;
1214                                var.add  = add;
1215                                work.op_stack.push( var );
1216                        }
1217                        else if( add )
1218                                work.add_stack.push(arg[1]);
1219                        else    work.sub_stack.push(arg[1]);
1220                        break;
1221
1222                        default:
1223                        CPPAD_ASSERT_UNKNOWN(false);
1224                }
1225        }
1226        // number of variables in this cummulative sum operator
1227        size_t n_add = work.add_stack.size();
1228        size_t n_sub = work.sub_stack.size();
1229        size_t old_arg, new_arg;
1230        rec->PutArg(n_add);                // arg[0]
1231        rec->PutArg(n_sub);                // arg[1]
1232        new_arg = rec->PutPar( sum_par );
1233        rec->PutArg(new_arg);              // arg[2]
1234        for(i = 0; i < n_add; i++)
1235        {       CPPAD_ASSERT_UNKNOWN( ! work.add_stack.empty() );
1236                old_arg = work.add_stack.top();
1237                new_arg = tape[old_arg].new_var;
1238                CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
1239                rec->PutArg(new_arg);      // arg[3+i]
1240                work.add_stack.pop();
1241        }
1242        for(i = 0; i < n_sub; i++)
1243        {       CPPAD_ASSERT_UNKNOWN( ! work.sub_stack.empty() );
1244                old_arg = work.sub_stack.top();
1245                new_arg = tape[old_arg].new_var;
1246                CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
1247                rec->PutArg(new_arg);      // arg[3 + arg[0] + i]
1248                work.sub_stack.pop();
1249        }
1250        rec->PutArg(n_add + n_sub);        // arg[3 + arg[0] + arg[1]]
1251
1252
1253        struct_size_pair ret;
1254        ret.i_op  = rec->num_op_rec();
1255        ret.i_var = rec->PutOp(CSumOp);
1256        CPPAD_ASSERT_UNKNOWN( new_arg < ret.i_var );
1257        return ret;
1258}
1259// ==========================================================================
1260/*!
1261Convert a player object to an optimized recorder object
1262
1263\tparam Base
1264base type for the operator; i.e., this operation was recorded
1265using AD< \a Base > and computations by this routine are done using type
1266\a Base.
1267
1268\param n
1269is the number of independent variables on the tape.
1270
1271\param dep_taddr
1272On input this vector contains the indices for each of the dependent
1273variable values in the operation sequence corresponding to \a play.
1274Upon return it contains the indices for the same variables but in
1275the operation sequence corresponding to \a rec.
1276
1277\param play
1278This is the operation sequence that we are optimizing.
1279It is essentially const, except for play back state which
1280changes while it plays back the operation seqeunce.
1281
1282\param rec
1283The input contents of this recording does not matter.
1284Upon return, it contains an optimized verison of the
1285operation sequence corresponding to \a play.
1286*/
1287
1288template <class Base>
1289void optimize_run(
1290        size_t                       n         ,
1291        CppAD::vector<size_t>&       dep_taddr ,
1292        player<Base>*                play      ,
1293        recorder<Base>*              rec       ) 
1294{
1295        // temporary indices
1296        size_t i, j, k;
1297
1298        // temporary variables
1299        OpCode        op;   // current operator
1300        const addr_t* arg;  // operator arguments
1301        size_t        i_var;  // index of first result for current operator
1302
1303        // range and domain dimensions for F
1304        size_t m = dep_taddr.size();
1305
1306        // number of variables in the player
1307        const size_t num_var = play->num_var_rec(); 
1308
1309# ifndef NDEBUG
1310        // number of parameters in the player
1311        const size_t num_par = play->num_par_rec();
1312# endif
1313
1314        // number of  VecAD indices
1315        size_t num_vecad_ind   = play->num_vec_ind_rec();
1316
1317        // number of VecAD vectors
1318        size_t num_vecad_vec   = play->num_vecad_vec_rec();
1319
1320        // -------------------------------------------------------------
1321        // data structure that maps variable index in original operation
1322        // sequence to corresponding operator information
1323        CppAD::vector<struct struct_old_variable> tape(num_var);
1324        // -------------------------------------------------------------
1325        // Determine how each variable is connected to the dependent variables
1326
1327        // initialize all variables has having no connections
1328        for(i = 0; i < num_var; i++)
1329                tape[i].connect_type = not_connected;
1330
1331        for(j = 0; j < m; j++)
1332        {       // mark dependent variables as having one or more connections
1333                tape[ dep_taddr[j] ].connect_type = yes_connected;
1334        }
1335
1336        // vecad_connect contains a value for each VecAD object.
1337        // vecad maps a VecAD index (which corresponds to the beginning of the
1338        // VecAD object) to the vecad_connect falg for the VecAD object.
1339        CppAD::vector<enum_connect_type>   vecad_connect(num_vecad_vec);
1340        CppAD::vector<size_t> vecad(num_vecad_ind);
1341        j = 0;
1342        for(i = 0; i < num_vecad_vec; i++)
1343        {       vecad_connect[i] = not_connected;
1344                // length of this VecAD
1345                size_t length = play->GetVecInd(j);
1346                // set to proper index for this VecAD
1347                vecad[j] = i; 
1348                for(k = 1; k <= length; k++)
1349                        vecad[j+k] = num_vecad_vec; // invalid index
1350                // start of next VecAD
1351                j       += length + 1;
1352        }
1353        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
1354
1355        // work space used by UserOp.
1356        typedef std::set<size_t> size_set;
1357        vector<size_set> user_r_set;   // set sparsity pattern for result
1358        vector<size_set> user_s_set;   // set sparisty pattern for argument
1359        vector<bool>     user_r_bool;  // bool sparsity pattern for result
1360        vector<bool>     user_s_bool;  // bool sparisty pattern for argument
1361        //
1362        size_t user_q     = 0;       // column dimension for sparsity patterns
1363        size_t user_index = 0;       // indentifier for this user_atomic operation
1364        size_t user_id    = 0;       // user identifier for this call to operator
1365        size_t user_i     = 0;       // index in result vector
1366        size_t user_j     = 0;       // index in argument vector
1367        size_t user_m     = 0;       // size of result vector
1368        size_t user_n     = 0;       // size of arugment vector
1369        //
1370        atomic_base<Base>* user_atom = CPPAD_NULL; // current user atomic function
1371        bool               user_set  = true;       // use set sparsity (or bool)
1372
1373        // next expected operator in a UserOp sequence
1374        enum { user_start, user_arg, user_ret, user_end } user_state;
1375
1376        // During reverse mode, compute type of connection for each call to
1377        // a user atomic function.
1378        CppAD::vector<struct_user_info>    user_info;
1379        size_t                             user_curr = 0;
1380
1381        /// During reverse mode, information for each CSkip operation
1382        CppAD::vector<struct_cskip_info>   cskip_info;
1383
1384        // Initialize a reverse mode sweep through the operation sequence
1385        size_t i_op;
1386        play->reverse_start(op, arg, i_op, i_var);
1387        CPPAD_ASSERT_UNKNOWN( op == EndOp );
1388        size_t mask;
1389        user_state = user_end;
1390        while(op != BeginOp)
1391        {       // next op
1392                play->reverse_next(op, arg, i_op, i_var);
1393
1394                // Store the operator corresponding to each variable
1395                if( NumRes(op) > 0 )
1396                {       tape[i_var].op = op;
1397                        tape[i_var].arg = arg;
1398                }
1399# ifndef NDEBUG
1400                if( i_op <= n )
1401                {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
1402                }
1403                else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
1404# endif
1405                enum_connect_type connect_type      = tape[i_var].connect_type;
1406                std::set<class_cexp_pair>& cexp_set = tape[i_var].cexp_set;
1407                switch( op )
1408                {
1409                        // One variable corresponding to arg[0]
1410                        case AbsOp:
1411                        case AcosOp:
1412                        case AsinOp:
1413                        case AtanOp:
1414                        case CosOp:
1415                        case CoshOp:
1416                        case DivvpOp:
1417                        case ExpOp:
1418                        case LogOp:
1419                        case PowvpOp:
1420                        case SignOp:
1421                        case SinOp:
1422                        case SinhOp:
1423                        case SqrtOp:
1424                        case TanOp:
1425                        case TanhOp:
1426                        switch( connect_type )
1427                        {       case not_connected:
1428                                break;
1429       
1430                                case yes_connected:
1431                                case sum_connected:
1432                                case csum_connected: 
1433                                tape[arg[0]].connect_type = yes_connected;
1434                                break;
1435
1436                                case cexp_connected:
1437                                if( tape[arg[0]].connect_type == not_connected )
1438                                {       tape[arg[0]].connect_type = cexp_connected;
1439                                        tape[arg[0]].cexp_set     = cexp_set;
1440                                }
1441                                else if( tape[arg[0]].connect_type == cexp_connected )
1442                                {       tape[arg[0]].cexp_set = intersection(
1443                                                tape[arg[0]].cexp_set, cexp_set
1444                                        );
1445                                        if( tape[arg[0]].cexp_set.empty() )
1446                                                tape[arg[0]].connect_type = yes_connected;
1447                                }
1448                                else    tape[arg[0]].connect_type = yes_connected;
1449                                break;
1450
1451                                default:
1452                                CPPAD_ASSERT_UNKNOWN(false);
1453                        }
1454                        break; // --------------------------------------------
1455
1456                        // One variable corresponding to arg[1]
1457                        case DisOp:
1458                        case DivpvOp:
1459                        case MulpvOp:
1460                        case PowpvOp:
1461                        switch( connect_type )
1462                        {       case not_connected:
1463                                break;
1464       
1465                                case yes_connected:
1466                                case sum_connected:
1467                                case csum_connected: 
1468                                tape[arg[1]].connect_type = yes_connected;
1469                                break;
1470
1471                                case cexp_connected:
1472                                if( tape[arg[1]].connect_type == not_connected )
1473                                {       tape[arg[1]].connect_type = cexp_connected;
1474                                        tape[arg[1]].cexp_set     = cexp_set;
1475                                }
1476                                else if( tape[arg[1]].connect_type == cexp_connected )
1477                                {       tape[arg[1]].cexp_set = intersection(
1478                                                tape[arg[1]].cexp_set, cexp_set
1479                                        );
1480                                        if( tape[arg[1]].cexp_set.empty() )
1481                                                tape[arg[1]].connect_type = yes_connected;
1482                                }
1483                                else    tape[arg[1]].connect_type = yes_connected;
1484                                break;
1485
1486                                default:
1487                                CPPAD_ASSERT_UNKNOWN(false);
1488                        }
1489                        break; // --------------------------------------------
1490               
1491                        // Special case for SubvpOp
1492                        case SubvpOp:
1493                        switch( connect_type )
1494                        {       case not_connected:
1495                                break;
1496       
1497                                case yes_connected:
1498                                case sum_connected:
1499                                case csum_connected: 
1500                                if( tape[arg[0]].connect_type == not_connected )
1501                                        tape[arg[0]].connect_type = sum_connected;
1502                                else    tape[arg[0]].connect_type = yes_connected;
1503                                break;
1504
1505                                case cexp_connected:
1506                                if( tape[arg[0]].connect_type == not_connected )
1507                                {       tape[arg[0]].connect_type = cexp_connected;
1508                                        tape[arg[0]].cexp_set     = cexp_set;
1509                                }
1510                                else if( tape[arg[0]].connect_type == cexp_connected )
1511                                {       tape[arg[0]].cexp_set = intersection(
1512                                                tape[arg[0]].cexp_set, cexp_set
1513                                        );
1514                                        if( tape[arg[0]].cexp_set.empty() )
1515                                                tape[arg[0]].connect_type = yes_connected;
1516                                }
1517                                else    tape[arg[0]].connect_type = yes_connected;
1518                                break;
1519
1520                                default:
1521                                CPPAD_ASSERT_UNKNOWN(false);
1522                        }
1523                        if( connect_type == sum_connected )
1524                        {       // convert sum to csum connection for this variable
1525                                tape[i_var].connect_type = connect_type = csum_connected;
1526                        }
1527                        break; // --------------------------------------------
1528               
1529                        // Special case for AddpvOp and SubpvOp
1530                        case AddpvOp:
1531                        case SubpvOp:
1532                        switch( connect_type )
1533                        {       case not_connected:
1534                                break;
1535       
1536                                case yes_connected:
1537                                case sum_connected:
1538                                case csum_connected:
1539                                if( tape[arg[1]].connect_type == not_connected )
1540                                        tape[arg[1]].connect_type = sum_connected;
1541                                else    tape[arg[1]].connect_type = yes_connected;
1542                                break;
1543
1544                                case cexp_connected:
1545                                if( tape[arg[1]].connect_type == not_connected )
1546                                {       tape[arg[1]].connect_type = cexp_connected;
1547                                        tape[arg[1]].cexp_set     = cexp_set;
1548                                }
1549                                else if( tape[arg[1]].connect_type == cexp_connected )
1550                                {       tape[arg[1]].cexp_set = intersection(
1551                                                tape[arg[1]].cexp_set, cexp_set
1552                                        );
1553                                        if( tape[arg[1]].cexp_set.empty() )
1554                                                tape[arg[1]].connect_type = yes_connected;
1555                                }
1556                                else    tape[arg[1]].connect_type = yes_connected;
1557                                break;
1558
1559                                default:
1560                                CPPAD_ASSERT_UNKNOWN(false);
1561                        }
1562                        if( connect_type == sum_connected )
1563                        {       // convert sum to csum connection for this variable
1564                                tape[i_var].connect_type = connect_type = csum_connected;
1565                        }
1566                        break; // --------------------------------------------
1567
1568               
1569                        // Special case for AddvvOp and SubvvOp
1570                        case AddvvOp:
1571                        case SubvvOp:
1572                        for(i = 0; i < 2; i++) switch( connect_type )
1573                        {       case not_connected:
1574                                break;
1575
1576                                case yes_connected:
1577                                case sum_connected:
1578                                case csum_connected:
1579                                if( tape[arg[i]].connect_type == not_connected )
1580                                        tape[arg[i]].connect_type = sum_connected;
1581                                else    tape[arg[i]].connect_type = yes_connected;
1582                                break;
1583
1584                                case cexp_connected:
1585                                if( tape[arg[i]].connect_type == not_connected )
1586                                {       tape[arg[i]].connect_type = cexp_connected;
1587                                        tape[arg[i]].cexp_set     = cexp_set;
1588                                }
1589                                else if( tape[arg[i]].connect_type == cexp_connected )
1590                                {       tape[arg[i]].cexp_set = intersection(
1591                                                tape[arg[i]].cexp_set, cexp_set
1592                                        );
1593                                        if( tape[arg[i]].cexp_set.empty() )
1594                                                tape[arg[i]].connect_type = yes_connected;
1595                                }
1596                                else    tape[arg[i]].connect_type = yes_connected;
1597                                break;
1598
1599                                default:
1600                                CPPAD_ASSERT_UNKNOWN(false);
1601                        }
1602                        if( connect_type == sum_connected )
1603                        {       // convert sum to csum connection for this variable
1604                                tape[i_var].connect_type = connect_type = csum_connected;
1605                        }
1606                        break; // --------------------------------------------
1607
1608                        // Other binary operators
1609                        // where operands are arg[0], arg[1]
1610                        case DivvvOp:
1611                        case MulvvOp:
1612                        case PowvvOp:
1613                        for(i = 0; i < 2; i++) switch( connect_type )
1614                        {       case not_connected:
1615                                break;
1616
1617                                case yes_connected:
1618                                case sum_connected:
1619                                case csum_connected:
1620                                tape[arg[i]].connect_type = yes_connected;
1621                                break;
1622
1623                                case cexp_connected:
1624                                if( tape[arg[i]].connect_type == not_connected )
1625                                {       tape[arg[i]].connect_type = cexp_connected;
1626                                        tape[arg[i]].cexp_set     = cexp_set;
1627                                }
1628                                else if( tape[arg[i]].connect_type == cexp_connected )
1629                                {       tape[arg[i]].cexp_set = intersection(
1630                                                tape[arg[i]].cexp_set, cexp_set
1631                                        );
1632                                        if( tape[arg[i]].cexp_set.empty() )
1633                                                tape[arg[i]].connect_type = yes_connected;
1634                                }
1635                                else    tape[arg[i]].connect_type = yes_connected;
1636                                break;
1637
1638                                default:
1639                                CPPAD_ASSERT_UNKNOWN(false);
1640                        }
1641                        break; // --------------------------------------------
1642
1643                        // Conditional expression operators
1644                        case CExpOp:
1645                        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
1646                        if( connect_type != not_connected )
1647                        {       struct_cskip_info info;
1648                                info.cop        = CompareOp( arg[0] );
1649                                info.flag       = arg[1];
1650                                info.left       = arg[2];
1651                                info.right      = arg[3];
1652                                info.n_op_true  = 0;
1653                                info.n_op_false = 0;
1654                                //
1655                                size_t index    = 0;
1656                                if( arg[1] & 1 )
1657                                {       index = std::max(index, info.left);
1658                                        tape[info.left].connect_type = yes_connected;
1659                                }
1660                                if( arg[1] & 2 )
1661                                {       index = std::max(index, info.right);
1662                                        tape[info.right].connect_type = yes_connected;
1663                                }
1664                                CPPAD_ASSERT_UNKNOWN( index > 0 );
1665                                info.max_left_right = index;
1666                                //
1667                                index = cskip_info.size();
1668                                cskip_info.push_back(info);
1669                                //
1670                                if( arg[1] & 4 )
1671                                {       tape[arg[4]].connect_type = cexp_connected;
1672                                        tape[arg[4]].cexp_set     = cexp_set;
1673                                        tape[arg[4]].cexp_set.insert(
1674                                                class_cexp_pair(true, index)
1675                                        );
1676                                }
1677                                if( arg[1] & 8 )
1678                                {       tape[arg[5]].connect_type = cexp_connected;
1679                                        tape[arg[5]].cexp_set     = cexp_set;
1680                                        tape[arg[5]].cexp_set.insert(
1681                                                class_cexp_pair(false, index)
1682                                        );
1683                                }
1684                        }
1685                        break;  // --------------------------------------------
1686
1687                        // Operations where there is nothing to do
1688                        case ComOp:
1689                        case EndOp:
1690                        case ParOp:
1691                        case PriOp:
1692                        break;  // --------------------------------------------
1693
1694                        // Operators that never get removed
1695                        case BeginOp:
1696                        case InvOp:
1697                        tape[i_var].connect_type = yes_connected;
1698                        break;
1699
1700                        // Load using a parameter index
1701                        case LdpOp:
1702                        if( tape[i_var].connect_type != not_connected )
1703                        {
1704                                i                = vecad[ arg[0] - 1 ];
1705                                vecad_connect[i] = yes_connected;
1706                        }
1707                        break; // --------------------------------------------
1708
1709                        // Load using a variable index
1710                        case LdvOp:
1711                        if( tape[i_var].connect_type != not_connected )
1712                        {
1713                                i                    = vecad[ arg[0] - 1 ];
1714                                vecad_connect[i]     = yes_connected;
1715                                tape[arg[1]].connect_type = yes_connected;
1716                        }
1717                        break; // --------------------------------------------
1718
1719                        // Store a variable using a parameter index
1720                        case StpvOp:
1721                        i = vecad[ arg[0] - 1 ];
1722                        if( vecad_connect[i] != not_connected )
1723                                tape[arg[2]].connect_type = yes_connected;
1724                        break; // --------------------------------------------
1725
1726                        // Store a variable using a variable index
1727                        case StvvOp:
1728                        i = vecad[ arg[0] - 1 ];
1729                        if( vecad_connect[i] )
1730                        {       tape[arg[1]].connect_type = yes_connected;
1731                                tape[arg[2]].connect_type = yes_connected;
1732                        }
1733                        break; 
1734                        // ============================================================
1735                        case UserOp:
1736                        // start or end atomic operation sequence
1737                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
1738                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
1739                        if( user_state == user_end )
1740                        {       user_index = arg[0];
1741                                user_id    = arg[1];
1742                                user_n     = arg[2];
1743                                user_m     = arg[3];
1744                                user_q     = 1;
1745                                user_atom  = atomic_base<Base>::class_object(user_index);
1746                                user_set   = user_atom->sparsity() ==
1747                                        atomic_base<Base>::set_sparsity_enum;
1748                                //
1749                                if(user_s_set.size() != user_n )
1750                                        user_s_set.resize(user_n);
1751                                if(user_r_set.size() != user_m )
1752                                        user_r_set.resize(user_m);
1753                                //
1754                                // Note user_q is 1, but use it for clarity of code
1755                                if(user_s_bool.size() != user_n * user_q )
1756                                        user_s_bool.resize(user_n * user_q);
1757                                if(user_r_bool.size() != user_m * user_q )
1758                                        user_r_bool.resize(user_m * user_q);
1759                                //
1760                                user_j     = user_n;
1761                                user_i     = user_m;
1762                                user_state = user_ret;
1763                                //
1764                                struct_user_info info;
1765                                info.connect_type = not_connected;
1766                                info.op_end       = i_op + 1;
1767                                user_info.push_back(info);
1768                               
1769                        }
1770                        else
1771                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
1772                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
1773                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
1774                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
1775                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
1776                                user_state = user_end;
1777                                //
1778                                CPPAD_ASSERT_UNKNOWN( user_curr + 1 == user_info.size() );
1779                                user_info[user_curr].op_begin = i_op;
1780                                user_curr                     = user_info.size();
1781               }
1782                        break;
1783
1784                        case UsrapOp:
1785                        // parameter argument in an atomic operation sequence
1786                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
1787                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
1788                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1789                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
1790                        --user_j;
1791                        if( user_j == 0 )
1792                                user_state = user_start;
1793                        break;
1794
1795                        case UsravOp:
1796                        // variable argument in an atomic operation sequence
1797                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
1798                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
1799                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1800                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
1801                        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
1802                        --user_j;
1803                        if( user_set )
1804                        {       if( ! user_s_set[user_j].empty() )
1805                                        tape[arg[0]].connect_type = 
1806                                                user_info[user_curr].connect_type;
1807                        }
1808                        else
1809                        {       if( user_s_bool[user_j] )
1810                                        tape[arg[0]].connect_type = 
1811                                                user_info[user_curr].connect_type;
1812                        }
1813                        if( user_j == 0 )
1814                                user_state = user_start;
1815                        break;
1816
1817                        case UsrrvOp:
1818                        // variable result in an atomic operation sequence
1819                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
1820                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
1821                        --user_i;
1822                        user_r_set[user_i].clear();
1823                        user_r_bool[user_i] = false;
1824                        switch( connect_type )
1825                        {       case not_connected:
1826                                break;
1827
1828                                case yes_connected:
1829                                case sum_connected:
1830                                case csum_connected:
1831                                user_info[user_curr].connect_type = yes_connected;
1832                                user_r_set[user_i].insert(0);
1833                                user_r_bool[user_i] = true;
1834                                break;
1835
1836                                case cexp_connected:
1837                                if( user_info[user_curr].connect_type == not_connected )
1838                                {       user_info[user_curr].connect_type  = connect_type;
1839                                        user_info[user_curr].cexp_set      = cexp_set;
1840                                }
1841                                else if(user_info[user_curr].connect_type==cexp_connected)
1842                                {       user_info[user_curr].cexp_set = intersection(
1843                                                user_info[user_curr].cexp_set, cexp_set
1844                                        );
1845                                        if( user_info[user_curr].cexp_set.empty() )
1846                                                user_info[user_curr].connect_type = yes_connected;
1847                                }
1848                                else    user_info[user_curr].connect_type = yes_connected;
1849                                break;
1850
1851                                default:
1852                                CPPAD_ASSERT_UNKNOWN(false);
1853                        }
1854                        // drop into op = UsrrpOp code to handle case where user_i == 0
1855                        // for both UsrrvOp and UsrrpOp together.
1856
1857                        case UsrrpOp:
1858                        if( op == UsrrpOp )
1859                        {       // parameter result in an atomic operation sequence
1860                                CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
1861                                CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
1862                                CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1863                                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
1864                                --user_i;
1865                                user_r_set[user_i].clear();
1866                                user_r_bool[user_i] = false;
1867                        }
1868                        if( user_i == 0 )
1869                        {       // call users function for this operation
1870                                user_atom->set_id(user_id);
1871# ifdef NDEBUG
1872                                if( user_set )
1873                                {       user_atom->
1874                                                rev_sparse_jac(user_q, user_r_set, user_s_set);
1875                                }
1876                                else
1877                                {       user_atom->
1878                                                rev_sparse_jac(user_q, user_r_bool, user_s_bool);
1879                                }
1880# else
1881                                bool flag;
1882                                if( user_set )
1883                                {       flag = user_atom->
1884                                                rev_sparse_jac(user_q, user_r_set, user_s_set);
1885                                }
1886                                else
1887                                {       flag = user_atom->
1888                                                rev_sparse_jac(user_q, user_r_bool, user_s_bool);
1889                                }
1890                                if( ! flag )
1891                                {       std::string s =
1892                                                "Optimizing an ADFun object"
1893                                                " that contains the atomic function\n\t";
1894                                        s += user_atom->afun_name();
1895                                        s += "\nCurrent atomic_sparsity is set to";
1896                                        //
1897                                        if( user_set )
1898                                                s += " std::set\nand std::set";
1899                                        else    s += " bool\nand bool";
1900                                        //
1901                                        s += " version of rev_sparse_jac returned false";
1902                                        CPPAD_ASSERT_KNOWN(false, s.c_str() );
1903                                }
1904# endif
1905                                user_state = user_arg;
1906                        }
1907                        break;
1908                        // ============================================================
1909
1910                        // all cases should be handled above
1911                        default:
1912                        CPPAD_ASSERT_UNKNOWN(0);
1913                }
1914        }
1915        // values corresponding to BeginOp
1916        CPPAD_ASSERT_UNKNOWN( i_op == 0 && i_var == 0 && op == BeginOp );
1917        tape[i_var].op = op;
1918        // -------------------------------------------------------------
1919
1920        // Determine which variables can be conditionally skipped
1921        for(i = 0; i < num_var; i++)
1922        {       if( tape[i].connect_type == cexp_connected )
1923                {       std::set<class_cexp_pair>::const_iterator itr = 
1924                                tape[i].cexp_set.begin();
1925                        while( itr != tape[i].cexp_set.end() )
1926                        {       j = itr->index_;
1927                                if( itr->compare_ == true )
1928                                        cskip_info[j].skip_var_false.push_back(i);
1929                                else cskip_info[j].skip_var_true.push_back(i);
1930                                itr++;
1931                        }
1932                }
1933        }
1934        // Determine size of skip information in user_info
1935        for(i = 0; i < user_info.size(); i++)
1936        {       if( user_info[i].connect_type == cexp_connected )
1937                {       std::set<class_cexp_pair>::const_iterator itr = 
1938                                user_info[i].cexp_set.begin();
1939                        while( itr != user_info[i].cexp_set.end() )
1940                        {       j = itr->index_;
1941                                if( itr->compare_ == true )
1942                                        cskip_info[j].n_op_false = 
1943                                                user_info[i].op_end - user_info[i].op_begin;
1944                                else
1945                                        cskip_info[j].n_op_true = 
1946                                                user_info[i].op_end - user_info[i].op_begin;
1947                                itr++;
1948                        }
1949                }
1950        }
1951
1952        // Sort the conditional skip information by the maximum of the
1953        // index for the left and right comparision operands
1954        CppAD::vector<size_t> cskip_info_order( cskip_info.size() );
1955        {       CppAD::vector<size_t> keys( cskip_info.size() );
1956                for(i = 0; i < cskip_info.size(); i++)
1957                        keys[i] = std::max( cskip_info[i].left, cskip_info[i].right );
1958                CppAD::index_sort(keys, cskip_info_order);
1959        }
1960        size_t cskip_info_next = 0;
1961
1962
1963        // Initilaize table mapping hash code to variable index in tape
1964        // as pointing to the BeginOp at the beginning of the tape
1965        CppAD::vector<size_t>  hash_table_var(CPPAD_HASH_TABLE_SIZE);
1966        for(i = 0; i < CPPAD_HASH_TABLE_SIZE; i++)
1967                hash_table_var[i] = 0;
1968        CPPAD_ASSERT_UNKNOWN( tape[0].op == BeginOp );
1969
1970        // initialize mapping from old variable index to new
1971        // operator and variable index
1972        for(i = 0; i < num_var; i++)
1973        {       tape[i].new_op  = 0;       // invalid index (except for BeginOp)
1974                tape[i].new_var = num_var; // invalid index
1975        }
1976
1977        // Erase all information in the old recording
1978        rec->free();
1979
1980        // initialize mapping from old VecAD index to new VecAD index
1981        CppAD::vector<size_t> new_vecad_ind(num_vecad_ind);
1982        for(i = 0; i < num_vecad_ind; i++)
1983                new_vecad_ind[i] = num_vecad_ind; // invalid index
1984
1985        j = 0;     // index into the old set of indices
1986        for(i = 0; i < num_vecad_vec; i++)
1987        {       // length of this VecAD
1988                size_t length = play->GetVecInd(j);
1989                if( vecad_connect[i] != not_connected )
1990                {       // Put this VecAD vector in new recording
1991                        CPPAD_ASSERT_UNKNOWN(length < num_vecad_ind);
1992                        new_vecad_ind[j] = rec->PutVecInd(length);
1993                        for(k = 1; k <= length; k++) new_vecad_ind[j+k] =
1994                                rec->PutVecInd(
1995                                        rec->PutPar(
1996                                                play->GetPar( 
1997                                                        play->GetVecInd(j+k)
1998                        ) ) );
1999                }
2000                // start of next VecAD
2001                j       += length + 1;
2002        }
2003        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
2004
2005        // start playing the operations in the forward direction
2006        play->forward_start(op, arg, i_op, i_var);
2007        CPPAD_ASSERT_UNKNOWN( user_curr == user_info.size() );
2008
2009        // playing forward skips BeginOp at the beginning, but not EndOp at
2010        // the end.  Put BeginOp at beginning of recording
2011        CPPAD_ASSERT_UNKNOWN( op == BeginOp );
2012        CPPAD_ASSERT_NARG_NRES(BeginOp, 1, 1);
2013        tape[i_var].new_op  = rec->num_op_rec();
2014        tape[i_var].new_var = rec->PutOp(BeginOp);
2015        rec->PutArg(0);
2016
2017
2018        // temporary buffer for new argument values
2019        addr_t new_arg[6];
2020
2021        // temporary work space used by record_csum
2022        // (decalared here to avoid realloaction of memory)
2023        struct_csum_stacks csum_work;
2024
2025        // tempory used to hold a size_pair
2026        struct_size_pair size_pair;
2027
2028        user_state = user_start;
2029        while(op != EndOp)
2030        {       // next op
2031                play->forward_next(op, arg, i_op, i_var);
2032                CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
2033                CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
2034
2035                // determine if we should insert a conditional skip here
2036                bool skip = cskip_info_next < cskip_info.size();
2037                skip     &= op != BeginOp;
2038                skip     &= op != InvOp;
2039                skip     &= user_state == user_start;
2040                if( skip )
2041                {       j     = cskip_info_order[cskip_info_next];
2042                        if( NumRes(op) > 0 )
2043                                skip &= cskip_info[j].max_left_right < i_var;
2044                        else
2045                                skip &= cskip_info[j].max_left_right <= i_var;
2046                }
2047                if( skip )
2048                {       cskip_info_next++;
2049                        skip &= cskip_info[j].skip_var_true.size() > 0 ||
2050                                        cskip_info[j].skip_var_false.size() > 0;
2051                        if( skip )
2052                        {       struct_cskip_info info = cskip_info[j];
2053                                CPPAD_ASSERT_UNKNOWN( NumRes(CSkipOp) == 0 );
2054                                size_t n_true  = 
2055                                        info.skip_var_true.size() + info.n_op_true;
2056                                size_t n_false = 
2057                                        info.skip_var_false.size() + info.n_op_false;
2058                                size_t n_arg   = 7 + n_true + n_false; 
2059                                // reserve space for the arguments to this operator but
2060                                // delay setting them until we have all the new addresses
2061                                cskip_info[j].i_arg = rec->ReserveArg(n_arg);
2062                                CPPAD_ASSERT_UNKNOWN( cskip_info[j].i_arg > 0 );
2063                                rec->PutOp(CSkipOp);
2064                        }
2065                        else    cskip_info[j].i_arg = 0;
2066                }
2067
2068                // determine if we should keep this operation in the new
2069                // operation sequence
2070                bool keep;
2071                switch( op )
2072                {       case ComOp:
2073                        case PriOp:
2074                        keep = false;
2075                        break;
2076
2077                        case InvOp:
2078                        case EndOp:
2079                        keep = true;
2080                        break;
2081
2082                        case StppOp:
2083                        case StvpOp:
2084                        case StpvOp:
2085                        case StvvOp:
2086                        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
2087                        i = vecad[ arg[0] - 1 ];
2088                        keep = vecad_connect[i] != not_connected;
2089                        break;
2090
2091                        case AddpvOp:
2092                        case AddvvOp:
2093                        case SubpvOp:
2094                        case SubvpOp:
2095                        case SubvvOp:
2096                        keep  = tape[i_var].connect_type != not_connected;
2097                        keep &= tape[i_var].connect_type != csum_connected;
2098                        break; 
2099
2100                        case UserOp:
2101                        case UsrapOp:
2102                        case UsravOp:
2103                        case UsrrpOp:
2104                        case UsrrvOp:
2105                        keep = true;
2106                        break;
2107
2108                        default:
2109                        keep = tape[i_var].connect_type != not_connected;
2110                        break;
2111                }
2112
2113                size_t         match_var    = 0;
2114                unsigned short code         = 0;
2115                bool           replace_hash = false;
2116                if( keep ) switch( op )
2117                {
2118                        // Unary operator where operand is arg[0]
2119                        case AbsOp:
2120                        case AcosOp:
2121                        case AsinOp:
2122                        case AtanOp:
2123                        case CosOp:
2124                        case CoshOp:
2125                        case ExpOp:
2126                        case LogOp:
2127                        case SignOp:
2128                        case SinOp:
2129                        case SinhOp:
2130                        case SqrtOp:
2131                        case TanOp:
2132                        case TanhOp:
2133                        match_var = unary_match(
2134                                tape                ,  // inputs
2135                                i_var               ,
2136                                play->num_par_rec() ,
2137                                play->GetPar()      ,
2138                                hash_table_var      ,
2139                                code                  // outputs
2140                        );
2141                        if( match_var > 0 )
2142                                tape[i_var].new_var = match_var;
2143                        else
2144                        {
2145                                replace_hash = true;
2146                                new_arg[0]   = tape[ arg[0] ].new_var;
2147                                rec->PutArg( new_arg[0] );
2148                                tape[i_var].new_op  = rec->num_op_rec();
2149                                tape[i_var].new_var = i = rec->PutOp(op);
2150                                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
2151                        }
2152                        break;
2153                        // ---------------------------------------------------
2154                        // Binary operators where
2155                        // left is a variable and right is a parameter
2156                        case SubvpOp:
2157                        if( tape[arg[0]].connect_type == csum_connected )
2158                        {
2159                                // convert to a sequence of summation operators
2160                                size_pair = record_csum(
2161                                        tape                , // inputs
2162                                        i_var               ,
2163                                        play->num_par_rec() ,
2164                                        play->GetPar()      ,
2165                                        rec                 ,
2166                                        csum_work
2167                                );
2168                                tape[i_var].new_op  = size_pair.i_op;
2169                                tape[i_var].new_var = size_pair.i_var;
2170                                // abort rest of this case
2171                                break;
2172                        }
2173                        case DivvpOp:
2174                        case PowvpOp:
2175                        match_var = binary_match(
2176                                tape                ,  // inputs
2177                                i_var               ,
2178                                play->num_par_rec() ,
2179                                play->GetPar()      ,
2180                                hash_table_var      ,
2181                                code                  // outputs
2182                        );
2183                        if( match_var > 0 )
2184                                tape[i_var].new_var = match_var;
2185                        else
2186                        {       size_pair = record_vp(
2187                                        tape                , // inputs
2188                                        i_var               ,
2189                                        play->num_par_rec() ,
2190                                        play->GetPar()      ,
2191                                        rec                 ,
2192                                        op                  ,
2193                                        arg
2194                                );
2195                                tape[i_var].new_op  = size_pair.i_op;
2196                                tape[i_var].new_var = size_pair.i_var;
2197                                replace_hash = true;
2198                        }
2199                        break;
2200                        // ---------------------------------------------------
2201                        // Binary operators where
2202                        // left is an index and right is a variable
2203                        case DisOp:
2204                        match_var = binary_match(
2205                                tape                ,  // inputs
2206                                i_var               ,
2207                                play->num_par_rec() ,
2208                                play->GetPar()      ,
2209                                hash_table_var      ,
2210                                code                  // outputs
2211                        );
2212                        if( match_var > 0 )
2213                                tape[i_var].new_var = match_var;
2214                        else
2215                        {       new_arg[0] = arg[0];
2216                                new_arg[1] = tape[ arg[1] ].new_var;
2217                                rec->PutArg( new_arg[0], new_arg[1] ); 
2218                                tape[i_var].new_op  = rec->num_op_rec();
2219                                tape[i_var].new_var = rec->PutOp(op);
2220                                CPPAD_ASSERT_UNKNOWN( 
2221                                        size_t(new_arg[1]) < tape[i_var].new_var
2222                                );
2223                                replace_hash = true;
2224                        }
2225                        break;
2226
2227                        // ---------------------------------------------------
2228                        // Binary operators where
2229                        // left is a parameter and right is a variable
2230                        case SubpvOp:
2231                        case AddpvOp:
2232                        if( tape[arg[1]].connect_type == csum_connected )
2233                        {
2234                                // convert to a sequence of summation operators
2235                                size_pair = record_csum(
2236                                        tape                , // inputs
2237                                        i_var               ,
2238                                        play->num_par_rec() ,
2239                                        play->GetPar()      ,
2240                                        rec                 ,
2241                                        csum_work
2242                                );
2243                                tape[i_var].new_op  = size_pair.i_op;
2244                                tape[i_var].new_var = size_pair.i_var;
2245                                // abort rest of this case
2246                                break;
2247                        }
2248                        case DivpvOp:
2249                        case MulpvOp:
2250                        case PowpvOp:
2251                        match_var = binary_match(
2252                                tape                ,  // inputs
2253                                i_var               ,
2254                                play->num_par_rec() ,
2255                                play->GetPar()      ,
2256                                hash_table_var      ,
2257                                code                  // outputs
2258                        );
2259                        if( match_var > 0 )
2260                                tape[i_var].new_var = match_var;
2261                        else
2262                        {       size_pair = record_pv(
2263                                        tape                , // inputs
2264                                        i_var               ,
2265                                        play->num_par_rec() ,
2266                                        play->GetPar()      ,
2267                                        rec                 ,
2268                                        op                  ,
2269                                        arg
2270                                );
2271                                tape[i_var].new_op  = size_pair.i_op;
2272                                tape[i_var].new_var = size_pair.i_var;
2273                                replace_hash = true;
2274                        }
2275                        break;
2276                        // ---------------------------------------------------
2277                        // Binary operator where
2278                        // both operators are variables
2279                        case AddvvOp:
2280                        case SubvvOp:
2281                        if( (tape[arg[0]].connect_type == csum_connected) |
2282                            (tape[arg[1]].connect_type == csum_connected)
2283                        )
2284                        {
2285                                // convert to a sequence of summation operators
2286                                size_pair = record_csum(
2287                                        tape                , // inputs
2288                                        i_var               ,
2289                                        play->num_par_rec() ,
2290                                        play->GetPar()      ,
2291                                        rec                 ,
2292                                        csum_work
2293                                );
2294                                tape[i_var].new_op  = size_pair.i_op;
2295                                tape[i_var].new_var = size_pair.i_var;
2296                                // abort rest of this case
2297                                break;
2298                        }
2299                        case DivvvOp:
2300                        case MulvvOp:
2301                        case PowvvOp:
2302                        match_var = binary_match(
2303                                tape                ,  // inputs
2304                                i_var               ,
2305                                play->num_par_rec() ,
2306                                play->GetPar()      ,
2307                                hash_table_var      ,
2308                                code                  // outputs
2309                        );
2310                        if( match_var > 0 )
2311                                tape[i_var].new_var = match_var;
2312                        else
2313                        {       size_pair = record_vv(
2314                                        tape                , // inputs
2315                                        i_var               ,
2316                                        play->num_par_rec() ,
2317                                        play->GetPar()      ,
2318                                        rec                 ,
2319                                        op                  ,
2320                                        arg
2321                                );
2322                                tape[i_var].new_op  = size_pair.i_op;
2323                                tape[i_var].new_var = size_pair.i_var;
2324                                replace_hash = true;
2325                        }
2326                        break;
2327                        // ---------------------------------------------------
2328                        // Conditional expression operators
2329                        case CExpOp:
2330                        CPPAD_ASSERT_NARG_NRES(op, 6, 1);
2331                        new_arg[0] = arg[0];
2332                        new_arg[1] = arg[1];
2333                        mask = 1;
2334                        for(i = 2; i < 6; i++)
2335                        {       if( arg[1] & mask )
2336                                {       new_arg[i] = tape[arg[i]].new_var;
2337                                        CPPAD_ASSERT_UNKNOWN( 
2338                                                size_t(new_arg[i]) < num_var
2339                                        );
2340                                }
2341                                else    new_arg[i] = rec->PutPar( 
2342                                                play->GetPar( arg[i] )
2343                                );
2344                                mask = mask << 1;
2345                        }
2346                        rec->PutArg(
2347                                new_arg[0] ,
2348                                new_arg[1] ,
2349                                new_arg[2] ,
2350                                new_arg[3] ,
2351                                new_arg[4] ,
2352                                new_arg[5] 
2353                        );
2354                        tape[i_var].new_op  = rec->num_op_rec();
2355                        tape[i_var].new_var = rec->PutOp(op);
2356                        break;
2357                        // ---------------------------------------------------
2358                        // Operations with no arguments and no results
2359                        case EndOp:
2360                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
2361                        rec->PutOp(op);
2362                        break;
2363                        // ---------------------------------------------------
2364                        // Operations with no arguments and one result
2365                        case InvOp:
2366                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
2367                        tape[i_var].new_op  = rec->num_op_rec();
2368                        tape[i_var].new_var = rec->PutOp(op);
2369                        break;
2370                        // ---------------------------------------------------
2371                        // Operations with one argument that is a parameter
2372                        case ParOp:
2373                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
2374                        new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
2375
2376                        rec->PutArg( new_arg[0] );
2377                        tape[i_var].new_op  = rec->num_op_rec();
2378                        tape[i_var].new_var = rec->PutOp(op);
2379                        break;
2380                        // ---------------------------------------------------
2381                        // Load using a parameter index
2382                        case LdpOp:
2383                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
2384                        new_arg[0] = new_vecad_ind[ arg[0] ];
2385                        new_arg[1] = arg[1];
2386                        new_arg[2] = rec->num_load_op_rec();
2387                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2388                        rec->PutArg( 
2389                                new_arg[0], 
2390                                new_arg[1], 
2391                                new_arg[2]
2392                        );
2393                        tape[i_var].new_op  = rec->num_op_rec();
2394                        tape[i_var].new_var = rec->PutLoadOp(op);
2395                        break;
2396                        // ---------------------------------------------------
2397                        // Load using a variable index
2398                        case LdvOp:
2399                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
2400                        new_arg[0] = new_vecad_ind[ arg[0] ];
2401                        new_arg[1] = tape[arg[1]].new_var;
2402                        new_arg[2] = rec->num_load_op_rec();
2403                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2404                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
2405                        rec->PutArg( 
2406                                new_arg[0], 
2407                                new_arg[1], 
2408                                new_arg[2]
2409                        );
2410                        tape[i_var].new_var = rec->num_op_rec();
2411                        tape[i_var].new_var = rec->PutLoadOp(op);
2412                        break;
2413                        // ---------------------------------------------------
2414                        // Store a parameter using a parameter index
2415                        case StppOp:
2416                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
2417                        new_arg[0] = new_vecad_ind[ arg[0] ];
2418                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
2419                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
2420                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2421                        rec->PutArg(
2422                                new_arg[0], 
2423                                new_arg[1], 
2424                                new_arg[2]
2425                        );
2426                        rec->PutOp(op);
2427                        break;
2428                        // ---------------------------------------------------
2429                        // Store a parameter using a variable index
2430                        case StvpOp:
2431                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
2432                        new_arg[0] = new_vecad_ind[ arg[0] ];
2433                        new_arg[1] = tape[arg[1]].new_var;
2434                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
2435                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2436                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
2437                        rec->PutArg(
2438                                new_arg[0], 
2439                                new_arg[1], 
2440                                new_arg[2]
2441                        );
2442                        rec->PutOp(op);
2443                        break;
2444                        // ---------------------------------------------------
2445                        // Store a variable using a parameter index
2446                        case StpvOp:
2447                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
2448                        new_arg[0] = new_vecad_ind[ arg[0] ];
2449                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
2450                        new_arg[2] = tape[arg[2]].new_var;
2451                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2452                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
2453                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
2454                        rec->PutArg(
2455                                new_arg[0], 
2456                                new_arg[1], 
2457                                new_arg[2]
2458                        );
2459                        rec->PutOp(op);
2460                        break;
2461                        // ---------------------------------------------------
2462                        // Store a variable using a variable index
2463                        case StvvOp:
2464                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
2465                        new_arg[0] = new_vecad_ind[ arg[0] ];
2466                        new_arg[1] = tape[arg[1]].new_var;
2467                        new_arg[2] = tape[arg[2]].new_var;
2468                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2469                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
2470                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
2471                        rec->PutArg(
2472                                new_arg[0], 
2473                                new_arg[1], 
2474                                new_arg[2]
2475                        );
2476                        rec->PutOp(op);
2477                        break;
2478
2479                        // -----------------------------------------------------------
2480                        case UserOp:
2481                        CPPAD_ASSERT_NARG_NRES(op, 4, 0);
2482                        if( user_state == user_start )
2483                        {       user_state = user_arg;
2484                                CPPAD_ASSERT_UNKNOWN( user_curr > 0 );
2485                                user_curr--;
2486                                user_info[user_curr].op_begin = rec->num_op_rec();
2487                        }
2488                        else
2489                        {       user_state = user_start;
2490                                user_info[user_curr].op_end = rec->num_op_rec() + 1;
2491                        }
2492                        // user_index, user_id, user_n, user_m
2493                        if( user_info[user_curr].connect_type != not_connected )
2494                        {       rec->PutArg(arg[0], arg[1], arg[2], arg[3]);
2495                                rec->PutOp(UserOp);
2496                        }
2497                        break;
2498
2499                        case UsrapOp:
2500                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
2501                        if( user_info[user_curr].connect_type != not_connected )
2502                        {       new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
2503                                rec->PutArg(new_arg[0]);
2504                                rec->PutOp(UsrapOp);
2505                        }
2506                        break;
2507
2508                        case UsravOp:
2509                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
2510                        if( user_info[user_curr].connect_type != not_connected )
2511                        {       new_arg[0] = tape[arg[0]].new_var;
2512                                if( size_t(new_arg[0]) < num_var )
2513                                {       rec->PutArg(new_arg[0]);
2514                                        rec->PutOp(UsravOp);
2515                                }
2516                                else
2517                                {       // This argument does not affect the result and
2518                                        // has been optimized out so use nan in its place.
2519                                        new_arg[0] = rec->PutPar( nan(Base(0)) );
2520                                        rec->PutArg(new_arg[0]);
2521                                        rec->PutOp(UsrapOp);
2522                                }
2523                        }
2524                        break;
2525
2526                        case UsrrpOp:
2527                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
2528                        if( user_info[user_curr].connect_type != not_connected )
2529                        {       new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
2530                                rec->PutArg(new_arg[0]);
2531                                rec->PutOp(UsrrpOp);
2532                        }
2533                        break;
2534                       
2535                        case UsrrvOp:
2536                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
2537                        if( user_info[user_curr].connect_type != not_connected )
2538                        {       tape[i_var].new_op  = rec->num_op_rec();
2539                                tape[i_var].new_var = rec->PutOp(UsrrvOp);
2540                        }
2541                        break;
2542                        // ---------------------------------------------------
2543
2544                        // all cases should be handled above
2545                        default:
2546                        CPPAD_ASSERT_UNKNOWN(false);
2547
2548                }
2549                if( replace_hash )
2550                {       // The old variable index i_var corresponds to the
2551                        // new variable index tape[i_var].new_var. In addition
2552                        // this is the most recent variable that has this code.
2553                        hash_table_var[code] = i_var;
2554                }
2555
2556        }
2557        // modify the dependent variable vector to new indices
2558        for(i = 0; i < dep_taddr.size(); i++ )
2559        {       CPPAD_ASSERT_UNKNOWN( size_t(tape[dep_taddr[i]].new_var) < num_var );
2560                dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
2561        }
2562
2563# ifndef NDEBUG
2564        size_t num_new_op = rec->num_op_rec();
2565        for(i_var = 0; i_var < tape.size(); i_var++)
2566                CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op < num_new_op );
2567# endif
2568
2569        // Move skip information from user_info to cskip_info
2570        for(i = 0; i < user_info.size(); i++)
2571        {       if( user_info[i].connect_type == cexp_connected )
2572                {       std::set<class_cexp_pair>::const_iterator itr =
2573                                user_info[i].cexp_set.begin();
2574                        while( itr != user_info[i].cexp_set.end() )
2575                        {       j = itr->index_;
2576                                k = user_info[i].op_begin;
2577                                while(k < user_info[i].op_end)
2578                                {       if( itr->compare_ == true )
2579                                                cskip_info[j].skip_op_false.push_back(k++);
2580                                        else    cskip_info[j].skip_op_true.push_back(k++);
2581                                }
2582                                itr++;
2583                        }
2584                }
2585        }
2586
2587        // fill in the arguments for the CSkip operations
2588        CPPAD_ASSERT_UNKNOWN( cskip_info_next == cskip_info.size() );
2589        for(i = 0; i < cskip_info.size(); i++)
2590        {       struct_cskip_info info = cskip_info[i];
2591                if( info.i_arg > 0 )
2592                {       CPPAD_ASSERT_UNKNOWN( info.n_op_true==info.skip_op_true.size() );
2593                        CPPAD_ASSERT_UNKNOWN(info.n_op_false==info.skip_op_false.size());
2594                        size_t n_true  = 
2595                                info.skip_var_true.size() + info.skip_op_true.size();
2596                        size_t n_false = 
2597                                info.skip_var_false.size() + info.skip_op_false.size();
2598                        size_t i_arg   = info.i_arg;
2599                        rec->ReplaceArg(i_arg++, info.cop   );
2600                        rec->ReplaceArg(i_arg++, info.flag  );
2601                        rec->ReplaceArg(i_arg++, info.left  );
2602                        rec->ReplaceArg(i_arg++, info.right );
2603                        rec->ReplaceArg(i_arg++, n_true     );
2604                        rec->ReplaceArg(i_arg++, n_false    );
2605                        for(j = 0; j < info.skip_var_true.size(); j++)
2606                        {       i_var = info.skip_var_true[j];
2607                                CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op > 0 );
2608                                rec->ReplaceArg(i_arg++, tape[i_var].new_op );
2609                        } 
2610                        for(j = 0; j < info.skip_op_true.size(); j++)
2611                        {       i_op = info.skip_op_true[j];
2612                                rec->ReplaceArg(i_arg++, i_op);
2613                        } 
2614                        for(j = 0; j < info.skip_var_false.size(); j++)
2615                        {       i_var = info.skip_var_false[j];
2616                                CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op > 0 );
2617                                rec->ReplaceArg(i_arg++, tape[i_var].new_op );
2618                        } 
2619                        for(j = 0; j < info.skip_op_false.size(); j++)
2620                        {       i_op = info.skip_op_false[j];
2621                                rec->ReplaceArg(i_arg++, i_op);
2622                        } 
2623                        rec->ReplaceArg(i_arg++, n_true + n_false);
2624# ifndef NDEBUG
2625                        size_t n_arg   = 7 + n_true + n_false; 
2626                        CPPAD_ASSERT_UNKNOWN( info.i_arg + n_arg == i_arg );
2627# endif
2628                }
2629        }
2630}
2631
2632} // END_CPPAD_OPTIMIZE_NAMESPACE
2633
2634/*!
2635Optimize a player object operation sequence
2636
2637The operation sequence for this object is replaced by one with fewer operations
2638but the same funcition and derivative values.
2639
2640\tparam Base
2641base type for the operator; i.e., this operation was recorded
2642using AD< \a Base > and computations by this routine are done using type
2643\a Base.
2644
2645*/
2646template <class Base>
2647void ADFun<Base>::optimize(void)
2648{       // place to store the optimized version of the recording
2649        recorder<Base> rec;
2650
2651        // number of independent variables
2652        size_t n = ind_taddr_.size();
2653
2654# ifndef NDEBUG
2655        size_t i, j, m = dep_taddr_.size();
2656        CppAD::vector<Base> x(n), y(m), check(m);
2657        Base max_taylor(0);
2658        bool check_zero_order = num_order_taylor_ > 0;
2659        if( check_zero_order )
2660        {       // zero order coefficients for independent vars
2661                for(j = 0; j < n; j++)
2662                {       CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
2663                        CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
2664                        x[j] = taylor_[ ind_taddr_[j] * cap_order_taylor_ + 0];
2665                }
2666                // zero order coefficients for dependent vars
2667                for(i = 0; i < m; i++)
2668                {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_  );
2669                        y[i] = taylor_[ dep_taddr_[i] * cap_order_taylor_ + 0];
2670                }
2671                // maximum zero order coefficient not counting BeginOp at beginning
2672                // (which is correpsonds to uninitialized memory).
2673                for(i = 1; i < num_var_tape_; i++)
2674                {       if(  abs_geq(taylor_[i*cap_order_taylor_+0] , max_taylor) )
2675                                max_taylor = taylor_[i*cap_order_taylor_+0];
2676                }
2677        }
2678# endif
2679
2680        // create the optimized recording
2681        CppAD::optimize::optimize_run<Base>(n, dep_taddr_, &play_, &rec);
2682
2683        // number of variables in the recording
2684        num_var_tape_  = rec.num_var_rec();
2685
2686        // now replace the recording
2687        play_.get(rec);
2688
2689        // free memory allocated for sparse Jacobian calculation
2690        // (the results are no longer valid)
2691        for_jac_sparse_pack_.resize(0, 0);
2692        for_jac_sparse_set_.resize(0,0);
2693
2694        // free old Taylor coefficient memory
2695        taylor_.free();
2696        num_order_taylor_     = 0;
2697        cap_order_taylor_     = 0;
2698
2699        // resize and initilaize conditional skip vector
2700        // (must use player size because it now has the recoreder information)
2701        cskip_op_.erase();
2702        cskip_op_.extend( play_.num_op_rec() );
2703
2704# ifndef NDEBUG
2705        if( check_zero_order )
2706        {
2707                // zero order forward calculation using new operation sequence
2708                check = Forward(0, x);
2709
2710                // check results
2711                Base eps = 10. * epsilon<Base>();
2712                for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
2713                        abs_geq( eps * max_taylor , check[i] - y[i] ) ,
2714                        "Error during check of f.optimize()."
2715                );
2716
2717                // Erase memory that this calculation was done so NDEBUG gives
2718                // same final state for this object (from users perspective)
2719                num_order_taylor_     = 0;
2720        }
2721# endif
2722}
2723
2724} // END_CPPAD_NAMESPACE
2725# endif
Note: See TracBrowser for help on using the repository browser.