source: trunk/test_more/rev_sparse_hes.cpp @ 2354

Last change on this file since 2354 was 1588, checked in by bradbell, 11 years ago

trunk: merge in branches/multiple_result:1584-158

  • Property svn:keywords set to Id
File size: 11.1 KB
Line 
1/* $Id: rev_sparse_hes.cpp 1588 2009-11-24 20:26:48Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Common Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13
14# include <cppad/cppad.hpp>
15 
16namespace { // Begin empty namespace
17
18bool case_one()
19{       bool ok = true;
20        using namespace CppAD;
21
22        // dimension of the domain space
23        size_t n = 10; 
24
25        // dimension of the range space
26        size_t m = 2;
27
28        // temporary indices
29        size_t i, j;
30
31        // initialize check values to false
32        CPPAD_TEST_VECTOR<bool> Check(n * n);
33        for(j = 0; j < n * n; j++)
34                Check[j] = false;
35
36        // independent variable vector
37        CPPAD_TEST_VECTOR< AD<double> > X(n);
38        for(j = 0; j < n; j++)
39                X[j] = AD<double>(j);
40        Independent(X);
41
42        // accumulate sum here
43        AD<double> sum(0.);
44
45        // variable * variable
46        sum += X[2] * X[3];
47        Check[2 * n + 3] = Check[3 * n + 2] = true;
48
49        // variable / variable
50        sum += X[4] / X[5];
51        Check[4 * n + 5] = Check[5 * n + 4] = Check[5 * n + 5] = true;
52
53        // CondExpLt(variable, variable, variable, variable)
54        sum += CondExpLt(X[1], X[2], sin(X[6]), cos(X[7]) );
55        Check[6 * n + 6] = true;
56        Check[7 * n + 7] = true;
57       
58        // pow(variable, variable)
59        sum += pow(X[8], X[9]);
60        Check[8 * n + 8] = Check[8 * n + 9] = true;
61        Check[9 * n + 8] = Check[9 * n + 9] = true;
62
63        // dependent variable vector
64        CPPAD_TEST_VECTOR< AD<double> > Y(m);
65        Y[0] = sum;
66
67        // variable - variable
68        Y[1]  = X[0] - X[1];
69
70        // create function object F : X -> Y
71        ADFun<double> F(X, Y);
72
73        // ------------------------------------------------------------------
74        // sparsity pattern for the identity function U(x) = x
75        CPPAD_TEST_VECTOR<bool> Px(n * n);
76        for(i = 0; i < n; i++)
77        {       for(j = 0; j < n; j++)
78                        Px[ i * n + j ] = false;
79                Px[ i * n + i ] = true;
80        }
81
82        // compute sparsity pattern for Jacobian of F(U(x))
83        F.ForSparseJac(n, Px);
84
85        // compute sparsity pattern for Hessian of F_0 ( U(x) )
86        CPPAD_TEST_VECTOR<bool> Py(m);
87        Py[0] = true;
88        Py[1] = false;
89        CPPAD_TEST_VECTOR<bool> Pxx(n * n);
90        Pxx = F.RevSparseHes(n, Py);
91
92        // check values
93        for(j = 0; j < n * n; j++)
94                ok &= (Pxx[j] == Check[j]);
95
96        // compute sparsity pattern for Hessian of F_1 ( U(x) )
97        Py[0] = false;
98        Py[1] = true;
99        Pxx = F.RevSparseHes(n, Py);
100        for(j = 0; j < n * n; j++)
101                ok &= (! Pxx[j]);  // Hessian is identically zero
102
103        // ------------------------------------------------------------------
104        // sparsity pattern for the identity function U(x) = x
105        CPPAD_TEST_VECTOR< std::set<size_t> > Sx(n);
106        for(i = 0; i < n; i++)
107                Sx[i].insert(i);
108
109        // compute sparsity pattern for Jacobian of F(U(x))
110        F.ForSparseJac(n, Sx);
111
112        // compute sparsity pattern for Hessian of F_0 ( U(x) )
113        CPPAD_TEST_VECTOR< std::set<size_t> > Sy(1);
114        Sy[0].insert(0);
115        CPPAD_TEST_VECTOR< std::set<size_t> > Sxx(n);
116        Sxx = F.RevSparseHes(n, Sy);
117
118        // check values
119        for(i = 0; i < n; i++)
120        {       for(j = 0; j < n; j++)
121                {       bool found = Sxx[i].find(j) != Sxx[i].end();
122                        ok &= (found == Check[i * n + j]);
123                }
124        }
125
126        // compute sparsity pattern for Hessian of F_1 ( U(x) )
127        Sy[0].clear();
128        Sy[0].insert(1);
129        Sxx = F.RevSparseHes(n, Sy);
130        for(i = 0; i < n; i++)
131        {       for(j = 0; j < n; j++)
132                {       bool found = Sxx[i].find(j) != Sxx[i].end();
133                        ok &= ! found;
134                }
135        }
136
137        return ok;
138}
139
140bool case_two()
141{       bool ok = true;
142        using namespace CppAD;
143
144        // dimension of the domain space
145        size_t n = 4; 
146
147        // dimension of the range space
148        size_t m = 1;
149
150        // temporary indices
151        size_t i, j;
152
153        // initialize check values to false
154        CPPAD_TEST_VECTOR<bool> Check(n * n);
155        for(j = 0; j < n * n; j++)
156                Check[j] = false;
157
158        // independent variable vector
159        CPPAD_TEST_VECTOR< AD<double> > X(n);
160        for(j = 0; j < n; j++)
161                X[j] = AD<double>(j);
162        Independent(X);
163
164        // Test the case where dependent variable is a non-linear function
165        // of the result of a conditional expression.
166        CPPAD_TEST_VECTOR< AD<double> > Y(m);
167        Y[0] = CondExpLt(X[0], X[1], X[2], X[3]);
168        Y[0] = cos(Y[0]) + X[0] + X[1];
169
170        // Hessian with respect to x[0] and x[1] is zero.
171        // Hessian with respect to x[2] and x[3] is full
172        // (although we know that there are no cross terms, this is an
173        // inefficiency of the conditional expression operator).
174        Check[2 * n + 2] = Check[ 2 * n + 3 ] = true;
175        Check[3 * n + 2] = Check[ 3 * n + 3 ] = true;
176       
177        // create function object F : X -> Y
178        ADFun<double> F(X, Y);
179
180        // -----------------------------------------------------------------
181        // sparsity pattern for the identity function U(x) = x
182        CPPAD_TEST_VECTOR<bool> Px(n * n);
183        for(i = 0; i < n; i++)
184        {       for(j = 0; j < n; j++)
185                        Px[ i * n + j ] = false;
186                Px[ i * n + i ] = true;
187        }
188
189        // compute sparsity pattern for Jacobian of F(U(x))
190        F.ForSparseJac(n, Px);
191
192        // compute sparsity pattern for Hessian of F_0 ( U(x) )
193        CPPAD_TEST_VECTOR<bool> Py(m);
194        Py[0] = true;
195        CPPAD_TEST_VECTOR<bool> Pxx(n * n);
196        Pxx = F.RevSparseHes(n, Py);
197
198        // check values
199        for(j = 0; j < n * n; j++)
200                ok &= (Pxx[j] == Check[j]);
201
202        // ------------------------------------------------------------------
203        // sparsity pattern for the identity function U(x) = x
204        CPPAD_TEST_VECTOR< std::set<size_t> > Sx(n);
205        for(i = 0; i < n; i++)
206                Sx[i].insert(i);
207
208        // compute sparsity pattern for Jacobian of F(U(x))
209        F.ForSparseJac(n, Sx);
210
211        // compute sparsity pattern for Hessian of F_0 ( U(x) )
212        CPPAD_TEST_VECTOR< std::set<size_t> > Sy(1);
213        Sy[0].insert(0);
214        CPPAD_TEST_VECTOR< std::set<size_t> > Sxx(n);
215        Sxx = F.RevSparseHes(n, Sy);
216
217        // check values
218        for(i = 0; i < n; i++)
219        {       for(j = 0; j < n; j++)
220                {       bool found = Sxx[i].find(j) != Sxx[i].end();
221                        ok &= (found == Check[i * n + j]);
222                }
223        }
224
225        return ok;
226}
227
228bool case_three()
229{       bool ok = true;
230        using CppAD::AD;
231
232        // domain space vector
233        size_t n = 1; 
234        CPPAD_TEST_VECTOR< AD<double> > X(n);
235        X[0] = 0.; 
236
237        // declare independent variables and start recording
238        CppAD::Independent(X);
239
240        // range space vector
241        size_t m = 1;
242        CPPAD_TEST_VECTOR< AD<double> > Y(m);
243
244        // make sure reverse jacobian is propogating dependency to
245        // intermediate values (not just final ones).
246        Y[0] = X[0] * X[0] + 2;
247
248        // create f: X -> Y and stop tape recording
249        CppAD::ADFun<double> f(X, Y);
250
251        // ------------------------------------------------------------------
252        // sparsity pattern for the identity matrix
253        CppAD::vector<bool> r(n * n);
254        size_t i, j;
255        for(i = 0; i < n; i++)
256        {       for(j = 0; j < n; j++)
257                        r[ i * n + j ] = false;
258                r[ i * n + i ] = true;
259        }
260
261        // compute sparsity pattern for J(x) = F^{(1)} (x)
262        f.ForSparseJac(n, r);
263
264        // compute sparsity pattern for H(x) = F_0^{(2)} (x)
265        CppAD::vector<bool> s(m);
266        for(i = 0; i < m; i++)
267                s[i] = false;
268        s[0] = true;
269        CppAD::vector<bool> h(n * n);
270        h    = f.RevSparseHes(n, s);
271
272        // check values
273        ok  &= (h[ 0 * n + 0 ] == true);  // second partial w.r.t x[0], x[0]
274
275        // ------------------------------------------------------------------
276        // sparsity pattern for the identity function U(x) = x
277        CPPAD_TEST_VECTOR< std::set<size_t> > Sx(n);
278        for(i = 0; i < n; i++)
279                Sx[i].insert(i);
280
281        // compute sparsity pattern for Jacobian of F(U(x))
282        f.ForSparseJac(n, Sx);
283
284        // compute sparsity pattern for Hessian of F_0 ( U(x) )
285        CPPAD_TEST_VECTOR< std::set<size_t> > Sy(1);
286        Sy[0].insert(0);
287        CPPAD_TEST_VECTOR< std::set<size_t> > Sxx(n);
288        Sxx = f.RevSparseHes(n, Sy);
289
290        // check value
291        bool found = Sxx[0].find(0) != Sxx[0].end();
292        ok &= (found == true);
293
294        return ok;
295}
296
297bool case_four()
298{       bool ok = true;
299        using namespace CppAD;
300
301        // dimension of the domain space
302        size_t n = 3; 
303
304        // dimension of the range space
305        size_t m = 1;
306
307        // inialize the vector as zero
308        CppAD::VecAD<double> Z(n - 1);
309        size_t k;
310        for(k = 0; k < n-1; k++)
311                Z[k] = 0.;
312
313        // independent variable vector
314        CPPAD_TEST_VECTOR< AD<double> > X(n);
315        X[0] = 0.; 
316        X[1] = 1.;
317        X[2] = 2.;
318        Independent(X);
319
320        // VecAD vector z depends on both x[1] and x[2]
321        // (component indices do not matter because they can change).
322        Z[ X[0] ] = X[1] * X[2];
323        Z[ X[1] ] = 0.; 
324
325        // dependent variable vector
326        CPPAD_TEST_VECTOR< AD<double> > Y(m);
327
328        // check results vector
329        CPPAD_TEST_VECTOR< bool >       Check(n * n);
330
331        // y = z[j] where j might be zero or one.
332        Y[0]  =  Z[ X[1] ];
333
334        Check[0 * n + 0] = false; // partial w.r.t x[0], x[0]
335        Check[0 * n + 1] = false; // partial w.r.t x[0], x[1]
336        Check[0 * n + 2] = false; // partial w.r.t x[0], x[2]
337
338        Check[1 * n + 0] = false; // partial w.r.t x[1], x[0]
339        Check[1 * n + 1] = false; // partial w.r.t x[1], x[1]
340        Check[1 * n + 2] = true;  // partial w.r.t x[1], x[2]
341
342        Check[2 * n + 0] = false; // partial w.r.t x[2], x[0]
343        Check[2 * n + 1] = true;  // partial w.r.t x[2], x[1]
344        Check[2 * n + 2] = false; // partial w.r.t x[2], x[2]
345
346        // create function object F : X -> Y
347        ADFun<double> F(X, Y);
348
349        // -----------------------------------------------------
350        // compute the forward Jacobian sparsity pattern for F
351        CPPAD_TEST_VECTOR< bool > r(n * n);
352        size_t i, j;
353        for(i = 0; i < n; i++)
354        {       for(j = 0; j < n; j++)
355                        r[ i * n + j ] = false;
356                r[ i * n + i ] = true;
357        }
358        F.ForSparseJac(n, r);
359
360        // compute the reverse Hessian sparsity pattern for F
361        CPPAD_TEST_VECTOR< bool > s(m), h(n * n);
362        s[0] = 1.;
363        h = F.RevSparseHes(n, s);
364
365        // check values
366        for(i = 0; i < n; i++)
367        {       for(j = 0; j < n; j++)
368                        ok &= (h[i * n + j] == Check[i * n + j]);
369        }       
370
371        // ------------------------------------------------------------------
372        // sparsity pattern for the identity function U(x) = x
373        CPPAD_TEST_VECTOR< std::set<size_t> > Sx(n);
374        for(i = 0; i < n; i++)
375                Sx[i].insert(i);
376
377        // compute sparsity pattern for Jacobian of F(U(x))
378        F.ForSparseJac(n, Sx);
379
380        // compute sparsity pattern for Hessian of F_0 ( U(x) )
381        CPPAD_TEST_VECTOR< std::set<size_t> > Sy(1);
382        Sy[0].insert(0);
383        CPPAD_TEST_VECTOR< std::set<size_t> > Sxx(n);
384        Sxx = F.RevSparseHes(n, Sy);
385
386        // check values
387        for(i = 0; i < n; i++)
388        {       for(j = 0; j < n; j++)
389                {       bool found = Sxx[i].find(j) != Sxx[i].end();
390                        ok &= (found == Check[i * n + j]);
391                }
392        }
393
394        return ok;
395}
396
397bool case_five(void)
398{       bool ok = true;
399        using CppAD::AD;
400        size_t i, j, k;
401
402        size_t n = 2;
403        CPPAD_TEST_VECTOR< AD<double> > X(n); 
404        X[0] = 1.;
405        X[1] = 2.;
406        CppAD::Independent(X);
407
408        size_t m = 2;
409        CPPAD_TEST_VECTOR< AD<double> > Y(m);
410        Y[0] = pow(X[0], 2.);
411        Y[1] = pow(2., X[1]);
412
413        // create function object F : X -> Y
414        CppAD::ADFun<double> F(X, Y);
415
416        // sparsity pattern for the identity function U(x) = x
417        CPPAD_TEST_VECTOR<bool> Px(n * n);
418        for(i = 0; i < n; i++)
419                for(j = 0; j < n; j++)
420                        Px[ i * n + j ] = (i == j);
421
422        // compute sparsity pattern for Jacobian of F(U(x))
423        F.ForSparseJac(n, Px);
424
425        // compute sparsity pattern for Hessian of F_k ( U(x) )
426        CPPAD_TEST_VECTOR<bool> Py(m);
427        CPPAD_TEST_VECTOR<bool> Pxx(n * n);
428        for(k = 0; k < m; k++)
429        {       for(i = 0; i < m; i++)
430                        Py[i] = (i == k);
431                Pxx = F.RevSparseHes(n, Py);
432                // check values
433                for(i = 0; i < n; i++)
434                {       for(j = 0; j < n; j++)
435                        {       bool check = (i == k) & (j == k);
436                                ok        &= Pxx[i * n + j] == check;
437                        }
438                }
439        }
440        return ok;
441}
442
443} // End of empty namespace
444
445bool rev_sparse_hes(void)
446{       bool ok = true;
447
448        ok &= case_one();
449        ok &= case_two();
450        ok &= case_three();
451        ok &= case_four();
452        ok &= case_five();
453
454        return ok;
455}
Note: See TracBrowser for help on using the repository browser.