Changeset 1521 for trunk/test_more


Ignore:
Timestamp:
Sep 21, 2009 6:46:25 AM (11 years ago)
Author:
bradbell
Message:

trunk: Fix bug in VecAD Hessian sparsity, convert load and store.

rev_sparse_jac.cpp: add a VecAD test.
rev_sparse_hes.cpp: add a VecAD test that demonstrates bug.
whats_new_09.omh: user's view of the changes

Convert load and store sparse operations to use connection.
prototype_op.hpp, for_jac_sweep.hpp, rev_jac_sweep.hpp, rev_hes_sweep.hpp.

convert sparse operations to use connection.
load_op.hpp, store_op.hpp.

Location:
trunk/test_more
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/test_more/rev_sparse_hes.cpp

    r1474 r1521  
    4141
    4242        // accumulate sum here
    43         AD<double> sum;
     43        AD<double> sum(0.);
    4444
    4545        // variable * variable
     
    216216}
    217217
     218bool case_four(void)
     219{       bool ok = true;
     220        using namespace CppAD;
     221
     222        // dimension of the domain space
     223        size_t n = 3;
     224
     225        // dimension of the range space
     226        size_t m = 1;
     227
     228        // inialize the vector as zero
     229        CppAD::VecAD<double> Z(n - 1);
     230        size_t k;
     231        for(k = 0; k < n-1; k++)
     232                Z[k] = 0.;
     233
     234        // independent variable vector
     235        CPPAD_TEST_VECTOR< AD<double> > X(n);
     236        X[0] = 0.;
     237        X[1] = 1.;
     238        X[2] = 2.;
     239        Independent(X);
     240
     241        // VecAD vector z depends on both x[1] and x[2]
     242        // (component indices do not matter because they can change).
     243        Z[ X[0] ] = X[1] * X[2];
     244        Z[ X[1] ] = 0.;
     245
     246        // dependent variable vector
     247        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     248
     249        // check results vector
     250        CPPAD_TEST_VECTOR< bool >       Check(n * n);
     251
     252        // y = z[j] where j might be zero or one.
     253        Y[0]  =  Z[ X[1] ];
     254
     255        Check[0 * n + 0] = false; // partial w.r.t x[0], x[0]
     256        Check[0 * n + 1] = false; // partial w.r.t x[0], x[1]
     257        Check[0 * n + 2] = false; // partial w.r.t x[0], x[2]
     258
     259        Check[1 * n + 0] = false; // partial w.r.t x[1], x[0]
     260        Check[1 * n + 1] = false; // partial w.r.t x[1], x[1]
     261        Check[1 * n + 2] = true;  // partial w.r.t x[1], x[2]
     262
     263        Check[2 * n + 0] = false; // partial w.r.t x[2], x[0]
     264        Check[2 * n + 1] = true;  // partial w.r.t x[2], x[1]
     265        Check[2 * n + 2] = false; // partial w.r.t x[2], x[2]
     266
     267        // create function object F : X -> Y
     268        ADFun<double> F(X, Y);
     269
     270        // compute the forward Jacobian sparsity pattern for F
     271        CPPAD_TEST_VECTOR< bool > r(n * n);
     272        size_t i, j;
     273        for(i = 0; i < n; i++)
     274        {       for(j = 0; j < n; j++)
     275                        r[ i * n + j ] = false;
     276                r[ i * n + i ] = true;
     277        }
     278        F.ForSparseJac(n, r);
     279
     280        // compute the reverse Hessian sparsity pattern for F
     281        CPPAD_TEST_VECTOR< bool > s(m), h(n * n);
     282        s[0] = 1.;
     283        h = F.RevSparseHes(n, s);
     284
     285        // check values
     286        for(i = 0; i < n; i++)
     287        {       for(j = 0; j < n; j++)
     288                        ok &= (h[i * n + j] == Check[i * n + j]);
     289        }       
     290
     291        return ok;
     292}
     293
    218294} // End of empty namespace
    219295
     
    223299        ok &= case_two();
    224300        ok &= case_three();
    225 
    226         return ok;
    227 }
     301        ok &= case_four();
     302
     303        return ok;
     304}
  • trunk/test_more/rev_sparse_jac.cpp

    r1471 r1521  
    6868        index++;
    6969
    70 
    71 bool rev_sparse_jac(void)
     70namespace { // BEGIN empty namespace
     71
     72bool case_one(void)
    7273{       bool ok = true;
    7374        using namespace CppAD;
     
    169170        return ok;
    170171}
     172
     173bool case_two(void)
     174{       bool ok = true;
     175        using namespace CppAD;
     176
     177        // dimension of the domain space
     178        size_t n = 3;
     179
     180        // dimension of the range space
     181        size_t m = 3;
     182
     183        // inialize the vector as zero
     184        CppAD::VecAD<double> Z(n - 1);
     185        size_t k;
     186        for(k = 0; k < n-1; k++)
     187                Z[k] = 0.;
     188
     189        // independent variable vector
     190        CPPAD_TEST_VECTOR< AD<double> > X(n);
     191        X[0] = 0.;
     192        X[1] = 1.;
     193        X[2] = 2.;
     194        Independent(X);
     195
     196        // VecAD vector is going to depend on X[1] and X[2]
     197        Z[ X[0] ] = X[1];
     198        Z[ X[1] ] = X[2];
     199
     200        // dependent variable vector
     201        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     202
     203        // check results vector
     204        CPPAD_TEST_VECTOR< bool >       Check(m * n);
     205
     206        // initialize index into Y
     207        size_t index = 0;
     208
     209        // First component only depends on X[0];
     210        Y[index]             = X[0];
     211        Check[index * n + 0] = true;
     212        Check[index * n + 1] = false;
     213        Check[index * n + 2] = false;
     214        index++;
     215
     216        // Second component depends on the vector Z
     217        AD<double> zero(0);
     218        Y[index]             = Z[zero]; // Load by a parameter
     219        Check[index * n + 0] = false;
     220        Check[index * n + 1] = true;
     221        Check[index * n + 2] = true;
     222        index++;
     223
     224        // Third component depends on the vector Z
     225        Y[index]             = Z[ X[0] ]; // Load by a variable
     226        Check[index * n + 0] = false;
     227        Check[index * n + 1] = true;
     228        Check[index * n + 2] = true;
     229        index++;
     230
     231        // check final index
     232        assert( index == m );
     233
     234        // create function object F : X -> Y
     235        ADFun<double> F(X, Y);
     236
     237        // dependency matrix for the identity function S(y) = y
     238        CPPAD_TEST_VECTOR< bool > Py(m * m);
     239        size_t i, j;
     240        for(i = 0; i < m; i++)
     241        {       for(j = 0; j < m; j++)
     242                        Py[ i * m + j ] = false;
     243                Py[ i * m + i ] = true;
     244        }
     245
     246        // evaluate the dependency matrix for S [ F(x) ]
     247        CPPAD_TEST_VECTOR< bool > Px(m * n);
     248        Px = F.RevSparseJac(m, Py);
     249
     250        // check values
     251        for(i = 0; i < m; i++)
     252        {       for(j = 0; j < n; j++)
     253                        ok &= (Px[i * n + j] == Check[i * n + j]);
     254        }       
     255
     256        return ok;
     257}
     258
     259} // END empty namespace
     260
     261bool rev_sparse_jac(void)
     262{       bool ok = true;
     263
     264        ok &= case_one();
     265        ok &= case_two();
     266
     267        return ok;
     268}
Note: See TracChangeset for help on using the changeset viewer.