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

Last change on this file since 3507 was 3507, checked in by bradbell, 5 years ago

Must use old variable index to mark variables used as replacements for
other variables (that have been removed).

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