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

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

atomic_atomic.sh: This bug has been fixed.
opt_atomic.sh: new bug being demonstrated.
template.sh: template file for creating new bug report.
optimize.hpp: fix spelling in error message.

  • Property svn:keywords set to Id
File size: 79.7 KB
Line 
1/* $Id: optimize.hpp 3364 2014-09-27 03:09:47Z 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) & (op != InvOp);
2038                if( skip )
2039                {       j     = cskip_info_order[cskip_info_next];
2040                        if( NumRes(op) > 0 )
2041                                skip &= cskip_info[j].max_left_right < i_var;
2042                        else
2043                                skip &= cskip_info[j].max_left_right <= i_var;
2044                }
2045                if( skip )
2046                {       cskip_info_next++;
2047                        skip &= cskip_info[j].skip_var_true.size() > 0 ||
2048                                        cskip_info[j].skip_var_false.size() > 0;
2049                        if( skip )
2050                        {       struct_cskip_info info = cskip_info[j];
2051                                CPPAD_ASSERT_UNKNOWN( NumRes(CSkipOp) == 0 );
2052                                size_t n_true  = 
2053                                        info.skip_var_true.size() + info.n_op_true;
2054                                size_t n_false = 
2055                                        info.skip_var_false.size() + info.n_op_false;
2056                                size_t n_arg   = 7 + n_true + n_false; 
2057                                // reserve space for the arguments to this operator but
2058                                // delay setting them until we have all the new addresses
2059                                cskip_info[j].i_arg = rec->ReserveArg(n_arg);
2060                                CPPAD_ASSERT_UNKNOWN( cskip_info[j].i_arg > 0 );
2061                                rec->PutOp(CSkipOp);
2062                        }
2063                        else    cskip_info[j].i_arg = 0;
2064                }
2065
2066                // determine if we should keep this operation in the new
2067                // operation sequence
2068                bool keep;
2069                switch( op )
2070                {       case ComOp:
2071                        case PriOp:
2072                        keep = false;
2073                        break;
2074
2075                        case InvOp:
2076                        case EndOp:
2077                        keep = true;
2078                        break;
2079
2080                        case StppOp:
2081                        case StvpOp:
2082                        case StpvOp:
2083                        case StvvOp:
2084                        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
2085                        i = vecad[ arg[0] - 1 ];
2086                        keep = vecad_connect[i] != not_connected;
2087                        break;
2088
2089                        case AddpvOp:
2090                        case AddvvOp:
2091                        case SubpvOp:
2092                        case SubvpOp:
2093                        case SubvvOp:
2094                        keep  = tape[i_var].connect_type != not_connected;
2095                        keep &= tape[i_var].connect_type != csum_connected;
2096                        break; 
2097
2098                        case UserOp:
2099                        case UsrapOp:
2100                        case UsravOp:
2101                        case UsrrpOp:
2102                        case UsrrvOp:
2103                        keep = true;
2104                        break;
2105
2106                        default:
2107                        keep = tape[i_var].connect_type != not_connected;
2108                        break;
2109                }
2110
2111                size_t         match_var    = 0;
2112                unsigned short code         = 0;
2113                bool           replace_hash = false;
2114                if( keep ) switch( op )
2115                {
2116                        // Unary operator where operand is arg[0]
2117                        case AbsOp:
2118                        case AcosOp:
2119                        case AsinOp:
2120                        case AtanOp:
2121                        case CosOp:
2122                        case CoshOp:
2123                        case ExpOp:
2124                        case LogOp:
2125                        case SignOp:
2126                        case SinOp:
2127                        case SinhOp:
2128                        case SqrtOp:
2129                        case TanOp:
2130                        case TanhOp:
2131                        match_var = unary_match(
2132                                tape                ,  // inputs
2133                                i_var               ,
2134                                play->num_par_rec() ,
2135                                play->GetPar()      ,
2136                                hash_table_var      ,
2137                                code                  // outputs
2138                        );
2139                        if( match_var > 0 )
2140                                tape[i_var].new_var = match_var;
2141                        else
2142                        {
2143                                replace_hash = true;
2144                                new_arg[0]   = tape[ arg[0] ].new_var;
2145                                rec->PutArg( new_arg[0] );
2146                                tape[i_var].new_op  = rec->num_op_rec();
2147                                tape[i_var].new_var = i = rec->PutOp(op);
2148                                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
2149                        }
2150                        break;
2151                        // ---------------------------------------------------
2152                        // Binary operators where
2153                        // left is a variable and right is a parameter
2154                        case SubvpOp:
2155                        if( tape[arg[0]].connect_type == csum_connected )
2156                        {
2157                                // convert to a sequence of summation operators
2158                                size_pair = record_csum(
2159                                        tape                , // inputs
2160                                        i_var               ,
2161                                        play->num_par_rec() ,
2162                                        play->GetPar()      ,
2163                                        rec                 ,
2164                                        csum_work
2165                                );
2166                                tape[i_var].new_op  = size_pair.i_op;
2167                                tape[i_var].new_var = size_pair.i_var;
2168                                // abort rest of this case
2169                                break;
2170                        }
2171                        case DivvpOp:
2172                        case PowvpOp:
2173                        match_var = binary_match(
2174                                tape                ,  // inputs
2175                                i_var               ,
2176                                play->num_par_rec() ,
2177                                play->GetPar()      ,
2178                                hash_table_var      ,
2179                                code                  // outputs
2180                        );
2181                        if( match_var > 0 )
2182                                tape[i_var].new_var = match_var;
2183                        else
2184                        {       size_pair = record_vp(
2185                                        tape                , // inputs
2186                                        i_var               ,
2187                                        play->num_par_rec() ,
2188                                        play->GetPar()      ,
2189                                        rec                 ,
2190                                        op                  ,
2191                                        arg
2192                                );
2193                                tape[i_var].new_op  = size_pair.i_op;
2194                                tape[i_var].new_var = size_pair.i_var;
2195                                replace_hash = true;
2196                        }
2197                        break;
2198                        // ---------------------------------------------------
2199                        // Binary operators where
2200                        // left is an index and right is a variable
2201                        case DisOp:
2202                        match_var = binary_match(
2203                                tape                ,  // inputs
2204                                i_var               ,
2205                                play->num_par_rec() ,
2206                                play->GetPar()      ,
2207                                hash_table_var      ,
2208                                code                  // outputs
2209                        );
2210                        if( match_var > 0 )
2211                                tape[i_var].new_var = match_var;
2212                        else
2213                        {       new_arg[0] = arg[0];
2214                                new_arg[1] = tape[ arg[1] ].new_var;
2215                                rec->PutArg( new_arg[0], new_arg[1] ); 
2216                                tape[i_var].new_op  = rec->num_op_rec();
2217                                tape[i_var].new_var = rec->PutOp(op);
2218                                CPPAD_ASSERT_UNKNOWN( 
2219                                        size_t(new_arg[1]) < tape[i_var].new_var
2220                                );
2221                                replace_hash = true;
2222                        }
2223                        break;
2224
2225                        // ---------------------------------------------------
2226                        // Binary operators where
2227                        // left is a parameter and right is a variable
2228                        case SubpvOp:
2229                        case AddpvOp:
2230                        if( tape[arg[1]].connect_type == csum_connected )
2231                        {
2232                                // convert to a sequence of summation operators
2233                                size_pair = record_csum(
2234                                        tape                , // inputs
2235                                        i_var               ,
2236                                        play->num_par_rec() ,
2237                                        play->GetPar()      ,
2238                                        rec                 ,
2239                                        csum_work
2240                                );
2241                                tape[i_var].new_op  = size_pair.i_op;
2242                                tape[i_var].new_var = size_pair.i_var;
2243                                // abort rest of this case
2244                                break;
2245                        }
2246                        case DivpvOp:
2247                        case MulpvOp:
2248                        case PowpvOp:
2249                        match_var = binary_match(
2250                                tape                ,  // inputs
2251                                i_var               ,
2252                                play->num_par_rec() ,
2253                                play->GetPar()      ,
2254                                hash_table_var      ,
2255                                code                  // outputs
2256                        );
2257                        if( match_var > 0 )
2258                                tape[i_var].new_var = match_var;
2259                        else
2260                        {       size_pair = record_pv(
2261                                        tape                , // inputs
2262                                        i_var               ,
2263                                        play->num_par_rec() ,
2264                                        play->GetPar()      ,
2265                                        rec                 ,
2266                                        op                  ,
2267                                        arg
2268                                );
2269                                tape[i_var].new_op  = size_pair.i_op;
2270                                tape[i_var].new_var = size_pair.i_var;
2271                                replace_hash = true;
2272                        }
2273                        break;
2274                        // ---------------------------------------------------
2275                        // Binary operator where
2276                        // both operators are variables
2277                        case AddvvOp:
2278                        case SubvvOp:
2279                        if( (tape[arg[0]].connect_type == csum_connected) |
2280                            (tape[arg[1]].connect_type == csum_connected)
2281                        )
2282                        {
2283                                // convert to a sequence of summation operators
2284                                size_pair = record_csum(
2285                                        tape                , // inputs
2286                                        i_var               ,
2287                                        play->num_par_rec() ,
2288                                        play->GetPar()      ,
2289                                        rec                 ,
2290                                        csum_work
2291                                );
2292                                tape[i_var].new_op  = size_pair.i_op;
2293                                tape[i_var].new_var = size_pair.i_var;
2294                                // abort rest of this case
2295                                break;
2296                        }
2297                        case DivvvOp:
2298                        case MulvvOp:
2299                        case PowvvOp:
2300                        match_var = binary_match(
2301                                tape                ,  // inputs
2302                                i_var               ,
2303                                play->num_par_rec() ,
2304                                play->GetPar()      ,
2305                                hash_table_var      ,
2306                                code                  // outputs
2307                        );
2308                        if( match_var > 0 )
2309                                tape[i_var].new_var = match_var;
2310                        else
2311                        {       size_pair = record_vv(
2312                                        tape                , // inputs
2313                                        i_var               ,
2314                                        play->num_par_rec() ,
2315                                        play->GetPar()      ,
2316                                        rec                 ,
2317                                        op                  ,
2318                                        arg
2319                                );
2320                                tape[i_var].new_op  = size_pair.i_op;
2321                                tape[i_var].new_var = size_pair.i_var;
2322                                replace_hash = true;
2323                        }
2324                        break;
2325                        // ---------------------------------------------------
2326                        // Conditional expression operators
2327                        case CExpOp:
2328                        CPPAD_ASSERT_NARG_NRES(op, 6, 1);
2329                        new_arg[0] = arg[0];
2330                        new_arg[1] = arg[1];
2331                        mask = 1;
2332                        for(i = 2; i < 6; i++)
2333                        {       if( arg[1] & mask )
2334                                {       new_arg[i] = tape[arg[i]].new_var;
2335                                        CPPAD_ASSERT_UNKNOWN( 
2336                                                size_t(new_arg[i]) < num_var
2337                                        );
2338                                }
2339                                else    new_arg[i] = rec->PutPar( 
2340                                                play->GetPar( arg[i] )
2341                                );
2342                                mask = mask << 1;
2343                        }
2344                        rec->PutArg(
2345                                new_arg[0] ,
2346                                new_arg[1] ,
2347                                new_arg[2] ,
2348                                new_arg[3] ,
2349                                new_arg[4] ,
2350                                new_arg[5] 
2351                        );
2352                        tape[i_var].new_op  = rec->num_op_rec();
2353                        tape[i_var].new_var = rec->PutOp(op);
2354                        break;
2355                        // ---------------------------------------------------
2356                        // Operations with no arguments and no results
2357                        case EndOp:
2358                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
2359                        rec->PutOp(op);
2360                        break;
2361                        // ---------------------------------------------------
2362                        // Operations with no arguments and one result
2363                        case InvOp:
2364                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
2365                        tape[i_var].new_op  = rec->num_op_rec();
2366                        tape[i_var].new_var = rec->PutOp(op);
2367                        break;
2368                        // ---------------------------------------------------
2369                        // Operations with one argument that is a parameter
2370                        case ParOp:
2371                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
2372                        new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
2373
2374                        rec->PutArg( new_arg[0] );
2375                        tape[i_var].new_op  = rec->num_op_rec();
2376                        tape[i_var].new_var = rec->PutOp(op);
2377                        break;
2378                        // ---------------------------------------------------
2379                        // Load using a parameter index
2380                        case LdpOp:
2381                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
2382                        new_arg[0] = new_vecad_ind[ arg[0] ];
2383                        new_arg[1] = arg[1];
2384                        new_arg[2] = rec->num_load_op_rec();
2385                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2386                        rec->PutArg( 
2387                                new_arg[0], 
2388                                new_arg[1], 
2389                                new_arg[2]
2390                        );
2391                        tape[i_var].new_op  = rec->num_op_rec();
2392                        tape[i_var].new_var = rec->PutLoadOp(op);
2393                        break;
2394                        // ---------------------------------------------------
2395                        // Load using a variable index
2396                        case LdvOp:
2397                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
2398                        new_arg[0] = new_vecad_ind[ arg[0] ];
2399                        new_arg[1] = tape[arg[1]].new_var;
2400                        new_arg[2] = rec->num_load_op_rec();
2401                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2402                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
2403                        rec->PutArg( 
2404                                new_arg[0], 
2405                                new_arg[1], 
2406                                new_arg[2]
2407                        );
2408                        tape[i_var].new_var = rec->num_op_rec();
2409                        tape[i_var].new_var = rec->PutLoadOp(op);
2410                        break;
2411                        // ---------------------------------------------------
2412                        // Store a parameter using a parameter index
2413                        case StppOp:
2414                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
2415                        new_arg[0] = new_vecad_ind[ arg[0] ];
2416                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
2417                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
2418                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2419                        rec->PutArg(
2420                                new_arg[0], 
2421                                new_arg[1], 
2422                                new_arg[2]
2423                        );
2424                        rec->PutOp(op);
2425                        break;
2426                        // ---------------------------------------------------
2427                        // Store a parameter using a variable index
2428                        case StvpOp:
2429                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
2430                        new_arg[0] = new_vecad_ind[ arg[0] ];
2431                        new_arg[1] = tape[arg[1]].new_var;
2432                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
2433                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2434                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
2435                        rec->PutArg(
2436                                new_arg[0], 
2437                                new_arg[1], 
2438                                new_arg[2]
2439                        );
2440                        rec->PutOp(op);
2441                        break;
2442                        // ---------------------------------------------------
2443                        // Store a variable using a parameter index
2444                        case StpvOp:
2445                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
2446                        new_arg[0] = new_vecad_ind[ arg[0] ];
2447                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
2448                        new_arg[2] = tape[arg[2]].new_var;
2449                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2450                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
2451                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
2452                        rec->PutArg(
2453                                new_arg[0], 
2454                                new_arg[1], 
2455                                new_arg[2]
2456                        );
2457                        rec->PutOp(op);
2458                        break;
2459                        // ---------------------------------------------------
2460                        // Store a variable using a variable index
2461                        case StvvOp:
2462                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
2463                        new_arg[0] = new_vecad_ind[ arg[0] ];
2464                        new_arg[1] = tape[arg[1]].new_var;
2465                        new_arg[2] = tape[arg[2]].new_var;
2466                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
2467                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
2468                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
2469                        rec->PutArg(
2470                                new_arg[0], 
2471                                new_arg[1], 
2472                                new_arg[2]
2473                        );
2474                        rec->PutOp(op);
2475                        break;
2476
2477                        // -----------------------------------------------------------
2478                        case UserOp:
2479                        CPPAD_ASSERT_NARG_NRES(op, 4, 0);
2480                        if( user_state == user_start )
2481                        {       user_state = user_arg;
2482                                CPPAD_ASSERT_UNKNOWN( user_curr > 0 );
2483                                user_curr--;
2484                                user_info[user_curr].op_begin = rec->num_op_rec();
2485                        }
2486                        else
2487                        {       user_state = user_start;
2488                                user_info[user_curr].op_end = rec->num_op_rec() + 1;
2489                        }
2490                        // user_index, user_id, user_n, user_m
2491                        if( user_info[user_curr].connect_type != not_connected )
2492                        {       rec->PutArg(arg[0], arg[1], arg[2], arg[3]);
2493                                rec->PutOp(UserOp);
2494                        }
2495                        break;
2496
2497                        case UsrapOp:
2498                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
2499                        if( user_info[user_curr].connect_type != not_connected )
2500                        {       new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
2501                                rec->PutArg(new_arg[0]);
2502                                rec->PutOp(UsrapOp);
2503                        }
2504                        break;
2505
2506                        case UsravOp:
2507                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
2508                        if( user_info[user_curr].connect_type != not_connected )
2509                        {       new_arg[0] = tape[arg[0]].new_var;
2510                                if( size_t(new_arg[0]) < num_var )
2511                                {       rec->PutArg(new_arg[0]);
2512                                        rec->PutOp(UsravOp);
2513                                }
2514                                else
2515                                {       // This argument does not affect the result and
2516                                        // has been optimized out so use nan in its place.
2517                                        new_arg[0] = rec->PutPar( nan(Base(0)) );
2518                                        rec->PutArg(new_arg[0]);
2519                                        rec->PutOp(UsrapOp);
2520                                }
2521                        }
2522                        break;
2523
2524                        case UsrrpOp:
2525                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
2526                        if( user_info[user_curr].connect_type != not_connected )
2527                        {       new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
2528                                rec->PutArg(new_arg[0]);
2529                                rec->PutOp(UsrrpOp);
2530                        }
2531                        break;
2532                       
2533                        case UsrrvOp:
2534                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
2535                        if( user_info[user_curr].connect_type != not_connected )
2536                        {       tape[i_var].new_op  = rec->num_op_rec();
2537                                tape[i_var].new_var = rec->PutOp(UsrrvOp);
2538                        }
2539                        break;
2540                        // ---------------------------------------------------
2541
2542                        // all cases should be handled above
2543                        default:
2544                        CPPAD_ASSERT_UNKNOWN(false);
2545
2546                }
2547                if( replace_hash )
2548                {       // The old variable index i_var corresponds to the
2549                        // new variable index tape[i_var].new_var. In addition
2550                        // this is the most recent variable that has this code.
2551                        hash_table_var[code] = i_var;
2552                }
2553
2554        }
2555        // modify the dependent variable vector to new indices
2556        for(i = 0; i < dep_taddr.size(); i++ )
2557        {       CPPAD_ASSERT_UNKNOWN( size_t(tape[dep_taddr[i]].new_var) < num_var );
2558                dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
2559        }
2560
2561# ifndef NDEBUG
2562        size_t num_new_op = rec->num_op_rec();
2563        for(i_var = 0; i_var < tape.size(); i_var++)
2564                CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op < num_new_op );
2565# endif
2566
2567        // Move skip information from user_info to cskip_info
2568        for(i = 0; i < user_info.size(); i++)
2569        {       if( user_info[i].connect_type == cexp_connected )
2570                {       std::set<class_cexp_pair>::const_iterator itr =
2571                                user_info[i].cexp_set.begin();
2572                        while( itr != user_info[i].cexp_set.end() )
2573                        {       j = itr->index_;
2574                                k = user_info[i].op_begin;
2575                                while(k < user_info[i].op_end)
2576                                {       if( itr->compare_ == true )
2577                                                cskip_info[j].skip_op_false.push_back(k++);
2578                                        else    cskip_info[j].skip_op_true.push_back(k++);
2579                                }
2580                                itr++;
2581                        }
2582                }
2583        }
2584
2585        // fill in the arguments for the CSkip operations
2586        CPPAD_ASSERT_UNKNOWN( cskip_info_next == cskip_info.size() );
2587        for(i = 0; i < cskip_info.size(); i++)
2588        {       struct_cskip_info info = cskip_info[i];
2589                if( info.i_arg > 0 )
2590                {       CPPAD_ASSERT_UNKNOWN( info.n_op_true==info.skip_op_true.size() );
2591                        CPPAD_ASSERT_UNKNOWN(info.n_op_false==info.skip_op_false.size());
2592                        size_t n_true  = 
2593                                info.skip_var_true.size() + info.skip_op_true.size();
2594                        size_t n_false = 
2595                                info.skip_var_false.size() + info.skip_op_false.size();
2596                        size_t i_arg   = info.i_arg;
2597                        rec->ReplaceArg(i_arg++, info.cop   );
2598                        rec->ReplaceArg(i_arg++, info.flag  );
2599                        rec->ReplaceArg(i_arg++, info.left  );
2600                        rec->ReplaceArg(i_arg++, info.right );
2601                        rec->ReplaceArg(i_arg++, n_true     );
2602                        rec->ReplaceArg(i_arg++, n_false    );
2603                        for(j = 0; j < info.skip_var_true.size(); j++)
2604                        {       i_var = info.skip_var_true[j];
2605                                CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op > 0 );
2606                                rec->ReplaceArg(i_arg++, tape[i_var].new_op );
2607                        } 
2608                        for(j = 0; j < info.skip_op_true.size(); j++)
2609                        {       i_op = info.skip_op_true[j];
2610                                rec->ReplaceArg(i_arg++, i_op);
2611                        } 
2612                        for(j = 0; j < info.skip_var_false.size(); j++)
2613                        {       i_var = info.skip_var_false[j];
2614                                CPPAD_ASSERT_UNKNOWN( tape[i_var].new_op > 0 );
2615                                rec->ReplaceArg(i_arg++, tape[i_var].new_op );
2616                        } 
2617                        for(j = 0; j < info.skip_op_false.size(); j++)
2618                        {       i_op = info.skip_op_false[j];
2619                                rec->ReplaceArg(i_arg++, i_op);
2620                        } 
2621                        rec->ReplaceArg(i_arg++, n_true + n_false);
2622# ifndef NDEBUG
2623                        size_t n_arg   = 7 + n_true + n_false; 
2624                        CPPAD_ASSERT_UNKNOWN( info.i_arg + n_arg == i_arg );
2625# endif
2626                }
2627        }
2628}
2629
2630} // END_CPPAD_OPTIMIZE_NAMESPACE
2631
2632/*!
2633Optimize a player object operation sequence
2634
2635The operation sequence for this object is replaced by one with fewer operations
2636but the same funcition and derivative values.
2637
2638\tparam Base
2639base type for the operator; i.e., this operation was recorded
2640using AD< \a Base > and computations by this routine are done using type
2641\a Base.
2642
2643*/
2644template <class Base>
2645void ADFun<Base>::optimize(void)
2646{       // place to store the optimized version of the recording
2647        recorder<Base> rec;
2648
2649        // number of independent variables
2650        size_t n = ind_taddr_.size();
2651
2652# ifndef NDEBUG
2653        size_t i, j, m = dep_taddr_.size();
2654        CppAD::vector<Base> x(n), y(m), check(m);
2655        Base max_taylor(0);
2656        bool check_zero_order = num_order_taylor_ > 0;
2657        if( check_zero_order )
2658        {       // zero order coefficients for independent vars
2659                for(j = 0; j < n; j++)
2660                {       CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
2661                        CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
2662                        x[j] = taylor_[ ind_taddr_[j] * cap_order_taylor_ + 0];
2663                }
2664                // zero order coefficients for dependent vars
2665                for(i = 0; i < m; i++)
2666                {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_  );
2667                        y[i] = taylor_[ dep_taddr_[i] * cap_order_taylor_ + 0];
2668                }
2669                // maximum zero order coefficient not counting BeginOp at beginning
2670                // (which is correpsonds to uninitialized memory).
2671                for(i = 1; i < num_var_tape_; i++)
2672                {       if(  abs_geq(taylor_[i*cap_order_taylor_+0] , max_taylor) )
2673                                max_taylor = taylor_[i*cap_order_taylor_+0];
2674                }
2675        }
2676# endif
2677
2678        // create the optimized recording
2679        CppAD::optimize::optimize_run<Base>(n, dep_taddr_, &play_, &rec);
2680
2681        // number of variables in the recording
2682        num_var_tape_  = rec.num_var_rec();
2683
2684        // now replace the recording
2685        play_.get(rec);
2686
2687        // free memory allocated for sparse Jacobian calculation
2688        // (the results are no longer valid)
2689        for_jac_sparse_pack_.resize(0, 0);
2690        for_jac_sparse_set_.resize(0,0);
2691
2692        // free old Taylor coefficient memory
2693        taylor_.free();
2694        num_order_taylor_     = 0;
2695        cap_order_taylor_     = 0;
2696
2697        // resize and initilaize conditional skip vector
2698        // (must use player size because it now has the recoreder information)
2699        cskip_op_.erase();
2700        cskip_op_.extend( play_.num_op_rec() );
2701
2702# ifndef NDEBUG
2703        if( check_zero_order )
2704        {
2705                // zero order forward calculation using new operation sequence
2706                check = Forward(0, x);
2707
2708                // check results
2709                Base eps = 10. * epsilon<Base>();
2710                for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
2711                        abs_geq( eps * max_taylor , check[i] - y[i] ) ,
2712                        "Error during check of f.optimize()."
2713                );
2714
2715                // Erase memory that this calculation was done so NDEBUG gives
2716                // same final state for this object (from users perspective)
2717                num_order_taylor_     = 0;
2718        }
2719# endif
2720}
2721
2722} // END_CPPAD_NAMESPACE
2723# endif
Note: See TracBrowser for help on using the repository browser.