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

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

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 071875a4beba3363e5fa9752426aec4762cd1caa
end hash code: 0bef506513a519e1073c6279d5c4cba9e5c3b180

commit 0bef506513a519e1073c6279d5c4cba9e5c3b180
Author: Brad Bell <bradbell@…>
Date: Thu May 7 12:14:32 2015 -0700

Add the acosh function (as an atomic operation when defined by compiler).

commit b3264fa17b2f65b65800423a0e243c9c3ccfe06a
Author: Brad Bell <bradbell@…>
Date: Wed May 6 20:25:38 2015 -0700

CMakeLists.txt: Change so test only check for compliation.

commit dcbac4d4f20cc383f2bd9edb02036659df40b791
Author: Brad Bell <bradbell@…>
Date: Wed May 6 15:06:28 2015 -0700

asinh.cpp: check higher orders, relax accuracy on test.

commit 5f8881993fedd18cccc3c74831133a8f8a9d17b0
Author: Brad Bell <bradbell@…>
Date: Wed May 6 14:36:18 2015 -0700

Change Acos to acos.
acos.cpp: remove trailing white space.

commit e828fa1f7c4c3848c727f14b1b7a8030071ee705
Author: Brad Bell <bradbell@…>
Date: Wed May 6 12:07:35 2015 -0700

Change Acos to acos.
acos.cpp: remove redundant index commands, remove trailing with space.

commit 3d16e5b9fe1bdafa4ad01d1d466bb72b792650fa
Author: Brad Bell <bradbell@…>
Date: Wed May 6 11:30:49 2015 -0700

op_code.hpp: Minor edits to AcosOp? commnets.

commit 58beaaad149b4ac29fae44589d7f8900bf8f4c40
Author: Brad Bell <bradbell@…>
Date: Wed May 6 10:51:43 2015 -0700

for_jac_sweep.hpp: Add missing AsinhOp? case.

commit 623c134870c522ae5e80bcf0f89d230902594c80
Author: Brad Bell <bradbell@…>
Date: Wed May 6 10:27:39 2015 -0700

Fix comment about AsinhOp? operator.

commit 226b14f6f4810f5abf1ca247aae541963efaf4d6
Author: Brad Bell <bradbell@…>
Date: Wed May 6 10:14:08 2015 -0700

Add derivative of F to make order zero case clearer.
acos_reverse.omh: Fix some sign errors.
asin_reverse.omh: Fix typo.
acos_forward.omh: Simplify by distributing minus sign.

commit 4682f4ee73e33b600b180086576e986f636a24dc
Author: Brad Bell <bradbell@…>
Date: Wed May 6 08:15:50 2015 -0700

acos_forward.omh: fix sign that depends on acos versus acosh.

commit 906ae10adf019ddda7f57dd165aab08fc55289c4
Author: Brad Bell <bradbell@…>
Date: Wed May 6 07:09:47 2015 -0700

  1. Fix inclusion of some temporary files in package (e.g., git_commit.sh).
  2. Simplify and improve using git ls-files and ls bin/check_*.
  3. Remove trailing white space.

commit 5096f4706a547bd76caa3766aa2c62802ef7f0bf
Author: Brad Bell <bradbell@…>
Date: Wed May 6 06:41:20 2015 -0700

Combine base type documentation for erf, asinh
(will add more functions to this list list).

commit b3535db5ad95bee90672abcaa686032d23bce2fc
Author: Brad Bell <bradbell@…>
Date: Tue May 5 18:01:11 2015 -0700

  1. Change Arc Cosine/Sine? to Inverse Cosine/Sine?.
  2. Change arcsin-> asin and arccos->acos.
  3. Remove index commands that are duplicates of words in titles.


acos_reverse.omh: Add acosh case to this page.

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