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

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

trunk: refactor and improve hash coding of operations.

whats_new_09.omh: user's view of the changes.
optimize.hpp: refactor unary operator optimization into a separate function.
hash_code.hpp: hash code parameter value instead of its argument index.

  • Property svn:keywords set to Id
File size: 33.7 KB
Line 
1/* $Id: optimize.hpp 1599 2009-12-04 14:01:50Z 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\param new_arg
239The input value of \c new_arg[0] does not matter.
240The output value of \c new_arg[0] is the argument for this unary
241operator in the new operation sequence.
242
243\return
244If the return value is zero,
245no match was found.
246If the return value is greater than zero,
247it is the index of a new variable that can be used to replace the
248old variable.
249
250\par Check Assertions
251\li
252*/
253template <class Base>
254inline size_t optimize_unary_match(
255        size_t                                         old_res        ,
256        OpCode                                         op             ,
257        const size_t*                                  old_arg        ,
258        size_t                                         npar           ,
259        const Base*                                    par            ,
260        const CppAD::vector<size_t>&                   hash_table_var ,
261        const CppAD::vector<struct optimize_variable>& tape           ,
262        unsigned short&                                code           ,
263        size_t*                                        new_arg        )
264{       
265        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
266        CPPAD_ASSERT_UNKNOWN( NumRes(op) > 0  );
267        CPPAD_ASSERT_UNKNOWN( old_arg[0] < old_res );
268        new_arg[0] = tape[old_arg[0]].new_var;
269        CPPAD_ASSERT_UNKNOWN( new_arg[0] < old_res );
270        code = hash_code(
271                op                  , 
272                new_arg             ,
273                npar                ,
274                par
275        );
276        size_t  i               = hash_table_var[code];
277        CPPAD_ASSERT_UNKNOWN( i < old_res );
278        if( op == tape[i].op )
279        {       size_t k = tape[i].arg[0];
280                CPPAD_ASSERT_UNKNOWN( k < i );
281                if (new_arg[0] == tape[k].new_var );
282                        return tape[i].new_var;
283        }
284        return 0;
285} 
286
287/*!
288Convert a player object to an optimized recorder object
289
290\tparam Base
291base type for the operator; i.e., this operation was recorded
292using AD< \a Base > and computations by this routine are done using type
293\a Base.
294
295\param n
296is the number of independent variables on the tape.
297
298\param dep_taddr
299On input this vector contains the indices for each of the dependent
300variable values in the operation sequence corresponding to \a play.
301Upon return it contains the indices for the same variables but in
302the operation sequence corresponding to \a rec.
303
304\param play
305This is the operation sequence that we are optimizing.
306It is essentially const, except for play back state which
307changes while it plays back the operation seqeunce.
308
309\param rec
310The input contents of this recording does not matter.
311Upon return, it contains an optimized verison of the
312operation sequence corresponding to \a play.
313*/
314
315template <class Base>
316void optimize(
317        size_t                       n         ,
318        CppAD::vector<size_t>&       dep_taddr ,
319        player<Base>*                play      ,
320        recorder<Base>*              rec       ) 
321{
322        // temporary indices
323        size_t i, j, k;
324
325        // temporary variables
326        OpCode        op;   // current operator
327        const size_t *arg;  // operator arguments
328        size_t        i_var;  // index of first result for current operator
329
330        // range and domain dimensions for F
331        size_t m = dep_taddr.size();
332
333        // number of variables in the player
334        const size_t num_var = play->num_rec_var(); 
335
336        // number of  VecAD indices
337        size_t num_vecad_ind   = play->num_rec_vecad_ind();
338
339        // number of VecAD vectors
340        size_t num_vecad_vec   = play->num_rec_vecad_vec();
341
342        // -------------------------------------------------------------
343        // data structure that maps variable index in original operation
344        // sequence to corresponding operator information
345        CppAD::vector<struct optimize_variable> tape(num_var);
346        // -------------------------------------------------------------
347        // Determine parrent value for each variable
348
349        // initialize all variables has having no parrent
350        for(i = 0; i < num_var; i++)
351                tape[i].parrent = 0;
352
353        for(j = 0; j < m; j++)
354        {       // mark dependent variables as having multiple parents
355                tape[ dep_taddr[j] ].parrent = num_var;
356        }
357
358        // vecad_parrent contains a flag for each VecAD object.
359        // vecad maps a VecAD index (which corresponds to the beginning of the
360        // VecAD object) to the vecad_parrent falg for the VecAD object.
361# if CPPAD_OPTIMIZE_TRACE
362        CppAD::vector<size_t> vecad_offset(num_vecad_vec);
363# endif
364        CppAD::vector<bool>   vecad_parrent(num_vecad_vec);
365        CppAD::vector<size_t> vecad(num_vecad_ind);
366        j = 0;
367        for(i = 0; i < num_vecad_vec; i++)
368        {       vecad_parrent[i] = false;
369# if CPPAD_OPTIMIZE_TRACE
370                // offset for this VecAD object
371                vecad_offset[i]  = j + 1;
372# endif
373                // length of this VecAD
374                size_t length = play->GetVecInd(j);
375                // set to proper index for this VecAD
376                vecad[j] = i; 
377                for(k = 1; k <= length; k++)
378                        vecad[j+k] = num_vecad_vec; // invalid index
379                // start of next VecAD
380                j       += length + 1;
381        }
382        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
383
384        // Initialize a reverse mode sweep through the operation sequence
385        size_t i_op;
386        play->start_reverse(op, arg, i_op, i_var);
387        CPPAD_ASSERT_UNKNOWN( op == EndOp );
388        size_t mask;
389# if CPPAD_OPTIMIZE_TRACE
390        std::cout << std::endl;
391# endif
392        while(op != BeginOp)
393        {       // next op
394                play->next_reverse(op, arg, i_op, i_var);
395                // This if is not necessary becasue last assignment
396                // with this value of i_var will have NumRes(op) > 0
397                if( NumRes(op) > 0 )
398                {       tape[i_var].op = op;
399                        tape[i_var].arg = arg;
400                }
401# ifndef NDEBUG
402                if( i_op <= n )
403                {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
404                }
405                else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
406# endif
407
408# if CPPAD_OPTIMIZE_TRACE
409                printOp(
410                        std::cout, 
411                        play,
412                        i_var,
413                        op, 
414                        arg,
415                        0, 
416                        (bool *) CPPAD_NULL,
417                        1, 
418                        & connected[i_var]
419                );
420# endif
421
422                switch( op )
423                {
424                        // Unary operator where operand is arg[0]
425                        case AbsOp:
426                        case AddvpOp:
427                        case AcosOp:
428                        case AsinOp:
429                        case AtanOp:
430                        case CosOp:
431                        case CoshOp:
432                        case DisOp:
433                        case DivvpOp:
434                        case ExpOp:
435                        case LogOp:
436                        case MulvpOp:
437                        case PowvpOp:
438                        case SinOp:
439                        case SinhOp:
440                        case SqrtOp:
441                        case SubvpOp:
442                        if( tape[i_var].parrent > 0 )
443                        {       if( tape[arg[0]].parrent == 0 )
444                                        tape[arg[0]].parrent = i_var;
445                                else    tape[arg[0]].parrent = num_var;
446                        }
447                        break;
448
449                        // Unary operator where operand is arg[1]
450                        case AddpvOp:
451                        case DivpvOp:
452                        case MulpvOp:
453                        case SubpvOp:
454                        case PowpvOp:
455                        case PrivOp:
456                        if( tape[i_var].parrent > 0 )
457                        {       if( tape[arg[1]].parrent == 0 )
458                                        tape[arg[1]].parrent = i_var;
459                                else    tape[arg[1]].parrent = num_var;
460                        }
461                        break;
462
463                        // Binary operator where operands are arg[0], arg[1]
464                        case AddvvOp:
465                        case DivvvOp:
466                        case MulvvOp:
467                        case PowvvOp:
468                        case SubvvOp:
469                        if( tape[i_var].parrent > 0 )
470                        {       if( tape[arg[0]].parrent == 0 )
471                                        tape[arg[0]].parrent = i_var;
472                                else    tape[arg[0]].parrent = num_var;
473                                if( tape[arg[1]].parrent == 0 )
474                                        tape[arg[1]].parrent = i_var;
475                                else    tape[arg[1]].parrent = num_var;
476                        }
477                        break;
478
479                        // Conditional expression operators
480                        case CExpOp:
481                        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
482                        if( tape[i_var].parrent > 0 )
483                        {       mask = 1;
484                                for(i = 2; i < 6; i++) if( arg[1] & mask )
485                                {       CPPAD_ASSERT_UNKNOWN( arg[i] < i_var );
486                                        if( tape[arg[i]].parrent == 0 )
487                                                tape[arg[i]].parrent = i_var;
488                                        else    tape[arg[i]].parrent = num_var;
489                                }
490                                mask = mask << 1;
491                        }
492
493                        // Operations where there is noting to do
494                        case BeginOp:
495                        case ComOp:
496                        case EndOp:
497                        case InvOp:
498                        case ParOp:
499                        case PripOp:
500                        case StppOp:
501                        break; 
502
503                        // Load using a parameter index
504                        case LdpOp:
505                        if( tape[i_var].parrent > 0 )
506                        {       i = vecad[ arg[0] - 1 ];
507                                vecad_parrent[i] = true;
508                        }
509                        break;
510
511                        // Load using a variable index
512                        case LdvOp:
513                        if( tape[i_var].parrent > 0 )
514                        {       i = vecad[ arg[0] - 1 ];
515                                vecad_parrent[i] = true;
516                                if( tape[arg[1]].parrent == 0 )
517                                        tape[arg[1]].parrent = i_var;
518                                else    tape[arg[1]].parrent = num_var;
519                        }
520                        break;
521
522                        // Store a variable using a parameter index
523                        case StpvOp:
524                        i = vecad[ arg[0] - 1 ];
525                        if( vecad_parrent[i] )
526                        {       if( tape[arg[2]].parrent == 0 )
527                                        tape[arg[2]].parrent = i_var;
528                                else    tape[arg[2]].parrent = num_var;
529                        }
530                        break;
531
532                        // Store a variable using a variable index
533                        case StvvOp:
534                        i = vecad[ arg[0] - 1 ];
535                        if( vecad_parrent[i] )
536                        {       if( tape[arg[1]].parrent == 0 )
537                                        tape[arg[1]].parrent = i_var;
538                                else    tape[arg[1]].parrent = num_var;
539                                if( tape[arg[2]].parrent == 0 )
540                                        tape[arg[2]].parrent = i_var;
541                                else    tape[arg[2]].parrent = num_var;
542                        }
543                        break;
544
545                        // all cases should be handled above
546                        default:
547                        CPPAD_ASSERT_UNKNOWN(0);
548                }
549        }
550        // values corresponding to BeginOp
551        CPPAD_ASSERT_UNKNOWN( i_op == 0 && i_var == 0 && op == BeginOp );
552        tape[i_var].op = op;
553# if CPPAD_OPTIMIZE_TRACE
554        std::cout << "VecAD information:" << std::endl;
555        for(i = 0; i < num_vecad_vec; i++) 
556        {       std::cout << "offset  = " << vecad_offset[i];
557                std::cout << ", parrent = " << vecad_parrent[i] << std::endl;
558        }
559# endif
560        // -------------------------------------------------------------
561
562        // Erase all information in the recording
563        rec->Erase();
564
565        // Initilaize table mapping hash code to variable index in tape
566        // as pointing to the BeginOp at the beginning of the tape
567        CppAD::vector<size_t>  hash_table_var(CPPAD_HASH_TABLE_SIZE);
568        for(i = 0; i < CPPAD_HASH_TABLE_SIZE; i++)
569                hash_table_var[i] = 0;
570        CPPAD_ASSERT_UNKNOWN( tape[0].op == BeginOp );
571
572        // initialize mapping from old variable index to new variable index
573        for(i = 0; i < num_var; i++)
574                tape[i].new_var = num_var; // invalid index
575       
576
577        // initialize mapping from old VecAD index to new VecAD index
578        CppAD::vector<size_t> new_vecad_ind(num_vecad_ind);
579        for(i = 0; i < num_vecad_ind; i++)
580                new_vecad_ind[i] = num_vecad_ind; // invalid index
581
582        j = 0;     // index into the old set of indices
583        for(i = 0; i < num_vecad_vec; i++)
584        {       // length of this VecAD
585                size_t length = play->GetVecInd(j);
586                if( vecad_parrent[i] )
587                {       // Put this VecAD vector in new recording
588                        CPPAD_ASSERT_UNKNOWN(length < num_vecad_ind);
589                        new_vecad_ind[j] = rec->PutVecInd(length);
590                        for(k = 1; k <= length; k++) new_vecad_ind[j+k] =
591                                rec->PutVecInd(
592                                        rec->PutPar(
593                                                play->GetPar( 
594                                                        play->GetVecInd(j+k)
595                        ) ) );
596                }
597                // start of next VecAD
598                j       += length + 1;
599        }
600        CPPAD_ASSERT_UNKNOWN( j == num_vecad_ind );
601
602        // start playing the operations in the forward direction
603        play->start_forward(op, arg, i_op, i_var);
604
605        // playing forward skips BeginOp at the beginning, but not EndOp at
606        // the end.  Put BeginOp at beginning of recording
607        CPPAD_ASSERT_UNKNOWN( op == BeginOp );
608        CPPAD_ASSERT_NARG_NRES(BeginOp, 0, 1);
609        tape[i_var].new_var = rec->PutOp(BeginOp);
610
611        // temporary buffer for new argument values
612        size_t new_arg[6];
613
614        while(op != EndOp)
615        {       // next op
616                play->next_forward(op, arg, i_op, i_var);
617                CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
618                CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
619
620                // determine if we should keep this operator
621                bool keep;
622                switch( op )
623                {       case ComOp:
624                        case PripOp:
625                        case PrivOp:
626                        keep = false;
627                        break;
628
629                        case InvOp:
630                        case EndOp:
631                        keep = true;
632                        break;
633
634                        case StppOp:
635                        case StvpOp:
636                        case StpvOp:
637                        case StvvOp:
638                        CPPAD_ASSERT_UNKNOWN( NumRes(op) == 0 );
639                        i = vecad[ arg[0] - 1 ];
640                        keep = vecad_parrent[i];
641                        break;
642
643                        default:
644                        keep = (tape[i_var].parrent > 0);
645                        break;
646                }
647
648                // start check if get a match in the hash table
649                unsigned short code      = hash_code(
650                        EndOp               , 
651                        arg                 ,
652                        play->num_rec_par() ,
653                        play->GetPar()
654                );
655                bool           match     = false;
656                size_t         match_var;
657                size_t         hash_var;
658                const size_t*  hash_arg;
659                Base           hash_par;
660                Base           par;
661                if( keep ) switch( op )
662                {
663                        // Unary operator where operand is arg[0]
664                        case AbsOp:
665                        case AcosOp:
666                        case AsinOp:
667                        case AtanOp:
668                        case CosOp:
669                        case CoshOp:
670                        case DisOp:
671                        case ExpOp:
672                        case LogOp:
673                        case SinOp:
674                        case SinhOp:
675                        case SqrtOp:
676                        match_var = optimize_unary_match(
677                                i_var               ,  // inputs
678                                op                  ,
679                                arg                 ,
680                                play->num_rec_par() ,
681                                play->GetPar()      ,
682                                hash_table_var      ,
683                                tape                , 
684                                code                , // outputs
685                                new_arg
686                        );
687                        if( match_var > 0 )
688                                tape[i_var].new_var = match_var;
689                        else
690                        {       rec->PutArg( new_arg[0] );
691                                tape[i_var].new_var = rec->PutOp(op);
692                        }
693                        break;
694                        // ---------------------------------------------------
695                        // Non-commutative binary operators where
696                        // left is a variable and right is a parameter
697                        case DivvpOp:
698                        case PowvpOp:
699                        case SubvpOp:
700                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
701                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
702                        new_arg[0] = tape[arg[0]].new_var;
703                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_var );
704                        new_arg[1] = arg[1];
705                        code      = hash_code(
706                                op                  , 
707                                new_arg             ,
708                                play->num_rec_par() ,
709                                play->GetPar()
710                        );
711                        hash_var   = hash_table_var[code];
712                        hash_arg   = tape[hash_var].arg;
713                        //
714                        par        = play->GetPar( arg[1] );
715                        if( op == tape[hash_var].op )
716                        {       match   = (
717                                        new_arg[0] == tape[hash_arg[0]].new_var
718                                );
719                                hash_par= play->GetPar( hash_arg[1] );
720                                match  &= IdenticalEqualPar(par, hash_par); 
721                        }
722                        if( match )
723                                tape[i_var].new_var = tape[hash_var].new_var;
724                        else
725                        {
726                                new_arg[1] = rec->PutPar(par);
727                                rec->PutArg( new_arg[0], new_arg[1] );
728                                tape[i_var].new_var = rec->PutOp(op);
729                        }
730                        break;
731                        // ---------------------------------------------------
732                        // Non-commutative binary operators where
733                        // left is a parameter and right is a variable
734                        case DivpvOp:
735                        case PowpvOp:
736                        case SubpvOp:
737                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
738                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
739                        new_arg[1] = tape[arg[1]].new_var;
740                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
741                        new_arg[0] = arg[0];
742                        code       = hash_code(
743                                op                  , 
744                                new_arg             ,
745                                play->num_rec_par() ,
746                                play->GetPar()
747                        );
748                        hash_var   = hash_table_var[code];
749                        hash_arg   = tape[hash_var].arg;
750                        //
751                        par        = play->GetPar( arg[0] );
752                        if( op == tape[hash_var].op )
753                        {       match   = (
754                                        new_arg[1] == tape[hash_arg[1]].new_var
755                                );
756                                hash_par= play->GetPar( hash_arg[0] );
757                                match  &= IdenticalEqualPar(par, hash_par); 
758                        }
759                        if( match )
760                                tape[i_var].new_var = tape[hash_var].new_var;
761                        else
762                        {
763                                new_arg[0] = rec->PutPar(par);
764                                rec->PutArg( new_arg[0], new_arg[1] );
765                                tape[i_var].new_var = rec->PutOp(op);
766                        }
767                        break;
768                        // ---------------------------------------------------
769                        // Non-commutative binary operator where
770                        // both operators are variables
771                        case DivvvOp:
772                        case PowvvOp:
773                        case SubvvOp:
774                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
775                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
776                        new_arg[0] = tape[arg[0]].new_var;
777                        new_arg[1] = tape[arg[1]].new_var;
778                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_var );
779                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
780                        code       = hash_code(
781                                op                  , 
782                                new_arg             ,
783                                play->num_rec_par() ,
784                                play->GetPar()
785                        );
786                        hash_var   = hash_table_var[code];
787                        hash_arg   = tape[hash_var].arg;
788                        //
789                        if( op == tape[hash_var].op )
790                        {       match  = (
791                                        new_arg[0] == tape[hash_arg[0]].new_var
792                                );
793                                match &= (
794                                        new_arg[1] == tape[hash_arg[1]].new_var
795                                );
796                        }
797                        if( match )
798                                tape[i_var].new_var = tape[hash_var].new_var;
799                        else
800                        {       rec->PutArg( new_arg[0], new_arg[1] );
801                                tape[i_var].new_var = rec->PutOp(op);
802                        }
803                        break;
804                        // ---------------------------------------------------
805                        // Commutative binary operators where
806                        // left is a variable and right is a parameter
807                        case AddvpOp:
808                        case MulvpOp:
809                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
810                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
811                        new_arg[0] = tape[arg[0]].new_var;
812                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_var );
813                        new_arg[1] = arg[1];
814                        code       = hash_code(
815                                op                  , 
816                                new_arg             ,
817                                play->num_rec_par() ,
818                                play->GetPar()
819                        );
820                        hash_var   = hash_table_var[code];
821                        hash_arg   = tape[hash_var].arg;
822                        //
823                        par        = play->GetPar( arg[1] );
824                        if( op == tape[hash_var].op )
825                        {       match   = (
826                                        new_arg[0] == tape[hash_arg[0]].new_var
827                                );
828                                hash_par= play->GetPar( hash_arg[1] );
829                                match  &= IdenticalEqualPar(par, hash_par); 
830                        }
831                        if(! match )
832                        {       OpCode tmp_op = AddpvOp;
833                                if(op == MulvpOp) 
834                                        tmp_op = MulpvOp;
835                                size_t tmp_arg[2];
836                                tmp_arg[0] = new_arg[1];
837                                tmp_arg[1] = new_arg[0];
838                                unsigned short tmp_code = hash_code(
839                                        tmp_op              , 
840                                        tmp_arg             ,
841                                        play->num_rec_par() ,
842                                        play->GetPar()
843                                );
844                                hash_var = hash_table_var[tmp_code];
845                                if( tmp_op == tape[hash_var].op )
846                                {
847                                hash_arg = tape[hash_var].arg;
848                                match    = (
849                                        new_arg[0] == tape[hash_arg[1]].new_var
850                                );
851                                hash_par = play->GetPar( hash_arg[0] );
852                                match   &= IdenticalEqualPar(par, hash_par);
853                                }
854                        }
855                        if( match )
856                                tape[i_var].new_var = tape[hash_var].new_var;
857                        else
858                        {
859                                new_arg[1] = rec->PutPar(par);
860                                rec->PutArg( new_arg[0], new_arg[1] );
861                                tape[i_var].new_var = rec->PutOp(op);
862                        }
863                        break;
864                        // ---------------------------------------------------
865                        // Commutative binary operators where
866                        // left is a parameter and right is a variable
867                        case AddpvOp:
868                        case MulpvOp:
869                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
870                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
871                        new_arg[1] = tape[arg[1]].new_var;
872                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
873                        new_arg[0] = arg[0];
874                        code       = hash_code(
875                                op                  , 
876                                new_arg             ,
877                                play->num_rec_par() ,
878                                play->GetPar()
879                        );
880                        hash_var   = hash_table_var[code];
881                        hash_arg   = tape[hash_var].arg;
882                        //
883                        par        = play->GetPar( arg[0] );
884                        if( op == tape[hash_var].op )
885                        {       match   = (
886                                        new_arg[1] == tape[hash_arg[1]].new_var
887                                );
888                                hash_par= play->GetPar( hash_arg[0] );
889                                match  &= IdenticalEqualPar(par, hash_par); 
890                        }
891                        if(! match )
892                        {       OpCode tmp_op = AddvpOp;
893                                if(op == MulpvOp) 
894                                        tmp_op = MulvpOp;
895                                size_t tmp_arg[2];
896                                tmp_arg[0] = new_arg[1];
897                                tmp_arg[1] = new_arg[0];
898                                unsigned short tmp_code = hash_code( 
899                                        tmp_op              , 
900                                        tmp_arg             ,
901                                        play->num_rec_par() ,
902                                        play->GetPar()
903                                );
904                                hash_var = hash_table_var[tmp_code];
905                                if( tmp_op == tape[hash_var].op )
906                                {
907                                hash_arg = tape[hash_var].arg;
908                                match    = (
909                                        new_arg[1] == tape[hash_arg[0]].new_var
910                                );
911                                hash_par = play->GetPar( hash_arg[1] );
912                                match   &= IdenticalEqualPar(par, hash_par);
913                                }
914                        }
915                        if( match )
916                                tape[i_var].new_var = tape[hash_var].new_var;
917                        else
918                        {
919                                new_arg[0] = rec->PutPar(par);
920                                rec->PutArg( new_arg[0], new_arg[1] );
921                                tape[i_var].new_var = rec->PutOp(op);
922                        }
923                        break;
924                        // ---------------------------------------------------
925                        // Commutative binary operator where
926                        // both operators are variables
927                        case AddvvOp:
928                        case MulvvOp:
929                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 2 );
930                        CPPAD_ASSERT_UNKNOWN( NumRes(op) >  0 );
931                        new_arg[0] = tape[arg[0]].new_var;
932                        new_arg[1] = tape[arg[1]].new_var;
933                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_var );
934                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
935                        code       = hash_code(
936                                op                  , 
937                                new_arg             ,
938                                play->num_rec_par() ,
939                                play->GetPar()
940                        );
941                        hash_var   = hash_table_var[code];
942                        hash_arg   = tape[hash_var].arg;
943                        //
944                        if( op == tape[hash_var].op )
945                        {       match   = (
946                                        new_arg[0] == tape[hash_arg[0]].new_var
947                                );
948                                match  &= (
949                                        new_arg[1] == tape[hash_arg[1]].new_var
950                                );
951                        }
952                        if(! match )
953                        {       size_t tmp_arg[2];
954                                tmp_arg[0] = new_arg[1];
955                                tmp_arg[1] = new_arg[0];
956                                unsigned short tmp_code = hash_code(
957                                        op                  , 
958                                        tmp_arg             ,
959                                        play->num_rec_par() ,
960                                        play->GetPar()
961                                );
962                                hash_var = hash_table_var[tmp_code];
963                                if( op == tape[hash_var].op )
964                                {
965                                hash_arg = tape[hash_var].arg;
966                                match    = (
967                                        new_arg[0] == tape[hash_arg[1]].new_var
968                                );
969                                match   &= (
970                                        new_arg[1] == tape[hash_arg[0]].new_var
971                                );
972                                }
973                        }
974                        if( match )
975                                tape[i_var].new_var = tape[hash_var].new_var;
976                        else
977                        {
978                                rec->PutArg( new_arg[0], new_arg[1] );
979                                tape[i_var].new_var = rec->PutOp(op);
980                        }
981                        break;
982                        // ---------------------------------------------------
983                        // Conditional expression operators
984                        case CExpOp:
985                        CPPAD_ASSERT_NARG_NRES(op, 6, 1);
986                        new_arg[0] = arg[0];
987                        new_arg[1] = arg[1];
988                        mask = 1;
989                        for(i = 2; i < 6; i++)
990                        {       if( arg[1] & mask )
991                                {       new_arg[i] = tape[arg[i]].new_var;
992                                        CPPAD_ASSERT_UNKNOWN( 
993                                                new_arg[i] < num_var
994                                        );
995                                }
996                                else    new_arg[i] = rec->PutPar( 
997                                                play->GetPar( arg[i] )
998                                );
999                                mask = mask << 1;
1000                        }
1001                        rec->PutArg(
1002                                new_arg[0] ,
1003                                new_arg[1] ,
1004                                new_arg[2] ,
1005                                new_arg[3] ,
1006                                new_arg[4] ,
1007                                new_arg[5] 
1008                        );
1009                        tape[i_var].new_var = rec->PutOp(op);
1010                        break;
1011                        // ---------------------------------------------------
1012                        // Operations with no arguments and no results
1013                        case EndOp:
1014                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
1015                        rec->PutOp(op);
1016                        break;
1017                        // ---------------------------------------------------
1018                        // Operations with no arguments and one result
1019                        case InvOp:
1020                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
1021                        tape[i_var].new_var = rec->PutOp(op);
1022                        break;
1023                        // ---------------------------------------------------
1024                        // Operations with one argument that is a parameter
1025                        case ParOp:
1026                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
1027                        new_arg[0] = rec->PutPar( play->GetPar(arg[0] ) );
1028
1029                        rec->PutArg( new_arg[0] );
1030                        tape[i_var].new_var = rec->PutOp(op);
1031                        break;
1032                        // ---------------------------------------------------
1033                        // Load using a parameter index
1034                        case LdpOp:
1035                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1036                        new_arg[0] = new_vecad_ind[ arg[0] ];
1037                        new_arg[1] = arg[1];
1038                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1039                        rec->PutArg( 
1040                                new_arg[0], 
1041                                new_arg[1], 
1042                                0
1043                        );
1044                        tape[i_var].new_var = rec->PutOp(op);
1045                        break;
1046                        // ---------------------------------------------------
1047                        // Load using a variable index
1048                        case LdvOp:
1049                        CPPAD_ASSERT_NARG_NRES(op, 3, 1);
1050                        new_arg[0] = new_vecad_ind[ arg[0] ];
1051                        new_arg[1] = tape[arg[1]].new_var;
1052                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1053                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1054                        rec->PutArg( 
1055                                new_arg[0], 
1056                                new_arg[1], 
1057                                0
1058                        );
1059                        tape[i_var].new_var = rec->PutOp(op);
1060                        break;
1061                        // ---------------------------------------------------
1062                        // Store a parameter using a parameter index
1063                        case StppOp:
1064                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1065                        new_arg[0] = new_vecad_ind[ arg[0] ];
1066                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1067                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1068                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1069                        rec->PutArg(
1070                                new_arg[0], 
1071                                new_arg[1], 
1072                                new_arg[2]
1073                        );
1074                        rec->PutOp(op);
1075                        break;
1076                        // ---------------------------------------------------
1077                        // Store a parameter using a variable index
1078                        case StvpOp:
1079                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1080                        new_arg[0] = new_vecad_ind[ arg[0] ];
1081                        new_arg[1] = tape[arg[1]].new_var;
1082                        new_arg[2] = rec->PutPar( play->GetPar(arg[2]) );
1083                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1084                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1085                        rec->PutArg(
1086                                new_arg[0], 
1087                                new_arg[1], 
1088                                new_arg[2]
1089                        );
1090                        rec->PutOp(op);
1091                        break;
1092                        // ---------------------------------------------------
1093                        // Store a variable using a parameter index
1094                        case StpvOp:
1095                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1096                        new_arg[0] = new_vecad_ind[ arg[0] ];
1097                        new_arg[1] = rec->PutPar( play->GetPar(arg[1]) );
1098                        new_arg[2] = tape[arg[2]].new_var;
1099                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1100                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1101                        CPPAD_ASSERT_UNKNOWN( new_arg[2] < num_var );
1102                        rec->PutArg(
1103                                new_arg[0], 
1104                                new_arg[1], 
1105                                new_arg[2]
1106                        );
1107                        rec->PutOp(op);
1108                        break;
1109                        // ---------------------------------------------------
1110                        // Store a variable using a variable index
1111                        case StvvOp:
1112                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
1113                        new_arg[0] = new_vecad_ind[ arg[0] ];
1114                        new_arg[1] = tape[arg[1]].new_var;
1115                        new_arg[2] = tape[arg[2]].new_var;
1116                        CPPAD_ASSERT_UNKNOWN( new_arg[0] < num_vecad_ind );
1117                        CPPAD_ASSERT_UNKNOWN( new_arg[1] < num_var );
1118                        CPPAD_ASSERT_UNKNOWN( new_arg[2] < num_var );
1119                        rec->PutArg(
1120                                new_arg[0], 
1121                                new_arg[1], 
1122                                new_arg[2]
1123                        );
1124                        rec->PutOp(op);
1125                        break;
1126                        // ---------------------------------------------------
1127                        // all cases should be handled above
1128                        default:
1129                        CPPAD_ASSERT_UNKNOWN(false);
1130
1131                }
1132                if( keep & (! match) & (NumRes(op) > 0) )
1133                {       // put most recent match for this code in hash table
1134                        hash_table_var[code] = i_var;
1135                }
1136
1137        }
1138        // modify the dependent variable vector to new indices
1139        for(i = 0; i < dep_taddr.size(); i++ )
1140        {       CPPAD_ASSERT_UNKNOWN( tape[ dep_taddr[i] ].new_var < num_var );
1141                dep_taddr[i] = tape[ dep_taddr[i] ].new_var;
1142        }
1143}
1144
1145/*!
1146Optimize a player object operation sequence
1147
1148The operation sequence for this object is replaced by one with fewer operations
1149but the same funcition and derivative values.
1150
1151\tparam Base
1152base type for the operator; i.e., this operation was recorded
1153using AD< \a Base > and computations by this routine are done using type
1154\a Base.
1155
1156*/
1157template <class Base>
1158void ADFun<Base>::optimize(void)
1159{       // place to store the optimized version of the recording
1160        recorder<Base> rec;
1161
1162        // number of independent variables
1163        size_t n = ind_taddr_.size();
1164
1165# ifndef NDEBUG
1166        size_t i, j, m = dep_taddr_.size();
1167        CppAD::vector<Base> x(n), y(m), check(m);
1168        bool check_zero_order = taylor_per_var_ > 0;
1169        if( check_zero_order )
1170        {       // zero order coefficients for independent vars
1171                for(j = 0; j < n; j++)
1172                {       CPPAD_ASSERT_UNKNOWN( play_.GetOp(j+1) == InvOp );
1173                        CPPAD_ASSERT_UNKNOWN( ind_taddr_[j]    == j+1   );
1174                        x[j] = taylor_[ ind_taddr_[j] * taylor_col_dim_ + 0];
1175                }
1176                // zero order coefficients for dependent vars
1177                for(i = 0; i < m; i++)
1178                {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < total_num_var_ );
1179                        y[i] = taylor_[ dep_taddr_[i] * taylor_col_dim_ + 0];
1180                }
1181        }
1182# endif
1183
1184        // create the optimized recording
1185        CppAD::optimize<Base>(n, dep_taddr_, &play_, &rec);
1186
1187        // now replace the recording
1188        play_ = rec;
1189
1190        // number of variables in the recording
1191        total_num_var_ = rec.num_rec_var();
1192
1193        // free memory allocated for sparse Jacobian calculation
1194        // (the results are no longer valid)
1195        for_jac_sparse_pack_.resize(0, 0);
1196        for_jac_sparse_set_.resize(0,0);
1197
1198        // free old Taylor coefficient memory
1199        if( taylor_ != CPPAD_NULL )
1200                CPPAD_TRACK_DEL_VEC(taylor_);
1201        taylor_         = CPPAD_NULL;
1202        taylor_per_var_ = 0;
1203        taylor_col_dim_ = 0;
1204
1205# ifndef NDEBUG
1206        if( check_zero_order )
1207        {
1208                // zero order forward calculation using new operation sequence
1209                check = Forward(0, x);
1210
1211                // check results
1212                for(i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 
1213                        check[i] == y[i] ,
1214                        "Error during check of f.optimize()."
1215                );
1216
1217                // Erase memory that this calculation was done so NDEBUG gives
1218                // same final state for this object (from users perspective)
1219                taylor_per_var_ = 0;
1220        }
1221# endif
1222}
1223
1224CPPAD_END_NAMESPACE
1225# endif
Note: See TracBrowser for help on using the repository browser.