source: trunk/cppad/local/cond_op.hpp @ 3301

Last change on this file since 3301 was 3301, checked in by bradbell, 6 years ago

merge in multiple forward direcitons from branches/forward_dir

  • Property svn:keywords set to Id
File size: 33.5 KB
Line 
1/* $Id: cond_op.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_COND_OP_INCLUDED
3# define CPPAD_COND_OP_INCLUDED
4/* --------------------------------------------------------------------------
5CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
6
7CppAD is distributed under multiple licenses. This distribution is under
8the terms of the
9                    Eclipse Public License Version 1.0.
10
11A copy of this license is included in the COPYING file of this distribution.
12Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
13-------------------------------------------------------------------------- */
14
15namespace CppAD { // BEGIN_CPPAD_NAMESPACE
16/*!
17\file cond_op.hpp
18Forward, reverse, and sparse operations for conditional expressions.
19*/
20
21/*!
22Shared documentation for conditional expressions (not called).
23
24<!-- define conditional_exp_op -->
25The C++ source code coresponding to this operation is
26\verbatim
27        z = CondExpRel(y_0, y_1, y_2, y_3)
28\endverbatim
29where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
30
31\tparam Base
32base type for the operator; i.e., this operation was recorded
33using AD< \a Base > and computations by this routine are done using type
34\a Base.
35
36\param i_z
37is the AD variable index corresponding to the variable z.
38
39\param arg
40\n
41\a arg[0]
42is static cast to size_t from the enum type
43\verbatim
44        enum CompareOp {
45                CompareLt,
46                CompareLe,
47                CompareEq,
48                CompareGe,
49                CompareGt,
50                CompareNe
51        }
52\endverbatim
53for this operation.
54Note that arg[0] cannot be equal to CompareNe.
55\n
56\n
57\a arg[1] & 1
58\n
59If this is zero, y_0 is a parameter. Otherwise it is a variable.
60\n
61\n
62\a arg[1] & 2
63\n
64If this is zero, y_1 is a parameter. Otherwise it is a variable.
65\n
66\n
67\a arg[1] & 4
68\n
69If this is zero, y_2 is a parameter. Otherwise it is a variable.
70\n
71\n
72\a arg[1] & 8
73\n
74If this is zero, y_3 is a parameter. Otherwise it is a variable.
75\n
76\n
77\a arg[2 + j ] for j = 0, 1, 2, 3
78\n
79is the index corresponding to y_j.
80
81\param num_par
82is the total number of values in the vector \a parameter.
83
84\param parameter
85For j = 0, 1, 2, 3,
86if y_j is a parameter, \a parameter [ arg[2 + j] ] is its value.
87
88\param cap_order
89number of columns in the matrix containing the Taylor coefficients.
90
91\par Checked Assertions
92\li NumArg(CExpOp) == 6
93\li NumRes(CExpOp) == 1
94\li arg[0] < static_cast<size_t> ( CompareNe )
95\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
96\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
97\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
98<!-- end conditional_exp_op -->
99*/
100template <class Base>
101inline void conditional_exp_op(
102        size_t         i_z         ,
103        const addr_t*  arg         , 
104        size_t         num_par     ,
105        const Base*    parameter   ,
106        size_t         cap_order   )
107{       // This routine is only for documentation, it should never be used
108        CPPAD_ASSERT_UNKNOWN( false );
109}
110
111/*!
112Shared documentation for conditional expression sparse operations (not called).
113
114<!-- define sparse_conditional_exp_op -->
115The C++ source code coresponding to this operation is
116\verbatim
117        z = CondExpRel(y_0, y_1, y_2, y_3)
118\endverbatim
119where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
120
121\tparam Vector_set
122is the type used for vectors of sets. It can be either
123\c sparse_pack, \c sparse_set, or \c sparse_list.
124
125\param i_z
126is the AD variable index corresponding to the variable z.
127
128\param arg
129\n
130\a arg[0]
131is static cast to size_t from the enum type
132\verbatim
133        enum CompareOp {
134                CompareLt,
135                CompareLe,
136                CompareEq,
137                CompareGe,
138                CompareGt,
139                CompareNe
140        }
141\endverbatim
142for this operation.
143Note that arg[0] cannot be equal to CompareNe.
144\n
145\n
146\a arg[1] & 1
147\n
148If this is zero, y_0 is a parameter. Otherwise it is a variable.
149\n
150\n
151\a arg[1] & 2
152\n
153If this is zero, y_1 is a parameter. Otherwise it is a variable.
154\n
155\n
156\a arg[1] & 4
157\n
158If this is zero, y_2 is a parameter. Otherwise it is a variable.
159\n
160\n
161\a arg[1] & 8
162\n
163If this is zero, y_3 is a parameter. Otherwise it is a variable.
164\n
165\n
166\a arg[2 + j ] for j = 0, 1, 2, 3
167\n
168is the index corresponding to y_j.
169
170\param num_par
171is the total number of values in the vector \a parameter.
172
173\par Checked Assertions
174\li NumArg(CExpOp) == 6
175\li NumRes(CExpOp) == 1
176\li arg[0] < static_cast<size_t> ( CompareNe )
177\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
178\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
179\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
180<!-- end sparse_conditional_exp_op -->
181*/
182template <class Vector_set>
183inline void sparse_conditional_exp_op(
184        size_t         i_z           ,
185        const addr_t*  arg           , 
186        size_t         num_par       )
187{       // This routine is only for documentation, it should never be used
188        CPPAD_ASSERT_UNKNOWN( false );
189}
190
191/*!
192Compute forward mode Taylor coefficients for op = CExpOp.
193
194<!-- replace conditional_exp_op -->
195The C++ source code coresponding to this operation is
196\verbatim
197        z = CondExpRel(y_0, y_1, y_2, y_3)
198\endverbatim
199where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
200
201\tparam Base
202base type for the operator; i.e., this operation was recorded
203using AD< \a Base > and computations by this routine are done using type
204\a Base.
205
206\param i_z
207is the AD variable index corresponding to the variable z.
208
209\param arg
210\n
211\a arg[0]
212is static cast to size_t from the enum type
213\verbatim
214        enum CompareOp {
215                CompareLt,
216                CompareLe,
217                CompareEq,
218                CompareGe,
219                CompareGt,
220                CompareNe
221        }
222\endverbatim
223for this operation.
224Note that arg[0] cannot be equal to CompareNe.
225\n
226\n
227\a arg[1] & 1
228\n
229If this is zero, y_0 is a parameter. Otherwise it is a variable.
230\n
231\n
232\a arg[1] & 2
233\n
234If this is zero, y_1 is a parameter. Otherwise it is a variable.
235\n
236\n
237\a arg[1] & 4
238\n
239If this is zero, y_2 is a parameter. Otherwise it is a variable.
240\n
241\n
242\a arg[1] & 8
243\n
244If this is zero, y_3 is a parameter. Otherwise it is a variable.
245\n
246\n
247\a arg[2 + j ] for j = 0, 1, 2, 3
248\n
249is the index corresponding to y_j.
250
251\param num_par
252is the total number of values in the vector \a parameter.
253
254\param parameter
255For j = 0, 1, 2, 3,
256if y_j is a parameter, \a parameter [ arg[2 + j] ] is its value.
257
258\param cap_order
259number of columns in the matrix containing the Taylor coefficients.
260
261\par Checked Assertions
262\li NumArg(CExpOp) == 6
263\li NumRes(CExpOp) == 1
264\li arg[0] < static_cast<size_t> ( CompareNe )
265\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
266\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
267\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
268<!-- end conditional_exp_op -->
269
270\param p
271is the lowest order of the Taylor coefficient of z that we are computing.
272
273\param q
274is the highest order of the Taylor coefficient of z that we are computing.
275
276\param taylor
277\b Input:
278For j = 0, 1, 2, 3 and k = 0 , ... , q,
279if y_j is a variable then
280<code>taylor [ arg[2+j] * cap_order + k ]</code>
281is the k-th order Taylor coefficient corresponding to y_j.
282\n
283\b Input: <code>taylor [ i_z * cap_order + k ]</code>
284for k = 0 , ... , p-1,
285is the k-th order Taylor coefficient corresponding to z.
286\n
287\b Output: <code>taylor [ i_z * cap_order + k ]</code>
288for k = p , ... , q,
289is the k-th order Taylor coefficient corresponding to z.
290
291*/
292template <class Base>
293inline void forward_cond_op(
294        size_t         p           ,
295        size_t         q           ,
296        size_t         i_z         ,
297        const addr_t*  arg         , 
298        size_t         num_par     ,
299        const Base*    parameter   ,
300        size_t         cap_order   ,
301        Base*          taylor      )
302{       Base y_0, y_1, y_2, y_3;
303        Base zero(0);
304        Base* z = taylor + i_z * cap_order;
305
306        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < static_cast<size_t> (CompareNe) );
307        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
308        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
309        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
310
311        if( arg[1] & 1 )
312        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
313                y_0 = taylor[ arg[2] * cap_order + 0 ];
314        }
315        else
316        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
317                y_0 = parameter[ arg[2] ];
318        }
319        if( arg[1] & 2 )
320        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < i_z );
321                y_1 = taylor[ arg[3] * cap_order + 0 ];
322        }
323        else
324        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < num_par );
325                y_1 = parameter[ arg[3] ];
326        }
327        if( p == 0 )
328        {       if( arg[1] & 4 )
329                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
330                        y_2 = taylor[ arg[4] * cap_order + 0 ];
331                }
332                else
333                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < num_par );
334                        y_2 = parameter[ arg[4] ];
335                }
336                if( arg[1] & 8 )
337                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
338                        y_3 = taylor[ arg[5] * cap_order + 0 ];
339                }
340                else
341                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < num_par );
342                        y_3 = parameter[ arg[5] ];
343                }
344                z[0] = CondExpOp(
345                        CompareOp( arg[0] ),
346                        y_0,
347                        y_1,
348                        y_2,
349                        y_3
350                );
351                p++;
352        }
353        for(size_t d = p; d <= q; d++)
354        {       if( arg[1] & 4 )
355                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
356                        y_2 = taylor[ arg[4] * cap_order + d];
357                }
358                else    y_2 = zero;
359                if( arg[1] & 8 )
360                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
361                        y_3 = taylor[ arg[5] * cap_order + d];
362                }
363                else    y_3 = zero;
364                z[d] = CondExpOp(
365                        CompareOp( arg[0] ),
366                        y_0,
367                        y_1,
368                        y_2,
369                        y_3
370                );
371        }
372        return;
373}
374
375/*!
376Multiple directions forward mode Taylor coefficients for op = CExpOp.
377
378<!-- replace conditional_exp_op -->
379The C++ source code coresponding to this operation is
380\verbatim
381        z = CondExpRel(y_0, y_1, y_2, y_3)
382\endverbatim
383where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
384
385\tparam Base
386base type for the operator; i.e., this operation was recorded
387using AD< \a Base > and computations by this routine are done using type
388\a Base.
389
390\param i_z
391is the AD variable index corresponding to the variable z.
392
393\param arg
394\n
395\a arg[0]
396is static cast to size_t from the enum type
397\verbatim
398        enum CompareOp {
399                CompareLt,
400                CompareLe,
401                CompareEq,
402                CompareGe,
403                CompareGt,
404                CompareNe
405        }
406\endverbatim
407for this operation.
408Note that arg[0] cannot be equal to CompareNe.
409\n
410\n
411\a arg[1] & 1
412\n
413If this is zero, y_0 is a parameter. Otherwise it is a variable.
414\n
415\n
416\a arg[1] & 2
417\n
418If this is zero, y_1 is a parameter. Otherwise it is a variable.
419\n
420\n
421\a arg[1] & 4
422\n
423If this is zero, y_2 is a parameter. Otherwise it is a variable.
424\n
425\n
426\a arg[1] & 8
427\n
428If this is zero, y_3 is a parameter. Otherwise it is a variable.
429\n
430\n
431\a arg[2 + j ] for j = 0, 1, 2, 3
432\n
433is the index corresponding to y_j.
434
435\param num_par
436is the total number of values in the vector \a parameter.
437
438\param parameter
439For j = 0, 1, 2, 3,
440if y_j is a parameter, \a parameter [ arg[2 + j] ] is its value.
441
442\param cap_order
443number of columns in the matrix containing the Taylor coefficients.
444
445\par Checked Assertions
446\li NumArg(CExpOp) == 6
447\li NumRes(CExpOp) == 1
448\li arg[0] < static_cast<size_t> ( CompareNe )
449\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
450\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
451\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
452<!-- end conditional_exp_op -->
453
454\param q
455is order of the Taylor coefficient of z that we are computing.
456
457\param r
458is the number of Taylor coefficient directions that we are computing.
459
460\par tpv
461We use the notation
462<code>tpv = (cap_order-1) * r + 1</code>
463which is the number of Taylor coefficients per variable
464
465\param taylor
466\b Input:
467For j = 0, 1, 2, 3, k = 1, ..., q,
468if y_j is a variable then
469<code>taylor [ arg[2+j] * tpv + 0 ]</code>
470is the zero order Taylor coefficient corresponding to y_j and
471<code>taylor [ arg[2+j] * tpv + (k-1)*r+1+ell</code> is its
472k-th order Taylor coefficient in the ell-th direction.
473\n
474\b Input:
475For j = 0, 1, 2, 3, k = 1, ..., q-1,
476<code>taylor [ i_z * tpv + 0 ]</code> 
477is the zero order Taylor coefficient corresponding to z and
478<code>taylor [ i_z * tpv + (k-1)*r+1+ell</code> is its
479k-th order Taylor coefficient in the ell-th direction.
480\n
481\b Output: <code>taylor [ i_z * tpv + (q-1)*r+1+ell ]</code>
482is the q-th order Taylor coefficient corresponding to z
483in the ell-th direction.
484*/
485template <class Base>
486inline void forward_cond_op_dir(
487        size_t         q           ,
488        size_t         r           ,
489        size_t         i_z         ,
490        const addr_t*  arg         , 
491        size_t         num_par     ,
492        const Base*    parameter   ,
493        size_t         cap_order   ,
494        Base*          taylor      )
495{       Base y_0, y_1, y_2, y_3;
496        Base zero(0);
497        size_t num_taylor_per_var = (cap_order-1) * r + 1;
498        Base* z = taylor + i_z * num_taylor_per_var;
499
500        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < static_cast<size_t> (CompareNe) );
501        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
502        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
503        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
504        CPPAD_ASSERT_UNKNOWN( 0 < q );
505        CPPAD_ASSERT_UNKNOWN( q < cap_order );
506
507        if( arg[1] & 1 )
508        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
509                y_0 = taylor[ arg[2] * num_taylor_per_var + 0 ];
510        }
511        else
512        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
513                y_0 = parameter[ arg[2] ];
514        }
515        if( arg[1] & 2 )
516        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < i_z );
517                y_1 = taylor[ arg[3] * num_taylor_per_var + 0 ];
518        }
519        else
520        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < num_par );
521                y_1 = parameter[ arg[3] ];
522        }
523        size_t m = (q-1) * r + 1;
524        for(size_t ell = 0; ell < r; ell++)
525        {       if( arg[1] & 4 )
526                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
527                        y_2 = taylor[ arg[4] * num_taylor_per_var + m + ell];
528                }
529                else    y_2 = zero;
530                if( arg[1] & 8 )
531                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
532                        y_3 = taylor[ arg[5] * num_taylor_per_var + m + ell];
533                }
534                else    y_3 = zero;
535                z[m+ell] = CondExpOp(
536                        CompareOp( arg[0] ),
537                        y_0,
538                        y_1,
539                        y_2,
540                        y_3
541                );
542        }
543        return;
544}
545
546/*!
547Compute zero order forward mode Taylor coefficients for op = CExpOp.
548
549<!-- replace conditional_exp_op -->
550The C++ source code coresponding to this operation is
551\verbatim
552        z = CondExpRel(y_0, y_1, y_2, y_3)
553\endverbatim
554where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
555
556\tparam Base
557base type for the operator; i.e., this operation was recorded
558using AD< \a Base > and computations by this routine are done using type
559\a Base.
560
561\param i_z
562is the AD variable index corresponding to the variable z.
563
564\param arg
565\n
566\a arg[0]
567is static cast to size_t from the enum type
568\verbatim
569        enum CompareOp {
570                CompareLt,
571                CompareLe,
572                CompareEq,
573                CompareGe,
574                CompareGt,
575                CompareNe
576        }
577\endverbatim
578for this operation.
579Note that arg[0] cannot be equal to CompareNe.
580\n
581\n
582\a arg[1] & 1
583\n
584If this is zero, y_0 is a parameter. Otherwise it is a variable.
585\n
586\n
587\a arg[1] & 2
588\n
589If this is zero, y_1 is a parameter. Otherwise it is a variable.
590\n
591\n
592\a arg[1] & 4
593\n
594If this is zero, y_2 is a parameter. Otherwise it is a variable.
595\n
596\n
597\a arg[1] & 8
598\n
599If this is zero, y_3 is a parameter. Otherwise it is a variable.
600\n
601\n
602\a arg[2 + j ] for j = 0, 1, 2, 3
603\n
604is the index corresponding to y_j.
605
606\param num_par
607is the total number of values in the vector \a parameter.
608
609\param parameter
610For j = 0, 1, 2, 3,
611if y_j is a parameter, \a parameter [ arg[2 + j] ] is its value.
612
613\param cap_order
614number of columns in the matrix containing the Taylor coefficients.
615
616\par Checked Assertions
617\li NumArg(CExpOp) == 6
618\li NumRes(CExpOp) == 1
619\li arg[0] < static_cast<size_t> ( CompareNe )
620\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
621\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
622\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
623<!-- end conditional_exp_op -->
624
625\param taylor
626\b Input:
627For j = 0, 1, 2, 3,
628if y_j is a variable then
629\a taylor [ \a arg[2+j] * cap_order + 0 ]
630is the zero order Taylor coefficient corresponding to y_j.
631\n
632\b Output: \a taylor [ \a i_z * \a cap_order + 0 ]
633is the zero order Taylor coefficient corresponding to z.
634*/
635template <class Base>
636inline void forward_cond_op_0(
637        size_t         i_z         ,
638        const addr_t*  arg         , 
639        size_t         num_par     ,
640        const Base*    parameter   ,
641        size_t         cap_order   ,
642        Base*          taylor      )
643{       Base y_0, y_1, y_2, y_3;
644        Base* z;
645
646        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < static_cast<size_t> (CompareNe) );
647        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
648        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
649        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
650
651        if( arg[1] & 1 )
652        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
653                y_0 = taylor[ arg[2] * cap_order + 0 ];
654        }
655        else
656        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
657                y_0 = parameter[ arg[2] ];
658        }
659        if( arg[1] & 2 )
660        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < i_z );
661                y_1 = taylor[ arg[3] * cap_order + 0 ];
662        }
663        else
664        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < num_par );
665                y_1 = parameter[ arg[3] ];
666        }
667        if( arg[1] & 4 )
668        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
669                y_2 = taylor[ arg[4] * cap_order + 0 ];
670        }
671        else
672        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < num_par );
673                y_2 = parameter[ arg[4] ];
674        }
675        if( arg[1] & 8 )
676        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
677                y_3 = taylor[ arg[5] * cap_order + 0 ];
678        }
679        else
680        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < num_par );
681                y_3 = parameter[ arg[5] ];
682        }
683        z = taylor + i_z * cap_order;
684        z[0] = CondExpOp(
685                CompareOp( arg[0] ),
686                y_0,
687                y_1,
688                y_2,
689                y_3
690        );
691        return;
692}
693
694/*!
695Compute reverse mode Taylor coefficients for op = CExpOp.
696
697This routine is given the partial derivatives of a function
698G( z , y , x , w , ... )
699and it uses them to compute the partial derivatives of
700\verbatim
701        H( y , x , w , u , ... ) = G[ z(y) , y , x , w , u , ... ]
702\endverbatim
703where y above represents y_0, y_1, y_2, y_3.
704
705<!-- replace conditional_exp_op -->
706The C++ source code coresponding to this operation is
707\verbatim
708        z = CondExpRel(y_0, y_1, y_2, y_3)
709\endverbatim
710where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
711
712\tparam Base
713base type for the operator; i.e., this operation was recorded
714using AD< \a Base > and computations by this routine are done using type
715\a Base.
716
717\param i_z
718is the AD variable index corresponding to the variable z.
719
720\param arg
721\n
722\a arg[0]
723is static cast to size_t from the enum type
724\verbatim
725        enum CompareOp {
726                CompareLt,
727                CompareLe,
728                CompareEq,
729                CompareGe,
730                CompareGt,
731                CompareNe
732        }
733\endverbatim
734for this operation.
735Note that arg[0] cannot be equal to CompareNe.
736\n
737\n
738\a arg[1] & 1
739\n
740If this is zero, y_0 is a parameter. Otherwise it is a variable.
741\n
742\n
743\a arg[1] & 2
744\n
745If this is zero, y_1 is a parameter. Otherwise it is a variable.
746\n
747\n
748\a arg[1] & 4
749\n
750If this is zero, y_2 is a parameter. Otherwise it is a variable.
751\n
752\n
753\a arg[1] & 8
754\n
755If this is zero, y_3 is a parameter. Otherwise it is a variable.
756\n
757\n
758\a arg[2 + j ] for j = 0, 1, 2, 3
759\n
760is the index corresponding to y_j.
761
762\param num_par
763is the total number of values in the vector \a parameter.
764
765\param parameter
766For j = 0, 1, 2, 3,
767if y_j is a parameter, \a parameter [ arg[2 + j] ] is its value.
768
769\param cap_order
770number of columns in the matrix containing the Taylor coefficients.
771
772\par Checked Assertions
773\li NumArg(CExpOp) == 6
774\li NumRes(CExpOp) == 1
775\li arg[0] < static_cast<size_t> ( CompareNe )
776\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
777\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
778\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
779<!-- end conditional_exp_op -->
780
781\param d
782is the order of the Taylor coefficient of z that we are  computing.
783
784\param taylor
785\b Input:
786For j = 0, 1, 2, 3 and k = 0 , ... , \a d,
787if y_j is a variable then
788\a taylor [ \a arg[2+j] * cap_order + k ]
789is the k-th order Taylor coefficient corresponding to y_j.
790\n
791\a taylor [ \a i_z * \a cap_order + k ]
792for k = 0 , ... , \a d
793is the k-th order Taylor coefficient corresponding to z.
794
795\param nc_partial
796number of columns in the matrix containing the Taylor coefficients.
797
798\param partial
799\b Input:
800For j = 0, 1, 2, 3 and k = 0 , ... , \a d,
801if y_j is a variable then
802\a partial [ \a arg[2+j] * nc_partial + k ]
803is the partial derivative of G( z , y , x , w , u , ... )
804with respect to the k-th order Taylor coefficient corresponding to y_j.
805\n
806\b Input: \a partial [ \a i_z * \a cap_order + k ]
807for k = 0 , ... , \a d
808is the partial derivative of G( z , y , x , w , u , ... )
809with respect to the k-th order Taylor coefficient corresponding to z.
810\n
811\b Output:
812For j = 0, 1, 2, 3 and k = 0 , ... , \a d,
813if y_j is a variable then
814\a partial [ \a arg[2+j] * nc_partial + k ]
815is the partial derivative of H( y , x , w , u , ... )
816with respect to the k-th order Taylor coefficient corresponding to y_j.
817
818*/
819template <class Base>
820inline void reverse_cond_op(
821        size_t         d           ,
822        size_t         i_z         ,
823        const addr_t*  arg         , 
824        size_t         num_par     ,
825        const Base*    parameter   ,
826        size_t         cap_order   ,
827        const Base*    taylor      ,
828        size_t         nc_partial  ,
829        Base*          partial     )
830{       Base y_0, y_1;
831        Base zero(0);
832        Base* pz;
833        Base* py_2;
834        Base* py_3;
835
836        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < static_cast<size_t> (CompareNe) );
837        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
838        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
839        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
840
841        pz = partial + i_z * nc_partial + 0;
842        if( arg[1] & 1 )
843        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
844                y_0 = taylor[ arg[2] * cap_order + 0 ];
845        }
846        else
847        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
848                y_0 = parameter[ arg[2] ];
849        }
850        if( arg[1] & 2 )
851        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < i_z );
852                y_1 = taylor[ arg[3] * cap_order + 0 ];
853        }
854        else
855        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < num_par );
856                y_1 = parameter[ arg[3] ];
857        }
858        if( arg[1] & 4 )
859        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
860                py_2 = partial + arg[4] * nc_partial;
861                size_t j = d + 1;
862                while(j--)
863                {       py_2[j] += CondExpOp(
864                                CompareOp( arg[0] ),
865                                y_0,
866                                y_1,
867                                pz[j],
868                                zero
869                        );
870                }
871        }
872        if( arg[1] & 8 )
873        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
874                py_3 = partial + arg[5] * nc_partial;
875                size_t j = d + 1;
876                while(j--)
877                {       py_3[j] += CondExpOp(
878                                CompareOp( arg[0] ),
879                                y_0,
880                                y_1,
881                                zero,
882                                pz[j]
883                        );
884                }
885        }
886        return;
887}
888
889/*!
890Compute forward Jacobian sparsity patterns for op = CExpOp.
891
892<!-- replace sparse_conditional_exp_op -->
893The C++ source code coresponding to this operation is
894\verbatim
895        z = CondExpRel(y_0, y_1, y_2, y_3)
896\endverbatim
897where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
898
899\tparam Vector_set
900is the type used for vectors of sets. It can be either
901\c sparse_pack, \c sparse_set, or \c sparse_list.
902
903\param i_z
904is the AD variable index corresponding to the variable z.
905
906\param arg
907\n
908\a arg[0]
909is static cast to size_t from the enum type
910\verbatim
911        enum CompareOp {
912                CompareLt,
913                CompareLe,
914                CompareEq,
915                CompareGe,
916                CompareGt,
917                CompareNe
918        }
919\endverbatim
920for this operation.
921Note that arg[0] cannot be equal to CompareNe.
922\n
923\n
924\a arg[1] & 1
925\n
926If this is zero, y_0 is a parameter. Otherwise it is a variable.
927\n
928\n
929\a arg[1] & 2
930\n
931If this is zero, y_1 is a parameter. Otherwise it is a variable.
932\n
933\n
934\a arg[1] & 4
935\n
936If this is zero, y_2 is a parameter. Otherwise it is a variable.
937\n
938\n
939\a arg[1] & 8
940\n
941If this is zero, y_3 is a parameter. Otherwise it is a variable.
942\n
943\n
944\a arg[2 + j ] for j = 0, 1, 2, 3
945\n
946is the index corresponding to y_j.
947
948\param num_par
949is the total number of values in the vector \a parameter.
950
951\par Checked Assertions
952\li NumArg(CExpOp) == 6
953\li NumRes(CExpOp) == 1
954\li arg[0] < static_cast<size_t> ( CompareNe )
955\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
956\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
957\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
958<!-- end sparse_conditional_exp_op -->
959
960\param sparsity
961\b Input:
962if y_2 is a variable, the set with index t is
963the sparsity pattern corresponding to y_2.
964This identifies which of the independent variables the variable y_2
965depends on.
966\n
967\b Input:
968if y_3 is a variable, the set with index t is
969the sparsity pattern corresponding to y_3.
970This identifies which of the independent variables the variable y_3
971depends on.
972\n
973\b Output:
974The set with index T is
975the sparsity pattern corresponding to z.
976This identifies which of the independent variables the variable z
977depends on.
978*/
979template <class Vector_set>
980inline void forward_sparse_jacobian_cond_op(
981        size_t             i_z           ,
982        const addr_t*      arg           , 
983        size_t             num_par       ,
984        Vector_set&        sparsity      )
985{
986        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < static_cast<size_t> (CompareNe) );
987        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
988        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
989        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
990
991# ifndef NDEBUG
992        if( arg[1] & 1 )
993        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
994        }
995        else
996        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
997        }
998        if( arg[1] & 2 )
999        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < i_z );
1000        }
1001        else
1002        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < num_par );
1003        }
1004# endif
1005        if( arg[1] & 4 )
1006        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
1007                if( arg[1] & 8 )
1008                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
1009                        sparsity.binary_union(i_z, arg[4], arg[5], sparsity);
1010                }
1011                else
1012                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < num_par );
1013                        sparsity.assignment(i_z, arg[4], sparsity);
1014                }
1015        }       
1016        else
1017        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < num_par );
1018                if( arg[1] & 8 )
1019                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
1020                        sparsity.assignment(i_z, arg[5], sparsity);
1021                }
1022                else
1023                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < num_par );
1024                        sparsity.clear(i_z);
1025                }
1026        }
1027        return;
1028}
1029
1030/*!
1031Compute reverse Jacobian sparsity patterns for op = CExpOp.
1032
1033This routine is given the sparsity patterns
1034for a function G(z, y, x, ... )
1035and it uses them to compute the sparsity patterns for
1036\verbatim
1037        H( y, x, w , u , ... ) = G[ z(x,y) , y , x , w , u , ... ]
1038\endverbatim
1039where y represents the combination of y_0, y_1, y_2, and y_3.
1040
1041<!-- replace sparse_conditional_exp_op -->
1042The C++ source code coresponding to this operation is
1043\verbatim
1044        z = CondExpRel(y_0, y_1, y_2, y_3)
1045\endverbatim
1046where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
1047
1048\tparam Vector_set
1049is the type used for vectors of sets. It can be either
1050\c sparse_pack, \c sparse_set, or \c sparse_list.
1051
1052\param i_z
1053is the AD variable index corresponding to the variable z.
1054
1055\param arg
1056\n
1057\a arg[0]
1058is static cast to size_t from the enum type
1059\verbatim
1060        enum CompareOp {
1061                CompareLt,
1062                CompareLe,
1063                CompareEq,
1064                CompareGe,
1065                CompareGt,
1066                CompareNe
1067        }
1068\endverbatim
1069for this operation.
1070Note that arg[0] cannot be equal to CompareNe.
1071\n
1072\n
1073\a arg[1] & 1
1074\n
1075If this is zero, y_0 is a parameter. Otherwise it is a variable.
1076\n
1077\n
1078\a arg[1] & 2
1079\n
1080If this is zero, y_1 is a parameter. Otherwise it is a variable.
1081\n
1082\n
1083\a arg[1] & 4
1084\n
1085If this is zero, y_2 is a parameter. Otherwise it is a variable.
1086\n
1087\n
1088\a arg[1] & 8
1089\n
1090If this is zero, y_3 is a parameter. Otherwise it is a variable.
1091\n
1092\n
1093\a arg[2 + j ] for j = 0, 1, 2, 3
1094\n
1095is the index corresponding to y_j.
1096
1097\param num_par
1098is the total number of values in the vector \a parameter.
1099
1100\par Checked Assertions
1101\li NumArg(CExpOp) == 6
1102\li NumRes(CExpOp) == 1
1103\li arg[0] < static_cast<size_t> ( CompareNe )
1104\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
1105\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
1106\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
1107<!-- end sparse_conditional_exp_op -->
1108
1109\param nz_compare
1110Are the derivatives with respect to left and right of the expression below
1111considered to be non-zero:
1112\code
1113        CondExpRel(left, right, if_true, if_false)
1114\endcode
1115This is used by the optimizer to obtain the correct dependency relations.
1116
1117
1118\param sparsity
1119if y_2 is a variable, the set with index t is
1120the sparsity pattern corresponding to y_2.
1121This identifies which of the dependent variables depend on the variable y_2.
1122On input, this pattern corresponds to the function G.
1123On ouput, it corresponds to the function H.
1124\n
1125\n
1126if y_3 is a variable, the set with index t is
1127the sparsity pattern corresponding to y_3.
1128This identifies which of the dependent variables depeond on the variable y_3.
1129On input, this pattern corresponds to the function G.
1130On ouput, it corresponds to the function H.
1131\n
1132\b Output:
1133The set with index T is
1134the sparsity pattern corresponding to z.
1135This identifies which of the dependent variables depend on the variable z.
1136On input and output, this pattern corresponds to the function G.
1137*/
1138template <class Vector_set>
1139inline void reverse_sparse_jacobian_cond_op(
1140        bool                nz_compare    ,
1141        size_t              i_z           ,
1142        const addr_t*       arg           , 
1143        size_t              num_par       ,
1144        Vector_set&         sparsity      )
1145{       
1146        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < static_cast<size_t> (CompareNe) );
1147        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
1148        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
1149        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
1150
1151# ifndef NDEBUG
1152        if( arg[1] & 1 )
1153        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
1154        }
1155        else
1156        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
1157        }
1158        if( arg[1] & 2 )
1159        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < i_z );
1160        }
1161        else
1162        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < num_par );
1163        }
1164        if( ! ( arg[1] & 4 ) )
1165        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < num_par );
1166        }
1167        if( ! ( arg[1] & 8 ) )
1168        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < num_par );
1169        }
1170# endif
1171        if( nz_compare )
1172        {       if( arg[1] & 1 )
1173                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
1174                        sparsity.binary_union(arg[2], arg[2], i_z, sparsity);
1175                }
1176                if( arg[1] & 2 )
1177                {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < i_z );
1178                        sparsity.binary_union(arg[3], arg[3], i_z, sparsity);
1179                }
1180        }
1181        // --------------------------------------------------------------------
1182        if( arg[1] & 4 )
1183        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
1184                sparsity.binary_union(arg[4], arg[4], i_z, sparsity);
1185        }
1186        if( arg[1] & 8 )
1187        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
1188                sparsity.binary_union(arg[5], arg[5], i_z, sparsity);
1189        }
1190        return;
1191}
1192
1193/*!
1194Compute reverse Hessian sparsity patterns for op = CExpOp.
1195
1196This routine is given the sparsity patterns
1197for a function G(z, y, x, ... )
1198and it uses them to compute the sparsity patterns for
1199\verbatim
1200        H( y, x, w , u , ... ) = G[ z(x,y) , y , x , w , u , ... ]
1201\endverbatim
1202where y represents the combination of y_0, y_1, y_2, and y_3.
1203
1204<!-- replace sparse_conditional_exp_op -->
1205The C++ source code coresponding to this operation is
1206\verbatim
1207        z = CondExpRel(y_0, y_1, y_2, y_3)
1208\endverbatim
1209where Rel is one of the following: Lt, Le, Eq, Ge, Gt.
1210
1211\tparam Vector_set
1212is the type used for vectors of sets. It can be either
1213\c sparse_pack, \c sparse_set, or \c sparse_list.
1214
1215\param i_z
1216is the AD variable index corresponding to the variable z.
1217
1218\param arg
1219\n
1220\a arg[0]
1221is static cast to size_t from the enum type
1222\verbatim
1223        enum CompareOp {
1224                CompareLt,
1225                CompareLe,
1226                CompareEq,
1227                CompareGe,
1228                CompareGt,
1229                CompareNe
1230        }
1231\endverbatim
1232for this operation.
1233Note that arg[0] cannot be equal to CompareNe.
1234\n
1235\n
1236\a arg[1] & 1
1237\n
1238If this is zero, y_0 is a parameter. Otherwise it is a variable.
1239\n
1240\n
1241\a arg[1] & 2
1242\n
1243If this is zero, y_1 is a parameter. Otherwise it is a variable.
1244\n
1245\n
1246\a arg[1] & 4
1247\n
1248If this is zero, y_2 is a parameter. Otherwise it is a variable.
1249\n
1250\n
1251\a arg[1] & 8
1252\n
1253If this is zero, y_3 is a parameter. Otherwise it is a variable.
1254\n
1255\n
1256\a arg[2 + j ] for j = 0, 1, 2, 3
1257\n
1258is the index corresponding to y_j.
1259
1260\param num_par
1261is the total number of values in the vector \a parameter.
1262
1263\par Checked Assertions
1264\li NumArg(CExpOp) == 6
1265\li NumRes(CExpOp) == 1
1266\li arg[0] < static_cast<size_t> ( CompareNe )
1267\li arg[1] != 0; i.e., not all of y_0, y_1, y_2, y_3 are parameters.
1268\li For j = 0, 1, 2, 3 if y_j is a parameter, arg[2+j] < num_par.
1269\li For j = 0, 1, 2, 3 if y_j is a variable, arg[2+j] < iz.
1270<!-- end sparse_conditional_exp_op -->
1271
1272
1273\param jac_reverse
1274\a jac_reverse[i_z]
1275is false (true) if the Jacobian of G with respect to z is always zero
1276(may be non-zero).
1277\n
1278\n
1279\a jac_reverse[ arg[4] ]
1280If y_2 is a variable,
1281\a jac_reverse[ arg[4] ]
1282is false (true) if the Jacobian with respect to y_2 is always zero
1283(may be non-zero).
1284On input, it corresponds to the function G,
1285and on output it corresponds to the function H.
1286\n
1287\n
1288\a jac_reverse[ arg[5] ]
1289If y_3 is a variable,
1290\a jac_reverse[ arg[5] ]
1291is false (true) if the Jacobian with respect to y_3 is always zero
1292(may be non-zero).
1293On input, it corresponds to the function G,
1294and on output it corresponds to the function H.
1295
1296\param hes_sparsity
1297The set with index \a i_z in \a hes_sparsity
1298is the Hessian sparsity pattern for the function G
1299where one of the partials is with respect to z.
1300\n
1301\n
1302If y_2 is a variable,
1303the set with index \a arg[4] in \a hes_sparsity
1304is the Hessian sparsity pattern
1305where one of the partials is with respect to y_2.
1306On input, this pattern corresponds to the function G.
1307On output, this pattern corresponds to the function H.
1308\n
1309\n
1310If y_3 is a variable,
1311the set with index \a arg[5] in \a hes_sparsity
1312is the Hessian sparsity pattern
1313where one of the partials is with respect to y_3.
1314On input, this pattern corresponds to the function G.
1315On output, this pattern corresponds to the function H.
1316*/
1317template <class Vector_set>
1318inline void reverse_sparse_hessian_cond_op(
1319        size_t               i_z           ,
1320        const addr_t*        arg           , 
1321        size_t               num_par       ,
1322        bool*                jac_reverse   ,
1323        Vector_set&          hes_sparsity  )
1324{       
1325
1326        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < static_cast<size_t> (CompareNe) );
1327        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
1328        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
1329        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
1330
1331# ifndef NDEBUG
1332        if( arg[1] & 1 )
1333        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < i_z );
1334        }
1335        else
1336        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par );
1337        }
1338        if( arg[1] & 2 )
1339        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < i_z );
1340        }
1341        else
1342        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[3]) < num_par );
1343        }
1344        if( ! ( arg[1] & 4 ) )
1345        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < num_par );
1346        }
1347        if( ! ( arg[1] & 8 ) )
1348        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < num_par );
1349        }
1350# endif
1351        if( arg[1] & 4 )
1352        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[4]) < i_z );
1353
1354                hes_sparsity.binary_union(arg[4], arg[4], i_z, hes_sparsity);
1355                jac_reverse[ arg[4] ] |= jac_reverse[i_z];
1356        }
1357        if( arg[1] & 8 )
1358        {       CPPAD_ASSERT_UNKNOWN( size_t(arg[5]) < i_z );
1359
1360                hes_sparsity.binary_union(arg[5], arg[5], i_z, hes_sparsity);
1361                jac_reverse[ arg[5] ] |= jac_reverse[i_z];
1362        }
1363        return;
1364}
1365
1366} // END_CPPAD_NAMESPACE
1367# endif
Note: See TracBrowser for help on using the repository browser.