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

Last change on this file since 2954 was 2954, checked in by bradbell, 7 years ago

optimize.hpp: prepare for having a conneciton index and type.

  • Property svn:keywords set to Id
File size: 60.4 KB
Line 
1/* $Id: optimize.hpp 2954 2013-10-16 04:32:12Z 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 operation that connects this variable to
158        /// the independent variables and it is one of the following:
159        ///  AddvvOp, AddpvOp, SubpvOp, SubvpOp, or SubvvOp.
160        sum_connected      ,
161
162        /// Satisfies sum_connected and in addition
163        /// the parrent operation is one of the following:
164        ///  AddvvOp, AddpvOp, SubpvOp, SubvpOp, or SubvvOp.
165        csum_connected
166
167};
168
169
170/*!
171Structure used by \c optimize to hold information about one variable.
172in the old operation seqeunce.
173*/
174struct optimize_old_variable {
175        /// Operator for which this variable is the result, \c NumRes(op) > 0.
176        /// Set by the reverse sweep at beginning of optimization.
177        OpCode              op;       
178
179        /// Pointer to first argument (child) for this operator.
180        /// Set by the reverse sweep at beginning of optimization.
181        const addr_t*       arg;
182
183        /// How is this variable connected to the independent variables
184        optimize_connection_type connect_type; 
185
186        /// Set during forward sweep to the index in the
187        /// new operation sequence corresponding to this old varable.
188        addr_t new_var;
189};
190
191/*!
192Structures used by \c optimize_record_csum
193to hold information about one variable.
194*/
195struct optimize_csum_variable {
196        /// Operator for which this variable is the result, \c NumRes(op) > 0.
197        OpCode              op;       
198
199        /// Pointer to first argument (child) for this operator.
200        /// Set by the reverse sweep at beginning of optimization.
201        const addr_t*       arg;
202
203        /// Is this variable added to the summation
204        /// (if not it is subtracted)
205        bool                add;
206};
207
208/*!
209Structure used to pass work space from \c optimize to \c optimize_record_csum
210(so that stacks do not start from zero size every time).
211*/
212struct optimize_csum_stacks {
213        /// stack of operations in the cummulative summation
214        std::stack<struct optimize_csum_variable>   op_stack;
215        /// stack of variables to be added
216        std::stack<size_t >                         add_stack;
217        /// stack of variables to be subtracted
218        std::stack<size_t >                         sub_stack;
219};
220
221
222/*!
223Shared documentation for optimization helper functions (not called).
224
225<!-- define optimize_prototype -->
226\param tape
227is a vector that maps a variable index, in the old operation sequence,
228to an <tt>optimize_old_variable</tt> information record.
229Note that the index for this vector must be greater than or equal zero and
230less than <tt>tape.size()</tt>.
231
232\li <tt>tape[i].op</tt>
233is the operator in the old operation sequence
234corresponding to the old variable index \c i.
235Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
236
237\li <tt>tape[i].arg</tt>
238for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
239is the j-th the argument, in the old operation sequence,
240corresponding to the old variable index \c i.
241Assertion: <tt>tape[i].arg[j] < i</tt>.
242
243\li <tt>tape[i].new_var</tt>
244Suppose
245<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
246and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
247It follows that <tt>tape[k].new_var</tt>
248has alread been set to the variable in the new operation sequence
249corresponding to the old variable index \c k.
250This means that the \c new_var value has been set
251for all the possible arguments that come before \a current.
252
253\param current
254is the index in the old operation sequence for
255the variable corresponding to the result for the current operator.
256Assertions:
257<tt>
258current < tape.size(),
259NumRes( tape[current].op ) > 0.
260</tt>
261
262\param npar
263is the number of parameters corresponding to this operation sequence.
264
265\param par
266is a vector of length \a npar containing the parameters
267for this operation sequence; i.e.,
268given a parameter index \c i, the corresponding parameter value is
269<tt>par[i]</tt>.
270<!-- end optimize_prototype -->
271*/
272template <class Base>
273void optimize_prototype(
274        const CppAD::vector<struct optimize_old_variable>& tape           ,
275        size_t                                             current        ,
276        size_t                                             npar           ,
277        const Base*                                        par            )
278{       CPPAD_ASSERT_UNKNOWN(false); }
279
280/*!
281Check a unary operator for a complete match with a previous operator.
282
283A complete match means that the result of the previous operator
284can be used inplace of the result for current operator.
285
286<!-- replace optimize_prototype -->
287\param tape
288is a vector that maps a variable index, in the old operation sequence,
289to an <tt>optimize_old_variable</tt> information record.
290Note that the index for this vector must be greater than or equal zero and
291less than <tt>tape.size()</tt>.
292
293\li <tt>tape[i].op</tt>
294is the operator in the old operation sequence
295corresponding to the old variable index \c i.
296Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
297
298\li <tt>tape[i].arg</tt>
299for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
300is the j-th the argument, in the old operation sequence,
301corresponding to the old variable index \c i.
302Assertion: <tt>tape[i].arg[j] < i</tt>.
303
304\li <tt>tape[i].new_var</tt>
305Suppose
306<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
307and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
308It follows that <tt>tape[k].new_var</tt>
309has alread been set to the variable in the new operation sequence
310corresponding to the old variable index \c k.
311This means that the \c new_var value has been set
312for all the possible arguments that come before \a current.
313
314\param current
315is the index in the old operation sequence for
316the variable corresponding to the result for the current operator.
317Assertions:
318<tt>
319current < tape.size(),
320NumRes( tape[current].op ) > 0.
321</tt>
322
323\param npar
324is the number of parameters corresponding to this operation sequence.
325
326\param par
327is a vector of length \a npar containing the parameters
328for this operation sequence; i.e.,
329given a parameter index \c i, the corresponding parameter value is
330<tt>par[i]</tt>.
331<!-- end optimize_prototype -->
332
333\param hash_table_var
334is a vector with size CPPAD_HASH_TABLE_SIZE
335that maps a hash code to the corresponding
336variable index in the old operation sequence.
337All the values in this table must be less than \a current.
338
339\param code
340The input value of code does not matter.
341The output value of code is the hash code corresponding to
342this operation in the new operation sequence.
343
344\return
345If the return value is zero,
346no match was found.
347If the return value is greater than zero,
348it is the index of a new variable that can be used to replace the
349old variable.
350
351\par Restrictions:
352NumArg( tape[current].op ) == 1
353*/
354template <class Base>
355size_t optimize_unary_match(
356        const CppAD::vector<struct optimize_old_variable>& tape           ,
357        size_t                                             current        ,
358        size_t                                             npar           ,
359        const Base*                                        par            ,
360        const CppAD::vector<size_t>&                       hash_table_var ,
361        unsigned short&                                    code           )
362{       const addr_t* arg = tape[current].arg;
363        OpCode        op  = tape[current].op;
364        addr_t new_arg[1];
365       
366        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
367        CPPAD_ASSERT_UNKNOWN( NumRes(op) > 0  );
368        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
369        new_arg[0] = tape[arg[0]].new_var;
370        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < current );
371        code = hash_code(
372                op                  , 
373                new_arg             ,
374                npar                ,
375                par
376        );
377        size_t  i               = hash_table_var[code];
378        CPPAD_ASSERT_UNKNOWN( i < current );
379        if( op == tape[i].op )
380        {       size_t k = tape[i].arg[0];
381                CPPAD_ASSERT_UNKNOWN( k < i );
382                if (new_arg[0] == tape[k].new_var )
383                        return tape[i].new_var;
384        }
385        return 0;
386} 
387
388/*!
389Check a binary operator for a complete match with a previous operator,
390
391<!-- replace optimize_prototype -->
392\param tape
393is a vector that maps a variable index, in the old operation sequence,
394to an <tt>optimize_old_variable</tt> information record.
395Note that the index for this vector must be greater than or equal zero and
396less than <tt>tape.size()</tt>.
397
398\li <tt>tape[i].op</tt>
399is the operator in the old operation sequence
400corresponding to the old variable index \c i.
401Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
402
403\li <tt>tape[i].arg</tt>
404for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
405is the j-th the argument, in the old operation sequence,
406corresponding to the old variable index \c i.
407Assertion: <tt>tape[i].arg[j] < i</tt>.
408
409\li <tt>tape[i].new_var</tt>
410Suppose
411<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
412and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
413It follows that <tt>tape[k].new_var</tt>
414has alread been set to the variable in the new operation sequence
415corresponding to the old variable index \c k.
416This means that the \c new_var value has been set
417for all the possible arguments that come before \a current.
418
419\param current
420is the index in the old operation sequence for
421the variable corresponding to the result for the current operator.
422Assertions:
423<tt>
424current < tape.size(),
425NumRes( tape[current].op ) > 0.
426</tt>
427
428\param npar
429is the number of parameters corresponding to this operation sequence.
430
431\param par
432is a vector of length \a npar containing the parameters
433for this operation sequence; i.e.,
434given a parameter index \c i, the corresponding parameter value is
435<tt>par[i]</tt>.
436<!-- end optimize_prototype -->
437
438\param hash_table_var
439is a vector with size CPPAD_HASH_TABLE_SIZE
440that maps a hash code to the corresponding
441variable index in the old operation sequence.
442All the values in this table must be less than \a current.
443
444\param code
445The input value of code does not matter.
446The output value of code is the hash code corresponding to
447this operation in the new operation sequence.
448
449\return
450If the return value is zero,
451no match was found.
452If the return value is greater than zero,
453it is the index of a new variable that can be used to replace the
454old variable.
455
456
457\par Restrictions:
458The binary operator must be an addition, subtraction, multiplication, division
459or power operator.  NumArg( tape[current].op ) == 1.
460*/
461template <class Base>
462inline size_t optimize_binary_match(
463        const CppAD::vector<struct optimize_old_variable>& tape           ,
464        size_t                                             current        ,
465        size_t                                             npar           ,
466        const Base*                                        par            ,
467        const CppAD::vector<size_t>&                       hash_table_var ,
468        unsigned short&                                    code           )
469{       OpCode        op         = tape[current].op;
470        const addr_t* arg        = tape[current].arg;
471        addr_t        new_arg[2];
472        bool          parameter[2];
473
474        // initialize return value
475        size_t  match_var = 0;
476
477        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
478        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
479        switch(op)
480        {       // parameter op variable ----------------------------------
481                case AddpvOp:
482                case MulpvOp:
483                case DivpvOp:
484                case PowpvOp:
485                case SubpvOp:
486                // arg[0]
487                parameter[0] = true;
488                new_arg[0]   = arg[0];
489                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar );
490                // arg[1]
491                parameter[1] = false;
492                new_arg[1]   = tape[arg[1]].new_var;
493                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
494                break;
495
496                // variable op parameter -----------------------------------
497                case DivvpOp:
498                case PowvpOp:
499                case SubvpOp:
500                // arg[0]
501                parameter[0] = false;
502                new_arg[0]   = tape[arg[0]].new_var;
503                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
504                // arg[1]
505                parameter[1] = true;
506                new_arg[1]   = arg[1];
507                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar );
508                break;
509
510                // variable op variable -----------------------------------
511                case AddvvOp:
512                case MulvvOp:
513                case DivvvOp:
514                case PowvvOp:
515                case SubvvOp:
516                // arg[0]
517                parameter[0] = false;
518                new_arg[0]   = tape[arg[0]].new_var;
519                CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
520                // arg[1]
521                parameter[1] = false;
522                new_arg[1]   = tape[arg[1]].new_var;
523                CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
524                break;
525
526                // must be one of the cases above
527                default:
528                CPPAD_ASSERT_UNKNOWN(false);
529        }
530        code = hash_code(
531                op                  , 
532                new_arg             ,
533                npar                ,
534                par
535        );
536        size_t  i  = hash_table_var[code];
537        CPPAD_ASSERT_UNKNOWN( i < current );
538        if( op == tape[i].op )
539        {       bool match = true;
540                size_t j;
541                for(j = 0; j < 2; j++)
542                {       size_t k = tape[i].arg[j];
543                        if( parameter[j] )
544                        {       CPPAD_ASSERT_UNKNOWN( k < npar );
545                                match &= IdenticalEqualPar(
546                                        par[ arg[j] ], par[k]
547                                );
548                        }
549                        else
550                        {       CPPAD_ASSERT_UNKNOWN( k < i );
551                                match &= (new_arg[j] == tape[k].new_var);
552                        }
553                }
554                if( match )
555                        match_var = tape[i].new_var;
556        }
557        if( (match_var > 0) | ( (op != AddvvOp) & (op != MulvvOp ) ) )
558                return match_var;
559
560        // check for match with argument order switched ----------------------
561        CPPAD_ASSERT_UNKNOWN( op == AddvvOp || op == MulvvOp );
562        i          = new_arg[0];
563        new_arg[0] = new_arg[1];
564        new_arg[1] = i;
565        unsigned short code_switch = hash_code(
566                op                  , 
567                new_arg             ,
568                npar                ,
569                par
570        );
571        i  = hash_table_var[code_switch];
572        CPPAD_ASSERT_UNKNOWN( i < current );
573        if( op == tape[i].op )
574        {       bool match = true;
575                size_t j;
576                for(j = 0; j < 2; j++)
577                {       size_t k = tape[i].arg[j];
578                        CPPAD_ASSERT_UNKNOWN( k < i );
579                        match &= (new_arg[j] == tape[k].new_var);
580                }
581                if( match )
582                        match_var = tape[i].new_var;
583        }
584        return match_var;
585} 
586
587/*!
588Record an operation of the form (parameter op variable).
589
590<!-- replace optimize_prototype -->
591\param tape
592is a vector that maps a variable index, in the old operation sequence,
593to an <tt>optimize_old_variable</tt> information record.
594Note that the index for this vector must be greater than or equal zero and
595less than <tt>tape.size()</tt>.
596
597\li <tt>tape[i].op</tt>
598is the operator in the old operation sequence
599corresponding to the old variable index \c i.
600Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
601
602\li <tt>tape[i].arg</tt>
603for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
604is the j-th the argument, in the old operation sequence,
605corresponding to the old variable index \c i.
606Assertion: <tt>tape[i].arg[j] < i</tt>.
607
608\li <tt>tape[i].new_var</tt>
609Suppose
610<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
611and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
612It follows that <tt>tape[k].new_var</tt>
613has alread been set to the variable in the new operation sequence
614corresponding to the old variable index \c k.
615This means that the \c new_var value has been set
616for all the possible arguments that come before \a current.
617
618\param current
619is the index in the old operation sequence for
620the variable corresponding to the result for the current operator.
621Assertions:
622<tt>
623current < tape.size(),
624NumRes( tape[current].op ) > 0.
625</tt>
626
627\param npar
628is the number of parameters corresponding to this operation sequence.
629
630\param par
631is a vector of length \a npar containing the parameters
632for this operation sequence; i.e.,
633given a parameter index \c i, the corresponding parameter value is
634<tt>par[i]</tt>.
635<!-- end optimize_prototype -->
636
637\param rec
638is the object that will record the operations.
639
640\param op
641is the operator that we are recording which must be one of the following:
642AddpvOp, DivpvOp, MulpvOp, PowvpOp, SubpvOp.
643 
644\param arg
645is the vector of arguments for this operator.
646
647\return
648the result value is the index corresponding to the current
649operation in the new operation sequence.
650*/
651template <class Base>
652size_t optimize_record_pv(
653        const CppAD::vector<struct optimize_old_variable>& tape           ,
654        size_t                                             current        ,
655        size_t                                             npar           ,
656        const Base*                                        par            ,
657        recorder<Base>*                                    rec            ,
658        OpCode                                             op             ,
659        const addr_t*                                      arg            )
660{
661# ifndef NDEBUG
662        switch(op)
663        {       case AddpvOp:
664                case DivpvOp:
665                case MulpvOp:
666                case PowpvOp:
667                case SubpvOp:
668                break;
669
670                default:
671                CPPAD_ASSERT_UNKNOWN(false);
672        }
673# endif
674        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar    );
675        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
676        addr_t new_arg[2];
677        new_arg[0]   = rec->PutPar( par[arg[0]] );
678        new_arg[1]   = tape[ arg[1] ].new_var;
679        rec->PutArg( new_arg[0], new_arg[1] );
680        size_t i     = rec->PutOp(op);
681        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
682        return i;
683}
684
685
686/*!
687Record an operation of the form (variable op parameter).
688
689<!-- replace optimize_prototype -->
690\param tape
691is a vector that maps a variable index, in the old operation sequence,
692to an <tt>optimize_old_variable</tt> information record.
693Note that the index for this vector must be greater than or equal zero and
694less than <tt>tape.size()</tt>.
695
696\li <tt>tape[i].op</tt>
697is the operator in the old operation sequence
698corresponding to the old variable index \c i.
699Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
700
701\li <tt>tape[i].arg</tt>
702for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
703is the j-th the argument, in the old operation sequence,
704corresponding to the old variable index \c i.
705Assertion: <tt>tape[i].arg[j] < i</tt>.
706
707\li <tt>tape[i].new_var</tt>
708Suppose
709<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
710and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
711It follows that <tt>tape[k].new_var</tt>
712has alread been set to the variable in the new operation sequence
713corresponding to the old variable index \c k.
714This means that the \c new_var value has been set
715for all the possible arguments that come before \a current.
716
717\param current
718is the index in the old operation sequence for
719the variable corresponding to the result for the current operator.
720Assertions:
721<tt>
722current < tape.size(),
723NumRes( tape[current].op ) > 0.
724</tt>
725
726\param npar
727is the number of parameters corresponding to this operation sequence.
728
729\param par
730is a vector of length \a npar containing the parameters
731for this operation sequence; i.e.,
732given a parameter index \c i, the corresponding parameter value is
733<tt>par[i]</tt>.
734<!-- end optimize_prototype -->
735
736\param rec
737is the object that will record the operations.
738
739\param op
740is the operator that we are recording which must be one of the following:
741DivvpOp, PowvpOp, SubvpOp.
742 
743\param arg
744is the vector of arguments for this operator.
745
746\return
747the result value is the index corresponding to the current
748operation in the new operation sequence.
749*/
750template <class Base>
751size_t optimize_record_vp(
752        const CppAD::vector<struct optimize_old_variable>& tape           ,
753        size_t                                             current        ,
754        size_t                                             npar           ,
755        const Base*                                        par            ,
756        recorder<Base>*                                    rec            ,
757        OpCode                                             op             ,
758        const addr_t*                                      arg            )
759{
760# ifndef NDEBUG
761        switch(op)
762        {       case DivvpOp:
763                case PowvpOp:
764                case SubvpOp:
765                break;
766
767                default:
768                CPPAD_ASSERT_UNKNOWN(false);
769        }
770# endif
771        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
772        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar    );
773        addr_t new_arg[2];
774        new_arg[0]   = tape[ arg[0] ].new_var;
775        new_arg[1]   = rec->PutPar( par[arg[1]] );
776        rec->PutArg( new_arg[0], new_arg[1] );
777        size_t i     = rec->PutOp(op);
778        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
779        return i;
780}
781
782/*!
783Record an operation of the form (variable op variable).
784
785<!-- replace optimize_prototype -->
786\param tape
787is a vector that maps a variable index, in the old operation sequence,
788to an <tt>optimize_old_variable</tt> information record.
789Note that the index for this vector must be greater than or equal zero and
790less than <tt>tape.size()</tt>.
791
792\li <tt>tape[i].op</tt>
793is the operator in the old operation sequence
794corresponding to the old variable index \c i.
795Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
796
797\li <tt>tape[i].arg</tt>
798for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
799is the j-th the argument, in the old operation sequence,
800corresponding to the old variable index \c i.
801Assertion: <tt>tape[i].arg[j] < i</tt>.
802
803\li <tt>tape[i].new_var</tt>
804Suppose
805<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
806and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
807It follows that <tt>tape[k].new_var</tt>
808has alread been set to the variable in the new operation sequence
809corresponding to the old variable index \c k.
810This means that the \c new_var value has been set
811for all the possible arguments that come before \a current.
812
813\param current
814is the index in the old operation sequence for
815the variable corresponding to the result for the current operator.
816Assertions:
817<tt>
818current < tape.size(),
819NumRes( tape[current].op ) > 0.
820</tt>
821
822\param npar
823is the number of parameters corresponding to this operation sequence.
824
825\param par
826is a vector of length \a npar containing the parameters
827for this operation sequence; i.e.,
828given a parameter index \c i, the corresponding parameter value is
829<tt>par[i]</tt>.
830<!-- end optimize_prototype -->
831
832\param rec
833is the object that will record the operations.
834
835\param op
836is the operator that we are recording which must be one of the following:
837AddvvOp, DivvvOp, MulvvOp, PowvpOp, SubvvOp.
838 
839\param arg
840is the vector of arguments for this operator.
841
842\return
843the result value is the index corresponding to the current
844operation in the new operation sequence.
845*/
846template <class Base>
847size_t optimize_record_vv(
848        const CppAD::vector<struct optimize_old_variable>& tape           ,
849        size_t                                             current        ,
850        size_t                                             npar           ,
851        const Base*                                        par            ,
852        recorder<Base>*                                    rec            ,
853        OpCode                                             op             ,
854        const addr_t*                                      arg            )
855{
856# ifndef NDEBUG
857        switch(op)
858        {       case AddvvOp:
859                case DivvvOp:
860                case MulvvOp:
861                case PowvvOp:
862                case SubvvOp:
863                break;
864
865                default:
866                CPPAD_ASSERT_UNKNOWN(false);
867        }
868# endif
869        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < current );
870        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < current );
871        addr_t new_arg[2];
872        new_arg[0]   = tape[ arg[0] ].new_var;
873        new_arg[1]   = tape[ arg[1] ].new_var;
874        rec->PutArg( new_arg[0], new_arg[1] );
875        size_t i     = rec->PutOp(op);
876        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
877        return i;
878}
879
880// ==========================================================================
881
882/*!
883Recording a cummulative cummulative summation starting at its highest parrent.
884
885<!-- replace optimize_prototype -->
886\param tape
887is a vector that maps a variable index, in the old operation sequence,
888to an <tt>optimize_old_variable</tt> information record.
889Note that the index for this vector must be greater than or equal zero and
890less than <tt>tape.size()</tt>.
891
892\li <tt>tape[i].op</tt>
893is the operator in the old operation sequence
894corresponding to the old variable index \c i.
895Assertion: <tt>NumRes(tape[i].op) > 0</tt>.
896
897\li <tt>tape[i].arg</tt>
898for <tt>j < NumArg( tape[i].op ), tape[i].arg[j]</tt>
899is the j-th the argument, in the old operation sequence,
900corresponding to the old variable index \c i.
901Assertion: <tt>tape[i].arg[j] < i</tt>.
902
903\li <tt>tape[i].new_var</tt>
904Suppose
905<tt>i <= current, j < NumArg( tape[i].op ), and k = tape[i].arg[j]</tt>,
906and \c j corresponds to a variable for operator <tt>tape[i].op</tt>.
907It follows that <tt>tape[k].new_var</tt>
908has alread been set to the variable in the new operation sequence
909corresponding to the old variable index \c k.
910This means that the \c new_var value has been set
911for all the possible arguments that come before \a current.
912
913\param current
914is the index in the old operation sequence for
915the variable corresponding to the result for the current operator.
916Assertions:
917<tt>
918current < tape.size(),
919NumRes( tape[current].op ) > 0.
920</tt>
921
922\param npar
923is the number of parameters corresponding to this operation sequence.
924
925\param par
926is a vector of length \a npar containing the parameters
927for this operation sequence; i.e.,
928given a parameter index \c i, the corresponding parameter value is
929<tt>par[i]</tt>.
930<!-- end optimize_prototype -->
931
932\param rec
933is the object that will record the operations.
934
935\param work
936Is used for computation. On input and output,
937<tt>work.op_stack.empty()</tt>,
938<tt>work.add_stack.empty()</tt>, and
939<tt>work.sub_stack.empty()</tt>,
940are all true true.
941These stacks are passed in so that elements can be allocated once
942and then the elements can be reused with calls to \c optimize_record_csum.
943
944\par Exception
945<tt>tape[i].new_var</tt>
946is not yet defined for any node \c i that is \c csum_connected
947to the \a current node.
948For example; suppose that index \c j corresponds to a variable
949in the current operator,
950<tt>i = tape[current].arg[j]</tt>,
951and
952<tt>tape[arg[j]].connect_type == csum_connected</tt>.
953It then follows that
954<tt>tape[i].new_var == tape.size()</tt>.
955
956\par Restriction:
957\li <tt>tape[current].op</tt>
958must be one of <tt>AddpvOp, AddvvOp, SubpvOp, SubvpOp, SubvvOp</tt>.
959
960\li <tt>tape[current].connect_type</tt> must be \c yes_connected.
961*/
962
963
964template <class Base>
965size_t optimize_record_csum(
966        const CppAD::vector<struct optimize_old_variable>& tape           ,
967        size_t                                             current        ,
968        size_t                                             npar           ,
969        const Base*                                        par            ,
970        recorder<Base>*                                    rec            ,
971        optimize_csum_stacks&                              work           )
972{
973       
974        CPPAD_ASSERT_UNKNOWN( work.op_stack.empty() );
975        CPPAD_ASSERT_UNKNOWN( work.add_stack.empty() );
976        CPPAD_ASSERT_UNKNOWN( work.sub_stack.empty() );
977        CPPAD_ASSERT_UNKNOWN( tape[current].connect_type == yes_connected );
978
979        size_t                        i;
980        OpCode                        op;
981        const addr_t*                 arg;
982        bool                          add;
983        struct optimize_csum_variable var;
984
985        var.op  = tape[current].op;
986        var.arg = tape[current].arg;
987        var.add = true; 
988        work.op_stack.push( var );
989        Base sum_par(0);
990        while( ! work.op_stack.empty() )
991        {       var     = work.op_stack.top();
992                work.op_stack.pop();
993                op      = var.op;
994                arg     = var.arg;
995                add     = var.add;
996                // process first argument to this operator
997                switch(op)
998                {       case AddpvOp:
999                        case SubpvOp:
1000                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < npar );
1001                        if( add )
1002                                sum_par += par[arg[0]];
1003                        else    sum_par -= par[arg[0]];
1004                        break;
1005
1006                        case AddvvOp:
1007                        case SubvpOp:
1008                        case SubvvOp:
1009                        if( tape[arg[0]].connect_type == csum_connected )
1010                        {       CPPAD_ASSERT_UNKNOWN(
1011                                        size_t(tape[arg[0]].new_var) == tape.size()
1012                                );
1013                                var.op  = tape[arg[0]].op;
1014                                var.arg = tape[arg[0]].arg;
1015                                var.add = add; 
1016                                work.op_stack.push( var );
1017                        }
1018                        else if( add )
1019                                work.add_stack.push(arg[0]);
1020                        else    work.sub_stack.push(arg[0]);
1021                        break;
1022
1023                        default:
1024                        CPPAD_ASSERT_UNKNOWN(false);
1025                }
1026                // process second argument to this operator
1027                switch(op)
1028                {
1029                        case SubvpOp:
1030                        CPPAD_ASSERT_UNKNOWN( size_t(arg[1]) < npar );
1031                        if( add )
1032                                sum_par -= par[arg[1]];
1033                        else    sum_par += par[arg[1]];
1034                        break;
1035
1036                        case SubvvOp:
1037                        case SubpvOp:
1038                        add = ! add;
1039
1040                        case AddvvOp:
1041                        case AddpvOp:
1042                        if( tape[arg[1]].connect_type == csum_connected )
1043                        {       CPPAD_ASSERT_UNKNOWN(
1044                                        size_t(tape[arg[1]].new_var) == tape.size()
1045                                );
1046                                var.op   = tape[arg[1]].op;
1047                                var.arg  = tape[arg[1]].arg;
1048                                var.add  = add;
1049                                work.op_stack.push( var );
1050                        }
1051                        else if( add )
1052                                work.add_stack.push(arg[1]);
1053                        else    work.sub_stack.push(arg[1]);
1054                        break;
1055
1056                        default:
1057                        CPPAD_ASSERT_UNKNOWN(false);
1058                }
1059        }
1060        // number of variables in this cummulative sum operator
1061        size_t n_add = work.add_stack.size();
1062        size_t n_sub = work.sub_stack.size();
1063        size_t old_arg, new_arg;
1064        rec->PutArg(n_add);                // arg[0]
1065        rec->PutArg(n_sub);                // arg[1]
1066        new_arg = rec->PutPar( sum_par );
1067        rec->PutArg(new_arg);              // arg[2]
1068        for(i = 0; i < n_add; i++)
1069        {       CPPAD_ASSERT_UNKNOWN( ! work.add_stack.empty() );
1070                old_arg = work.add_stack.top();
1071                new_arg = tape[old_arg].new_var;
1072                CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
1073                rec->PutArg(new_arg);      // arg[3+i]
1074                work.add_stack.pop();
1075        }
1076        for(i = 0; i < n_sub; i++)
1077        {       CPPAD_ASSERT_UNKNOWN( ! work.sub_stack.empty() );
1078                old_arg = work.sub_stack.top();
1079                new_arg = tape[old_arg].new_var;
1080                CPPAD_ASSERT_UNKNOWN( new_arg < tape.size() );
1081                rec->PutArg(new_arg);      // arg[3 + arg[0] + i]
1082                work.sub_stack.pop();
1083        }
1084        rec->PutArg(n_add + n_sub);        // arg[3 + arg[0] + arg[1]]
1085        i = rec->PutOp(CSumOp);
1086        CPPAD_ASSERT_UNKNOWN(new_arg < tape.size());
1087
1088        return i;
1089}
1090// ==========================================================================
1091/*!
1092Convert a player object to an optimized recorder object
1093
1094\tparam Base
1095base type for the operator; i.e., this operation was recorded
1096using AD< \a Base > and computations by this routine are done using type
1097\a Base.
1098
1099\param n
1100is the number of independent variables on the tape.
1101
1102\param dep_taddr
1103On input this vector contains the indices for each of the dependent
1104variable values in the operation sequence corresponding to \a play.
1105Upon return it contains the indices for the same variables but in
1106the operation sequence corresponding to \a rec.
1107
1108\param play
1109This is the operation sequence that we are optimizing.
1110It is essentially const, except for play back state which
1111changes while it plays back the operation seqeunce.
1112
1113\param rec
1114The input contents of this recording does not matter.
1115Upon return, it contains an optimized verison of the
1116operation sequence corresponding to \a play.
1117*/
1118
1119template <class Base>
1120void optimize(
1121        size_t                       n         ,
1122        CppAD::vector<size_t>&       dep_taddr ,
1123        player<Base>*                play      ,
1124        recorder<Base>*              rec       ) 
1125{
1126        // temporary indices
1127        size_t i, j, k;
1128
1129        // temporary variables
1130        OpCode        op;   // current operator
1131        const addr_t* arg;  // operator arguments
1132        size_t        i_var;  // index of first result for current operator
1133
1134        // range and domain dimensions for F
1135        size_t m = dep_taddr.size();
1136
1137        // number of variables in the player
1138        const size_t num_var = play->num_rec_var(); 
1139
1140# ifndef NDEBUG
1141        // number of parameters in the player
1142        const size_t num_par = play->num_rec_par();
1143# endif
1144
1145        // number of  VecAD indices
1146        size_t num_vecad_ind   = play->num_rec_vecad_ind();
1147
1148        // number of VecAD vectors
1149        size_t num_vecad_vec   = play->num_rec_vecad_vec();
1150
1151        // -------------------------------------------------------------
1152        // data structure that maps variable index in original operation
1153        // sequence to corresponding operator information
1154        CppAD::vector<struct optimize_old_variable> tape(num_var);
1155        // -------------------------------------------------------------
1156        // Determine how each variable is connected to the dependent variables
1157
1158        // initialize all variables has having no connections
1159        for(i = 0; i < num_var; i++)
1160                tape[i].connect_type = not_connected;
1161
1162        for(j = 0; j < m; j++)
1163        {       // mark dependent variables as having one or more connections
1164                tape[ dep_taddr[j] ].connect_type = yes_connected;
1165        }
1166
1167        // vecad_connect contains a value for each VecAD object.
1168        // vecad maps a VecAD index (which corresponds to the beginning of the
1169        // VecAD object) to the vecad_connect falg for the VecAD object.
1170        CppAD::vector<optimize_connection_type>   vecad_connect(num_vecad_vec);
1171        CppAD::vector<size_t> vecad(num_vecad_ind);
1172        j = 0;
1173        for(i = 0; i < num_vecad_vec; i++)
1174        {       vecad_connect[i] = not_connected;
1175                // length of this VecAD
1176                size_t length = play->GetVecInd(j);
1177                // set to proper index for this VecAD
1178                vecad[j] = i; 
1179                for(k = 1; k <= length; k++)
1180                        vecad[j+k] = num_vecad_vec; // invalid index
1181                // start of next VecAD
1182                j       += length + 1;
1183        }
1184        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
1185
1186        // work space used by UserOp.
1187        typedef std::set<size_t> size_set;
1188        size_t user_q     = 0;       // maximum set element plus one
1189        vector< size_set > user_r;   // sparsity pattern for the argument x
1190        vector< size_set > user_s;   // sparisty pattern for the result y
1191        size_t user_index = 0;       // indentifier for this user_atomic operation
1192        size_t user_id    = 0;       // user identifier for this call to operator
1193        size_t user_i     = 0;       // index in result vector
1194        size_t user_j     = 0;       // index in argument vector
1195        size_t user_m     = 0;       // size of result vector
1196        size_t user_n     = 0;       // size of arugment vector
1197        // next expected operator in a UserOp sequence
1198        enum { user_start, user_arg, user_ret, user_end } user_state;
1199
1200        // During reverse mode, push true if user operation is connected,
1201        // push false otherwise. During forward mode, use to determine if
1202        // we are keeping this operation and then pop.
1203        std::stack<bool> user_keep;
1204
1205        // Initialize a reverse mode sweep through the operation sequence
1206        size_t i_op;
1207        play->start_reverse(op, arg, i_op, i_var);
1208        CPPAD_ASSERT_UNKNOWN( op == EndOp );
1209        size_t mask;
1210        user_state = user_end;
1211        while(op != BeginOp)
1212        {       // next op
1213                play->next_reverse(op, arg, i_op, i_var);
1214                // This if is not necessary becasue last assignment
1215                // with this value of i_var will have NumRes(op) > 0
1216                if( NumRes(op) > 0 )
1217                {       tape[i_var].op = op;
1218                        tape[i_var].arg = arg;
1219                }
1220# ifndef NDEBUG
1221                if( i_op <= n )
1222                {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
1223                }
1224                else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
1225# endif
1226                switch( op )
1227                {
1228                        // Unary operator where operand is arg[0]
1229                        case AbsOp:
1230                        case AcosOp:
1231                        case AsinOp:
1232                        case AtanOp:
1233                        case CosOp:
1234                        case CoshOp:
1235                        case DivvpOp:
1236                        case ExpOp:
1237                        case LogOp:
1238                        case PowvpOp:
1239                        case SinOp:
1240                        case SinhOp:
1241                        case SqrtOp:
1242                        case TanOp:
1243                        case TanhOp:
1244                        if( tape[i_var].connect_type != not_connected )
1245                                tape[arg[0]].connect_type = yes_connected;
1246                        break; // --------------------------------------------
1247
1248                        // Unary operator where operand is arg[1]
1249                        case DisOp:
1250                        case DivpvOp:
1251                        case MulpvOp:
1252                        case PowpvOp:
1253                        if( tape[i_var].connect_type != not_connected )
1254                                tape[arg[1]].connect_type = yes_connected;
1255                        break; // --------------------------------------------
1256               
1257                        // Special case for SubvpOp
1258                        case SubvpOp:
1259                        if( tape[i_var].connect_type != not_connected )
1260                        {
1261                                if( tape[arg[0]].connect_type == not_connected )
1262                                        tape[arg[0]].connect_type = sum_connected;
1263                                else
1264                                        tape[arg[0]].connect_type = yes_connected;
1265                                if( tape[i_var].connect_type == sum_connected )
1266                                        tape[i_var].connect_type = csum_connected;
1267                        }
1268                        break; // --------------------------------------------
1269               
1270                        // Special case for AddpvOp and SubpvOp
1271                        case AddpvOp:
1272                        case SubpvOp:
1273                        if( tape[i_var].connect_type != not_connected )
1274                        {
1275                                if( tape[arg[1]].connect_type == not_connected )
1276                                        tape[arg[1]].connect_type = sum_connected;
1277                                else
1278                                        tape[arg[1]].connect_type = yes_connected;
1279                                if( tape[i_var].connect_type == sum_connected )
1280                                        tape[i_var].connect_type = csum_connected;
1281                        }
1282                        break; // --------------------------------------------
1283
1284               
1285                        // Special case for AddvvOp and SubvvOp
1286                        case AddvvOp:
1287                        case SubvvOp:
1288                        if( tape[i_var].connect_type != not_connected )
1289                        {
1290                                if( tape[arg[0]].connect_type == not_connected )
1291                                        tape[arg[0]].connect_type = sum_connected;
1292                                else
1293                                        tape[arg[0]].connect_type = yes_connected;
1294
1295                                if( tape[arg[1]].connect_type == not_connected )
1296                                        tape[arg[1]].connect_type = sum_connected;
1297                                else
1298                                        tape[arg[1]].connect_type = yes_connected;
1299                                if( tape[i_var].connect_type == sum_connected )
1300                                        tape[i_var].connect_type = csum_connected;
1301                        }
1302                        break; // --------------------------------------------
1303
1304                        // Other binary operators
1305                        // where operands are arg[0], arg[1]
1306                        case DivvvOp:
1307                        case MulvvOp:
1308                        case PowvvOp:
1309                        if( tape[i_var].connect_type != not_connected )
1310                        {
1311                                tape[arg[0]].connect_type = yes_connected;
1312                                tape[arg[1]].connect_type = yes_connected;
1313                        }
1314                        break; // --------------------------------------------
1315
1316                        // Conditional expression operators
1317                        case CExpOp:
1318                        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
1319                        if( tape[i_var].connect_type != not_connected )
1320                        {
1321                                mask = 1;
1322                                for(i = 2; i < 6; i++)
1323                                {       if( arg[1] & mask )
1324                                        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[i]) < i_var );
1325                                                tape[arg[i]].connect_type = yes_connected;
1326                                        }
1327                                        mask = mask << 1;
1328                                }
1329                        }
1330                        break;  // --------------------------------------------
1331
1332                        // Operations where there is noting to do
1333                        case BeginOp:
1334                        case ComOp:
1335                        case EndOp:
1336                        case InvOp:
1337                        case ParOp:
1338                        case PriOp:
1339                        break;  // --------------------------------------------
1340
1341                        // Load using a parameter index
1342                        case LdpOp:
1343                        if( tape[i_var].connect_type != not_connected )
1344                        {
1345                                i                = vecad[ arg[0] - 1 ];
1346                                vecad_connect[i] = yes_connected;
1347                        }
1348                        break; // --------------------------------------------
1349
1350                        // Load using a variable index
1351                        case LdvOp:
1352                        if( tape[i_var].connect_type != not_connected )
1353                        {
1354                                i                    = vecad[ arg[0] - 1 ];
1355                                vecad_connect[i]     = yes_connected;
1356                                tape[arg[1]].connect_type = yes_connected;
1357                        }
1358                        break; // --------------------------------------------
1359
1360                        // Store a variable using a parameter index
1361                        case StpvOp:
1362                        i = vecad[ arg[0] - 1 ];
1363                        if( vecad_connect[i] != not_connected )
1364                                tape[arg[2]].connect_type = yes_connected;
1365                        break; // --------------------------------------------
1366
1367                        // Store a variable using a variable index
1368                        case StvvOp:
1369                        i = vecad[ arg[0] - 1 ];
1370                        if( vecad_connect[i] )
1371                        {       tape[arg[1]].connect_type = yes_connected;
1372                                tape[arg[2]].connect_type = yes_connected;
1373                        }
1374                        break; 
1375                        // ============================================================
1376                        case UserOp:
1377                        // start or end atomic operation sequence
1378                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
1379                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
1380                        if( user_state == user_end )
1381                        {       user_index = arg[0];
1382                                user_id    = arg[1];
1383                                user_n     = arg[2];
1384                                user_m     = arg[3];
1385                                user_q     = 1;
1386                                if(user_s.size() != user_n )
1387                                        user_s.resize(user_n);
1388                                if(user_r.size() != user_m )
1389                                        user_r.resize(user_m);
1390                                user_j     = user_n;
1391                                user_i     = user_m;
1392                                user_state = user_ret;
1393                                user_keep.push(false);
1394                        }
1395                        else
1396                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
1397                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
1398                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
1399                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
1400                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
1401                                user_state = user_end;
1402               }
1403                        break;
1404
1405                        case UsrapOp:
1406                        // parameter argument in an atomic operation sequence
1407                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
1408                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
1409                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1410                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
1411                        --user_j;
1412                        if( user_j == 0 )
1413                                user_state = user_start;
1414                        break;
1415
1416                        case UsravOp:
1417                        // variable argument in an atomic operation sequence
1418                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
1419                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
1420                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1421                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
1422                        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
1423                        --user_j;
1424                        if( ! user_s[user_j].empty() )
1425                        {       tape[arg[0]].connect_type = yes_connected;
1426                                user_keep.top() = true;
1427                        }
1428                        if( user_j == 0 )
1429                                user_state = user_start;
1430                        break;
1431
1432                        case UsrrpOp:
1433                        // parameter result in an atomic operation sequence
1434                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
1435                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
1436                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
1437                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
1438                        --user_i;
1439                        user_r[user_i].clear();
1440                        if( user_i == 0 )
1441                        {       // call users function for this operation
1442                                atomic_base<Base>* atom = 
1443                                        atomic_base<Base>::class_object(user_index);
1444                                atom->set_id(user_id);
1445                                atom->rev_sparse_jac(
1446                                        user_q, user_r, user_s
1447                                );
1448                                user_state = user_arg;
1449                        }
1450                        break;
1451
1452                        case UsrrvOp:
1453                        // variable result in an atomic operation sequence
1454                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
1455                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
1456                        --user_i;
1457                        user_r[user_i].clear();
1458                        if( tape[i_var].connect_type != not_connected )
1459                                user_r[user_i].insert(0);
1460                        if( user_i == 0 )
1461                        {       // call users function for this operation
1462                                atomic_base<Base>* atom = 
1463                                        atomic_base<Base>::class_object(user_index);
1464                                atom->set_id(user_id);
1465                                atom->rev_sparse_jac(
1466                                        user_q, user_r, user_s
1467                                );
1468                                user_state = user_arg;
1469                        }
1470                        break;
1471                        // ============================================================
1472
1473                        // all cases should be handled above
1474                        default:
1475                        CPPAD_ASSERT_UNKNOWN(0);
1476                }
1477        }
1478        // values corresponding to BeginOp
1479        CPPAD_ASSERT_UNKNOWN( i_op == 0 && i_var == 0 && op == BeginOp );
1480        tape[i_var].op = op;
1481        // -------------------------------------------------------------
1482
1483        // Erase all information in the recording
1484        rec->free();
1485
1486        // Initilaize table mapping hash code to variable index in tape
1487        // as pointing to the BeginOp at the beginning of the tape
1488        CppAD::vector<size_t>  hash_table_var(CPPAD_HASH_TABLE_SIZE);
1489        for(i = 0; i < CPPAD_HASH_TABLE_SIZE; i++)
1490                hash_table_var[i] = 0;
1491        CPPAD_ASSERT_UNKNOWN( tape[0].op == BeginOp );
1492
1493        // initialize mapping from old variable index to new variable index
1494        for(i = 0; i < num_var; i++)
1495                tape[i].new_var = num_var; // invalid index
1496       
1497
1498        // initialize mapping from old VecAD index to new VecAD index
1499        CppAD::vector<size_t> new_vecad_ind(num_vecad_ind);
1500        for(i = 0; i < num_vecad_ind; i++)
1501                new_vecad_ind[i] = num_vecad_ind; // invalid index
1502
1503        j = 0;     // index into the old set of indices
1504        for(i = 0; i < num_vecad_vec; i++)
1505        {       // length of this VecAD
1506                size_t length = play->GetVecInd(j);
1507                if( vecad_connect[i] != not_connected )
1508                {       // Put this VecAD vector in new recording
1509                        CPPAD_ASSERT_UNKNOWN(length < num_vecad_ind);
1510                        new_vecad_ind[j] = rec->PutVecInd(length);
1511                        for(k = 1; k <= length; k++) new_vecad_ind[j+k] =
1512                                rec->PutVecInd(
1513                                        rec->PutPar(
1514                                                play->GetPar( 
1515                                                        play->GetVecInd(j+k)
1516                        ) ) );
1517                }
1518                // start of next VecAD
1519                j       += length + 1;
1520        }
1521        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
1522
1523        // start playing the operations in the forward direction
1524        play->start_forward(op, arg, i_op, i_var);
1525
1526        // playing forward skips BeginOp at the beginning, but not EndOp at
1527        // the end.  Put BeginOp at beginning of recording
1528        CPPAD_ASSERT_UNKNOWN( op == BeginOp );
1529        CPPAD_ASSERT_NARG_NRES(BeginOp, 0, 1);
1530        tape[i_var].new_var = rec->PutOp(BeginOp);
1531
1532        // temporary buffer for new argument values
1533        addr_t new_arg[6];
1534
1535        // temporary work space used by optimize_record_csum
1536        // (decalared here to avoid realloaction of memory)
1537        optimize_csum_stacks csum_work;
1538
1539        user_state = user_start;
1540        while(op != EndOp)
1541        {       // next op
1542                play->next_forward(op, arg, i_op, i_var);
1543                CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
1544                CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
1545
1546                // determine if we should keep this operation in the new
1547                // operation sequence
1548                bool keep;
1549                switch( op )
1550                {       case ComOp:
1551                        case PriOp:
1552                        keep = false;
1553                        break;
1554
1555                        case InvOp:
1556                        case EndOp:
1557                        keep = true;
1558                        break;
1559
1560                        case StppOp:
1561                        case StvpOp:
1562                        case StpvOp:
1563                        case StvvOp:
1564                        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
1565                        i = vecad[ arg[0] - 1 ];
1566                        keep = vecad_connect[i] != not_connected;
1567                        break;
1568
1569                        case AddpvOp:
1570                        case AddvvOp:
1571                        case SubpvOp:
1572                        case SubvpOp:
1573                        case SubvvOp:
1574                        keep  = tape[i_var].connect_type != not_connected;
1575                        keep &= tape[i_var].connect_type != csum_connected;
1576                        break; 
1577
1578                        case UserOp:
1579                        case UsrapOp:
1580                        case UsravOp:
1581                        case UsrrpOp:
1582                        case UsrrvOp:
1583                        keep = user_keep.top();
1584                        break;
1585
1586                        default:
1587                        keep = tape[i_var].connect_type != not_connected;
1588                        break;
1589                }
1590
1591                size_t         match_var    = 0;
1592                unsigned short code         = 0;
1593                bool           replace_hash = false;
1594                if( keep ) switch( op )
1595                {
1596                        // Unary operator where operand is arg[0]
1597                        case AbsOp:
1598                        case AcosOp:
1599                        case AsinOp:
1600                        case AtanOp:
1601                        case CosOp:
1602                        case CoshOp:
1603                        case ExpOp:
1604                        case LogOp:
1605                        case SinOp:
1606                        case SinhOp:
1607                        case SqrtOp:
1608                        case TanOp:
1609                        case TanhOp:
1610                        match_var = optimize_unary_match(
1611                                tape                ,  // inputs
1612                                i_var               ,
1613                                play->num_rec_par() ,
1614                                play->GetPar()      ,
1615                                hash_table_var      ,
1616                                code                  // outputs
1617                        );
1618                        if( match_var > 0 )
1619                                tape[i_var].new_var = match_var;
1620                        else
1621                        {
1622                                replace_hash = true;
1623                                new_arg[0]   = tape[ arg[0] ].new_var;
1624                                rec->PutArg( new_arg[0] );
1625                                i                   = rec->PutOp(op);
1626                                tape[i_var].new_var = i;
1627                                CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < i );
1628                        }
1629                        break;
1630                        // ---------------------------------------------------
1631                        // Binary operators where
1632                        // left is a variable and right is a parameter
1633                        case SubvpOp:
1634                        if( tape[arg[0]].connect_type == csum_connected )
1635                        {
1636                                // convert to a sequence of summation operators
1637                                tape[i_var].new_var = optimize_record_csum(
1638                                        tape                , // inputs
1639                                        i_var               ,
1640                                        play->num_rec_par() ,
1641                                        play->GetPar()      ,
1642                                        rec                 ,
1643                                        csum_work
1644                                );
1645                                // abort rest of this case
1646                                break;
1647                        }
1648                        case DivvpOp:
1649                        case PowvpOp:
1650                        match_var = optimize_binary_match(
1651                                tape                ,  // inputs
1652                                i_var               ,
1653                                play->num_rec_par() ,
1654                                play->GetPar()      ,
1655                                hash_table_var      ,
1656                                code                  // outputs
1657                        );
1658                        if( match_var > 0 )
1659                                tape[i_var].new_var = match_var;
1660                        else
1661                        {       tape[i_var].new_var = optimize_record_vp(
1662                                        tape                , // inputs
1663                                        i_var               ,
1664                                        play->num_rec_par() ,
1665                                        play->GetPar()      ,
1666                                        rec                 ,
1667                                        op                  ,
1668                                        arg
1669                                );
1670                                replace_hash = true;
1671                        }
1672                        break;
1673                        // ---------------------------------------------------
1674                        // Binary operators where
1675                        // left is a parameter and right is a variable
1676                        case SubpvOp:
1677                        case AddpvOp:
1678                        if( tape[arg[1]].connect_type == csum_connected )
1679                        {
1680                                // convert to a sequence of summation operators
1681                                tape[i_var].new_var = optimize_record_csum(
1682                                        tape                , // inputs
1683                                        i_var               ,
1684                                        play->num_rec_par() ,
1685                                        play->GetPar()      ,
1686                                        rec                 ,
1687                                        csum_work
1688                                );
1689                                // abort rest of this case
1690                                break;
1691                        }
1692                        case DivpvOp:
1693                        case MulpvOp:
1694                        case PowpvOp:
1695                        match_var = optimize_binary_match(
1696                                tape                ,  // inputs
1697                                i_var               ,
1698                                play->num_rec_par() ,
1699                                play->GetPar()      ,
1700                                hash_table_var      ,
1701                                code                  // outputs
1702                        );
1703                        if( match_var > 0 )
1704                                tape[i_var].new_var = match_var;
1705                        else
1706                        {       tape[i_var].new_var = optimize_record_pv(
1707                                        tape                , // inputs
1708                                        i_var               ,
1709                                        play->num_rec_par() ,
1710                                        play->GetPar()      ,
1711                                        rec                 ,
1712                                        op                  ,
1713                                        arg
1714                                );
1715                                replace_hash = true;
1716                        }
1717                        break;
1718                        // ---------------------------------------------------
1719                        // Binary operator where
1720                        // both operators are variables
1721                        case AddvvOp:
1722                        case SubvvOp:
1723                        if( (tape[arg[0]].connect_type == csum_connected) |
1724                            (tape[arg[1]].connect_type == csum_connected)
1725                        )
1726                        {
1727                                // convert to a sequence of summation operators
1728                                tape[i_var].new_var = optimize_record_csum(
1729                                        tape                , // inputs
1730                                        i_var               ,
1731                                        play->num_rec_par() ,
1732                                        play->GetPar()      ,
1733                                        rec                 ,
1734                                        csum_work
1735                                );
1736                                // abort rest of this case
1737                                break;
1738                        }
1739                        case DivvvOp:
1740                        case MulvvOp:
1741                        case PowvvOp:
1742                        match_var = optimize_binary_match(
1743                                tape                ,  // inputs
1744                                i_var               ,
1745                                play->num_rec_par() ,
1746                                play->GetPar()      ,
1747                                hash_table_var      ,
1748                                code                  // outputs
1749                        );
1750                        if( match_var > 0 )
1751                                tape[i_var].new_var = match_var;
1752                        else
1753                        {       tape[i_var].new_var = optimize_record_vv(
1754                                        tape                , // inputs
1755                                        i_var               ,
1756                                        play->num_rec_par() ,
1757                                        play->GetPar()      ,
1758                                        rec                 ,
1759                                        op                  ,
1760                                        arg
1761                                );
1762                                replace_hash = true;
1763                        }
1764                        break;
1765                        // ---------------------------------------------------
1766                        // Conditional expression operators
1767                        case CExpOp:
1768                        CPPAD_ASSERT_NARG_NRES(op, 6, 1);
1769                        new_arg[0] = arg[0];
1770                        new_arg[1] = arg[1];
1771                        mask = 1;
1772                        for(i = 2; i < 6; i++)
1773                        {       if( arg[1] & mask )
1774                                {       new_arg[i] = tape[arg[i]].new_var;
1775                                        CPPAD_ASSERT_UNKNOWN( 
1776                                                size_t(new_arg[i]) < num_var
1777                                        );
1778                                }
1779                                else    new_arg[i] = rec->PutPar( 
1780                                                play->GetPar( arg[i] )
1781                                );
1782                                mask = mask << 1;
1783                        }
1784                        rec->PutArg(
1785                                new_arg[0] ,
1786                                new_arg[1] ,
1787                                new_arg[2] ,
1788                                new_arg[3] ,
1789                                new_arg[4] ,
1790                                new_arg[5] 
1791                        );
1792                        tape[i_var].new_var = rec->PutOp(op);
1793                        break;
1794                        // ---------------------------------------------------
1795                        // Operations with no arguments and no results
1796                        case EndOp:
1797                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
1798                        rec->PutOp(op);
1799                        break;
1800                        // ---------------------------------------------------
1801                        // Operations with no arguments and one result
1802                        case InvOp:
1803                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
1804                        tape[i_var].new_var = rec->PutOp(op);
1805                        break;
1806                        // ---------------------------------------------------
1807                        // Operations with one argument that is a parameter
1808                        case ParOp:
1809                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
1810                        new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
1811
1812                        rec->PutArg( new_arg[0] );
1813                        tape[i_var].new_var = rec->PutOp(op);
1814                        break;
1815                        // ---------------------------------------------------
1816                        // Load using a parameter index
1817                        case LdpOp:
1818                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1819                        new_arg[0] = new_vecad_ind[ arg[0] ];
1820                        new_arg[1] = arg[1];
1821                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1822                        rec->PutArg( 
1823                                new_arg[0], 
1824                                new_arg[1], 
1825                                0
1826                        );
1827                        tape[i_var].new_var = rec->PutOp(op);
1828                        break;
1829                        // ---------------------------------------------------
1830                        // Load using a variable index
1831                        case LdvOp:
1832                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1833                        new_arg[0] = new_vecad_ind[ arg[0] ];
1834                        new_arg[1] = tape[arg[1]].new_var;
1835                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1836                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1837                        rec->PutArg( 
1838                                new_arg[0], 
1839                                new_arg[1], 
1840                                0
1841                        );
1842                        tape[i_var].new_var = rec->PutOp(op);
1843                        break;
1844                        // ---------------------------------------------------
1845                        // Store a parameter using a parameter index
1846                        case StppOp:
1847                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1848                        new_arg[0] = new_vecad_ind[ arg[0] ];
1849                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1850                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1851                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1852                        rec->PutArg(
1853                                new_arg[0], 
1854                                new_arg[1], 
1855                                new_arg[2]
1856                        );
1857                        rec->PutOp(op);
1858                        break;
1859                        // ---------------------------------------------------
1860                        // Store a parameter using a variable index
1861                        case StvpOp:
1862                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1863                        new_arg[0] = new_vecad_ind[ arg[0] ];
1864                        new_arg[1] = tape[arg[1]].new_var;
1865                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1866                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1867                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1868                        rec->PutArg(
1869                                new_arg[0], 
1870                                new_arg[1], 
1871                                new_arg[2]
1872                        );
1873                        rec->PutOp(op);
1874                        break;
1875                        // ---------------------------------------------------
1876                        // Store a variable using a parameter index
1877                        case StpvOp:
1878                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1879                        new_arg[0] = new_vecad_ind[ arg[0] ];
1880                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1881                        new_arg[2] = tape[arg[2]].new_var;
1882                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1883                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1884                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
1885                        rec->PutArg(
1886                                new_arg[0], 
1887                                new_arg[1], 
1888                                new_arg[2]
1889                        );
1890                        rec->PutOp(op);
1891                        break;
1892                        // ---------------------------------------------------
1893                        // Store a variable using a variable index
1894                        case StvvOp:
1895                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1896                        new_arg[0] = new_vecad_ind[ arg[0] ];
1897                        new_arg[1] = tape[arg[1]].new_var;
1898                        new_arg[2] = tape[arg[2]].new_var;
1899                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[0]) < num_vecad_ind );
1900                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[1]) < num_var );
1901                        CPPAD_ASSERT_UNKNOWN( size_t(new_arg[2]) < num_var );
1902                        rec->PutArg(
1903                                new_arg[0], 
1904                                new_arg[1], 
1905                                new_arg[2]
1906                        );
1907                        rec->PutOp(op);
1908                        break;
1909
1910                        // -----------------------------------------------------------
1911                        case UserOp:
1912                        CPPAD_ASSERT_NARG_NRES(op, 4, 0);
1913                        if( user_state == user_start )
1914                                user_state = user_arg;
1915                        else
1916                        {       user_state = user_start;
1917                                user_keep.pop();       
1918                        }
1919                        // user_index, user_id, user_n, user_m
1920                        rec->PutArg(arg[0], arg[1], arg[2], arg[3]);
1921                        rec->PutOp(UserOp);
1922                        break;
1923
1924                        case UsrapOp:
1925                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1926                        new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
1927                        rec->PutArg(new_arg[0]);
1928                        rec->PutOp(UsrapOp);
1929                        break;
1930
1931                        case UsravOp:
1932                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1933                        new_arg[0] = tape[arg[0]].new_var;
1934                        if( size_t(new_arg[0]) < num_var )
1935                        {       rec->PutArg(new_arg[0]);
1936                                rec->PutOp(UsravOp);
1937                        }
1938                        else
1939                        {       // This argument does not affect the result and
1940                                // has been optimized out so use nan in its place.
1941                                new_arg[0] = rec->PutPar( nan(Base(0)) );
1942                                rec->PutArg(new_arg[0]);
1943                                rec->PutOp(UsrapOp);
1944                        }
1945                        break;
1946
1947                        case UsrrpOp:
1948                        CPPAD_ASSERT_NARG_NRES(op, 1, 0);
1949                        new_arg[0] = rec->PutPar( play->GetPar(arg[0]) );
1950                        rec->PutArg(new_arg[0]);
1951                        rec->PutOp(UsrrpOp);
1952                        break;
1953                       
1954                        case UsrrvOp:
1955                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
1956                        tape[i_var].new_var = rec->PutOp(UsrrvOp);
1957                        break;
1958                        // ---------------------------------------------------
1959
1960                        // all cases should be handled above
1961                        default:
1962                        CPPAD_ASSERT_UNKNOWN(false);
1963
1964                }
1965                if( replace_hash )
1966                {       // The old variable index i_var corresponds to the
1967                        // new variable index tape[i_var].new_var. In addition
1968                        // this is the most recent variable that has this code.
1969                        hash_table_var[code] = i_var;
1970                }
1971
1972        }
1973        // modify the dependent variable vector to new indices
1974        for(i = 0; i < dep_taddr.size(); i++ )
1975        {       CPPAD_ASSERT_UNKNOWN( size_t(tape[ dep_taddr[i] ].new_var) < num_var );
1976                dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
1977        }
1978}
1979
1980/*!
1981Optimize a player object operation sequence
1982
1983The operation sequence for this object is replaced by one with fewer operations
1984but the same funcition and derivative values.
1985
1986\tparam Base
1987base type for the operator; i.e., this operation was recorded
1988using AD< \a Base > and computations by this routine are done using type
1989\a Base.
1990
1991*/
1992template <class Base>
1993void ADFun<Base>::optimize(void)
1994{       // place to store the optimized version of the recording
1995        recorder<Base> rec;
1996
1997        // number of independent variables
1998        size_t n = ind_taddr_.size();
1999
2000# ifndef NDEBUG
2001        size_t i, j, m = dep_taddr_.size();
2002        CppAD::vector<Base> x(n), y(m), check(m);
2003        bool check_zero_order = taylor_per_var_ > 0;
2004        Base max_taylor(0);
2005        if( check_zero_order )
2006        {       // zero order coefficients for independent vars
2007                for(j = 0; j < n; j++)
2008                {       CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
2009                        CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
2010                        x[j] = taylor_[ ind_taddr_[j] * taylor_col_dim_ + 0];
2011                }
2012                // zero order coefficients for dependent vars
2013                for(i = 0; i < m; i++)
2014                {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
2015                        y[i] = taylor_[ dep_taddr_[i] * taylor_col_dim_ + 0];
2016                }
2017                // maximum zero order coefficient not counting BeginOp at beginning
2018                // (which is correpsonds to uninitialized memory).
2019                for(i = 1; i < total_num_var_; i++)
2020                {       if(  abs_geq(taylor_[i*taylor_col_dim_+0] , max_taylor) )
2021                                max_taylor = taylor_[i*taylor_col_dim_+0];
2022                }
2023        }
2024# endif
2025
2026        // create the optimized recording
2027        CppAD::optimize<Base>(n, dep_taddr_, &play_, &rec);
2028
2029        // number of variables in the recording
2030        total_num_var_ = rec.num_rec_var();
2031
2032        // now replace the recording
2033        play_.get(rec);
2034
2035        // free memory allocated for sparse Jacobian calculation
2036        // (the results are no longer valid)
2037        for_jac_sparse_pack_.resize(0, 0);
2038        for_jac_sparse_set_.resize(0,0);
2039
2040        // free old Taylor coefficient memory
2041        taylor_.free();
2042        taylor_per_var_ = 0;
2043        taylor_col_dim_ = 0;
2044
2045# ifndef NDEBUG
2046        if( check_zero_order )
2047        {
2048                // zero order forward calculation using new operation sequence
2049                check = Forward(0, x);
2050
2051                // check results
2052                Base eps = 10. * epsilon<Base>();
2053                for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
2054                        abs_geq( eps * max_taylor , check[i] - y[i] ) ,
2055                        "Error during check of f.optimize()."
2056                );
2057
2058                // Erase memory that this calculation was done so NDEBUG gives
2059                // same final state for this object (from users perspective)
2060                taylor_per_var_ = 0;
2061        }
2062# endif
2063}
2064
2065/*! \} */
2066} // END_CPPAD_NAMESPACE
2067# endif
Note: See TracBrowser for help on using the repository browser.