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

Change comment so end of group is included in doxygen processing
(do not know why works without this command at end).

• Property svn:keywords set to Id
File size: 16.8 KB
Line
1/* $Id: for_jac_sweep.hpp 2625 2012-12-23 14:34:12Z bradbell$ */
4/* --------------------------------------------------------------------------
6
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.
13-------------------------------------------------------------------------- */
14# include <set>
16
18/*!
19\defgroup for_jac_sweep_hpp for_jac_sweep.hpp
20\{
21\file for_jac_sweep.hpp
22Compute Forward mode Jacobian sparsity patterns.
23*/
24
25/*!
27This value is either zero or one.
28Zero is the normal operational value.
29If it is one, a trace of every for_jac_sweep computation is printed.
30*/
32
33/*!
34Given the sparsity pattern for the independent variables,
35ForJacSweep computes the sparsity pattern for all the other variables.
36
37\tparam Base
38base type for the operator; i.e., this operation sequence was recorded
39using AD< \a Base > and computations by this routine are done using type
40\a Base.
41
42\tparam Vector_set
43is the type used for vectors of sets. It can be either
44\c sparse_pack, \c sparse_set, or \c sparse_list.
45
46\param n
47is the number of independent variables on the tape.
48
49\param numvar
50is the total number of variables on the tape; i.e.,
51\a play->num_rec_var().
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
67\b Input: For j = 1 , ... , \a n,
68the sparsity pattern for the independent variable with index (j-1)
69corresponds to the set with index j in \a var_sparsity.
70\n
71\n
72\b Output: For i = \a n + 1 , ... , \a numvar - 1,
73the sparsity pattern for the variable with index i on the tape
74corresponds to the set with index i in \a var_sparsity.
75
76\par Checked Assertions:
77\li numvar == var_sparsity.n_set()
78\li numvar == play->num_rec_var()
79*/
80
81template <class Base, class Vector_set>
82void ForJacSweep(
83        size_t                n            ,
84        size_t                numvar       ,
85        player<Base>*         play         ,
86        Vector_set&           var_sparsity )
87{
88        OpCode           op;
89        size_t         i_op;
90        size_t        i_var;
91
92        const addr_t*   arg = 0;
93
94        size_t            i, j, k;
95
96        // check numvar argument
97        CPPAD_ASSERT_UNKNOWN( play->num_rec_var()  == numvar );
98        CPPAD_ASSERT_UNKNOWN( var_sparsity.n_set() == numvar );
99
100        // length of the parameter vector (used by CppAD assert macros)
101        const size_t num_par = play->num_rec_par();
102
103        // cum_sparsity accumulates sparsity pattern a cummulative sum
104        size_t limit = var_sparsity.end();
105
107        // to all the other variables.
115        if( num_vecad_vec > 0 )
116        {       size_t length;
118                j             = 0;
119                for(i = 0; i < num_vecad_vec; i++)
120                {       // length of this VecAD
121                        length   = play->GetVecInd(j);
122                        // set to proper index for this VecAD
124                        for(k = 1; k <= length; k++)
126                        // start of next VecAD
127                        j       += length + 1;
128                }
130        }
131
132        // work space used by UserOp.
133        typedef std::set<size_t> size_set;
134        const size_t user_q = limit; // maximum element plus one
135        size_set::iterator set_itr;  // iterator for a standard set
136        size_set::iterator set_end;  // end of iterator sequence
137        vector< size_set > user_r;   // sparsity pattern for the argument x
138        vector< size_set > user_s;   // sparisty pattern for the result y
139        size_t user_index = 0;       // indentifier for this user_atomic operation
140        size_t user_id    = 0;       // user identifier for this call to operator
141        size_t user_i     = 0;       // index in result vector
142        size_t user_j     = 0;       // index in argument vector
143        size_t user_m     = 0;       // size of result vector
144        size_t user_n     = 0;       // size of arugment vector
145        // next expected operator in a UserOp sequence
146        enum { user_start, user_arg, user_ret, user_end } user_state = user_start;
147
149        std::cout << std::endl;
151# endif
152
153        // skip the BeginOp at the beginning of the recording
154        play->start_forward(op, arg, i_op, i_var);
155        CPPAD_ASSERT_UNKNOWN( op == BeginOp );
156        while(op != EndOp)
157        {
158                // this op
159                play->next_forward(op, arg, i_op, i_var);
160                CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );
161                CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );
162
163                // rest of information depends on the case
164                switch( op )
165                {
166                        case AbsOp:
168                        forward_sparse_jacobian_unary_op(
169                                i_var, arg[0], var_sparsity
170                        );
171                        break;
172                        // -------------------------------------------------
173
176                        forward_sparse_jacobian_binary_op(
177                                i_var, arg, var_sparsity
178                        );
179                        break;
180                        // -------------------------------------------------
181
184                        forward_sparse_jacobian_unary_op(
185                                i_var, arg[1], var_sparsity
186                        );
187                        break;
188                        // -------------------------------------------------
189
190                        case AcosOp:
191                        // sqrt(1 - x * x), acos(x)
193                        forward_sparse_jacobian_unary_op(
194                                i_var, arg[0], var_sparsity
195                        );
196                        break;
197                        // -------------------------------------------------
198
199                        case AsinOp:
200                        // sqrt(1 - x * x), asin(x)
202                        forward_sparse_jacobian_unary_op(
203                                i_var, arg[0], var_sparsity
204                        );
205                        break;
206                        // -------------------------------------------------
207
208                        case AtanOp:
209                        // 1 + x * x, atan(x)
211                        forward_sparse_jacobian_unary_op(
212                                i_var, arg[0], var_sparsity
213                        );
214                        break;
215                        // -------------------------------------------------
216
217                        case CSumOp:
218                        // CSumOp has a variable number of arguments and
219                        // next_forward thinks it one has one argument.
220                        // we must inform next_forward of this special case.
221                        play->forward_csum(op, arg, i_op, i_var);
222                        forward_sparse_jacobian_csum_op(
223                                i_var, arg, var_sparsity
224                        );
225                        break;
226                        // -------------------------------------------------
227
228                        case CExpOp:
229                        forward_sparse_jacobian_cond_op(
230                                i_var, arg, num_par, var_sparsity
231                        );
232                        break;
233                        // ---------------------------------------------------
234
235                        case ComOp:
237                        CPPAD_ASSERT_UNKNOWN( arg[1] > 1 );
238                        break;
239                        // --------------------------------------------------
240
241                        case CosOp:
242                        // sin(x), cos(x)
244                        forward_sparse_jacobian_unary_op(
245                                i_var, arg[0], var_sparsity
246                        );
247                        break;
248                        // ---------------------------------------------------
249
250                        case CoshOp:
251                        // sinh(x), cosh(x)
253                        forward_sparse_jacobian_unary_op(
254                                i_var, arg[0], var_sparsity
255                        );
256                        break;
257                        // -------------------------------------------------
258
259                        case DisOp:
261                        var_sparsity.clear(i_var);
262                        break;
263                        // -------------------------------------------------
264
265                        case DivvvOp:
267                        forward_sparse_jacobian_binary_op(
268                                i_var, arg, var_sparsity
269                        );
270                        break;
271                        // -------------------------------------------------
272
273                        case DivpvOp:
275                        forward_sparse_jacobian_unary_op(
276                                i_var, arg[1], var_sparsity
277                        );
278                        break;
279                        // -------------------------------------------------
280
281                        case DivvpOp:
283                        forward_sparse_jacobian_unary_op(
284                                i_var, arg[0], var_sparsity
285                        );
286                        break;
287                        // -------------------------------------------------
288
289                        case EndOp:
291                        break;
292                        // -------------------------------------------------
293
294                        case ExpOp:
296                        forward_sparse_jacobian_unary_op(
297                                i_var, arg[0], var_sparsity
298                        );
299                        break;
300                        // -------------------------------------------------
301
302                        case InvOp:
304                        // sparsity pattern is already defined
305                        break;
306                        // -------------------------------------------------
307
308                        case LdpOp:
310                                op,
311                                i_var,
312                                arg,
315                                var_sparsity,
317                        );
318                        break;
319                        // -------------------------------------------------
320
321                        case LdvOp:
323                                op,
324                                i_var,
325                                arg,
328                                var_sparsity,
330                        );
331                        break;
332                        // -------------------------------------------------
333
334                        case LogOp:
336                        forward_sparse_jacobian_unary_op(
337                                i_var, arg[0], var_sparsity
338                        );
339                        break;
340                        // -------------------------------------------------
341
342                        case MulvvOp:
344                        forward_sparse_jacobian_binary_op(
345                                i_var, arg, var_sparsity
346                        );
347                        break;
348                        // -------------------------------------------------
349
350                        case MulpvOp:
352                        forward_sparse_jacobian_unary_op(
353                                i_var, arg[1], var_sparsity
354                        );
355                        break;
356                        // -------------------------------------------------
357
358                        case ParOp:
360                        var_sparsity.clear(i_var);
361                        break;
362                        // -------------------------------------------------
363
364                        case PowvpOp:
366                        forward_sparse_jacobian_unary_op(
367                                i_var, arg[0], var_sparsity
368                        );
369                        break;
370                        // -------------------------------------------------
371
372                        case PowpvOp:
374                        forward_sparse_jacobian_unary_op(
375                                i_var, arg[1], var_sparsity
376                        );
377                        break;
378                        // -------------------------------------------------
379
380                        case PowvvOp:
382                        forward_sparse_jacobian_binary_op(
383                                i_var, arg, var_sparsity
384                        );
385                        break;
386                        // -------------------------------------------------
387
388                        case PriOp:
390                        break;
391                        // -------------------------------------------------
392
393                        case SignOp:
395                        forward_sparse_jacobian_unary_op(
396                                i_var, arg[0], var_sparsity
397                        );
398                        break;
399                        // -------------------------------------------------
400
401                        case SinOp:
402                        // cos(x), sin(x)
404                        forward_sparse_jacobian_unary_op(
405                                i_var, arg[0], var_sparsity
406                        );
407                        break;
408                        // -------------------------------------------------
409
410                        case SinhOp:
411                        // cosh(x), sinh(x)
413                        forward_sparse_jacobian_unary_op(
414                                i_var, arg[0], var_sparsity
415                        );
416                        break;
417                        // -------------------------------------------------
418
419                        case SqrtOp:
421                        forward_sparse_jacobian_unary_op(
422                                i_var, arg[0], var_sparsity
423                        );
424                        break;
425                        // -------------------------------------------------
426
427                        case StppOp:
429                        // storing a parameter does not affect vector sparsity
430                        break;
431                        // -------------------------------------------------
432
433                        case StpvOp:
434                        forward_sparse_store_op(
435                                op,
436                                arg,
439                                var_sparsity,
441                        );
442                        break;
443                        // -------------------------------------------------
444
445                        case StvpOp:
447                        // storing a parameter does not affect vector sparsity
448                        break;
449                        // -------------------------------------------------
450
451                        case StvvOp:
452                        forward_sparse_store_op(
453                                op,
454                                arg,
457                                var_sparsity,
459                        );
460                        break;
461                        // -------------------------------------------------
462
463                        case SubvvOp:
465                        forward_sparse_jacobian_binary_op(
466                                i_var, arg, var_sparsity
467                        );
468                        break;
469                        // -------------------------------------------------
470
471                        case SubpvOp:
473                        forward_sparse_jacobian_unary_op(
474                                i_var, arg[1], var_sparsity
475                        );
476                        break;
477                        // -------------------------------------------------
478
479                        case SubvpOp:
481                        forward_sparse_jacobian_unary_op(
482                                i_var, arg[0], var_sparsity
483                        );
484                        break;
485                        // -------------------------------------------------
486
487                        case TanOp:
488                        // tan(x)^2, tan(x)
490                        forward_sparse_jacobian_unary_op(
491                                i_var, arg[0], var_sparsity
492                        );
493                        break;
494                        // -------------------------------------------------
495
496                        case TanhOp:
497                        // tanh(x)^2, tanh(x)
499                        forward_sparse_jacobian_unary_op(
500                                i_var, arg[0], var_sparsity
501                        );
502                        break;
503                        // -------------------------------------------------
504
505                        case UserOp:
506                        // start or end an atomic operation sequence
507                        CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
508                        CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
509                        if( user_state == user_start )
510                        {       user_index = arg[0];
511                                user_id    = arg[1];
512                                user_n     = arg[2];
513                                user_m     = arg[3];
514                                if(user_r.size() < user_n )
515                                        user_r.resize(user_n);
516                                if(user_s.size() < user_m )
517                                        user_s.resize(user_m);
518                                user_j     = 0;
519                                user_i     = 0;
520                                user_state = user_arg;
521                        }
522                        else
523                        {       CPPAD_ASSERT_UNKNOWN( user_state == user_end );
524                                CPPAD_ASSERT_UNKNOWN( user_index == size_t(arg[0]) );
525                                CPPAD_ASSERT_UNKNOWN( user_id    == size_t(arg[1]) );
526                                CPPAD_ASSERT_UNKNOWN( user_n     == size_t(arg[2]) );
527                                CPPAD_ASSERT_UNKNOWN( user_m     == size_t(arg[3]) );
528                                user_state = user_start;
529                        }
530                        break;
531
532                        case UsrapOp:
533                        // parameter argument in an atomic operation sequence
534                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
535                        CPPAD_ASSERT_UNKNOWN( user_j < user_n );
536                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) < num_par );
537                        // parameters have an empty sparsity pattern
538                        user_r[user_j].clear();
539                        ++user_j;
540                        if( user_j == user_n )
541                        {       // call users function for this operation
542                                user_atomic<Base>::for_jac_sparse(user_index, user_id,
543                                        user_n, user_m, user_q, user_r, user_s
544                                );
545                                user_state = user_ret;
546                        }
547                        break;
548
549                        case UsravOp:
550                        // variable argument in an atomic operation sequence
551                        CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
552                        CPPAD_ASSERT_UNKNOWN( user_j < user_n );
553                        CPPAD_ASSERT_UNKNOWN( size_t(arg[0]) <= i_var );
554                        // set user_r[user_j] to sparsity pattern for variable arg[0]
555                        user_r[user_j].clear();
556                        var_sparsity.begin(arg[0]);
557                        i = var_sparsity.next_element();
558                        while( i < user_q )
559                        {       user_r[user_j].insert(i);
560                                i = var_sparsity.next_element();
561                        }
562                        ++user_j;
563                        if( user_j == user_n )
564                        {       // call users function for this operation
565                                user_atomic<Base>::for_jac_sparse(user_index, user_id,
566                                        user_n, user_m, user_q, user_r, user_s
567                                );
568                                user_state = user_ret;
569                        }
570                        break;
571
572                        case UsrrpOp:
573                        // parameter result in an atomic operation sequence
574                        CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
575                        CPPAD_ASSERT_UNKNOWN( user_i < user_m );
576                        user_i++;
577                        if( user_i == user_m )
578                                user_state = user_end;
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( user_i < user_m );
585                        // It might be faster if we add set union to var_sparsity
586                        // where one of the sets is not in var_sparsity
587                        set_itr = user_s[user_i].begin();
588                        set_end = user_s[user_i].end();
589                        while( set_itr != set_end )
591                        user_i++;
592                        if( user_i == user_m )
593                                user_state = user_end;
594                        break;
595                        // -------------------------------------------------
596
597                        default:
599                }
601                // value for this variable
602                for(j = 0; j < limit; j++)
603                        z_value[j] = false;
604                var_sparsity.begin(i_var);
605                j = var_sparsity.next_element();
606                while( j < limit )
607                {       z_value[j] = true;
608                        j = var_sparsity.next_element();
609                }
610                printOp(
611                        std::cout,
612                        play,
613                        i_var,
614                        op,
615                        arg,
616                        1,
617                        &z_value,
618                        0,
620                );
621        }
622        std::cout << std::endl;
623# else
624        }
625# endif
626        CPPAD_ASSERT_UNKNOWN( i_var + 1 == play->num_rec_var() );
627
628        return;
629}
630
631/*! \} */