source: trunk/cppad/local/for_jac_sweep.hpp @ 2625

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).

ad_assign.hpp: fix doxygen brief description.
ad_ctor.hpp: fix doxygen brief description.

  • 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 $ */
2# ifndef CPPAD_FOR_JAC_SWEEP_INCLUDED
3# define CPPAD_FOR_JAC_SWEEP_INCLUDED
4/* --------------------------------------------------------------------------
5CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 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# include <set>
15# include <cppad/local/pod_vector.hpp>
16
17CPPAD_BEGIN_NAMESPACE
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/*!
26\def CPPAD_FOR_JAC_SWEEP_TRACE
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*/
31# define CPPAD_FOR_JAC_SWEEP_TRACE 0
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
106        // vecad_sparsity contains a sparsity pattern from each VecAD object
107        // to all the other variables.
108        // vecad_ind maps a VecAD index (the beginning of the
109        // VecAD object) to its from index in vecad_sparsity
110        size_t num_vecad_ind   = play->num_rec_vecad_ind();
111        size_t num_vecad_vec   = play->num_rec_vecad_vec();
112        Vector_set  vecad_sparsity;
113        vecad_sparsity.resize(num_vecad_vec, limit);
114        pod_vector<size_t> vecad_ind;
115        if( num_vecad_vec > 0 )
116        {       size_t length;
117                vecad_ind.extend(num_vecad_ind);
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
123                        vecad_ind[j] = i; 
124                        for(k = 1; k <= length; k++)
125                                vecad_ind[j+k] = num_vecad_vec; // invalid index
126                        // start of next VecAD
127                        j       += length + 1;
128                }
129                CPPAD_ASSERT_UNKNOWN( j == play->num_rec_vecad_ind() );
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
148# if CPPAD_FOR_JAC_SWEEP_TRACE
149        std::cout << std::endl;
150        CppAD::vector<bool> z_value(limit);
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:
167                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
168                        forward_sparse_jacobian_unary_op(
169                                i_var, arg[0], var_sparsity
170                        );
171                        break;
172                        // -------------------------------------------------
173
174                        case AddvvOp:
175                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
176                        forward_sparse_jacobian_binary_op(
177                                i_var, arg, var_sparsity
178                        );
179                        break;
180                        // -------------------------------------------------
181
182                        case AddpvOp:
183                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
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)
192                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
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)
201                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
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)
210                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
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:
236                        CPPAD_ASSERT_NARG_NRES(op, 4, 0);
237                        CPPAD_ASSERT_UNKNOWN( arg[1] > 1 );
238                        break;
239                        // --------------------------------------------------
240
241                        case CosOp:
242                        // sin(x), cos(x)
243                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
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)
252                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
253                        forward_sparse_jacobian_unary_op(
254                                i_var, arg[0], var_sparsity
255                        );
256                        break;
257                        // -------------------------------------------------
258
259                        case DisOp:
260                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
261                        var_sparsity.clear(i_var);
262                        break;
263                        // -------------------------------------------------
264
265                        case DivvvOp:
266                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
267                        forward_sparse_jacobian_binary_op(
268                                i_var, arg, var_sparsity
269                        );
270                        break;
271                        // -------------------------------------------------
272
273                        case DivpvOp:
274                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
275                        forward_sparse_jacobian_unary_op(
276                                i_var, arg[1], var_sparsity
277                        );
278                        break;
279                        // -------------------------------------------------
280
281                        case DivvpOp:
282                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
283                        forward_sparse_jacobian_unary_op(
284                                i_var, arg[0], var_sparsity
285                        );
286                        break;
287                        // -------------------------------------------------
288
289                        case EndOp:
290                        CPPAD_ASSERT_NARG_NRES(op, 0, 0);
291                        break;
292                        // -------------------------------------------------
293
294                        case ExpOp:
295                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
296                        forward_sparse_jacobian_unary_op(
297                                i_var, arg[0], var_sparsity
298                        );
299                        break;
300                        // -------------------------------------------------
301
302                        case InvOp:
303                        CPPAD_ASSERT_NARG_NRES(op, 0, 1);
304                        // sparsity pattern is already defined
305                        break;
306                        // -------------------------------------------------
307
308                        case LdpOp:
309                        forward_sparse_load_op(
310                                op,
311                                i_var,
312                                arg,
313                                num_vecad_ind,
314                                vecad_ind.data(),
315                                var_sparsity,
316                                vecad_sparsity
317                        );
318                        break;
319                        // -------------------------------------------------
320
321                        case LdvOp:
322                        forward_sparse_load_op(
323                                op,
324                                i_var,
325                                arg,
326                                num_vecad_ind,
327                                vecad_ind.data(),
328                                var_sparsity,
329                                vecad_sparsity
330                        );
331                        break;
332                        // -------------------------------------------------
333
334                        case LogOp:
335                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
336                        forward_sparse_jacobian_unary_op(
337                                i_var, arg[0], var_sparsity
338                        );
339                        break;
340                        // -------------------------------------------------
341
342                        case MulvvOp:
343                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
344                        forward_sparse_jacobian_binary_op(
345                                i_var, arg, var_sparsity
346                        );
347                        break;
348                        // -------------------------------------------------
349
350                        case MulpvOp:
351                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
352                        forward_sparse_jacobian_unary_op(
353                                i_var, arg[1], var_sparsity
354                        );
355                        break;
356                        // -------------------------------------------------
357
358                        case ParOp:
359                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
360                        var_sparsity.clear(i_var);
361                        break;
362                        // -------------------------------------------------
363
364                        case PowvpOp:
365                        CPPAD_ASSERT_NARG_NRES(op, 2, 3);
366                        forward_sparse_jacobian_unary_op(
367                                i_var, arg[0], var_sparsity
368                        );
369                        break;
370                        // -------------------------------------------------
371
372                        case PowpvOp:
373                        CPPAD_ASSERT_NARG_NRES(op, 2, 3);
374                        forward_sparse_jacobian_unary_op(
375                                i_var, arg[1], var_sparsity
376                        );
377                        break;
378                        // -------------------------------------------------
379
380                        case PowvvOp:
381                        CPPAD_ASSERT_NARG_NRES(op, 2, 3);
382                        forward_sparse_jacobian_binary_op(
383                                i_var, arg, var_sparsity
384                        );
385                        break;
386                        // -------------------------------------------------
387
388                        case PriOp:
389                        CPPAD_ASSERT_NARG_NRES(op, 5, 0);
390                        break;
391                        // -------------------------------------------------
392
393                        case SignOp:
394                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
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)
403                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
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)
412                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
413                        forward_sparse_jacobian_unary_op(
414                                i_var, arg[0], var_sparsity
415                        );
416                        break;
417                        // -------------------------------------------------
418
419                        case SqrtOp:
420                        CPPAD_ASSERT_NARG_NRES(op, 1, 1);
421                        forward_sparse_jacobian_unary_op(
422                                i_var, arg[0], var_sparsity
423                        );
424                        break;
425                        // -------------------------------------------------
426
427                        case StppOp:
428                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
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,
437                                num_vecad_ind,
438                                vecad_ind.data(),
439                                var_sparsity,
440                                vecad_sparsity
441                        );
442                        break;
443                        // -------------------------------------------------
444
445                        case StvpOp:
446                        CPPAD_ASSERT_NARG_NRES(op, 3, 0);
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,
455                                num_vecad_ind,
456                                vecad_ind.data(),
457                                var_sparsity,
458                                vecad_sparsity
459                        );
460                        break;
461                        // -------------------------------------------------
462
463                        case SubvvOp:
464                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
465                        forward_sparse_jacobian_binary_op(
466                                i_var, arg, var_sparsity
467                        );
468                        break;
469                        // -------------------------------------------------
470
471                        case SubpvOp:
472                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
473                        forward_sparse_jacobian_unary_op(
474                                i_var, arg[1], var_sparsity
475                        );
476                        break;
477                        // -------------------------------------------------
478
479                        case SubvpOp:
480                        CPPAD_ASSERT_NARG_NRES(op, 2, 1);
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)
489                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
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)
498                        CPPAD_ASSERT_NARG_NRES(op, 1, 2);
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 )
590                                var_sparsity.add_element(i_var, *set_itr++);
591                        user_i++;
592                        if( user_i == user_m )
593                                user_state = user_end;
594                        break;
595                        // -------------------------------------------------
596
597                        default:
598                        CPPAD_ASSERT_UNKNOWN(0);
599                }
600# if CPPAD_FOR_JAC_SWEEP_TRACE
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,
619                        (CppAD::vector<bool> *) CPPAD_NULL
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/*! \} */
632CPPAD_END_NAMESPACE
633
634// preprocessor symbols that are local to this file
635# undef CPPAD_FOR_JAC_SWEEP_TRACE
636
637# endif
Note: See TracBrowser for help on using the repository browser.