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

Last change on this file since 2979 was 2979, checked in by bradbell, 7 years ago
  1. Add operator index to trace output.
  2. Check that new operator values have been set and are valid.
  3. Reinitialize cskip_op_ after optimization.

optimize.hpp: Add test of entirely removing an atomic call.

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