Ignore:
Timestamp:
May 24, 2012 12:30:25 PM (8 years ago)
Author:
bradbell
Message:

Merge in changes from branches/sparse.

  1. Update version number to today.
  2. Add sparse return user API interfaces to SparseJacobian? and SparseHessian?.
  3. Include new interface and imporve sparse Hessian and Jacobian examples.
  4. More testing of sparse Hessian and Jacobians.
  5. Add is_element to sparse_set and sparse_pack.

.: merge information.
search.sh: search cppad/speed directory, drop makefile.in from results.
svn_merge.sh: parameters for this merge.
parallel_ad.hpp: initilize is_element statics.
sparse_evaluate.hpp: fix typo in documentation.
glossary.omh: convert to modern omhelp, improve AD Function entry.
link_sparse_hessian.cpp: fix typo in documentation.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/test_more/sparse_hessian.cpp

    r1558 r2401  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    1717# include <cppad/cppad.hpp>
    1818namespace { // ---------------------------------------------------------
     19
     20bool rc_tridiagonal(void)
     21{       bool ok = true;
     22        using CppAD::AD;
     23        using CppAD::NearEqual;
     24        size_t i, j, k, ell;
     25        double eps = 10. * CppAD::epsilon<double>();
     26
     27        size_t n = 12; // must be greater than or equal 3; see n_sweep.
     28        size_t m = n - 1;
     29        CPPAD_TEST_VECTOR< AD<double> > a_x(n), a_y(m);
     30        CPPAD_TEST_VECTOR<double> x(n), check(n * n), w(m);
     31        for(j = 0; j < n; j++)
     32                a_x[j] = x[j] = double(j+1);
     33
     34        // declare independent variables and start taping
     35        CppAD::Independent(a_x);
     36
     37        for(ell = 0; ell < n * n; ell++)
     38                check[ell] = 0.0;
     39
     40        for(i = 0; i < m; i++)
     41        {       AD<double> diff = a_x[i+1] - a_x[i];   
     42                a_y[i] = 0.5 * diff * diff / double(i+2);
     43                w[i] = double(i+1);
     44                ell         = i * n + i;
     45                check[ell] += w[i] / double(i+2);
     46                ell         = (i+1) * n + i+1;
     47                check[ell] += w[i] / double(i+2);
     48                ell         = (i+1) * n + i;
     49                check[ell] -= w[i] / double(i+2);
     50                ell         = i * n + i+1;
     51                check[ell] -= w[i] / double(i+2);
     52        }
     53
     54        // create f: x -> y
     55        CppAD::ADFun<double> f(a_x, a_y);
     56                 
     57        // determine the sparsity pattern p for Hessian of w^T f
     58        typedef CppAD::vector< std::set<size_t> > VectorSet;
     59        VectorSet p_r(n);
     60        for(j = 0; j < n; j++)
     61                p_r[j].insert(j);
     62        f.ForSparseJac(n, p_r);
     63        //
     64        VectorSet p_s(1);
     65        for(i = 0; i < m; i++)
     66                if( w[i] != 0 )
     67                        p_s[0].insert(i);
     68        VectorSet p_h = f.RevSparseHes(n, p_s);
     69
     70        // requres the upper triangle of the Hessian
     71        size_t K = 2 * n - 1;
     72        CPPAD_TEST_VECTOR<size_t> r(K), c(K);
     73        CPPAD_TEST_VECTOR<double> hes(K);
     74        k = 0;
     75        for(i = 0; i < n; i++)
     76        {       r[k] = i;
     77                c[k] = i;
     78                k++;
     79                if( i < n-1 )
     80                {       r[k] = i;
     81                        c[k] = i+1;
     82                        k++;
     83                }
     84        }
     85        ok &= k == K;
     86
     87        // test computing sparse Hessian
     88        CppAD::sparse_hessian_work work;
     89        size_t n_sweep = f.SparseHessian(x, w, p_h, r, c, hes, work);
     90        ok &= n_sweep == 3;
     91        for(k = 0; k < K; k++)
     92        {       ell = r[k] * n + c[k];
     93                ok &=  NearEqual(check[ell], hes[k], eps, eps);
     94                ell = c[k] * n + r[k];
     95                ok &=  NearEqual(check[ell], hes[k], eps, eps);
     96        }
     97
     98        return ok;
     99}
     100
     101
    19102template <class VectorBase, class VectorBool>
    20103bool bool_case()
     
    63146
    64147        // determine the sparsity pattern p for Hessian of w^T F
    65         VectorBool r(n * n);
    66         for(j = 0; j < n; j++)
    67         {       for(k = 0; k < n; k++)
    68                         r[j * n + k] = false;
    69                 r[j * n + j] = true;
    70         }
    71         f.ForSparseJac(n, r);
    72         //
    73         VectorBool s(m);
    74         for(i = 0; i < m; i++)
    75                 s[i] = w[i] != 0;
    76         VectorBool p = f.RevSparseHes(n, s);
     148        VectorBool r(n * n);
     149        for(j = 0; j < n; j++)
     150        {       for(k = 0; k < n; k++)
     151                r[j * n + k] = false;
     152                r[j * n + j] = true;
     153        }
     154        f.ForSparseJac(n, r);
     155        //
     156        VectorBool s(m);
     157        for(i = 0; i < m; i++)
     158                s[i] = w[i] != 0;
     159        VectorBool p = f.RevSparseHes(n, s);
    77160
    78161        // test passing sparsity pattern
     
    129212
    130213        // determine the sparsity pattern p for Hessian of w^T F
    131         VectorSet r(n);
    132         for(j = 0; j < n; j++)
     214        VectorSet r(n);
     215        for(j = 0; j < n; j++)
    133216                r[j].insert(j);
    134         f.ForSparseJac(n, r);
    135         //
    136         VectorSet s(1);
    137         for(i = 0; i < m; i++)
    138                 if( w[i] != 0 )
     217        f.ForSparseJac(n, r);
     218        //
     219        VectorSet s(1);
     220        for(i = 0; i < m; i++)
     221                if( w[i] != 0 )
    139222                        s[0].insert(i);
    140         VectorSet p = f.RevSparseHes(n, s);
     223        VectorSet p = f.RevSparseHes(n, s);
    141224
    142225        // test passing sparsity pattern
     
    152235bool sparse_hessian(void)
    153236{       bool ok = true;
     237
     238        ok &= rc_tridiagonal();
     239        // ---------------------------------------------------------------
    154240        // vector of bool cases
    155241        ok &= bool_case< CppAD::vector  <double>, CppAD::vectorBool   >();
Note: See TracChangeset for help on using the changeset viewer.