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

Last change on this file since 1600 was 1600, checked in by bradbell, 11 years ago

trunk: In optimize, factor out (parameter op variable) case to a function.

  • Property svn:keywords set to Id
File size: 36.9 KB
Line 
1/* $Id: optimize.hpp 1600 2009-12-04 16:54:16Z bradbell $ */
2# ifndef CPPAD_OPTIMIZE_INCLUDED
3# define CPPAD_OPTIMIZE_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Common 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$$
24
25$section Optimize the Tape Corresponding to an ADFun object$$
26
27$index optimize$$
28$index tape, optimize$$
29$index sequence, optimize operations$$
30$index operations, optimize sequence$$
31$index speed, optimize$$
32$index memory, optimize$$
33
34$head Syntax$$
35$icode%f%.optimize()%$$
36
37
38$head Purpose$$
39The operation sequence corresponding to an $cref/ADFun/$$ object can
40be very large and involve many operations; see the
41size functions in $cref/seq_property/$$.
42This enables one to reduce the number of operations
43and thereby reduce the time and the memory required to
44compute function and derivative values using an $code ADFun$$ object.
45
46$head f$$
47The object $icode f$$ has prototype
48$icode%
49        ADFun<%Base%> %f%
50%$$
51
52$head Improvements$$
53You can see the reduction in number of variables in the operation sequence
54by calling the function $cref/f.size_var()/seq_property/size_var/$$
55before and after the optimization procedure.
56Given that the optimization procedure takes time,
57it may be faster to skip this optimize procedure and just compute
58derivatives using the original operation sequence.
59
60$subhead Testing$$
61You can run the CppAD $cref/speed/speed_main/$$ tests and see
62the corresponding changes in number of variables and execution time.
63This will require that you specify the
64$cref/--with-Speed/InstallUnix/--with-Speed/$$ option when you
65configure CppAD.
66In this case, you can run the speed tests with optimization
67using the command
68$codep
69        speed/cppad/cppad speed 123 optimize
70$$
71and without optimization using the command
72$codep
73        speed/cppad/cppad speed 123
74$$
75Note that $code 123$$ is used for a random number seed and can be
76replaced by any integer in the commands above.
77
78$head Efficiency$$
79The $code optimize$$ member function
80may greatly reduce the number of variables
81in the operation sequence; see $cref/size_var/seq_property/size_var/$$.
82If a $cref/zero order forward/ForwardZero/$$ calculation is done during
83the construction of f, it will require more memory
84and time that after the optimization procedure.
85In addition, it will need to be redone.
86For this reason, it is more efficient to use
87$codei%
88        ADFun<%Base%> %f%;
89        %f%.Dependent(%x%, %y%);
90        %f%.optimize();
91%$$
92instead of
93$codei%
94        ADFun<%Base%> %f%(%x%, %y%)
95        %f%.optimize();
96%$$
97See the discussion about
98$cref/sequence constructors/FunConstruct/Sequence Constructor/$$.
99
100$head Comparison Operators$$
101Any comparison operators that are in the tape are removed by this operation.
102Hence the return value of $cref/CompareChange/$$ will always be zero
103for an optimized tape (even if $code NDEBUG$$ is not defined).
104
105$head Checking Optimization$$
106$index NDEBUG$$
107If $code NDEBUG$$ is not defined
108and $cref/f.size_taylor()/size_taylor/$$ is greater than zero,
109a $cref/ForwardZero/$$ calculation is done using the optimized version
110of $icode f$$ and the results are checked to see that they are
111the same as before.
112If they are not the same, the
113$cref/ErrorHandler/$$ is called with a known error message
114related to $icode%f%.optimize()%$$.
115
116$head Example$$
117$children%
118        example/optimize.cpp
119%$$
120The file
121$xref/optimize.cpp/$$
122contains an example and test of this operation.
123It returns true if it succeeds and false otherwise.
124
125$end
126-----------------------------------------------------------------------------
127*/
128
129/*!
130\file optimize.hpp
131Routines for optimizing a tape
132*/
133
134/*!
135\def CPPAD_OPTIMIZE_TRACE
136This value is either zero or one.
137Zero is the normal operation value.
138If it is one, a trace of the reverse sweep is printed.
139This sweep determines which variables are connected to the
140dependent variables.
141*/
142# define CPPAD_OPTIMIZE_TRACE 0
143
144CPPAD_BEGIN_NAMESPACE
145/*!
146Structure used by \c optimize to hold information about one variable.
147This \c struct would be local inside of \c optimize,
148but the current C++ standard does not support local template parameters.
149*/
150struct optimize_variable {
151        /// Operator for which this variable is the result, \c NumRes(op) > 0.
152        /// Set by the reverse sweep at beginning of optimization.
153        OpCode         op;       
154        /// Pointer to first argument (child) for this operator.
155        /// Set by the reverse sweep at beginning of optimization.
156        const size_t*  arg;
157        /*!
158        Information about the parrents for this variable.
159
160        \li
161        If \c parrent == 0, this variable has no parrents.
162        \li
163        If \c 0 < parrent < num_var, then \c parrent is the only parrent.
164        \li
165        If \c parrent == num_var, this variable has more than one parrent.
166
167        Set by the reverse sweep at beginning of optimization.
168        May be changed during the forward sweep when a variable is
169        optimized out and its \c parrent is set to zero.
170        */
171        size_t         parrent; 
172        /// Index of this variable in the optimized operation sequence.
173        /// Set by the forward sweep at end of optimization.
174        size_t         new_var; 
175};
176
177/*!
178Check a unary operator for a complete match with a previous operator.
179
180A complete match means that the result of the previous operator
181can be used inplace of the result for this operator.
182
183\param old_res
184is the index in the old operation sequence for
185the variable corresponding to the result for this unary operator.
186
187\param op
188is the unary operator that we are checking a match for.
189Assertion: NumArg(op) == 1 and NumRes(op) > 0.
190
191\param old_arg
192the value \c old_arg[0] is the argument for this unary operator
193in the old operation sequence.
194Assertion: old_arg[0] < old_res.
195
196\param npar
197is the number of paraemters corresponding to this operation sequence.
198
199\param par
200is a vector of length \a npar containing the parameters
201for this operation sequence; i.e.,
202given a parameter index \c i, the corresponding parameter value is
203\a par[i].
204
205\param hash_table_var
206is a vector with size CPPAD_HASH_TABLE_SIZE
207that maps a hash code to the corresponding
208variable index in the old operation sequence.
209All the values in this table must be less than \c old_res.
210
211\param tape
212is a vector that maps a variable index, in the old operation sequence,
213to the corresponding \c optimze_variable information.
214For indices \c i less than \c old_res,
215the following must be true of \c tape[i]:
216
217\li tape[i].op
218is the operator in the old operation sequence
219corresponding to the old variable index \c i.
220Assertion: NumRes(tape[i].op) > 0.
221
222\li tape[i].arg
223For j < NumArg( tape[i].op ),
224tape[i].arg[j] is the j-th the argument, in the old operation sequence,
225corresponding to the old variable index \c i.
226Assertion: tape[i].arg[j] < i.
227
228\li \c tape[i].new_var
229is the variable in the new operation sequence
230corresponding to the old variable index \c i.
231Assertion: tape[i].new_var < old_res.
232
233\param code
234The input value of \c code does not matter.
235The output value of \c code is the hash code corresponding to
236this operation in the new operation sequence.
237
238\return
239If the return value is zero,
240no match was found.
241If the return value is greater than zero,
242it is the index of a new variable that can be used to replace the
243old variable.
244
245\par Check Assertions
246\li
247*/
248template <class Base>
249inline size_t optimize_unary_match(
250        size_t                                         old_res        ,
251        OpCode                                         op             ,
252        const size_t*                                  old_arg        ,
253        size_t                                         npar           ,
254        const Base*                                    par            ,
255        const CppAD::vector<size_t>&                   hash_table_var ,
256        const CppAD::vector<struct optimize_variable>& tape           ,
257        unsigned short&                                code           )
258{       size_t new_arg[1];
259       
260        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
261        CPPAD_ASSERT_UNKNOWN( NumRes(op) > 0  );
262        CPPAD_ASSERT_UNKNOWN( old_arg[0] < old_res );
263        new_arg[0] = tape[old_arg[0]].new_var;
264        CPPAD_ASSERT_UNKNOWN( new_arg[0] < old_res );
265        code = hash_code(
266                op                  , 
267                new_arg             ,
268                npar                ,
269                par
270        );
271        size_t  i               = hash_table_var[code];
272        CPPAD_ASSERT_UNKNOWN( i < old_res );
273        if( op == tape[i].op )
274        {       size_t k = tape[i].arg[0];
275                CPPAD_ASSERT_UNKNOWN( k < i );
276                if (new_arg[0] == tape[k].new_var );
277                        return tape[i].new_var;
278        }
279        return 0;
280} 
281
282/*!
283Check a (parameter op variable) for a complete match with a previous operator.
284
285A complete match means that the result of the previous operator
286can be used inplace of the result for this operator.
287
288\param old_res
289is the index in the old operation sequence for
290the variable corresponding to the result for this unary operator.
291
292\param op
293is the unary operator that we are checking a match for.
294The value of op be the same for a match; i.e., the commutative case
295of an AddpvOp matching an AddvpOp or MulpvOp matching a MulvpOp
296is not checked for.
297Assertion: NumArg(op) == 2 and NumRes(op) > 0.
298
299\param old_arg
300The value \c old_arg[0] is the index of the
301first operand (parameter) in the \a par vector.
302Assertion: old_arg[0] < npar.
303The value \c old_arg[1] is the index of the second operand (variable)
304in the old operation sequence.
305Assertion: old_arg[1] < old_res.
306
307\param npar
308is the number of paraemters corresponding to this operation sequence.
309
310\param par
311is a vector of length \a npar containing the parameters
312for this operation sequence; i.e.,
313given a parameter index \c i, the corresponding parameter value is
314\a par[i].
315
316\param hash_table_var
317is a vector with size CPPAD_HASH_TABLE_SIZE
318that maps a hash code to the corresponding
319variable index in the old operation sequence.
320All the values in this table must be less than \c old_res.
321
322\param tape
323is a vector that maps a variable index, in the old operation sequence,
324to the corresponding \c optimze_variable information.
325For indices \c i less than \c old_res,
326the following must be true of \c tape[i]:
327
328\li tape[i].op
329is the operator in the old operation sequence
330corresponding to the old variable index \c i.
331Assertion: NumRes(tape[i].op) > 0.
332
333\li tape[i].arg
334For j < NumArg( tape[i].op ),
335tape[i].arg[j] is the j-th the operand, in the old operation sequence,
336corresponding to the old variable index \c i.
337Assertion: tape[i].arg[j] < i.
338
339\li \c tape[i].new_var
340is the variable in the new operation sequence
341corresponding to the old variable index \c i.
342Assertion: tape[i].new_var < old_res.
343
344\param code
345The input value of \c code does not matter.
346The output value of \c code is the hash code corresponding to
347this operation in the new operation sequence.
348
349\return
350If the return value is zero,
351no match was found.
352If the return value is greater than zero,
353it is the index of a new variable that can be used to replace the
354old variable.
355
356\par Check Assertions
357\li
358*/
359template <class Base>
360inline size_t optimize_par_op_var_match(
361        size_t                                         old_res        ,
362        OpCode                                         op             ,
363        const size_t*                                  old_arg        ,
364        size_t                                         npar           ,
365        const Base*                                    par            ,
366        const CppAD::vector<size_t>&                   hash_table_var ,
367        const CppAD::vector<struct optimize_variable>& tape           ,
368        unsigned short&                                code           )
369{       size_t new_arg[2];
370       
371        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
372        CPPAD_ASSERT_UNKNOWN( NumRes(op) > 0  );
373        CPPAD_ASSERT_UNKNOWN( old_arg[0] < npar );
374        CPPAD_ASSERT_UNKNOWN( old_arg[1] < old_res );
375        new_arg[0] = old_arg[0];
376        new_arg[1] = tape[old_arg[1]].new_var;
377        CPPAD_ASSERT_UNKNOWN( new_arg[1] < old_res );
378        code = hash_code(
379                op                  , 
380                new_arg             ,
381                npar                ,
382                par
383        );
384        size_t  i  = hash_table_var[code];
385        CPPAD_ASSERT_UNKNOWN( i < old_res );
386        if( op == tape[i].op )
387        {       size_t k = tape[i].arg[0];
388                CPPAD_ASSERT_UNKNOWN( k < npar );
389                // check if parameter values match
390                if( IdenticalEqualPar( par[ old_arg[0] ] , par[k] ) )
391                {       k = tape[i].arg[1];
392                        CPPAD_ASSERT_UNKNOWN( k < i );
393                        // check if same new variable
394                        if (new_arg[1] == tape[k].new_var );
395                                return tape[i].new_var;
396                }
397        }
398        return 0;
399} 
400
401
402/*!
403Convert a player object to an optimized recorder object
404
405\tparam Base
406base type for the operator; i.e., this operation was recorded
407using AD< \a Base > and computations by this routine are done using type
408\a Base.
409
410\param n
411is the number of independent variables on the tape.
412
413\param dep_taddr
414On input this vector contains the indices for each of the dependent
415variable values in the operation sequence corresponding to \a play.
416Upon return it contains the indices for the same variables but in
417the operation sequence corresponding to \a rec.
418
419\param play
420This is the operation sequence that we are optimizing.
421It is essentially const, except for play back state which
422changes while it plays back the operation seqeunce.
423
424\param rec
425The input contents of this recording does not matter.
426Upon return, it contains an optimized verison of the
427operation sequence corresponding to \a play.
428*/
429
430template <class Base>
431void optimize(
432        size_t                       n         ,
433        CppAD::vector<size_t>&       dep_taddr ,
434        player<Base>*                play      ,
435        recorder<Base>*              rec       ) 
436{
437        // temporary indices
438        size_t i, j, k;
439
440        // temporary variables
441        OpCode        op;   // current operator
442        const size_t *arg;  // operator arguments
443        size_t        i_var;  // index of first result for current operator
444
445        // range and domain dimensions for F
446        size_t m = dep_taddr.size();
447
448        // number of variables in the player
449        const size_t num_var = play->num_rec_var(); 
450
451        // number of  VecAD indices
452        size_t num_vecad_ind   = play->num_rec_vecad_ind();
453
454        // number of VecAD vectors
455        size_t num_vecad_vec   = play->num_rec_vecad_vec();
456
457        // -------------------------------------------------------------
458        // data structure that maps variable index in original operation
459        // sequence to corresponding operator information
460        CppAD::vector<struct optimize_variable> tape(num_var);
461        // -------------------------------------------------------------
462        // Determine parrent value for each variable
463
464        // initialize all variables has having no parrent
465        for(i = 0; i < num_var; i++)
466                tape[i].parrent = 0;
467
468        for(j = 0; j < m; j++)
469        {       // mark dependent variables as having multiple parents
470                tape[ dep_taddr[j] ].parrent = num_var;
471        }
472
473        // vecad_parrent contains a flag for each VecAD object.
474        // vecad maps a VecAD index (which corresponds to the beginning of the
475        // VecAD object) to the vecad_parrent falg for the VecAD object.
476# if CPPAD_OPTIMIZE_TRACE
477        CppAD::vector<size_t> vecad_offset(num_vecad_vec);
478# endif
479        CppAD::vector<bool>   vecad_parrent(num_vecad_vec);
480        CppAD::vector<size_t> vecad(num_vecad_ind);
481        j = 0;
482        for(i = 0; i < num_vecad_vec; i++)
483        {       vecad_parrent[i] = false;
484# if CPPAD_OPTIMIZE_TRACE
485                // offset for this VecAD object
486                vecad_offset[i]  = j + 1;
487# endif
488                // length of this VecAD
489                size_t length = play->GetVecInd(j);
490                // set to proper index for this VecAD
491                vecad[j] = i; 
492                for(k = 1; k <= length; k++)
493                        vecad[j+k] = num_vecad_vec; // invalid index
494                // start of next VecAD
495                j       += length + 1;
496        }
497        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
498
499        // Initialize a reverse mode sweep through the operation sequence
500        size_t i_op;
501        play->start_reverse(op, arg, i_op, i_var);
502        CPPAD_ASSERT_UNKNOWN( op == EndOp );
503        size_t mask;
504# if CPPAD_OPTIMIZE_TRACE
505        std::cout << std::endl;
506# endif
507        while(op != BeginOp)
508        {       // next op
509                play->next_reverse(op, arg, i_op, i_var);
510                // This if is not necessary becasue last assignment
511                // with this value of i_var will have NumRes(op) > 0
512                if( NumRes(op) > 0 )
513                {       tape[i_var].op = op;
514                        tape[i_var].arg = arg;
515                }
516# ifndef NDEBUG
517                if( i_op <= n )
518                {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
519                }
520                else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
521# endif
522
523# if CPPAD_OPTIMIZE_TRACE
524                printOp(
525                        std::cout, 
526                        play,
527                        i_var,
528                        op, 
529                        arg,
530                        0, 
531                        (bool *) CPPAD_NULL,
532                        1, 
533                        & connected[i_var]
534                );
535# endif
536
537                switch( op )
538                {
539                        // Unary operator where operand is arg[0]
540                        case AbsOp:
541                        case AddvpOp:
542                        case AcosOp:
543                        case AsinOp:
544                        case AtanOp:
545                        case CosOp:
546                        case CoshOp:
547                        case DisOp:
548                        case DivvpOp:
549                        case ExpOp:
550                        case LogOp:
551                        case MulvpOp:
552                        case PowvpOp:
553                        case SinOp:
554                        case SinhOp:
555                        case SqrtOp:
556                        case SubvpOp:
557                        if( tape[i_var].parrent > 0 )
558                        {       if( tape[arg[0]].parrent == 0 )
559                                        tape[arg[0]].parrent = i_var;
560                                else    tape[arg[0]].parrent = num_var;
561                        }
562                        break;
563
564                        // Unary operator where operand is arg[1]
565                        case AddpvOp:
566                        case DivpvOp:
567                        case MulpvOp:
568                        case SubpvOp:
569                        case PowpvOp:
570                        case PrivOp:
571                        if( tape[i_var].parrent > 0 )
572                        {       if( tape[arg[1]].parrent == 0 )
573                                        tape[arg[1]].parrent = i_var;
574                                else    tape[arg[1]].parrent = num_var;
575                        }
576                        break;
577
578                        // Binary operator where operands are arg[0], arg[1]
579                        case AddvvOp:
580                        case DivvvOp:
581                        case MulvvOp:
582                        case PowvvOp:
583                        case SubvvOp:
584                        if( tape[i_var].parrent > 0 )
585                        {       if( tape[arg[0]].parrent == 0 )
586                                        tape[arg[0]].parrent = i_var;
587                                else    tape[arg[0]].parrent = num_var;
588                                if( tape[arg[1]].parrent == 0 )
589                                        tape[arg[1]].parrent = i_var;
590                                else    tape[arg[1]].parrent = num_var;
591                        }
592                        break;
593
594                        // Conditional expression operators
595                        case CExpOp:
596                        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
597                        if( tape[i_var].parrent > 0 )
598                        {       mask = 1;
599                                for(i = 2; i < 6; i++) if( arg[1] & mask )
600                                {       CPPAD_ASSERT_UNKNOWN( arg[i] < i_var );
601                                        if( tape[arg[i]].parrent == 0 )
602                                                tape[arg[i]].parrent = i_var;
603                                        else    tape[arg[i]].parrent = num_var;
604                                }
605                                mask = mask << 1;
606                        }
607
608                        // Operations where there is noting to do
609                        case BeginOp:
610                        case ComOp:
611                        case EndOp:
612                        case InvOp:
613                        case ParOp:
614                        case PripOp:
615                        case StppOp:
616                        break; 
617
618                        // Load using a parameter index
619                        case LdpOp:
620                        if( tape[i_var].parrent > 0 )
621                        {       i = vecad[ arg[0] - 1 ];
622                                vecad_parrent[i] = true;
623                        }
624                        break;
625
626                        // Load using a variable index
627                        case LdvOp:
628                        if( tape[i_var].parrent > 0 )
629                        {       i = vecad[ arg[0] - 1 ];
630                                vecad_parrent[i] = true;
631                                if( tape[arg[1]].parrent == 0 )
632                                        tape[arg[1]].parrent = i_var;
633                                else    tape[arg[1]].parrent = num_var;
634                        }
635                        break;
636
637                        // Store a variable using a parameter index
638                        case StpvOp:
639                        i = vecad[ arg[0] - 1 ];
640                        if( vecad_parrent[i] )
641                        {       if( tape[arg[2]].parrent == 0 )
642                                        tape[arg[2]].parrent = i_var;
643                                else    tape[arg[2]].parrent = num_var;
644                        }
645                        break;
646
647                        // Store a variable using a variable index
648                        case StvvOp:
649                        i = vecad[ arg[0] - 1 ];
650                        if( vecad_parrent[i] )
651                        {       if( tape[arg[1]].parrent == 0 )
652                                        tape[arg[1]].parrent = i_var;
653                                else    tape[arg[1]].parrent = num_var;
654                                if( tape[arg[2]].parrent == 0 )
655                                        tape[arg[2]].parrent = i_var;
656                                else    tape[arg[2]].parrent = num_var;
657                        }
658                        break;
659
660                        // all cases should be handled above
661                        default:
662                        CPPAD_ASSERT_UNKNOWN(0);
663                }
664        }
665        // values corresponding to BeginOp
666        CPPAD_ASSERT_UNKNOWN( i_op == 0 && i_var == 0 && op == BeginOp );
667        tape[i_var].op = op;
668# if CPPAD_OPTIMIZE_TRACE
669        std::cout << "VecAD information:" << std::endl;
670        for(i = 0; i < num_vecad_vec; i++) 
671        {       std::cout << "offset  = " << vecad_offset[i];
672                std::cout << ", parrent = " << vecad_parrent[i] << std::endl;
673        }
674# endif
675        // -------------------------------------------------------------
676
677        // Erase all information in the recording
678        rec->Erase();
679
680        // Initilaize table mapping hash code to variable index in tape
681        // as pointing to the BeginOp at the beginning of the tape
682        CppAD::vector<size_t>  hash_table_var(CPPAD_HASH_TABLE_SIZE);
683        for(i = 0; i < CPPAD_HASH_TABLE_SIZE; i++)
684                hash_table_var[i] = 0;
685        CPPAD_ASSERT_UNKNOWN( tape[0].op == BeginOp );
686
687        // initialize mapping from old variable index to new variable index
688        for(i = 0; i < num_var; i++)
689                tape[i].new_var = num_var; // invalid index
690       
691
692        // initialize mapping from old VecAD index to new VecAD index
693        CppAD::vector<size_t> new_vecad_ind(num_vecad_ind);
694        for(i = 0; i < num_vecad_ind; i++)
695                new_vecad_ind[i] = num_vecad_ind; // invalid index
696
697        j = 0;     // index into the old set of indices
698        for(i = 0; i < num_vecad_vec; i++)
699        {       // length of this VecAD
700                size_t length = play->GetVecInd(j);
701                if( vecad_parrent[i] )
702                {       // Put this VecAD vector in new recording
703                        CPPAD_ASSERT_UNKNOWN(length < num_vecad_ind);
704                        new_vecad_ind[j] = rec->PutVecInd(length);
705                        for(k = 1; k <= length; k++) new_vecad_ind[j+k] =
706                                rec->PutVecInd(
707                                        rec->PutPar(
708                                                play->GetPar( 
709                                                        play->GetVecInd(j+k)
710                        ) ) );
711                }
712                // start of next VecAD
713                j       += length + 1;
714        }
715        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
716
717        // start playing the operations in the forward direction
718        play->start_forward(op, arg, i_op, i_var);
719
720        // playing forward skips BeginOp at the beginning, but not EndOp at
721        // the end.  Put BeginOp at beginning of recording
722        CPPAD_ASSERT_UNKNOWN( op == BeginOp );
723        CPPAD_ASSERT_NARG_NRES(BeginOp, 0, 1);
724        tape[i_var].new_var = rec->PutOp(BeginOp);
725
726        // temporary buffer for new argument values
727        size_t new_arg[6];
728
729        while(op != EndOp)
730        {       // next op
731                play->next_forward(op, arg, i_op, i_var);
732                CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
733                CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
734
735                // determine if we should keep this operator
736                bool keep;
737                switch( op )
738                {       case ComOp:
739                        case PripOp:
740                        case PrivOp:
741                        keep = false;
742                        break;
743
744                        case InvOp:
745                        case EndOp:
746                        keep = true;
747                        break;
748
749                        case StppOp:
750                        case StvpOp:
751                        case StpvOp:
752                        case StvvOp:
753                        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
754                        i = vecad[ arg[0] - 1 ];
755                        keep = vecad_parrent[i];
756                        break;
757
758                        default:
759                        keep = (tape[i_var].parrent > 0);
760                        break;
761                }
762
763                // start check if get a match in the hash table
764                unsigned short code      = hash_code(
765                        EndOp               , 
766                        arg                 ,
767                        play->num_rec_par() ,
768                        play->GetPar()
769                );
770                bool           match     = false;
771                size_t         match_var;
772                size_t         hash_var;
773                const size_t*  hash_arg;
774                Base           hash_par;
775                Base           par;
776                if( keep ) switch( op )
777                {
778                        // Unary operator where operand is arg[0]
779                        case AbsOp:
780                        case AcosOp:
781                        case AsinOp:
782                        case AtanOp:
783                        case CosOp:
784                        case CoshOp:
785                        case DisOp:
786                        case ExpOp:
787                        case LogOp:
788                        case SinOp:
789                        case SinhOp:
790                        case SqrtOp:
791                        match_var = optimize_unary_match(
792                                i_var               ,  // inputs
793                                op                  ,
794                                arg                 ,
795                                play->num_rec_par() ,
796                                play->GetPar()      ,
797                                hash_table_var      ,
798                                tape                , 
799                                code                  // outputs
800                        );
801                        if( match_var > 0 )
802                                tape[i_var].new_var = match_var;
803                        else
804                        {       new_arg[0] = tape[ arg[0] ].new_var;
805                                rec->PutArg( new_arg[0] );
806                                i                   = rec->PutOp(op);
807                                tape[i_var].new_var = i;
808                                CPPAD_ASSERT_UNKNOWN( new_arg[0] < i );
809                        }
810                        break;
811                        // ---------------------------------------------------
812                        // Non-commutative binary operators where
813                        // left is a variable and right is a parameter
814                        case DivvpOp:
815                        case PowvpOp:
816                        case SubvpOp:
817                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
818                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
819                        new_arg[0] = tape[arg[0]].new_var;
820                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_var );
821                        new_arg[1] = arg[1];
822                        code      = hash_code(
823                                op                  , 
824                                new_arg             ,
825                                play->num_rec_par() ,
826                                play->GetPar()
827                        );
828                        hash_var   = hash_table_var[code];
829                        hash_arg   = tape[hash_var].arg;
830                        //
831                        par        = play->GetPar( arg[1] );
832                        if( op == tape[hash_var].op )
833                        {       match   = (
834                                        new_arg[0] == tape[hash_arg[0]].new_var
835                                );
836                                hash_par= play->GetPar( hash_arg[1] );
837                                match  &= IdenticalEqualPar(par, hash_par); 
838                        }
839                        if( match )
840                                tape[i_var].new_var = tape[hash_var].new_var;
841                        else
842                        {
843                                new_arg[1] = rec->PutPar(par);
844                                rec->PutArg( new_arg[0], new_arg[1] );
845                                tape[i_var].new_var = rec->PutOp(op);
846                        }
847                        break;
848                        // ---------------------------------------------------
849                        // Non-commutative binary operators where
850                        // left is a parameter and right is a variable
851                        case DivpvOp:
852                        case PowpvOp:
853                        case SubpvOp:
854                        match_var = optimize_par_op_var_match(
855                                i_var               ,  // inputs
856                                op                  ,
857                                arg                 ,
858                                play->num_rec_par() ,
859                                play->GetPar()      ,
860                                hash_table_var      ,
861                                tape                , 
862                                code                  // outputs
863                        );
864                        if( match_var > 0 )
865                                tape[i_var].new_var = match_var;
866                        else
867                        {       new_arg[0] = rec->PutPar(
868                                        play->GetPar( arg[0] )
869                                );
870                                new_arg[1] = tape[ arg[1] ].new_var;
871                                rec->PutArg( new_arg[0], new_arg[1] );
872                                i                   = rec->PutOp(op);
873                                tape[i_var].new_var = i;
874                                CPPAD_ASSERT_UNKNOWN( new_arg[1] < i );
875                        }
876                        break;
877                        // ---------------------------------------------------
878                        // Non-commutative binary operator where
879                        // both operators are variables
880                        case DivvvOp:
881                        case PowvvOp:
882                        case SubvvOp:
883                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
884                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
885                        new_arg[0] = tape[arg[0]].new_var;
886                        new_arg[1] = tape[arg[1]].new_var;
887                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_var );
888                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
889                        code       = hash_code(
890                                op                  , 
891                                new_arg             ,
892                                play->num_rec_par() ,
893                                play->GetPar()
894                        );
895                        hash_var   = hash_table_var[code];
896                        hash_arg   = tape[hash_var].arg;
897                        //
898                        if( op == tape[hash_var].op )
899                        {       match  = (
900                                        new_arg[0] == tape[hash_arg[0]].new_var
901                                );
902                                match &= (
903                                        new_arg[1] == tape[hash_arg[1]].new_var
904                                );
905                        }
906                        if( match )
907                                tape[i_var].new_var = tape[hash_var].new_var;
908                        else
909                        {       rec->PutArg( new_arg[0], new_arg[1] );
910                                tape[i_var].new_var = rec->PutOp(op);
911                        }
912                        break;
913                        // ---------------------------------------------------
914                        // Commutative binary operators where
915                        // left is a variable and right is a parameter
916                        case AddvpOp:
917                        case MulvpOp:
918                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
919                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
920                        new_arg[0] = tape[arg[0]].new_var;
921                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_var );
922                        new_arg[1] = arg[1];
923                        code       = hash_code(
924                                op                  , 
925                                new_arg             ,
926                                play->num_rec_par() ,
927                                play->GetPar()
928                        );
929                        hash_var   = hash_table_var[code];
930                        hash_arg   = tape[hash_var].arg;
931                        //
932                        par        = play->GetPar( arg[1] );
933                        if( op == tape[hash_var].op )
934                        {       match   = (
935                                        new_arg[0] == tape[hash_arg[0]].new_var
936                                );
937                                hash_par= play->GetPar( hash_arg[1] );
938                                match  &= IdenticalEqualPar(par, hash_par); 
939                        }
940                        if(! match )
941                        {       OpCode tmp_op = AddpvOp;
942                                if(op == MulvpOp) 
943                                        tmp_op = MulpvOp;
944                                size_t tmp_arg[2];
945                                tmp_arg[0] = arg[1];
946                                tmp_arg[1] = arg[0];
947                                unsigned short tmp_code;
948                                match_var = optimize_par_op_var_match(
949                                        i_var               ,  // inputs
950                                        tmp_op              ,
951                                        tmp_arg             ,
952                                        play->num_rec_par() ,
953                                        play->GetPar()      ,
954                                        hash_table_var      ,
955                                        tape                , 
956                                        tmp_code              // outputs
957                                );
958                                hash_var = hash_table_var[tmp_code];
959                                match    = match_var > 0;
960                        }
961                        if( match )
962                                tape[i_var].new_var = tape[hash_var].new_var;
963                        else
964                        {
965                                new_arg[1] = rec->PutPar(par);
966                                rec->PutArg( new_arg[0], new_arg[1] );
967                                tape[i_var].new_var = rec->PutOp(op);
968                        }
969                        break;
970                        // ---------------------------------------------------
971                        // Commutative binary operators where
972                        // left is a parameter and right is a variable
973                        case AddpvOp:
974                        case MulpvOp:
975                        match_var = optimize_par_op_var_match(
976                                i_var               ,  // inputs
977                                op                  ,
978                                arg                 ,
979                                play->num_rec_par() ,
980                                play->GetPar()      ,
981                                hash_table_var      ,
982                                tape                , 
983                                code                  // outputs
984                        );
985                        if( match_var == 0 )
986                        {       OpCode tmp_op = AddvpOp;
987                                if(op == MulpvOp) 
988                                        tmp_op = MulvpOp;
989                                size_t tmp_arg[2];
990                                tmp_arg[0] = tape[arg[1]].new_var;
991                                tmp_arg[1] = arg[0];
992                                unsigned short tmp_code = hash_code( 
993                                        tmp_op              , 
994                                        tmp_arg             ,
995                                        play->num_rec_par() ,
996                                        play->GetPar()
997                                );
998                                hash_var = hash_table_var[tmp_code];
999                                if( tmp_op == tape[hash_var].op )
1000                                {
1001                                hash_arg = tape[hash_var].arg;
1002                                match    = (
1003                                        tmp_arg[0] == tape[hash_arg[0]].new_var
1004                                );
1005                                par      = play->GetPar( arg[0] );
1006                                hash_par = play->GetPar( hash_arg[1] );
1007                                match   &= IdenticalEqualPar(par, hash_par);
1008                                }
1009                                if( match )
1010                                        match_var = tape[hash_var].new_var;
1011                                       
1012                        }
1013                        if( match_var > 0 )
1014                                tape[i_var].new_var = match_var;
1015                        else
1016                        {
1017                                new_arg[1] = tape[arg[1]].new_var;
1018                                new_arg[0] = rec->PutPar(
1019                                        play->GetPar( arg[0] )
1020                                );
1021                                rec->PutArg( new_arg[0], new_arg[1] );
1022                                i                   = rec->PutOp(op);
1023                                tape[i_var].new_var = i;
1024                                CPPAD_ASSERT_UNKNOWN( new_arg[1] < i );
1025                        }
1026                        break;
1027                        // ---------------------------------------------------
1028                        // Commutative binary operator where
1029                        // both operators are variables
1030                        case AddvvOp:
1031                        case MulvvOp:
1032                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
1033                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
1034                        new_arg[0] = tape[arg[0]].new_var;
1035                        new_arg[1] = tape[arg[1]].new_var;
1036                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_var );
1037                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1038                        code       = hash_code(
1039                                op                  , 
1040                                new_arg             ,
1041                                play->num_rec_par() ,
1042                                play->GetPar()
1043                        );
1044                        hash_var   = hash_table_var[code];
1045                        hash_arg   = tape[hash_var].arg;
1046                        //
1047                        if( op == tape[hash_var].op )
1048                        {       match   = (
1049                                        new_arg[0] == tape[hash_arg[0]].new_var
1050                                );
1051                                match  &= (
1052                                        new_arg[1] == tape[hash_arg[1]].new_var
1053                                );
1054                        }
1055                        if(! match )
1056                        {       size_t tmp_arg[2];
1057                                tmp_arg[0] = new_arg[1];
1058                                tmp_arg[1] = new_arg[0];
1059                                unsigned short tmp_code = hash_code(
1060                                        op                  , 
1061                                        tmp_arg             ,
1062                                        play->num_rec_par() ,
1063                                        play->GetPar()
1064                                );
1065                                hash_var = hash_table_var[tmp_code];
1066                                if( op == tape[hash_var].op )
1067                                {
1068                                hash_arg = tape[hash_var].arg;
1069                                match    = (
1070                                        new_arg[0] == tape[hash_arg[1]].new_var
1071                                );
1072                                match   &= (
1073                                        new_arg[1] == tape[hash_arg[0]].new_var
1074                                );
1075                                }
1076                        }
1077                        if( match )
1078                                tape[i_var].new_var = tape[hash_var].new_var;
1079                        else
1080                        {
1081                                rec->PutArg( new_arg[0], new_arg[1] );
1082                                tape[i_var].new_var = rec->PutOp(op);
1083                        }
1084                        break;
1085                        // ---------------------------------------------------
1086                        // Conditional expression operators
1087                        case CExpOp:
1088                        CPPAD_ASSERT_NARG_NRES(op, 6, 1);
1089                        new_arg[0] = arg[0];
1090                        new_arg[1] = arg[1];
1091                        mask = 1;
1092                        for(i = 2; i < 6; i++)
1093                        {       if( arg[1] & mask )
1094                                {       new_arg[i] = tape[arg[i]].new_var;
1095                                        CPPAD_ASSERT_UNKNOWN( 
1096                                                new_arg[i] < num_var
1097                                        );
1098                                }
1099                                else    new_arg[i] = rec->PutPar( 
1100                                                play->GetPar( arg[i] )
1101                                );
1102                                mask = mask << 1;
1103                        }
1104                        rec->PutArg(
1105                                new_arg[0] ,
1106                                new_arg[1] ,
1107                                new_arg[2] ,
1108                                new_arg[3] ,
1109                                new_arg[4] ,
1110                                new_arg[5] 
1111                        );
1112                        tape[i_var].new_var = rec->PutOp(op);
1113                        break;
1114                        // ---------------------------------------------------
1115                        // Operations with no arguments and no results
1116                        case EndOp:
1117                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
1118                        rec->PutOp(op);
1119                        break;
1120                        // ---------------------------------------------------
1121                        // Operations with no arguments and one result
1122                        case InvOp:
1123                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
1124                        tape[i_var].new_var = rec->PutOp(op);
1125                        break;
1126                        // ---------------------------------------------------
1127                        // Operations with one argument that is a parameter
1128                        case ParOp:
1129                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
1130                        new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
1131
1132                        rec->PutArg( new_arg[0] );
1133                        tape[i_var].new_var = rec->PutOp(op);
1134                        break;
1135                        // ---------------------------------------------------
1136                        // Load using a parameter index
1137                        case LdpOp:
1138                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1139                        new_arg[0] = new_vecad_ind[ arg[0] ];
1140                        new_arg[1] = arg[1];
1141                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1142                        rec->PutArg( 
1143                                new_arg[0], 
1144                                new_arg[1], 
1145                                0
1146                        );
1147                        tape[i_var].new_var = rec->PutOp(op);
1148                        break;
1149                        // ---------------------------------------------------
1150                        // Load using a variable index
1151                        case LdvOp:
1152                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1153                        new_arg[0] = new_vecad_ind[ arg[0] ];
1154                        new_arg[1] = tape[arg[1]].new_var;
1155                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1156                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1157                        rec->PutArg( 
1158                                new_arg[0], 
1159                                new_arg[1], 
1160                                0
1161                        );
1162                        tape[i_var].new_var = rec->PutOp(op);
1163                        break;
1164                        // ---------------------------------------------------
1165                        // Store a parameter using a parameter index
1166                        case StppOp:
1167                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1168                        new_arg[0] = new_vecad_ind[ arg[0] ];
1169                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1170                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1171                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1172                        rec->PutArg(
1173                                new_arg[0], 
1174                                new_arg[1], 
1175                                new_arg[2]
1176                        );
1177                        rec->PutOp(op);
1178                        break;
1179                        // ---------------------------------------------------
1180                        // Store a parameter using a variable index
1181                        case StvpOp:
1182                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1183                        new_arg[0] = new_vecad_ind[ arg[0] ];
1184                        new_arg[1] = tape[arg[1]].new_var;
1185                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1186                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1187                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1188                        rec->PutArg(
1189                                new_arg[0], 
1190                                new_arg[1], 
1191                                new_arg[2]
1192                        );
1193                        rec->PutOp(op);
1194                        break;
1195                        // ---------------------------------------------------
1196                        // Store a variable using a parameter index
1197                        case StpvOp:
1198                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1199                        new_arg[0] = new_vecad_ind[ arg[0] ];
1200                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1201                        new_arg[2] = tape[arg[2]].new_var;
1202                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1203                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1204                        CPPAD_ASSERT_UNKNOWN( new_arg[2] < num_var );
1205                        rec->PutArg(
1206                                new_arg[0], 
1207                                new_arg[1], 
1208                                new_arg[2]
1209                        );
1210                        rec->PutOp(op);
1211                        break;
1212                        // ---------------------------------------------------
1213                        // Store a variable using a variable index
1214                        case StvvOp:
1215                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1216                        new_arg[0] = new_vecad_ind[ arg[0] ];
1217                        new_arg[1] = tape[arg[1]].new_var;
1218                        new_arg[2] = tape[arg[2]].new_var;
1219                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1220                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1221                        CPPAD_ASSERT_UNKNOWN( new_arg[2] < num_var );
1222                        rec->PutArg(
1223                                new_arg[0], 
1224                                new_arg[1], 
1225                                new_arg[2]
1226                        );
1227                        rec->PutOp(op);
1228                        break;
1229                        // ---------------------------------------------------
1230                        // all cases should be handled above
1231                        default:
1232                        CPPAD_ASSERT_UNKNOWN(false);
1233
1234                }
1235                if( keep & (! match) & (NumRes(op) > 0) )
1236                {       // put most recent match for this code in hash table
1237                        hash_table_var[code] = i_var;
1238                }
1239
1240        }
1241        // modify the dependent variable vector to new indices
1242        for(i = 0; i < dep_taddr.size(); i++ )
1243        {       CPPAD_ASSERT_UNKNOWN( tape[ dep_taddr[i] ].new_var < num_var );
1244                dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
1245        }
1246}
1247
1248/*!
1249Optimize a player object operation sequence
1250
1251The operation sequence for this object is replaced by one with fewer operations
1252but the same funcition and derivative values.
1253
1254\tparam Base
1255base type for the operator; i.e., this operation was recorded
1256using AD< \a Base > and computations by this routine are done using type
1257\a Base.
1258
1259*/
1260template <class Base>
1261void ADFun<Base>::optimize(void)
1262{       // place to store the optimized version of the recording
1263        recorder<Base> rec;
1264
1265        // number of independent variables
1266        size_t n = ind_taddr_.size();
1267
1268# ifndef NDEBUG
1269        size_t i, j, m = dep_taddr_.size();
1270        CppAD::vector<Base> x(n), y(m), check(m);
1271        bool check_zero_order = taylor_per_var_ > 0;
1272        if( check_zero_order )
1273        {       // zero order coefficients for independent vars
1274                for(j = 0; j < n; j++)
1275                {       CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
1276                        CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
1277                        x[j] = taylor_[ ind_taddr_[j] * taylor_col_dim_ + 0];
1278                }
1279                // zero order coefficients for dependent vars
1280                for(i = 0; i < m; i++)
1281                {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
1282                        y[i] = taylor_[ dep_taddr_[i] * taylor_col_dim_ + 0];
1283                }
1284        }
1285# endif
1286
1287        // create the optimized recording
1288        CppAD::optimize<Base>(n, dep_taddr_, &play_, &rec);
1289
1290        // now replace the recording
1291        play_ = rec;
1292
1293        // number of variables in the recording
1294        total_num_var_ = rec.num_rec_var();
1295
1296        // free memory allocated for sparse Jacobian calculation
1297        // (the results are no longer valid)
1298        for_jac_sparse_pack_.resize(0, 0);
1299        for_jac_sparse_set_.resize(0,0);
1300
1301        // free old Taylor coefficient memory
1302        if( taylor_ != CPPAD_NULL )
1303                CPPAD_TRACK_DEL_VEC(taylor_);
1304        taylor_         = CPPAD_NULL;
1305        taylor_per_var_ = 0;
1306        taylor_col_dim_ = 0;
1307
1308# ifndef NDEBUG
1309        if( check_zero_order )
1310        {
1311                // zero order forward calculation using new operation sequence
1312                check = Forward(0, x);
1313
1314                // check results
1315                for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
1316                        check[i] == y[i] ,
1317                        "Error during check of f.optimize()."
1318                );
1319
1320                // Erase memory that this calculation was done so NDEBUG gives
1321                // same final state for this object (from users perspective)
1322                taylor_per_var_ = 0;
1323        }
1324# endif
1325}
1326
1327CPPAD_END_NAMESPACE
1328# endif
Note: See TracBrowser for help on using the repository browser.