Changeset 1608 for trunk/test_more


Ignore:
Timestamp:
Dec 12, 2009 2:10:19 PM (10 years ago)
Author:
bradbell
Message:

trunk: Optimization conversion of an addition sequence to one variable.

optimize.cpp: Test addition of new operators for this purpose.
whats_new_09.omh: Changes from the user's point of view.
optimize.hpp: Implement detection and conversion of addition sequences.

Add the new operators CAddOp, CSubOp, and CSumOp:
for_jac_sweep.hpp, op_code.hpp, forward0sweep.hpp, forward_sweep.hpp,
rev_jac_sweep.hpp, rev_hes_sweep.hpp, reverse_sweep.hpp

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/test_more/optimize.cpp

    r1607 r1608  
    318318
    319319                // result vector
    320                 y[0] = a1 + a2;
    321                 y[1] = b1 + b2;
    322                 y[2] = c1 + c2;
    323                 y[3] = d1 + d2;
    324                 y[4] = e1 + e2;
    325                 y[5] = f1 + f2;
    326                 y[6] = g1 + g2;
     320                y[0] = a1 * a2;
     321                y[1] = b1 * b2;
     322                y[2] = c1 * c2;
     323                y[3] = d1 * d2;
     324                y[4] = e1 * e2;
     325                y[5] = f1 * f2;
     326                y[6] = g1 * g2;
    327327                original += 7;
    328328                opt      += 7;
     
    509509                        ok &= ( y[i] == Value( Y[i] ) );
    510510
     511                return ok;
     512        }
     513        // -------------------------------------------------------------------
     514        bool cummulative_sum(void)
     515        {       // test conversion of a sequence of additions and subtraction
     516                // to a cummulative summation sequence.
     517                bool ok = true;
     518                using CppAD::AD;
     519                size_t i, j;
     520
     521                // domain space vector
     522                size_t n  = 7;
     523                CPPAD_TEST_VECTOR< AD<double> > X(n);
     524                for(j = 0; j < n; j++)
     525                        X[j] = double(j + 2);
     526
     527                size_t n_original = 1 + n;
     528                size_t n_optimize = 1 + n;
     529
     530                // range space vector
     531                size_t m = 2;
     532                CPPAD_TEST_VECTOR< AD<double> > Y(m);
     533
     534                // declare independent variables and start tape recording
     535                CppAD::Independent(X);
     536
     537                // Operations inside of optimize_cadd
     538                Y[0] = 5. + (X[0] + X[1]) + (X[1] - X[2]) // Addvv, Subvv
     539                     + (X[2] - 1.) + (2. - X[3])   // Subvp, Subpv
     540                     + (X[4] + 3.) + (4. + X[5]);  // Addpv, Addpv (no Addvp)
     541                n_original += 12;
     542                n_optimize += 1;
     543
     544
     545                // Operations inside of optimize_csub
     546                Y[1] = 5. - (X[1] + X[2]) - (X[2] - X[3]) // Addvv, Subvv
     547                     - (X[3] - 1.) - (2. - X[4])   // Subvp, Subpv
     548                     - (X[5] + 3.) - (4. + X[6]);  // Addpv, Addpv (no Addvp)
     549                n_original += 12;
     550                n_optimize += 1;
     551
     552                CppAD::ADFun<double> F;
     553                F.Dependent(X, Y);
     554
     555                // check number of variables in original function
     556                ok &= (F.size_var() ==  n_original );
     557       
     558                CPPAD_TEST_VECTOR<double> x(n), y(m);
     559                for(j = 0; j < n; j++)
     560                        x[j] = double(j + 2);
     561
     562                y   = F.Forward(0, x);
     563                for(i = 0; i < m; i++)
     564                        ok &= ( y[i] == Value( Y[i] ) );
     565
     566                F.optimize();
     567
     568                // check number of variables  in optimized version
     569                ok &= (F.size_var() == n_optimize );
     570
     571                y   = F.Forward(0, x);
     572                for(i = 0; i < m; i++)
     573                        ok &= ( y[i] == Value( Y[i] ) );
     574
     575                return ok;
     576        }
     577        // -------------------------------------------------------------------
     578        bool forward_csum(void)
     579        {       bool ok = true;
     580       
     581                using namespace CppAD;
     582       
     583                // independent variable vector
     584                CPPAD_TEST_VECTOR< AD<double> > X(2);
     585                X[0] = 0.;
     586                X[1] = 1.;
     587                Independent(X);
     588       
     589                // compute sum of elements in X
     590                CPPAD_TEST_VECTOR< AD<double> > Y(1);
     591                Y[0] = X[0] + X[0] + X[1];
     592       
     593                // create function object F : X -> Y
     594                ADFun<double> F(X, Y);
     595       
     596                // now optimize the operation sequence
     597                F.optimize();
     598       
     599                // use zero order to evaluate F[ (3, 4) ]
     600                CPPAD_TEST_VECTOR<double>  x0( F.Domain() );
     601                CPPAD_TEST_VECTOR<double>  y0( F.Range() );
     602                x0[0]    = 3.;
     603                x0[1]    = 4.;
     604                y0       = F.Forward(0, x0);
     605                ok      &= NearEqual(y0[0] , x0[0]+x0[0]+x0[1], 1e-10, 1e-10);
     606       
     607                // evaluate derivative of F in X[0] direction
     608                CPPAD_TEST_VECTOR<double> x1( F.Domain() );
     609                CPPAD_TEST_VECTOR<double> y1( F.Range() );
     610                x1[0]    = 1.;
     611                x1[1]    = 0.;
     612                y1       = F.Forward(1, x1);
     613                ok      &= NearEqual(y1[0] , x1[0]+x1[0]+x1[1], 1e-10, 1e-10);
     614       
     615                // evaluate second derivative of F in X[0] direction
     616                CPPAD_TEST_VECTOR<double> x2( F.Domain() );
     617                CPPAD_TEST_VECTOR<double> y2( F.Range() );
     618                x2[0]       = 0.;
     619                x2[1]       = 0.;
     620                y2          = F.Forward(2, x2);
     621                double F_00 = 2. * y2[0];
     622                ok         &= NearEqual(F_00, 0., 1e-10, 1e-10);
     623       
     624                return ok;
     625        }
     626        // -------------------------------------------------------------------
     627        bool reverse_csum(void)
     628        {       bool ok = true;
     629       
     630                using namespace CppAD;
     631       
     632                // independent variable vector
     633                CPPAD_TEST_VECTOR< AD<double> > X(2);
     634                X[0] = 0.;
     635                X[1] = 1.;
     636                Independent(X);
     637       
     638                // compute sum of elements in X
     639                CPPAD_TEST_VECTOR< AD<double> > Y(1);
     640                Y[0] = X[0] - (X[0] - X[1]);
     641       
     642                // create function object F : X -> Y
     643                ADFun<double> F(X, Y);
     644       
     645                // now optimize the operation sequence
     646                F.optimize();
     647       
     648                // use zero order to evaluate F[ (3, 4) ]
     649                CPPAD_TEST_VECTOR<double>  x0( F.Domain() );
     650                CPPAD_TEST_VECTOR<double>  y0( F.Range() );
     651                x0[0]    = 3.;
     652                x0[1]    = 4.;
     653                y0       = F.Forward(0, x0);
     654                ok      &= NearEqual(y0[0] , x0[0]-x0[0]+x0[1], 1e-10, 1e-10);
     655       
     656                // evaluate derivative of F
     657                CPPAD_TEST_VECTOR<double> dF( F.Domain() );
     658                CPPAD_TEST_VECTOR<double> w( F.Range() );
     659                w[0]    = 1.;
     660                dF      = F.Reverse(1, w);
     661                ok     &= NearEqual(dF[0] , 0., 1e-10, 1e-10);
     662                ok     &= NearEqual(dF[1] , 1., 1e-10, 1e-10);
     663       
     664                return ok;
     665        }
     666        bool forward_sparse_jacobian_csum()
     667        {       bool ok = true;
     668                using namespace CppAD;
     669       
     670                // dimension of the domain space
     671                size_t n = 3;
     672       
     673                // dimension of the range space
     674                size_t m = 2;
     675       
     676                // independent variable vector
     677                CPPAD_TEST_VECTOR< AD<double> > X(n);
     678                X[0] = 2.;
     679                X[1] = 3.;
     680                Independent(X);
     681       
     682                // dependent variable vector
     683                CPPAD_TEST_VECTOR< AD<double> > Y(m);
     684       
     685                // check results vector
     686                CPPAD_TEST_VECTOR< bool >       Check(m * n);
     687       
     688                // initialize index into Y
     689                size_t index = 0;
     690       
     691                // Y[0]
     692                Y[index]             = X[0] + X[1] + 5.;
     693                Check[index * n + 0] = true;
     694                Check[index * n + 1] = true;
     695                Check[index * n + 2] = false;
     696                index++;
     697       
     698                // Y[1]
     699                Y[index]             = Y[0] - (X[1] + X[2]);
     700                Check[index * n + 0] = true;
     701                Check[index * n + 1] = true;
     702                Check[index * n + 2] = true;
     703                index++;
     704       
     705                // check final index
     706                assert( index == m );
     707       
     708                // create function object F : X -> Y
     709                ADFun<double> F(X, Y);
     710                F.optimize();
     711       
     712                // ---------------------------------------------------------
     713                // dependency matrix for the identity function
     714                CPPAD_TEST_VECTOR< std::set<size_t> > Sx(n);
     715                size_t i, j;
     716                for(i = 0; i < n; i++)
     717                {       assert( Sx[i].empty() );
     718                        Sx[i].insert(i);
     719                }
     720       
     721                // evaluate the dependency matrix for F(x)
     722                CPPAD_TEST_VECTOR< std::set<size_t> > Sy(m);
     723                Sy = F.ForSparseJac(n, Sx);
     724       
     725                // check values
     726                bool found;
     727                for(i = 0; i < m; i++)
     728                {       for(j = 0; j < n; j++)
     729                        {       found = Sy[i].find(j) != Sy[i].end();
     730                                ok &= (found == Check[i * n + j]);
     731                        }
     732                }       
     733       
     734                return ok;
     735        }
     736        bool reverse_sparse_jacobian_csum()
     737        {       bool ok = true;
     738                using namespace CppAD;
     739       
     740                // dimension of the domain space
     741                size_t n = 3;
     742       
     743                // dimension of the range space
     744                size_t m = 2;
     745       
     746                // independent variable vector
     747                CPPAD_TEST_VECTOR< AD<double> > X(n);
     748                X[0] = 2.;
     749                X[1] = 3.;
     750                Independent(X);
     751       
     752                // dependent variable vector
     753                CPPAD_TEST_VECTOR< AD<double> > Y(m);
     754       
     755                // check results vector
     756                CPPAD_TEST_VECTOR< bool >       Check(m * n);
     757       
     758                // initialize index into Y
     759                size_t index = 0;
     760       
     761                // Y[0]
     762                Y[index]             = X[0] + X[1] + 5.;
     763                Check[index * n + 0] = true;
     764                Check[index * n + 1] = true;
     765                Check[index * n + 2] = false;
     766                index++;
     767       
     768                // Y[1]
     769                Y[index]             = Y[0] - (X[1] + X[2]);
     770                Check[index * n + 0] = true;
     771                Check[index * n + 1] = true;
     772                Check[index * n + 2] = true;
     773                index++;
     774       
     775                // check final index
     776                assert( index == m );
     777       
     778                // create function object F : X -> Y
     779                ADFun<double> F(X, Y);
     780                F.optimize();
     781       
     782                // ----------------------------------------------------------
     783                // dependency matrix for the identity function
     784                CPPAD_TEST_VECTOR< bool > Py(m * m);
     785                size_t i, j;
     786                for(i = 0; i < m; i++)
     787                {       for(j = 0; j < m; j++)
     788                                Py[ i * m + j ] = (i == j);
     789                }
     790       
     791                // evaluate the dependency matrix for F(x)
     792                CPPAD_TEST_VECTOR< bool > Px(m * n);
     793                Px = F.RevSparseJac(m, Py);
     794       
     795                // check values
     796                for(i = 0; i < m; i++)
     797                {       for(j = 0; j < n; j++)
     798                                ok &= (Px[i * n + j] == Check[i * n + j]);
     799                }       
     800       
     801                return ok;
     802        }
     803        bool reverse_sparse_hessian_csum(void)
     804        {       bool ok = true;
     805                using CppAD::AD;
     806                size_t i, j;
     807       
     808                size_t n = 2;
     809                CPPAD_TEST_VECTOR< AD<double> > X(n);
     810                X[0] = 1.;
     811                X[1] = 2.;
     812                CppAD::Independent(X);
     813       
     814                size_t m = 1;
     815                CPPAD_TEST_VECTOR< AD<double> > Y(m);
     816                Y[0] = (2. + X[0] + X[1] + 3.) * X[0];
     817       
     818                CPPAD_TEST_VECTOR<bool> check(n * n);
     819                check[0 * n + 0] = true;  // partial w.r.t. x[0], x[0]
     820                check[0 * n + 1] = true;  //                x[0], x[1]
     821                check[1 * n + 0] = true;  //                x[1], x[0]
     822                check[1 * n + 1] = false; //                x[1], x[1]
     823       
     824                // create function object F : X -> Y
     825                CppAD::ADFun<double> F(X, Y);
     826                F.optimize();
     827       
     828                // sparsity pattern for the identity function U(x) = x
     829                CPPAD_TEST_VECTOR<bool> Px(n * n);
     830                for(i = 0; i < n; i++)
     831                        for(j = 0; j < n; j++)
     832                                Px[ i * n + j ] = (i == j);
     833       
     834                // compute sparsity pattern for Jacobian of F(U(x))
     835                CPPAD_TEST_VECTOR<bool> P_jac(m * n);
     836                P_jac = F.ForSparseJac(n, Px);
     837       
     838                // compute sparsity pattern for Hessian of F_k ( U(x) )
     839                CPPAD_TEST_VECTOR<bool> Py(m);
     840                CPPAD_TEST_VECTOR<bool> Pxx(n * n);
     841                Py[0] = true;
     842                Pxx = F.RevSparseHes(n, Py);
     843                // check values
     844                for(i = 0; i < n * n; i++)
     845                        ok &= (Pxx[i] == check[i]);
     846       
    511847                return ok;
    512848        }
     
    523859        ok     &= duplicate_two();
    524860        ok     &= duplicate_three();
     861        // convert sequence of additions to cummulative summation
     862        ok     &= cummulative_sum();
     863        ok     &= forward_csum();
     864        ok     &= reverse_csum();
     865        ok     &= forward_sparse_jacobian_csum();
     866        ok     &= reverse_sparse_jacobian_csum();
     867        ok     &= reverse_sparse_hessian_csum();
    525868        return ok;
    526869}
Note: See TracChangeset for help on using the changeset viewer.