source: trunk/cppad/local/rev_jac_sweep.hpp @ 2794

Last change on this file since 2794 was 2794, checked in by bradbell, 7 years ago
  1. Use CPPAD_NULL, intead of 0, for null pointer.

check_if_0.sh: Ignore subdirectories of new directories.
jenkins.sh: output logs when an error occurs.
acos_op.hpp: avoid use of q (will use it for an order index).
asin_op.hpp: avoid use of q (will use it for an order index).
forward_sweep.hpp: chnage d to p, use const in prototype.
div_op.hpp: minor edit of white space.
atom_usead_2.cpp: use ADFUN to compute variable/parameter information.

  • Property svn:keywords set to Id
File size: 16.9 KB
Line 
1/* $Id: rev_jac_sweep.hpp 2794 2013-05-02 08:20:30Z bradbell $ */
2# ifndef CPPAD_REV_JAC_SWEEP_INCLUDED
3# define CPPAD_REV_JAC_SWEEP_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16CPPAD_BEGIN_NAMESPACE
17/*!
18\defgroup rev_jac_sweep_hpp rev_jac_sweep.hpp
19\{
20\file rev_jac_sweep.hpp
21Compute Reverse mode Jacobian sparsity patterns.
22*/
23
24/*!
25\def CPPAD_REV_JAC_SWEEP_TRACE
26This value is either zero or one.
27Zero is the normal operational value.
28If it is one, a trace of every rev_jac_sweep computation is printed.
29*/
30# define CPPAD_REV_JAC_SWEEP_TRACE 0
31
32/*!
33Given the sparsity pattern for the dependent variables,
34RevJacSweep computes the sparsity pattern for all the independent variables.
35
36\tparam Base
37base type for the operator; i.e., this operation sequence was recorded
38using AD< \a Base > and computations by this routine are done using type
39\a Base.
40
41\tparam Vector_set
42is the type used for vectors of sets. It can be either
43\c sparse_pack, \c sparse_set, or \c sparse_list.
44
45\param n
46is the number of independent variables on the tape.
47
48\param numvar
49is the total number of variables on the tape; i.e.,
50\a play->num_rec_var().
51This is also the number of rows in the entire sparsity pattern \a RevJac.
52
53\param play
54The information stored in \a play
55is a recording of the operations corresponding to a function
56\f[
57        F : {\bf R}^n \rightarrow {\bf R}^m
58\f]
59where \f$ n \f$ is the number of independent variables
60and \f$ m \f$ is the number of dependent variables.
61The object \a play is effectly constant.
62It is not declared const because while playing back the tape
63the object \a play holds information about the currentl location
64with in the tape and this changes during playback.
65
66\param var_sparsity
67For i = 0 , ... , \a numvar - 1,
68(all the variables on the tape)
69the forward Jacobian sparsity pattern for variable i
70corresponds to the set with index i in \a var_sparsity.
71\b
72\b
73\b Input:
74For i = 0 , ... , \a numvar - 1,
75the forward Jacobian sparsity pattern for variable i is an input
76if i corresponds to a dependent variable.
77Otherwise the sparsity patten is empty.
78\n
79\n
80\b Output: For j = 1 , ... , \a n,
81the sparsity pattern for the dependent variable with index (j-1)
82is given by the set with index index j in \a var_sparsity.
83*/
84
85template <class Base, class Vector_set>
86void RevJacSweep(
87        size_t                n,
88        size_t                numvar,
89        player<Base>         *play,
90        Vector_set&           var_sparsity
91)
92{
93        OpCode           op;
94        size_t         i_op;
95        size_t        i_var;
96
97        const addr_t*   arg = CPPAD_NULL;
98
99        size_t            i, j, k;
100
101        // length of the parameter vector (used by CppAD assert macros)
102        const size_t num_par = play->num_rec_par();
103
104        // check numvar argument
105        CPPAD_ASSERT_UNKNOWN( numvar > 0 );
106        CPPAD_ASSERT_UNKNOWN( play->num_rec_var()   == numvar );
107        CPPAD_ASSERT_UNKNOWN( var_sparsity.n_set() == numvar );
108
109        // upper limit (exclusive) for elements in the set
110        size_t limit = var_sparsity.end();
111
112        // vecad_sparsity contains a sparsity pattern for each VecAD object.
113        // vecad_ind maps a VecAD index (beginning of the VecAD object)
114        // to the index of the corresponding set in vecad_sparsity.
115        size_t num_vecad_ind   = play->num_rec_vecad_ind();
116        size_t num_vecad_vec   = play->num_rec_vecad_vec();
117        Vector_set  vecad_sparsity;
118        vecad_sparsity.resize(num_vecad_vec, limit);
119        pod_vector<size_t> vecad_ind;
120        if( num_vecad_vec > 0 )
121        {       size_t length;
122                vecad_ind.extend(num_vecad_ind);
123                j             = 0;
124                for(i = 0; i < num_vecad_vec; i++)
125                {       // length of this VecAD
126                        length   = play->GetVecInd(j);
127                        // set to proper index for this VecAD
128                        vecad_ind[j] = i; 
129                        for(k = 1; k <= length; k++)
130                                vecad_ind[j+k] = num_vecad_vec; // invalid index
131                        // start of next VecAD
132                        j       += length + 1;
133                }
134                CPPAD_ASSERT_UNKNOWN( j == play->num_rec_vecad_ind() );
135        }
136
137        // work space used by UserOp.
138        typedef std::set<size_t> size_set;
139        const size_t user_q = limit; // maximum element plus one
140        size_set::iterator set_itr;  // iterator for a standard set
141        size_set::iterator set_end;  // end of iterator sequence
142        vector< size_set > user_r;   // sparsity pattern for the argument x
143        vector< size_set > user_s;   // sparisty pattern for the result y
144        size_t user_index = 0;       // indentifier for this user_atomic operation
145        size_t user_id    = 0;       // user identifier for this call to operator
146        size_t user_i     = 0;       // index in result vector
147        size_t user_j     = 0;       // index in argument vector
148        size_t user_m     = 0;       // size of result vector
149        size_t user_n     = 0;       // size of arugment vector
150        // next expected operator in a UserOp sequence
151        enum { user_start, user_arg, user_ret, user_end } user_state = user_end;
152
153        // Initialize
154        play->start_reverse(op, arg, i_op, i_var);
155        CPPAD_ASSERT_UNKNOWN( op == EndOp );
156# if CPPAD_REV_JAC_SWEEP_TRACE
157        std::cout << std::endl;
158        CppAD::vector<bool> z_value(limit);
159# endif
160        while(op != BeginOp )
161        {
162                // next op
163                play->next_reverse(op, arg, i_op, i_var);
164# ifndef NDEBUG
165                if( i_op <= n )
166                {       CPPAD_ASSERT_UNKNOWN((op == InvOp) | (op == BeginOp));
167                }
168                else    CPPAD_ASSERT_UNKNOWN((op != InvOp) & (op != BeginOp));
169# endif
170
171                // rest of information depends on the case
172                switch( op )
173                {
174                        case AbsOp:
175                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
176                        reverse_sparse_jacobian_unary_op(
177                                i_var, arg[0], var_sparsity
178                        );
179                        break;
180                        // -------------------------------------------------
181
182                        case AddvvOp:
183                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
184                        reverse_sparse_jacobian_binary_op(
185                                i_var, arg, var_sparsity
186                        );
187                        break;
188                        // -------------------------------------------------
189
190                        case AddpvOp:
191                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
192                        reverse_sparse_jacobian_unary_op(
193                                i_var, arg[1], var_sparsity
194                        );
195                        break;
196                        // -------------------------------------------------
197
198                        case AcosOp:
199                        // sqrt(1 - x * x), acos(x)
200                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
201                        reverse_sparse_jacobian_unary_op(
202                                i_var, arg[0], var_sparsity
203                        );
204                        break;
205                        // -------------------------------------------------
206
207                        case AsinOp:
208                        // sqrt(1 - x * x), asin(x)
209                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
210                        reverse_sparse_jacobian_unary_op(
211                                i_var, arg[0], var_sparsity
212                        );
213                        break;
214                        // -------------------------------------------------
215
216                        case AtanOp:
217                        // 1 + x * x, atan(x)
218                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
219                        reverse_sparse_jacobian_unary_op(
220                                i_var, arg[0], var_sparsity
221                        );
222                        break;
223                        // -------------------------------------------------
224
225                        case BeginOp:
226                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
227                        break;
228                        // -------------------------------------------------
229
230                        case CSumOp:
231                        // CSumOp has a variable number of arguments and
232                        // next_reverse thinks it one has one argument.
233                        // We must inform next_reverse of this special case.
234                        play->reverse_csum(op, arg, i_op, i_var);
235                        reverse_sparse_jacobian_csum_op(
236                                i_var, arg, var_sparsity
237                        );
238                        break;
239                        // -------------------------------------------------
240
241                        case CExpOp:
242                        reverse_sparse_jacobian_cond_op(
243                                i_var, arg, num_par, var_sparsity
244                        );
245                        break;
246                        // ---------------------------------------------------
247
248                        case ComOp:
249                        CPPAD_ASSERT_NARG_NRES(op, 4, 0);
250                        CPPAD_ASSERT_UNKNOWN( arg[1] > 1 );
251                        break;
252                        // --------------------------------------------------
253
254                        case CosOp:
255                        // sin(x), cos(x)
256                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
257                        reverse_sparse_jacobian_unary_op(
258                                i_var, arg[0], var_sparsity
259                        );
260                        break;
261                        // ---------------------------------------------------
262
263                        case CoshOp:
264                        // sinh(x), cosh(x)
265                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
266                        reverse_sparse_jacobian_unary_op(
267                                i_var, arg[0], var_sparsity
268                        );
269                        break;
270                        // -------------------------------------------------
271
272                        case DisOp:
273                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
274                        // derivative is identically zero
275                        break;
276                        // -------------------------------------------------
277
278                        case DivvvOp:
279                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
280                        reverse_sparse_jacobian_binary_op(
281                                i_var, arg, var_sparsity
282                        );
283                        break;
284                        // -------------------------------------------------
285
286                        case DivpvOp:
287                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
288                        reverse_sparse_jacobian_unary_op(
289                                i_var, arg[1], var_sparsity
290                        );
291                        break;
292                        // -------------------------------------------------
293
294                        case DivvpOp:
295                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
296                        reverse_sparse_jacobian_unary_op(
297                                i_var, arg[0], var_sparsity
298                        );
299                        break;
300                        // -------------------------------------------------
301
302                        case ExpOp:
303                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
304                        reverse_sparse_jacobian_unary_op(
305                                i_var, arg[0], var_sparsity
306                        );
307                        break;
308                        // -------------------------------------------------
309
310                        case InvOp:
311                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
312                        break;
313                        // -------------------------------------------------
314
315                        case LdpOp:
316                        reverse_sparse_jacobian_load_op(
317                                op,
318                                i_var,
319                                arg,
320                                num_vecad_ind,
321                                vecad_ind.data(),
322                                var_sparsity,
323                                vecad_sparsity
324                        );
325                        break;
326                        // -------------------------------------------------
327
328                        case LdvOp:
329                        reverse_sparse_jacobian_load_op(
330                                op,
331                                i_var,
332                                arg,
333                                num_vecad_ind,
334                                vecad_ind.data(),
335                                var_sparsity,
336                                vecad_sparsity
337                        );
338                        break;
339                        // -------------------------------------------------
340
341                        case LogOp:
342                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
343                        reverse_sparse_jacobian_unary_op(
344                                i_var, arg[0], var_sparsity
345                        );
346                        break;
347                        // -------------------------------------------------
348
349                        case MulvvOp:
350                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
351                        reverse_sparse_jacobian_binary_op(
352                                i_var, arg, var_sparsity
353                        );
354                        break;
355                        // -------------------------------------------------
356
357                        case MulpvOp:
358                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
359                        reverse_sparse_jacobian_unary_op(
360                                i_var, arg[1], var_sparsity
361                        );
362                        break;
363                        // -------------------------------------------------
364
365                        case ParOp:
366                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
367
368                        break;
369                        // -------------------------------------------------
370
371                        case PowvpOp:
372                        reverse_sparse_jacobian_unary_op(
373                                i_var, arg[0], var_sparsity
374                        );
375                        break;
376                        // -------------------------------------------------
377
378                        case PowpvOp:
379                        CPPAD_ASSERT_NARG_NRES(op, 2, 3);
380                        reverse_sparse_jacobian_unary_op(
381                                i_var, arg[1], var_sparsity
382                        );
383                        break;
384                        // -------------------------------------------------
385
386                        case PowvvOp:
387                        CPPAD_ASSERT_NARG_NRES(op, 2, 3);
388                        reverse_sparse_jacobian_binary_op(
389                                i_var, arg, var_sparsity
390                        );
391                        break;
392                        // -------------------------------------------------
393
394                        case PriOp:
395                        CPPAD_ASSERT_NARG_NRES(op, 5, 0);
396                        break;
397                        // -------------------------------------------------
398
399                        case SignOp:
400                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
401                        // derivative is identically zero
402                        break;
403                        // -------------------------------------------------
404
405                        case SinOp:
406                        // cos(x), sin(x)
407                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
408                        reverse_sparse_jacobian_unary_op(
409                                i_var, arg[0], var_sparsity
410                        );
411                        break;
412                        // -------------------------------------------------
413
414                        case SinhOp:
415                        // cosh(x), sinh(x)
416                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
417                        reverse_sparse_jacobian_unary_op(
418                                i_var, arg[0], var_sparsity
419                        );
420                        break;
421                        // -------------------------------------------------
422
423                        case SqrtOp:
424                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
425                        reverse_sparse_jacobian_unary_op(
426                                i_var, arg[0], var_sparsity
427                        );
428                        break;
429                        // -------------------------------------------------
430
431                        case StppOp:
432                        // sparsity cannot proagate through a parameter
433                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
434                        break;
435                        // -------------------------------------------------
436
437                        case StpvOp:
438                        reverse_sparse_jacobian_store_op(
439                                op,
440                                arg,
441                                num_vecad_ind,
442                                vecad_ind.data(),
443                                var_sparsity,
444                                vecad_sparsity
445                        );
446                        break;
447                        // -------------------------------------------------
448
449                        case StvpOp:
450                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
451                        break;
452                        // -------------------------------------------------
453
454                        case StvvOp:
455                        reverse_sparse_jacobian_store_op(
456                                op,
457                                arg,
458                                num_vecad_ind,
459                                vecad_ind.data(),
460                                var_sparsity,
461                                vecad_sparsity
462                        );
463                        break;
464                        // -------------------------------------------------
465
466                        case SubvvOp:
467                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
468                        reverse_sparse_jacobian_binary_op(
469                                i_var, arg, var_sparsity
470                        );
471                        break;
472                        // -------------------------------------------------
473
474                        case SubpvOp:
475                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
476                        reverse_sparse_jacobian_unary_op(
477                                i_var, arg[1], var_sparsity
478                        );
479                        break;
480                        // -------------------------------------------------
481
482                        case SubvpOp:
483                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
484                        reverse_sparse_jacobian_unary_op(
485                                i_var, arg[0], var_sparsity
486                        );
487                        break;
488                        // -------------------------------------------------
489
490                        case TanOp:
491                        // tan(x)^2, tan(x)
492                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
493                        reverse_sparse_jacobian_unary_op(
494                                i_var, arg[0], var_sparsity
495                        );
496                        break;
497                        // -------------------------------------------------
498
499                        case TanhOp:
500                        // tanh(x)^2, tanh(x)
501                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
502                        reverse_sparse_jacobian_unary_op(
503                                i_var, arg[0], var_sparsity
504                        );
505                        break;
506                        // -------------------------------------------------
507
508                        case UserOp:
509                        // start or end atomic operation sequence
510                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
511                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
512                        if( user_state == user_end )
513                        {       user_index = arg[0];
514                                user_id    = arg[1];
515                                user_n     = arg[2];
516                                user_m     = arg[3];
517                                if(user_r.size() < user_n )
518                                        user_r.resize(user_n);
519                                if(user_s.size() < user_m )
520                                        user_s.resize(user_m);
521                                user_j     = user_n;
522                                user_i     = user_m;
523                                user_state = user_ret;
524                        }
525                        else
526                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_start );
527                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
528                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
529                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
530                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
531                                user_state = user_end;
532               }
533                        break;
534
535                        case UsrapOp:
536                        // parameter argument in an atomic operation sequence
537                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
538                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
539                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
540                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
541                        --user_j;
542                        if( user_j == 0 )
543                                user_state = user_start;
544                        break;
545
546                        case UsravOp:
547                        // variable argument in an atomic operation sequence
548                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
549                        CPPAD_ASSERT_UNKNOWN( 0 < user_j && user_j <= user_n );
550                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
551                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
552                        CPPAD_ASSERT_UNKNOWN( 0 < arg[0] );
553                        --user_j;
554                        // It might be faster if we add set union to var_sparsity
555                        // where one of the sets is not in var_sparsity.
556                        set_itr = user_r[user_j].begin();
557                        set_end = user_r[user_j].end();
558                        while( set_itr != set_end )
559                                var_sparsity.add_element(arg[0], *set_itr++);   
560                        if( user_j == 0 )
561                                user_state = user_start;
562                        break;
563
564                        case UsrrpOp:
565                        // parameter result in an atomic operation sequence
566                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
567                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
568                        CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
569                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
570                        --user_i;
571                        user_s[user_i].clear();
572                        if( user_i == 0 )
573                        {       // call users function for this operation
574                                user_atomic<Base>::rev_jac_sparse(user_index, user_id,
575                                        user_n, user_m, user_q, user_r, user_s
576                                );
577                                user_state = user_arg;
578                        }
579                        break;
580
581                        case UsrrvOp:
582                        // variable result in an atomic operation sequence
583                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
584                        CPPAD_ASSERT_UNKNOWN( 0 < user_i && user_i <= user_m );
585                        --user_i;
586                        user_s[user_i].clear();
587                        var_sparsity.begin(i_var);
588                        i = var_sparsity.next_element();
589                        while( i < user_q )
590                        {       user_s[user_i].insert(i);
591                                i = var_sparsity.next_element();
592                        }
593                        if( user_i == 0 )
594                        {       // call users function for this operation
595                                user_atomic<Base>::rev_jac_sparse(user_index, user_id,
596                                        user_n, user_m, user_q, user_r, user_s
597                                );
598                                user_state = user_arg;
599                        }
600                        break;
601                        // -------------------------------------------------
602
603                        default:
604                        CPPAD_ASSERT_UNKNOWN(0);
605                }
606# if CPPAD_REV_JAC_SWEEP_TRACE
607                for(j = 0; j < limit; j++)
608                        z_value[j] = false;
609                var_sparsity.begin(i_var);
610                j = var_sparsity.next_element();
611                while( j < limit )
612                {       z_value[j] = true;
613                        j          = var_sparsity.next_element();
614                }
615                printOp(
616                        std::cout, 
617                        play,
618                        i_var,
619                        op, 
620                        arg,
621                        0, 
622                        (CppAD::vector<bool> *) CPPAD_NULL,
623                        1, 
624                        &z_value
625                );
626# endif
627        }
628        // values corresponding to BeginOp
629        CPPAD_ASSERT_UNKNOWN( i_op == 0 );
630        CPPAD_ASSERT_UNKNOWN( i_var == 0 );
631
632        return;
633}
634/*! \} */
635CPPAD_END_NAMESPACE
636
637// preprocessor symbols that are local to this file
638# undef CPPAD_REV_JAC_SWEEP_TRACE
639
640# endif
Note: See TracBrowser for help on using the repository browser.