Changeset 3721


Ignore:
Timestamp:
Sep 3, 2015 11:15:40 PM (4 years ago)
Author:
bradbell
Message:

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: a2ae3a72194537866d9c632b4fecd8623246c7cc
end hash code: 5967e213f6fd2be5adf4adef224b3c86f6532496

commit 5967e213f6fd2be5adf4adef224b3c86f6532496
Author: Brad Bell <bradbell@…>
Date: Thu Sep 3 20:14:10 2015 -0700

Test for and fix bug in reverse hessian sparsity calculation.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/cppad/local/rev_hes_sweep.hpp

    r3717 r3721  
    880880                                if( user_bool )
    881881                                        bool_u[user_i * user_q + j] = true;
    882                                 else
     882                                if( user_set )
    883883                                        set_u[user_i].insert(j);
    884884                                j = rev_hes_sparse.next_element();
  • trunk/example/atomic/checkpoint.cpp

    r3717 r3721  
    3131defined by
    3232$latex \[
    33         F(x) = \left( \begin{array}{c}
    34                 x_0 \cdot x_0 \cdot x_0 \cdot x_0 \cdot x_0
    35                 \\
    36                 x_1 \cdot x_1 \cdot x_1 \cdot x_1 \cdot x_1
    37         \end{array} \right)
    38         \; , \;
    39         G(y) = \left( \begin{array}{c}
     33        F(y) = \left( \begin{array}{c}
    4034                y_0 + y_0 + y_0
    4135                \\
    4236                y_1 + y_1 + y_1
     37        \end{array} \right)
     38        \; , \;
     39        G(x) = \left( \begin{array}{c}
     40                \cdot x_0 \cdot x_0 \cdot x_0
     41                \\
     42                \cdot x_1 \cdot x_1 \cdot x_1
    4343        \end{array} \right)
    4444\] $$
     
    5757namespace {
    5858        using CppAD::AD;
    59         typedef CPPAD_TESTVECTOR(AD<double>) ADVector;
     59        typedef CPPAD_TESTVECTOR(AD<double>)            ADVector;
     60        typedef CppAD::atomic_base<double>::option_enum option_enum;
    6061
    61         void f_algo(const ADVector& x, ADVector& y)
    62         {       y[0] = 1.0;
    63                 y[1] = 1.0;
    64                 for(size_t k = 0; k < 5; k++)
    65                 {       y[0] *= x[0];
    66                         y[1] *= x[1];
    67                 }
    68                 return;
    69         }
    70         void g_algo(const ADVector& y, ADVector& z)
     62        void f_algo(const ADVector& y, ADVector& z)
    7163        {       z[0] = 0.0;
    7264                z[1] = 0.0;
     
    7769                return;
    7870        }
     71        void g_algo(const ADVector& x, ADVector& y)
     72        {       y[0] = 1.0;
     73                y[1] = 1.0;
     74                for(size_t k = 0; k < 3; k++)
     75                {       y[0] *= x[0];
     76                        y[1] *= x[1];
     77                }
     78                return;
     79        }
     80        bool test_case(option_enum f_sparsity, option_enum g_sparsity)
     81        {       bool ok = true;
     82                using CppAD::checkpoint;
     83                using CppAD::ADFun;
     84                using CppAD::NearEqual;
     85                size_t i, j, k, n = 2, m = n;
     86                double eps = 10. * std::numeric_limits<double>::epsilon();
     87
     88                // checkpoint version of the function F(x)
     89                ADVector ax(n), ay(n), az(m);
     90                for(j = 0; j < n; j++)
     91                        ax[j] = double(j + 1);
     92                // could also use bool_sparsity_enum or set_sparsity_enum
     93                checkpoint<double> atom_f("atom_f", f_algo, ax, ay, f_sparsity);
     94                checkpoint<double> atom_g("atom_g", g_algo, ay, az, g_sparsity);
     95
     96                // Record a version of z = g[f(x)] without checkpointing
     97                Independent(ax);
     98                f_algo(ax, ay);
     99                g_algo(ay, az);
     100                ADFun<double> check_not(ax, az);
     101
     102                // Record a version of z = g[f(x)] with checkpointing
     103                Independent(ax);
     104                atom_f(ax, ay);
     105                atom_g(ay, az);
     106                ADFun<double> check_yes(ax, az);
     107
     108                // checkpointing should use fewer operations
     109                ok &= check_yes.size_var() < check_not.size_var();
     110
     111                // this does not really save space becasue f and g are only used once
     112                ok &= check_not.size_var() <=
     113                        check_yes.size_var() + atom_f.size_var() + atom_g.size_var();
     114
     115                // compare forward mode results for orders 0, 1, 2
     116                size_t q = 2;
     117                CPPAD_TESTVECTOR(double) x_q(n*(q+1)), z_not(m*(q+1)), z_yes(m*(q+1));
     118                for(j = 0; j < n; j++)
     119                {       for(k = 0; k <= q; k++)
     120                                x_q[ j * (q+1) + k ] = 1.0 / (q + 1 - k);
     121                }
     122                z_not = check_not.Forward(q, x_q);
     123                z_yes = check_yes.Forward(q, x_q);
     124                for(i = 0; i < m; i++)
     125                {       for(k = 0; k <= q; k++)
     126                        {       double zik_not = z_not[ i * (q+1) + k];
     127                                double zik_yes = z_yes[ i * (q+1) + k];
     128                                ok &= NearEqual(zik_not, zik_yes, eps, eps);
     129                        }
     130                }
     131
     132                // compare reverse mode results
     133                CPPAD_TESTVECTOR(double) w(m*(q+1)), dw_not(n*(q+1)), dw_yes(n*(q+1));
     134                dw_not = check_not.Reverse(q+1, w);
     135                dw_yes = check_yes.Reverse(q+1, w);
     136                for(j = 0; j < n; j++)
     137                {       for(k = 0; k <= q; k++)
     138                        {       double dwjk_not = dw_not[ j * (q+1) + k];
     139                                double dwjk_yes = dw_yes[ j * (q+1) + k];
     140                                ok &= NearEqual(dwjk_not, dwjk_yes, eps, eps);
     141                        }
     142                }
     143
     144                // compare forward mode Jacobian sparsity patterns
     145                CppAD::vector< std::set<size_t> > r(n), s_not(m), s_yes(m);
     146                for(j = 0; j < n; j++)
     147                        r[j].insert(j);
     148                s_not = check_not.ForSparseJac(n, r);
     149                s_yes = check_yes.ForSparseJac(n, r);
     150                for(i = 0; i < m; i++)
     151                        ok &= s_not[i] == s_yes[i];
     152
     153                // compare reverse mode Jacobian sparsity patterns
     154                CppAD::vector< std::set<size_t> > s(m), r_not(m), r_yes(m);
     155                for(i = 0; i < m; i++)
     156                        s[i].insert(i);
     157                r_not = check_not.RevSparseJac(m, s);
     158                r_yes = check_yes.RevSparseJac(m, s);
     159                for(i = 0; i < m; i++)
     160                        ok &= s_not[i] == s_yes[i];
     161
     162
     163                // compare reverse mode Hessian sparsity patterns
     164                CppAD::vector< std::set<size_t> > s_one(1), h_not(n), h_yes(n);
     165                for(i = 0; i < m; i++)
     166                        s_one[0].insert(i);
     167                h_not = check_not.RevSparseHes(n, s_one);
     168                h_yes = check_yes.RevSparseHes(n, s_one);
     169                for(i = 0; i < n; i++)
     170                        ok &= h_not[i] == h_yes[i];
     171
     172                return ok;
     173        }
    79174}
    80175
    81176bool checkpoint(void)
    82177{       bool ok = true;
    83         using CppAD::checkpoint;
    84         using CppAD::ADFun;
    85         using CppAD::NearEqual;
    86         size_t i, j, k, n = 2, m = n;
    87         double eps = 10. * std::numeric_limits<double>::epsilon();
    88178
    89         // checkpoint version of the function F(x)
    90         ADVector ax(n), ay(n), az(m);
    91         for(j = 0; j < n; j++)
    92                 ax[j] = double(j);
    93         // could also use bool_sparsity_enum or set_sparsity_enum
    94         CppAD::atomic_base<double>::option_enum sparsity =
    95                 CppAD::atomic_base<double>::pack_sparsity_enum;
    96         checkpoint<double> atom_f("atom_f", f_algo, ax, ay, sparsity);
    97         checkpoint<double> atom_g("atom_g", g_algo, ay, az, sparsity);
     179        // different types of sparsity
     180        option_enum pack_sparsity = CppAD::atomic_base<double>::pack_sparsity_enum;
     181        option_enum bool_sparsity = CppAD::atomic_base<double>::bool_sparsity_enum;
     182        option_enum set_sparsity  = CppAD::atomic_base<double>::set_sparsity_enum;
    98183
    99         // Record a version of z = g[f(x)] without checkpointing
    100         Independent(ax);
    101         f_algo(ax, ay);
    102         g_algo(ay, az);
    103         ADFun<double> check_not(ax, az);
    104 
    105         // Record a version of z = g[f(x)] with checkpointing
    106         Independent(ax);
    107         atom_f(ax, ay);
    108         atom_g(ay, az);
    109         ADFun<double> check_yes(ax, az);
    110 
    111         // checkpointing should use fewer operations
    112         ok &= check_yes.size_var() < check_not.size_var();
    113 
    114         // this does not really save space becasue f and g are only used once
    115         ok &= check_not.size_var() <=
    116                 check_yes.size_var() + atom_f.size_var() + atom_g.size_var();
    117 
    118         // compare forward mode results for orders 0, 1, 2
    119         size_t q = 2;
    120         CPPAD_TESTVECTOR(double) x_q(n*(q+1)), z_not(m*(q+1)), z_yes(m*(q+1));
    121         for(j = 0; j < n; j++)
    122         {       for(k = 0; k <= q; k++)
    123                         x_q[ j * (q+1) + k ] = 1.0 / (q + 1 - k);
    124         }
    125         z_not = check_not.Forward(q, x_q);
    126         z_yes = check_yes.Forward(q, x_q);
    127         for(i = 0; i < m; i++)
    128         {       for(k = 0; k <= q; k++)
    129                 {       double zik_not = z_not[ i * (q+1) + k];
    130                         double zik_yes = z_yes[ i * (q+1) + k];
    131                         ok &= NearEqual(zik_not, zik_yes, eps, eps);
    132                 }
    133         }
    134 
    135         // compare reverse mode results
    136         CPPAD_TESTVECTOR(double) w(m*(q+1)), dw_not(n*(q+1)), dw_yes(n*(q+1));
    137         dw_not = check_not.Reverse(q+1, w);
    138         dw_yes = check_yes.Reverse(q+1, w);
    139         for(j = 0; j < n; j++)
    140         {       for(k = 0; k <= q; k++)
    141                 {       double dwjk_not = dw_not[ j * (q+1) + k];
    142                         double dwjk_yes = dw_yes[ j * (q+1) + k];
    143                         ok &= NearEqual(dwjk_not, dwjk_yes, eps, eps);
    144                 }
    145         }
    146 
    147         // mix sparsity so test both cases
    148         atom_f.option( CppAD::atomic_base<double>::bool_sparsity_enum );
    149         atom_g.option( CppAD::atomic_base<double>::set_sparsity_enum );
    150 
    151         // compare forward mode Jacobian sparsity patterns
    152         CppAD::vector< std::set<size_t> > r(n), s_not(m), s_yes(m);
    153         for(j = 0; j < n; j++)
    154                 r[j].insert(j);
    155         s_not = check_not.ForSparseJac(n, r);
    156         s_yes = check_yes.ForSparseJac(n, r);
    157         for(i = 0; i < m; i++)
    158                 ok &= s_not[i] == s_yes[i];
    159 
    160         // compare reverse mode Jacobian sparsity patterns
    161         CppAD::vector< std::set<size_t> > s(m), r_not(m), r_yes(m);
    162         for(i = 0; i < m; i++)
    163                 s[i].insert(i);
    164         r_not = check_not.RevSparseJac(m, s);
    165         r_yes = check_yes.RevSparseJac(m, s);
    166         for(i = 0; i < m; i++)
    167                 ok &= s_not[i] == s_yes[i];
    168 
    169 
    170         // compare reverse mode Hessian sparsity patterns
    171         CppAD::vector< std::set<size_t> > s_one(1), h_not(n), h_yes(n);
    172         for(i = 0; i < m; i++)
    173                 s_one[0].insert(i);
    174         h_not = check_not.RevSparseHes(n, s_one);
    175         h_yes = check_yes.RevSparseHes(n, s_one);
    176         for(i = 0; i < n; i++)
    177                 ok &= h_not[i] == h_yes[i];
     184        // test some different cases
     185        ok &= test_case(pack_sparsity, pack_sparsity);
     186        ok &= test_case(pack_sparsity, bool_sparsity);
     187        ok &= test_case(bool_sparsity, set_sparsity);
     188        ok &= test_case(set_sparsity,  set_sparsity);
    178189
    179190        return ok;
  • trunk/omh/whats_new/whats_new_15.omh

    r3720 r3721  
    6565
    6666$head 09-03$$
     67$list number$$
    6768There was a bug in the $cref/vectorBool/CppAD_vector/vectorBool/$$
    6869$cref/assignment/CppAD_vector/Assignment/$$.
     
    7071it not allow a size zero vector to be assigned using a vector any other size.
    7172This has been fixed.
     73$lnext
     74The addition of the
     75$cref/pack/atomic_option/atomic_sparsity/pack_sparsity_enum/$$ option
     76on 08-31 introduced a bug in the calculation of $cref RevSparseHes$$.
     77The $cref checkpoint.cpp$$ example was changed to demonstrate this
     78problem and the bug was fixed.
     79$lend
    7280
    7381$head 09-02$$
Note: See TracChangeset for help on using the changeset viewer.