source: trunk/test_more/optimize.cpp @ 3607

Last change on this file since 3607 was 3607, checked in by bradbell, 5 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: f4baab0ea7c2679ba9351b0439b6efebfa77ba04
end hash code: 894d554a00ceec8a3545f05eca62708a7b5cb43d

commit 894d554a00ceec8a3545f05eca62708a7b5cb43d
Author: Brad Bell <bradbell@…>
Date: Tue Jan 20 09:06:10 2015 -0700

Fix copyright end date.


whats_new_15.omh: Add comment about date of deprecation.

commit 611e982000168db91aba22b763c14bb78d57da47
Author: Brad Bell <bradbell@…>
Date: Tue Jan 20 08:53:00 2015 -0700

Squashed commit from old/compare_op to master:


In addition, fix copyright end date for some files, and add note about
deprecated date in whats_new_15.omh.


commit 6e46df5c850ecd58d7a886db4043bc3f2d4579d1
Author: Brad Bell <bradbell@…>
Date: Tue Jan 20 08:16:57 2015 -0700


Always return f.compare_change_op_index() == 0 after f.optimize().


checkpoint.hpp: ignore comparison operators.
fun_construct.hpp: remove double initilaization of values.
compare_change.cpp: demonstrate more features of new interface.
whats_new_15.omh: prepare for merging this branch into master.
wish_list.omh: update wish list item to new compare_change interface.


commit 45315907c70e5b383d984fb9498b54a474001af0
Author: Brad Bell <bradbell@…>
Date: Tue Jan 20 05:04:37 2015 -0700


Use the new f.compare_change_count(0) option in speed tests.


commit bb6e72befd6d01f1fb62c43b9b19471ffaa7cc2c
Author: Brad Bell <bradbell@…>
Date: Tue Jan 20 04:51:16 2015 -0700


Move CompareChange? -> deprecated and plan -> compare_change.


forward0sweep.hpp: skip comparison operator when count is zero.
forward1sweep.hpp: skip comparison operator when count is zero.
compare_change.cpp: demonstrate count == 0 case.


commit 622a13c568c612d9dfe9ccd1a01f4ac5f74ba824
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 23:17:42 2015 -0700


Add example and test of new compare change user API.


ad_fun.hpp: fix f.compare_change_op_index.
compare_change.cpp: Change example to use new API.
compare_change.cpp: Move old example here (just test now).


commit ec4c1613eae8df56fbf31e7b8711ce69cc41df83
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 21:12:11 2015 -0700


Implement the compare change plan, still needs an example and test.
Also have the change from 'plan' to just plain documentation.


commit a81a46f27011bee08ba072551044dc9f4a99a906
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 17:49:05 2015 -0700


Change name of compare_change functions and partially implement this new plan.


commit 146faad48a700a56362e74f9c3a3c39144a79738
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 14:22:40 2015 -0700


Branch: compare_op
plan.omh: change name of count function.


commit 35d91d126765d1a0ab4bfe9e2b006bbf535cd648
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 13:19:07 2015 -0700


Add deprecation date for some more cases.


commit 5bb65a8c48fae4263b66fcd04520e10e66febc11
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 13:13:51 2015 -0700


Add date of deprecation for some more cases.


commit e95ee6b209601cd9a075d2e37c602e73c32fb6ab
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 12:58:44 2015 -0700


Add date of deprecation for some more cases.


commit 0ea84ccd87383edc62a6ae1711da104b12e8c444
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 12:47:01 2015 -0700


Add date of deprecation for some cases.


commit 17755e609ea8e03472b08dcc2fb5ad347eb723cb
Author: Brad Bell <bradbell@…>
Date: Mon Jan 19 08:20:45 2015 -0700


plan.omh: changs some names.


commit 29f369c06d4d0ee284c4c668d52d8461613066dc
Author: Brad Bell <bradbell@…>
Date: Fri Jan 16 06:39:45 2015 -0700


Document the plan for compare_change user API.


compare_change.omh: fix minor typo.
plan.txt: change to the omhelp file plan.omh.


commit a3a2f4dedd202a722812b6eb2714851b40726e6e
Author: Brad Bell <bradbell@…>
Date: Thu Jan 15 21:03:44 2015 -0700


new_branch.sh: remove unused variable.
push_git2svn.py: move repository from github/bradbell to coin-or/bradbell.


commit 3751a197ab2897e76616f9d9b0915148bd855356
Author: Brad Bell <bradbell@…>
Date: Thu Jan 15 20:56:17 2015 -0700


plan.txt: plan for this branches API will go here.


commit 76013ec2ad7baacdeab5e761812d542867910174
Author: Brad Bell <bradbell@…>
Date: Thu Jan 15 18:04:33 2015 -0700


Store the operator index for the first comparision change in the ADFun object.


commit 9caf25014079a60df5de17bcac76775daf8ee201
Author: Brad Bell <bradbell@…>
Date: Thu Jan 15 12:45:56 2015 -0700


Make compare_change a parameter (so will be easy to add compare_first).


commit 2246d22fe82b8909d432f82ab0d783ce3351a02f
Author: Brad Bell <bradbell@…>
Date: Thu Jan 15 09:12:40 2015 -0700


speed_branch.sh: fix directory (before cd).


commit b3910de86558a97749741bfb728e45c5a86d1c73
Author: Brad Bell <bradbell@…>
Date: Thu Jan 15 05:14:01 2015 -0700


search.sh: use git to get list of source files.
ad_fun.hpp: imporve doxygen doc for compare_change_.
ad_tape.hpp: remove RecordCompare? (no longer used).
atomic_base.hpp: minor edit to user documentation.


commit dd74f331386cadc9cc272c264296e575691aa3f8
Author: Brad Bell <bradbell@…>
Date: Thu Jan 15 04:12:34 2015 -0700


Change Leq..Op -> Le..Op and leq.._op -> le.._op.


commit ae729296323eb7f4f4a7c0e90a303a8d7f4ed42a
Author: Brad Bell <bradbell@…>
Date: Wed Jan 14 21:19:55 2015 -0700


comp_op.hpp: Add doxygen documentaiton for compare operators.
compare.hpp: avoid an extra branch.


commit b064d59a5ad01dff5c708cc8c02f628f58c863ec
Author: Brad Bell <bradbell@…>
Date: Wed Jan 14 16:11:25 2015 -0700


Remove ComOp?, replaced by specific cases Eq..Op, Ne..Op, Leq..Op, Lt..Op,
which should run faster.


forward0sweep.hpp: Remove out-dated comment about CompOp?.
forward1sweep.hpp: Remove out-dated comment about CompOp?.


commit 5bb0a70d1151e9086b88024050cea6cf11e83aa7
Author: Brad Bell <bradbell@…>
Date: Wed Jan 14 09:02:17 2015 -0700


Use Eq..Op, Ne..Op to implement == and !=.


commit 0ebeec61bc040a00c50db41ca5da31fb87194f93
Author: Brad Bell <bradbell@…>
Date: Wed Jan 14 07:19:26 2015 -0700


compare.hpp: Finish folding < <= > >= into Lt..Op and Leq..Op.


commit c949e5b72158b98cbab61c8aea98e76008e9c2f4
Author: Brad Bell <bradbell@…>
Date: Wed Jan 14 06:28:41 2015 -0700


  1. Change plan to fold all compare operations into Leq..Op and Lt..Op cases.
  2. Fix alphabetical order between Ld and Leq.


commit 6ffee88b68b682359d62bc75a8c2ba3e28d012ac
Author: Brad Bell <bradbell@…>
Date: Tue Jan 13 22:40:18 2015 -0700


Splite parameter <= variable and parameter > variable as separate compare
operators.


commit 0841014db4ead690d1c2358f5e09494030ae1e5f
Author: Brad Bell <bradbell@…>
Date: Tue Jan 13 21:12:57 2015 -0700


Attempting to recording and playback of compare operators faster:
Split variable <= variable and variable > variable out as separate cases.


compare.hpp: remove white space at end of lines.


commit cfd3afceebd4b3383b3042cbca98caf82ff77670
Author: Brad Bell <bradbell@…>
Date: Tue Jan 13 08:39:12 2015 -0700


speed_branch.sh: compare speed_cppad between two git branches.
speed_new.sh: correct list of options, remove unused variable.
wish_list.omh: correct discussion of effect of keeping compare operators.


commit 9ed03e1ee2c258ca38561ad083067e235d032e14
Author: Brad Bell <bradbell@…>
Date: Tue Jan 13 05:54:30 2015 -0700


Change test so it corresponds to optimization keeping compare operators.


commit 5c418d477d58984b094bba027ebb0794e759e557
Author: Brad Bell <bradbell@…>
Date: Tue Jan 13 05:12:50 2015 -0700


Include example and test of using CompareChange? with optimized tape.


commit c1b48edfa56ca096ce8c2db1dbadb2658cb13fe3
Author: Brad Bell <bradbell@…>
Date: Tue Jan 13 04:24:54 2015 -0700


Fix optimization so variables used in compare opertions are alwasy available.


commit 6fe607ca30db07b356fd3a9fe9779fa2dfd382d8
Author: Brad Bell <bradbell@…>
Date: Mon Jan 12 07:56:41 2015 -0700


Keep CompareChange? in optimized versions of tape and when NDEBUG is defined.

commit b4c0e51489cc0878499a331b4af4875b2781d8d8
Author: Brad Bell <bradbell@…>
Date: Sun Jan 11 17:01:59 2015 -0700

new_branch.sh: final procedure (only hand tested).

  • Property svn:keywords set to Id
File size: 42.2 KB
Line 
1/* $Id: optimize.cpp 3607 2015-01-20 16:20:41Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse 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// 2DO: Test that optimize.hpp use of base_atomic<Base>::rev_sparse_jac works.
13
14# include <limits>
15# include <cppad/cppad.hpp>
16
17
18namespace {
19        // note this enum type is not part of the API (but its values are)
20        CppAD::atomic_base<double>::option_enum atomic_sparsity_option;
21        //
22        // ----------------------------------------------------------------
23        // Test nested conditional expressions.
24        int nested_cond_exp(void)
25        {       bool ok = true;
26                using CppAD::AD;
27                using CppAD::vector;
28       
29                // independent variable vector
30                vector< AD<double> > ax(2), ay(1);
31                ax[0] = 1.0;
32                ax[1] = 2.0;
33                Independent(ax);
34       
35                // first conditional expression
36                AD<double> ac1 = CondExpLe(ax[0], ax[1], 2.0 * ax[0], 3.0 * ax[1] );
37       
38                // second conditional expression
39                AD<double> ac2 = CondExpGe(ax[0], ax[1], 4.0 * ax[0], 5.0 * ax[1] );
40       
41                // third conditional expression
42                AD<double> ac3 = CondExpLt(ax[0], ax[1], 6.0 * ac1, 7.0 * ac2 );
43       
44                // create function object f : ax -> ay
45                ay[0] = ac3;
46                CppAD::ADFun<double> f(ax, ay);
47       
48                // now optimize the operation sequence
49                f.optimize();
50       
51                // now zero order forward
52                vector<double> x(2), y(1);
53                for(size_t i = 0; i < 3; i++)
54                {       x[0] = 1.0 - double(i);
55                        x[1] = - x[0];
56                        y    = f.Forward(0, x);
57                        //
58                        // first conditional expression
59                        double c1;
60                        if( x[0] <= x[1] )
61                                c1 = 2.0 * x[0];
62                        else    c1 = 3.0 * x[1];
63                        //
64                        // second conditional expression
65                        double c2;
66                        if( x[0] >= x[1] )
67                                c2 = 4.0 * x[0];
68                        else    c2 = 5.0 * x[1];
69       
70                        // third conditional expression
71                        double c3;
72                        if( x[0] < x[1] )
73                                c3 = 6.0 * c1;
74                        else    c3 = 7.0 * c2;
75       
76                        ok &= y[0] == c3;
77                }
78                return ok;
79        }
80        // ----------------------------------------------------------------
81        // Test for bug where checkpoint function did not depend on
82        // the operands in the logical comparison because of the CondExp
83        // sparsity pattern.
84        void j_algo( 
85                const CppAD::vector< CppAD::AD<double> >& ax ,
86                      CppAD::vector< CppAD::AD<double> >& ay )
87        {       ay[0] = CondExpGt(ax[0], ax[1], ax[2], ax[3]); }
88       
89        bool atomic_cond_exp_sparsity(void)
90        {       bool ok = true;
91                using CppAD::AD;
92                using CppAD::vector;
93       
94                // Create a checkpoint version of the function g
95                vector< AD<double> > au(4), av(1);
96                for(size_t i = 0; i < 4; i++)
97                        au[i] = AD<double>(i);
98                CppAD::checkpoint<double> j_check("j_check", j_algo, au, av);
99       
100                // independent variable vector
101                vector< AD<double> > ax(2), ay(1);
102                ax[0] = 1.;
103                ax[1] = 1.;
104                Independent(ax);
105       
106                // call atomic function that does not get used
107                for(size_t i = 0; i < 4; i++) 
108                        au[i] = ax[0] + AD<double>(i + 1) * ax[1];
109                j_check(au, ay);
110       
111                // create function object f : ax -> ay
112                CppAD::ADFun<double> f(ax, ay);
113
114       
115                // now optimize the operation sequence
116                j_check.option( atomic_sparsity_option );
117                f.optimize();
118       
119                // check result where true case is used; i.e., au[0] > au[1]
120                vector<double> x(2), y(1);
121                x[0] = 1.;
122                x[1] = -1;
123                y    = f.Forward(0, x);
124                ok  &= y[0] == x[0] + double(3) * x[1];
125               
126       
127                // check result where false case is used; i.e., au[0] <= au[1]
128                x[0] = 1.;
129                x[1] = 1;
130                y    = f.Forward(0, x);
131                ok  &= y[0] == x[0] + double(4) * x[1];
132               
133                return ok;
134        }
135        // -------------------------------------------------------------------
136        // Test conditional optimizing out call to an atomic function call
137        void k_algo(
138                const CppAD::vector< CppAD::AD<double> >& x ,
139                      CppAD::vector< CppAD::AD<double> >& y )
140        {       y[0] = x[0] + x[1]; }
141
142        void h_algo(
143                const CppAD::vector< CppAD::AD<double> >& x ,
144                      CppAD::vector< CppAD::AD<double> >& y )
145        {       y[0] = x[0] - x[1]; }
146
147        bool atomic_cond_exp(void)
148        {       bool ok = true;
149                typedef CppAD::vector< CppAD::AD<double> > ADVector;
150
151                // Create a checkpoint version of the function g
152                ADVector ax(2), ag(1), ah(1), ay(1);
153                ax[0] = 0.;
154                ax[1] = 1.;
155                CppAD::checkpoint<double> k_check("k_check", k_algo, ax, ag);
156                CppAD::checkpoint<double> h_check("h_check", h_algo, ax, ah);
157
158                // independent variable vector
159                Independent(ax);
160
161                // atomic function calls that get conditionally used
162                k_check(ax, ag);
163                h_check(ax, ah);
164
165                // conditional expression
166                ay[0] = CondExpLt(ax[0], ax[1], ag[0], ah[0]); 
167       
168                // create function object f : ax -> ay
169                CppAD::ADFun<double> f;
170                f.Dependent(ax, ay);
171       
172                // use zero order to evaluate when condition is true
173                CppAD::vector<double>  x(2), dx(2);
174                CppAD::vector<double>  y(1), dy(1), w(1);
175                x[0] = 3.;
176                x[1] = 4.;
177                y    = f.Forward(0, x);
178                ok  &= y[0] == x[0] + x[1];
179
180                // before optimize
181                k_check.option( atomic_sparsity_option );
182                h_check.option( atomic_sparsity_option );
183                ok  &= f.number_skip() == 0;
184
185                // now optimize the operation sequence
186                f.optimize();
187
188                // optimized zero order forward when condition is false
189                x[0] = 4.;
190                x[1] = 3.;
191                y    = f.Forward(0, x);
192                ok  &= y[0] == x[0] - x[1];
193
194                // after optimize can skip either call to g or call to h
195                ok  &= f.number_skip() == 1;
196
197                // optimized first order forward
198                dx[0] = 2.;
199                dx[1] = 1.;
200                dy    = f.Forward(1, dx);
201                ok   &= dy[0] == dx[0] - dx[1];
202
203                // optimized first order reverse
204                w[0]  = 1.;
205                dx    = f.Reverse(1, w);
206                ok   &= dx[0] == 1.;
207                ok   &= dx[1] == -1.;
208       
209                return ok;
210        }
211        // -------------------------------------------------------------------
212        // Test of optimizing out arguments to an atomic function
213        void g_algo( 
214                const CppAD::vector< CppAD::AD<double> >& ax ,
215                      CppAD::vector< CppAD::AD<double> >& ay )
216        {       ay = ax; }
217
218        bool atomic_no_used(void)
219        {       bool ok = true;
220                using CppAD::AD;
221                using CppAD::vector;
222       
223                // Create a checkpoint version of the function g
224                vector< AD<double> > ax(2), ay(2), az(1);
225                ax[0] = 0.;
226                ax[1] = 1.;
227                CppAD::checkpoint<double> g_check("g_check", g_algo, ax, ay);
228       
229                // independent variable vector
230                Independent(ax);
231       
232                // call atomic function that does not get used
233                g_check(ax, ay);
234       
235                // conditional expression
236                az[0] = CondExpLt(ax[0], ax[1], ax[0] + ax[1], ax[0] - ax[1]); 
237               
238                // create function object f : ax -> az
239                CppAD::ADFun<double> f(ax, az);
240
241                // number of variables before optimization
242                // (include ay[0] and ay[1])
243                size_t n_before = f.size_var();
244               
245                // now optimize the operation sequence
246                g_check.option( atomic_sparsity_option );
247                f.optimize();
248
249                // number of variables after optimization
250                // (does not include ay[0] and ay[1])
251                size_t n_after = f.size_var();
252                ok            &= n_after + 2 == n_before;
253       
254                // check optimization works ok
255                vector<double> x(2), z(1);
256                x[0] = 4.;
257                x[1] = 3.;
258                z    = f.Forward(0, x);
259                ok  &= z[0] == x[0] - x[1];
260               
261                return ok;
262        }
263        bool atomic_arguments(void)
264        {       bool ok = true;
265                using CppAD::AD;
266                using CppAD::vector;
267                vector< AD<double> > au(2), aw(2), ax(2), ay(1);
268
269                // create atomic function corresponding to g_algo
270                au[0] = 1.0;
271                au[1] = 2.0;
272                CppAD::checkpoint<double> g_check("g_algo", g_algo, au, ax);
273
274                // start recording a new function
275                CppAD::Independent(ax);
276
277                // now use g_check during the recording
278                au[0] = ax[0] + ax[1]; // this argument requires a new variable
279                au[1] = ax[0] - ax[1]; // this argument also requires a new variable
280                g_check(au, aw);
281
282                // now create f(x) = x_0 - x_1
283                ay[0] = aw[0];
284                CppAD::ADFun<double> f(ax, ay);
285
286                // number of variables before optimization
287                size_t n_before = f.size_var();
288 
289                // now optimize f so that the calculation of au[1] is removed
290                g_check.option( atomic_sparsity_option );
291                f.optimize();
292
293                // check difference in number of variables
294                size_t n_after = f.size_var();
295                ok &= n_before == n_after + 1;
296
297                // now compute and check a forward mode calculation
298                vector<double> x(2), y(1);
299                x[0] = 5.0;
300                x[1] = 6.0;
301                y    = f.Forward(0, x);
302                ok  &= (y[0] == x[0] + x[1]); 
303
304                return ok;
305        }
306
307        // -------------------------------------------------------------------
308        // Test the reverse dependency analysis optimization
309        template <class Vector>
310        void depend_fun
311        (const Vector& x, Vector& y, size_t& original, size_t& opt)
312        {       typedef typename Vector::value_type Scalar;
313                Scalar not_used;
314                Scalar one(1), two(2), three(3), four(4);
315
316                // independent variable and phantom at beginning
317                original = 1 + x.size();
318                opt      = 1 + x.size();
319
320                // unary operator where operand is arg[0]
321                // (note that sin corresponds to two tape variables)
322                not_used = CppAD::abs(x[0]); 
323                y[0]     = sin(x[0]);
324                original += 3;
325                opt      += 2;
326
327                // binary operator where left operand is a variable
328                // and right operand is a parameter
329                not_used = not_used + 2.;
330                y[1]     = x[1] * 3.;
331                original += 2;
332                opt      += 1;
333
334                // binary operator where left operand is a parameter
335                // and right operation is a variable
336                not_used = 2. - not_used;
337                y[2]     = 3. / x[2];
338                original += 2;
339                opt      += 1;
340
341                // binary operator where both operands are variables
342                not_used  = x[3] - not_used;
343                y[3]      = x[3] / x[2];
344                original += 2;
345                opt      += 1;
346
347                // conditional expression that will be optimized out
348                not_used = CppAD::CondExpLt(x[0], x[1], x[2], x[3]) + not_used;
349                y[4]     = CppAD::CondExpLt(x[4], one, two, three);
350                original += 3;
351                opt      += 1;
352
353                // y[5] does not depend on the value of not_used.
354                // Make sure a parameter, corresponding to a dependent variable,
355                // is not optimized out of the operation sequence.
356                y[5]      = 0.0 * not_used;
357                original += 1;
358                opt      += 1;
359
360                // Wwe do not use the argument x[5], to
361                // make sure it is not optimized out.
362
363                return;
364        }
365
366        bool depend_one(void)
367        {       // Test all except for VecAD operations
368                bool ok = true;
369                using CppAD::AD;
370                size_t original;
371                size_t opt;
372                size_t i, j;
373       
374                // domain space vector
375                size_t n  = 6;
376                CppAD::vector< AD<double> > X(n);
377                for(j = 0; j < n; j++)
378                        X[j] = 1. / double(j + 1); 
379       
380                // declare independent variables and start tape recording
381                CppAD::Independent(X);
382       
383                // range space vector
384                size_t m = n;
385                CppAD::vector< AD<double> > Y(m);
386                depend_fun(X, Y, original, opt);
387       
388                // create f: X -> Y and stop tape recording
389                CppAD::ADFun<double> F;
390                F.Dependent(X, Y); 
391       
392                CppAD::vector<double> x(n), y(m), check(m);
393                for(j = 0; j < n; j++)
394                        x[j] = Value(X[j]);
395                y = F.Forward(0, x);
396                depend_fun(x, check, original, opt);
397                for(i = 0; i < m; i++)
398                        ok &= (y[i] == check[i]);
399       
400                // Check size before optimization
401                ok &= F.size_var() == original;
402       
403                // Optimize the operation sequence
404                F.optimize();
405       
406                // Check size after optimization
407                ok &= F.size_var() == opt;
408       
409                // check result now
410                // (should have already been checked if NDEBUG not defined)
411                y = F.Forward(0, x);
412                for(i = 0; i < m; i++)
413                        ok &= (y[i] == check[i]);
414       
415                return ok;
416        }
417
418        bool depend_two(void)
419        {       // Test VecAD operations
420                bool ok = true;
421                using CppAD::AD;
422                size_t i, j;
423
424                // domain space vector
425                size_t n  = 2;
426                CppAD::vector< AD<double> > X(n);
427                for(j = 0; j < n; j++)
428                        X[j] = double(j); 
429
430                // range space vector
431                size_t m = 3;
432                CppAD::vector< AD<double> > Y(m);
433
434                CppAD::VecAD<double> U(m);
435                CppAD::VecAD<double> V(n);
436                for(i = 0; i < m; i++)
437                        U[i] = 0;
438                for(j = 0; j < n; j++)
439                        V[j] = 0;
440       
441                // declare independent variables and start tape recording
442                CppAD::Independent(X);
443
444                // first vecad vector that is a variable
445                U[ X[0] ] = X[1];
446
447                // second vecad vector that is a variable
448                V[ X[0] ] = X[1];
449
450                // Make dependency for vecad vectors different that for
451                // variables because original code used worng dependency info.
452                // Y does not depend on the first variable in the tape; i.e.
453                // the one corresponding to the BeginOp. So make it depend
454                // on the first vecad vector in the tape.
455                for(i = 0; i < m; i++)
456                {       AD<double> I(i);
457                        Y[i] = U[I];
458                }
459       
460                // create f: X -> Y and stop tape recording
461                // Y[ X[0] ] = X[1] and other components of Y are zero.
462                CppAD::ADFun<double> F;
463                F.Dependent(X, Y); 
464
465                // Check number of VecAD vectors plus number of VecAD elements
466                ok &= (F.size_VecAD() == 2 + n + m); 
467       
468                CppAD::vector<double> x(n), y(m);
469                for(j = 0; j < n; j++)
470                        x[j] = double(j);
471
472                y = F.Forward(0, x);
473                for(i = 0; i < m; i++)
474                {       if( i != static_cast<size_t>(x[0]) )
475                                ok &= (y[i] == 0.);
476                        else    ok &= (y[i] == x[1]);
477                }
478
479                F.optimize();
480
481                // Check number of VecAD vectors plus number of VecAD elements
482                ok &= (F.size_VecAD() == 1 + m); 
483                y = F.Forward(0, x);
484                for(i = 0; i < m; i++)
485                {       if( i != static_cast<size_t>(x[0]) )
486                                ok &= (y[i] == 0.);
487                        else    ok &= (y[i] == x[1]);
488                }
489               
490                return ok;
491        }
492        bool depend_three(void)
493        {       // Power function is a special case for optimize
494                bool ok = true;
495                using CppAD::AD;
496                using CppAD::vector;
497
498                size_t n = 3;
499                size_t j;
500       
501                vector< AD<double> >    X(n), Y(n);
502                vector<double>          x(n), y(n); 
503       
504                for(j = 0; j < n; j++)
505                        X[j] = x[j] = double(j+2);
506       
507                CppAD::Independent(X);
508                       
509                Y[0] = pow(X[0], 2.0);
510                Y[1] = pow(2.0, X[1]);
511                Y[2] = pow(X[0], X[1]);
512       
513                CppAD::ADFun<double> F(X, Y);
514                F.optimize();
515                y = F.Forward(0, x);
516
517                // Use identically equal because the result of the operations
518                // have been stored as double and gaurd bits have been dropped.
519                // (This may not be true for some compiler in the future).
520                for(j = 0; j < n; j++)
521                        ok &= ( y[j] == Value(Y[j]) );
522
523                // check reverse mode derivative
524                vector<double>   w(n), dw(n); 
525                w[0] = 0.;
526                w[1] = 0.;
527                w[2] = 1.;
528                dw = F.Reverse(1, w);
529
530                double eps = 20. * std::numeric_limits<double>::epsilon();
531                double check = x[1] * pow( x[0], x[1] - 1. );
532                ok &= CppAD::NearEqual( dw[0], check, eps, eps );
533
534                check = log( x[0] ) * pow( x[0], x[1] );
535                ok &= CppAD::NearEqual( dw[1], check, eps, eps );
536
537                check = 0.;
538                ok &= CppAD::NearEqual( dw[2], check, eps, eps );
539       
540                return ok;
541        }
542        // ===================================================================
543        // Test duplicate operation analysis
544
545        template <class Vector>
546        void duplicate_fun
547        (const Vector& x, Vector& y, size_t& original, size_t& opt)
548        {       typedef typename Vector::value_type Scalar;
549                original = 0;
550                opt      = 0;
551
552                // unary operator where operand is arg[0] and one result
553                Scalar a1 = CppAD::exp(x[0]); 
554                original += 1;
555                opt      += 1;
556
557                // unary operator where operand is arg[0] and two results
558                Scalar b1 = CppAD::sin(x[1]);
559                original += 2;
560                opt      += 2;
561
562                // non-commutative binary operator where left is a variable
563                // and right is a parameter
564                Scalar c1 = x[2] - 3.;
565                original += 1;
566                opt      += 1;
567
568                // non-commutative binary operator where left is a parameter
569                // and right is a variable
570                Scalar d1 = 3. / x[3];
571                original += 1;
572                opt      += 1;
573
574                // non-commutative binary operator where left is a variable
575                // and right is a variable
576                Scalar e1 = pow(x[3], x[4]);
577                original += 3;
578                opt      += 3;
579
580                // commutative binary operator where  left is a variable
581                // and right is a parameter
582                Scalar f1 = x[5] * 5.;
583                original += 1;
584                opt      += 1;
585
586                // commutative binary operator where  left is a variable
587                // and right is a variable
588                Scalar g1 = x[5] + x[6];
589                original += 1;
590                opt      += 1;
591
592                // duplicate variables
593                Scalar a2 = CppAD::exp(x[0]);
594                Scalar b2 = CppAD::sin(x[1]);  // counts for 2 variables
595                Scalar c2 = x[2] - 3.;
596                Scalar d2 = 3. / x[3];
597                Scalar e2 = pow(x[3], x[4]);   // counts for 3 variables
598                Scalar f2 = 5. * x[5];
599                Scalar g2 = x[6] + x[5];
600                original += 10;
601
602                // result vector
603                y[0] = a1 * a2;
604                y[1] = b1 * b2;
605                y[2] = c1 * c2;
606                y[3] = d1 * d2;
607                y[4] = e1 * e2;
608                y[5] = f1 * f2;
609                y[6] = g1 * g2;
610                original += 7;
611                opt      += 7;
612
613                return;
614        }
615        bool duplicate_one(void)
616        {
617                bool ok = true;
618                using CppAD::AD;
619                size_t original;
620                size_t opt;
621                size_t i, j;
622       
623                // domain space vector
624                size_t n  = 7;
625                CppAD::vector< AD<double> > X(n);
626                for(j = 0; j < n; j++)
627                        X[j] = 1. / double(j + 1); 
628       
629                // declare independent variables and start tape recording
630                CppAD::Independent(X);
631       
632                // range space vector
633                size_t m = n;
634                CppAD::vector< AD<double> > Y(m);
635                duplicate_fun(X, Y, original, opt);
636       
637                // create f: X -> Y and stop tape recording
638                CppAD::ADFun<double> F;
639                F.Dependent(X, Y); 
640       
641                CppAD::vector<double> x(n), y(m), check(m);
642                for(j = 0; j < n; j++)
643                        x[j] = Value(X[j]);
644                y = F.Forward(0, x);
645                duplicate_fun(x, check, original, opt);
646                for(i = 0; i < m; i++)
647                        ok &= (y[i] == check[i]);
648       
649                // Check size before optimization
650                ok &= F.size_var() == (n + 1 + original);
651       
652                // Optimize the operation sequence
653                F.optimize();
654       
655                // Check size after optimization
656                ok &= F.size_var() == (n + 1 + opt);
657       
658                // check result now
659                // (should have already been checked if NDEBUG not defined)
660                y = F.Forward(0, x);
661                for(i = 0; i < m; i++)
662                        ok &= (y[i] == check[i]);
663       
664                return ok;
665        }
666        // -------------------------------------------------------------------
667        bool duplicate_two(void)
668        {       // test that duplicate expression removal is relative to
669                // new and not just old argument indices.
670                bool ok = true;
671                using CppAD::AD;
672                size_t i, j;
673
674                // domain space vector
675                size_t n  = 1;
676                CppAD::vector< AD<double> > X(n);
677                for(j = 0; j < n; j++)
678                        X[j] = double(j + 2); 
679
680                // range space vector
681                size_t m = 1;
682                CppAD::vector< AD<double> > Y(m);
683
684                // declare independent variables and start tape recording
685                CppAD::Independent(X);
686
687                // create a new variable
688                AD<double> A1 = X[0] - 2.;
689
690                // create a duplicate variable
691                AD<double> A2 = X[0] - 2.;
692
693                // create a new variable using first version of duplicate
694                AD<double> B1 = A1 / 2.;
695
696                // create a duplicate that can only be dectected using new
697                // argument indices
698                AD<double> B2 = A2 / 2.; 
699
700                // Make a new variable for result
701                // and make it depend on all the variables
702                Y[0] = B1 + B2;
703
704                // create f: X -> Y and stop tape recording
705                CppAD::ADFun<double> F;
706                F.Dependent(X, Y); 
707
708                // check number of variables in original function
709                ok &= (F.size_var() ==  1 + n + m + 4 ); 
710       
711                CppAD::vector<double> x(n), y(m);
712                for(j = 0; j < n; j++)
713                        x[j] = double(j + 2);
714
715                y   = F.Forward(0, x);
716                for(i = 0; i < m; i++)
717                        ok &= ( y[i] == Value( Y[i] ) );
718
719                F.optimize();
720
721                // check number of variables  in optimized version
722                ok &= (F.size_var() == 1 + n + m + 2 ); 
723
724                y   = F.Forward(0, x);
725                for(i = 0; i < m; i++)
726                        ok &= ( y[i] == Value( Y[i] ) );
727
728                return ok;
729        }
730        // -------------------------------------------------------------------
731        bool duplicate_three(void)
732        {       // test that duplicate expression removal is relative to
733                // new and not just old argument indices (commutative case).
734                bool ok = true;
735                using CppAD::AD;
736                size_t i, j;
737
738                // domain space vector
739                size_t n  = 1;
740                CppAD::vector< AD<double> > X(n);
741                for(j = 0; j < n; j++)
742                        X[j] = double(j + 2); 
743
744                // range space vector
745                size_t m = 1;
746                CppAD::vector< AD<double> > Y(m);
747
748                // declare independent variables and start tape recording
749                CppAD::Independent(X);
750
751                // create a new variable
752                AD<double> A1 = X[0] + 2.;
753
754                // create a duplicate variable
755                AD<double> A2 = 2. + X[0];
756
757                // create a new variable using first version of duplicate
758                AD<double> B1 = A1 * 2.;
759
760                // create a duplicate that can only be dectected using new
761                // argument indices
762                AD<double> B2 = 2. * A2; 
763
764                // Make a new variable for result
765                // and make it depend on all the variables
766                Y[0] = B1 + B2;
767
768                // create f: X -> Y and stop tape recording
769                CppAD::ADFun<double> F;
770                F.Dependent(X, Y); 
771
772                // check number of variables in original function
773                ok &= (F.size_var() ==  1 + n + m + 4 ); 
774       
775                CppAD::vector<double> x(n), y(m);
776                for(j = 0; j < n; j++)
777                        x[j] = double(j + 2);
778
779                y   = F.Forward(0, x);
780                for(i = 0; i < m; i++)
781                        ok &= ( y[i] == Value( Y[i] ) );
782
783                F.optimize();
784
785                // check number of variables  in optimized version
786                ok &= (F.size_var() == 1 + n + m + 2 ); 
787
788                y   = F.Forward(0, x);
789                for(i = 0; i < m; i++)
790                        ok &= ( y[i] == Value( Y[i] ) );
791
792                return ok;
793        }
794        // -------------------------------------------------------------------
795        bool duplicate_four(void)
796        {       // Check that unary expression matching not only checks hash code,
797                // and operator, but also operand (old bug that has been fixed).
798                bool ok = true;
799                using CppAD::AD;
800                size_t j;
801
802                // domain space vector
803                size_t n  = 1;
804                CppAD::vector< AD<double> > X(n);
805                X[0] = 1.;
806
807                // range space vector
808                size_t m = 1;
809                CppAD::vector< AD<double> > Y(m);
810
811                // declare independent variables and start tape recording
812                CppAD::Independent(X);
813
814                // check a huge number of same operation with different operands
815                size_t n_operations = std::min(
816                        size_t(CPPAD_HASH_TABLE_SIZE) + 5,
817                        size_t(std::numeric_limits<CPPAD_TAPE_ADDR_TYPE>::max()) - 5
818                );
819                Y[0] = X[0];
820                for(j = 0; j < n_operations; j++)
821                        Y[0] = abs(Y[0]);
822
823                // create f: X -> Y and stop tape recording
824                CppAD::ADFun<double> F;
825                F.Dependent(X, Y); 
826
827                // check number of variables in original function
828                ok &= (F.size_var() ==  1 + n + n_operations ); 
829       
830                CppAD::vector<double> x(n), y(m);
831                x[0] = 1.;
832
833                y   = F.Forward(0, x);
834                ok &= ( y[0] == Value( Y[0] ) );
835
836                F.optimize();
837
838                // check same number of variables in optimized version
839                ok &= (F.size_var() == 1 + n + n_operations ); 
840
841                y   = F.Forward(0, x);
842                ok &= ( y[0] == Value( Y[0] ) );
843
844                return ok;
845        }
846        // ====================================================================
847        bool cummulative_sum(void)
848        {       // test conversion of a sequence of additions and subtraction
849                // to a cummulative summation sequence.
850                bool ok = true;
851                using CppAD::AD;
852                size_t i, j;
853
854                // domain space vector
855                size_t n  = 7;
856                CppAD::vector< AD<double> > X(n);
857                for(j = 0; j < n; j++)
858                        X[j] = double(j + 2); 
859
860                size_t n_original = 1 + n;
861                size_t n_optimize = 1 + n;
862
863                // range space vector
864                size_t m = 2;
865                CppAD::vector< AD<double> > Y(m);
866
867                // declare independent variables and start tape recording
868                CppAD::Independent(X);
869
870                // Operations inside of optimize_cadd
871                Y[0] = 5. + (X[0] + X[1]) + (X[1] - X[2]) // Addvv, Subvv
872                     + (X[2] - 1.) + (2. - X[3])   // Subvp, Subpv
873                     + (X[4] + 3.) + (4. + X[5]);  // Addpv, Addpv (no Addvp)
874                n_original += 12;
875                n_optimize += 1;
876
877
878                // Operations inside of optimize_csub
879                Y[1] = 5. - (X[1] + X[2]) - (X[2] - X[3]) // Addvv, Subvv
880                     - (X[3] - 1.) - (2. - X[4])   // Subvp, Subpv
881                     - (X[5] + 3.) - (4. + X[6]);  // Addpv, Addpv (no Addvp)
882                n_original += 12;
883                n_optimize += 1;
884
885                CppAD::ADFun<double> F;
886                F.Dependent(X, Y); 
887
888                // check number of variables in original function
889                ok &= (F.size_var() ==  n_original ); 
890       
891                CppAD::vector<double> x(n), y(m);
892                for(j = 0; j < n; j++)
893                        x[j] = double(j + 2);
894
895                y   = F.Forward(0, x);
896                for(i = 0; i < m; i++)
897                        ok &= ( y[i] == Value( Y[i] ) );
898
899                F.optimize();
900
901                // check number of variables  in optimized version
902                ok &= (F.size_var() == n_optimize ); 
903
904                y   = F.Forward(0, x);
905                for(i = 0; i < m; i++)
906                        ok &= ( y[i] == Value( Y[i] ) );
907
908                return ok;
909        }
910        // -------------------------------------------------------------------
911        bool forward_csum(void)
912        {       bool ok = true;
913       
914                using namespace CppAD;
915       
916                // independent variable vector
917                CppAD::vector< AD<double> > X(2);
918                X[0] = 0.; 
919                X[1] = 1.;
920                Independent(X);
921       
922                // compute sum of elements in X
923                CppAD::vector< AD<double> > Y(1);
924                Y[0] = X[0] + X[0] + X[1];
925       
926                // create function object F : X -> Y
927                ADFun<double> F(X, Y);
928       
929                // now optimize the operation sequence
930                F.optimize();
931       
932                // use zero order to evaluate F[ (3, 4) ]
933                CppAD::vector<double>  x0( F.Domain() );
934                CppAD::vector<double>  y0( F.Range() );
935                x0[0]    = 3.;
936                x0[1]    = 4.;
937                y0       = F.Forward(0, x0);
938                ok      &= NearEqual(y0[0] , x0[0]+x0[0]+x0[1], 1e-10, 1e-10);
939       
940                // evaluate derivative of F in X[0] direction
941                CppAD::vector<double> x1( F.Domain() );
942                CppAD::vector<double> y1( F.Range() );
943                x1[0]    = 1.;
944                x1[1]    = 0.;
945                y1       = F.Forward(1, x1);
946                ok      &= NearEqual(y1[0] , x1[0]+x1[0]+x1[1], 1e-10, 1e-10);
947       
948                // evaluate second derivative of F in X[0] direction
949                CppAD::vector<double> x2( F.Domain() );
950                CppAD::vector<double> y2( F.Range() );
951                x2[0]       = 0.;
952                x2[1]       = 0.;
953                y2          = F.Forward(2, x2);
954                double F_00 = 2. * y2[0];
955                ok         &= NearEqual(F_00, 0., 1e-10, 1e-10);
956       
957                return ok;
958        }
959        // -------------------------------------------------------------------
960        bool reverse_csum(void)
961        {       bool ok = true;
962       
963                using namespace CppAD;
964       
965                // independent variable vector
966                CppAD::vector< AD<double> > X(2);
967                X[0] = 0.; 
968                X[1] = 1.;
969                Independent(X);
970       
971                // compute sum of elements in X
972                CppAD::vector< AD<double> > Y(1);
973                Y[0] = X[0] - (X[0] - X[1]);
974       
975                // create function object F : X -> Y
976                ADFun<double> F(X, Y);
977       
978                // now optimize the operation sequence
979                F.optimize();
980       
981                // use zero order to evaluate F[ (3, 4) ]
982                CppAD::vector<double>  x0( F.Domain() );
983                CppAD::vector<double>  y0( F.Range() );
984                x0[0]    = 3.;
985                x0[1]    = 4.;
986                y0       = F.Forward(0, x0);
987                ok      &= NearEqual(y0[0] , x0[0]-x0[0]+x0[1], 1e-10, 1e-10);
988       
989                // evaluate derivative of F
990                CppAD::vector<double> dF( F.Domain() );
991                CppAD::vector<double> w( F.Range() );
992                w[0]    = 1.;
993                dF      = F.Reverse(1, w);
994                ok     &= NearEqual(dF[0] , 0., 1e-10, 1e-10);
995                ok     &= NearEqual(dF[1] , 1., 1e-10, 1e-10);
996       
997                return ok;
998        }
999        bool forward_sparse_jacobian()
1000        {       bool ok = true;
1001                using namespace CppAD;
1002       
1003                // dimension of the domain space
1004                size_t n = 3; 
1005       
1006                // dimension of the range space
1007                size_t m = 3;
1008       
1009                // independent variable vector
1010                CppAD::vector< AD<double> > X(n);
1011                X[0] = 2.; 
1012                X[1] = 3.;
1013                X[2] = 4.;
1014                Independent(X);
1015       
1016                // dependent variable vector
1017                CppAD::vector< AD<double> > Y(m);
1018       
1019                // check results vector
1020                CppAD::vector< bool >       Check(m * n);
1021       
1022                // initialize index into Y
1023                size_t index = 0;
1024       
1025                // Y[0]
1026                Y[index]             = X[0] + X[1] + 5.;
1027                Check[index * n + 0] = true;
1028                Check[index * n + 1] = true;
1029                Check[index * n + 2] = false;
1030                index++;
1031       
1032                // Y[1]
1033                Y[index]             = Y[0] - (X[1] + X[2]);
1034                Check[index * n + 0] = true;
1035                Check[index * n + 1] = true;
1036                Check[index * n + 2] = true;
1037                index++;
1038
1039                // Y[2]
1040                // 2DO: There is a subtitle issue that has to do with using reverse
1041                // jacobian sparsity patterns during the optimization process.
1042                // We need an option to include X[0] in the sparsity pattern
1043                // so the optimizer can know it affects the results.
1044                Y[index]             = CondExpLe(X[0], X[1], X[1]+X[1], X[2]-X[2]);
1045                Check[index * n + 0] = false;
1046                Check[index * n + 1] = true;
1047                Check[index * n + 2] = true;
1048                index++;
1049       
1050                // check final index
1051                assert( index == m );
1052       
1053                // create function object F : X -> Y
1054                ADFun<double> F(X, Y);
1055                F.optimize();
1056       
1057                // ---------------------------------------------------------
1058                // dependency matrix for the identity function
1059                CppAD::vector< std::set<size_t> > Sx(n);
1060                size_t i, j;
1061                for(i = 0; i < n; i++)
1062                {       assert( Sx[i].empty() );
1063                        Sx[i].insert(i);
1064                }
1065       
1066                // evaluate the dependency matrix for F(x)
1067                CppAD::vector< std::set<size_t> > Sy(m);
1068                Sy = F.ForSparseJac(n, Sx);
1069       
1070                // check values
1071                bool found;
1072                for(i = 0; i < m; i++)
1073                {       for(j = 0; j < n; j++)
1074                        {       found = Sy[i].find(j) != Sy[i].end();
1075                                ok &= (found == Check[i * n + j]);
1076                        }
1077                }       
1078       
1079                return ok;
1080        }
1081        bool reverse_sparse_jacobian()
1082        {       bool ok = true;
1083                using namespace CppAD;
1084       
1085                // dimension of the domain space
1086                size_t n = 3; 
1087       
1088                // dimension of the range space
1089                size_t m = 3;
1090       
1091                // independent variable vector
1092                CppAD::vector< AD<double> > X(n);
1093                X[0] = 2.; 
1094                X[1] = 3.;
1095                X[2] = 4.;
1096                Independent(X);
1097       
1098                // dependent variable vector
1099                CppAD::vector< AD<double> > Y(m);
1100       
1101                // check results vector
1102                CppAD::vector< bool >       Check(m * n);
1103       
1104                // initialize index into Y
1105                size_t index = 0;
1106       
1107                // Y[0]
1108                Y[index]             = X[0] + X[1] + 5.;
1109                Check[index * n + 0] = true;
1110                Check[index * n + 1] = true;
1111                Check[index * n + 2] = false;
1112                index++;
1113       
1114                // Y[1]
1115                Y[index]             = Y[0] - (X[1] + X[2]);
1116                Check[index * n + 0] = true;
1117                Check[index * n + 1] = true;
1118                Check[index * n + 2] = true;
1119                index++;
1120
1121                // Y[2]
1122                Y[index]             = CondExpLe(X[0], X[1], X[1]+X[1], X[2]-X[2]);
1123                Check[index * n + 0] = false;
1124                Check[index * n + 1] = true;
1125                Check[index * n + 2] = true;
1126                index++;
1127       
1128                // check final index
1129                assert( index == m );
1130       
1131                // create function object F : X -> Y
1132                ADFun<double> F(X, Y);
1133                F.optimize();
1134       
1135                // ----------------------------------------------------------
1136                // dependency matrix for the identity function
1137                CppAD::vector< bool > Py(m * m);
1138                size_t i, j;
1139                for(i = 0; i < m; i++)
1140                {       for(j = 0; j < m; j++)
1141                                Py[ i * m + j ] = (i == j);
1142                }
1143       
1144                // evaluate the dependency matrix for F(x)
1145                CppAD::vector< bool > Px(m * n);
1146                Px = F.RevSparseJac(m, Py);
1147       
1148                // check values
1149                for(i = 0; i < m; i++)
1150                {       for(j = 0; j < n; j++)
1151                                ok &= (Px[i * n + j] == Check[i * n + j]);
1152                }       
1153       
1154                return ok;
1155        }
1156        bool reverse_sparse_hessian(void)
1157        {       bool ok = true;
1158                using CppAD::AD;
1159                size_t i, j;
1160       
1161                size_t n = 3;
1162                CppAD::vector< AD<double> > X(n); 
1163                X[0] = 1.;
1164                X[1] = 2.;
1165                X[2] = 3.;
1166                CppAD::Independent(X);
1167       
1168                size_t m = 1;
1169                CppAD::vector< AD<double> > Y(m);
1170                Y[0] = CondExpGe( X[0], X[1], 
1171                        X[0] + (2. + X[1] + 3.) * X[1], 
1172                        X[0] + (2. + X[2] + 3.) * X[1]
1173                );
1174       
1175                CppAD::vector<bool> check(n * n);
1176                check[0 * n + 0] = false; // partial w.r.t. x[0], x[0]
1177                check[0 * n + 1] = false; //                x[0], x[1]
1178                check[0 * n + 2] = false; //                x[0], x[2]
1179
1180                check[1 * n + 0] = false; // partial w.r.t. x[1], x[0]
1181                check[1 * n + 1] = true;  //                x[1], x[1]
1182                check[1 * n + 2] = true;  //                x[1], x[2]
1183
1184                check[2 * n + 0] = false; // partial w.r.t. x[2], x[0]
1185                check[2 * n + 1] = true;  //                x[2], x[1]
1186                check[2 * n + 2] = false; //                x[2], x[2]
1187       
1188                // create function object F : X -> Y
1189                CppAD::ADFun<double> F(X, Y);
1190                F.optimize();
1191       
1192                // sparsity pattern for the identity function U(x) = x
1193                CppAD::vector<bool> Px(n * n);
1194                for(i = 0; i < n; i++)
1195                        for(j = 0; j < n; j++)
1196                                Px[ i * n + j ] = (i == j);
1197       
1198                // compute sparsity pattern for Jacobian of F(U(x))
1199                CppAD::vector<bool> P_jac(m * n);
1200                P_jac = F.ForSparseJac(n, Px);
1201       
1202                // compute sparsity pattern for Hessian of F_k ( U(x) )
1203                CppAD::vector<bool> Py(m);
1204                CppAD::vector<bool> Pxx(n * n);
1205                Py[0] = true;
1206                Pxx = F.RevSparseHes(n, Py);
1207                // check values
1208                for(i = 0; i < n * n; i++)
1209                        ok &= (Pxx[i] == check[i]);
1210       
1211                return ok;
1212        }
1213        // check that CondExp properly detects dependencies
1214        bool cond_exp_depend(void)
1215        {       bool ok = true;
1216                using CppAD::AD;
1217
1218                AD<double> zero(0.), one(1.), two(2.), three(3.);
1219       
1220                size_t n = 4;
1221                CppAD::vector< AD<double> > X(n); 
1222                X[0] = zero;
1223                X[1] = one;
1224                X[2] = two;
1225                X[3] = three;
1226                CppAD::Independent(X);
1227       
1228                size_t m = 4;
1229                CppAD::vector< AD<double> > Y(m);
1230                Y[0] = CondExpLt(X[0] + .5,  one,  two, three);
1231                Y[1] = CondExpLt(zero, X[1] + .5,  two, three);
1232                Y[2] = CondExpLt(zero,  one, X[2] + .5, three);
1233                Y[3] = CondExpLt(zero,  one,  two,  X[3] + .5);
1234
1235                CppAD::ADFun<double> f(X, Y);
1236                f.optimize();
1237
1238                CppAD::vector<double> x(n), y(m);
1239                size_t i;
1240                for(i = 0; i < n; i++)
1241                        x[i] = double(n - i);
1242                y    = f.Forward(0, x);
1243
1244                if( x[0] + .5 < 1. )
1245                        ok  &= y[0] == 2.;
1246                else    ok  &= y[0] == 3.;
1247                if( 0. < x[1] + .5 )
1248                        ok  &= y[1] == 2.;
1249                else    ok  &= y[1] == 3.;
1250                ok  &= y[2] == x[2] + .5;;
1251                ok  &= y[3] == 2.;
1252
1253                return ok;
1254        }
1255        // -------------------------------------------------------------------
1256        void my_union(
1257                std::set<size_t>&         result  ,
1258                const std::set<size_t>&   left    ,
1259                const std::set<size_t>&   right   )
1260        {       std::set<size_t> temp;
1261                std::set_union(
1262                        left.begin()              ,
1263                        left.end()                ,
1264                        right.begin()             ,
1265                        right.end()               ,
1266                        std::inserter(temp, temp.begin())
1267                );
1268                result.swap(temp);
1269        }
1270
1271        bool old_atomic_forward(
1272                size_t                         id ,
1273                size_t                          k , 
1274                size_t                          n ,
1275                size_t                          m ,
1276                const CppAD::vector<bool>&     vx ,
1277                CppAD::vector<bool>&           vy ,
1278                const CppAD::vector<double>&   tx , 
1279                CppAD::vector<double>&         ty )
1280        {       assert(n == 3 && m == 2);
1281                if( k > 0 ) 
1282                        return false;
1283
1284                // y[0] = x[0] + x[1]
1285                ty[0] = tx[0] + tx[1];
1286
1287                // y[1] = x[1] + x[2]
1288                ty[1] = tx[1] + tx[2];
1289               
1290                if( vy.size() > 0 )
1291                {       vy[0] = (vx[0] | vx[1]);
1292                        vy[1] = (vx[1] | vx[2]);
1293                }
1294                return true; 
1295        }
1296
1297        bool old_atomic_reverse(
1298                size_t                         id ,
1299                size_t                          k , 
1300                size_t                          n , 
1301                size_t                          m , 
1302                const CppAD::vector<double>&   tx , 
1303                const CppAD::vector<double>&   ty ,
1304                CppAD::vector<double>&         px ,
1305                const CppAD::vector<double>&   py )
1306        {       return false; }
1307
1308        bool old_atomic_for_jac_sparse(
1309                size_t                                  id ,
1310                size_t                                   n ,
1311                size_t                                   m ,
1312                size_t                                   q ,
1313                const CppAD::vector< std::set<size_t> >& r ,
1314                CppAD::vector< std::set<size_t>  >&      s )
1315        {       return false; }
1316
1317        bool old_atomic_rev_jac_sparse(
1318                size_t                                  id ,
1319                size_t                                   n ,
1320                size_t                                   m ,
1321                size_t                                   q ,
1322                CppAD::vector< std::set<size_t> >&       r ,
1323                const CppAD::vector< std::set<size_t> >& s )
1324        {       assert(n == 3 && m == 2);
1325                r[0].clear();
1326                r[1].clear();
1327                r[2].clear();
1328                // y[0] = x[0] + x[1]
1329                my_union(r[0], r[0], s[0]);
1330                my_union(r[1], r[1], s[0]);
1331                // y[1] = x[1] + x[2]
1332                my_union(r[1], r[1], s[1]);
1333                my_union(r[2], r[2], s[1]);
1334
1335                return true; 
1336        }
1337
1338        bool old_atomic_rev_hes_sparse(
1339                size_t                                  id ,
1340                size_t                                   n ,
1341                size_t                                   m ,
1342                size_t                                   q ,
1343                const CppAD::vector< std::set<size_t> >& r ,
1344                const CppAD::vector<bool>&               s ,
1345                CppAD::vector<bool>&                     t ,
1346                const CppAD::vector< std::set<size_t> >& u ,
1347                CppAD::vector< std::set<size_t> >&       v )
1348        {       return false; }
1349
1350        CPPAD_USER_ATOMIC(
1351                my_old_atomic             ,
1352                CppAD::vector              ,
1353                double                     ,
1354                old_atomic_forward        ,
1355                old_atomic_reverse        ,
1356                old_atomic_for_jac_sparse ,
1357                old_atomic_rev_jac_sparse ,
1358                old_atomic_rev_hes_sparse
1359        )
1360
1361        bool old_atomic_test(void)
1362        {       bool ok = true;
1363
1364                using CppAD::AD;
1365                size_t j;
1366                size_t n = 3;
1367                size_t m = 2;
1368                CppAD::vector< AD<double> > ax(n), ay(m), az(m);
1369                for(j = 0; j < n; j++)
1370                        ax[j] = AD<double>(j + 1);
1371                CppAD::Independent(ax);
1372
1373                size_t id = 0;
1374                // first call should stay in the tape
1375                my_old_atomic(id++, ax, ay);
1376                // second call will not get used
1377                my_old_atomic(id++, ax, az);
1378                // create function
1379                CppAD::ADFun<double> g(ax, ay);
1380                // should have 1 + n + m + m varaibles
1381                ok &= g.size_var() == (1 + n + m + m);
1382                g.optimize();
1383                // should have 1 + n + m varaibles
1384                ok &= g.size_var() == (1 + n + m);
1385
1386                // now test that the optimized function gives same results
1387                CppAD::vector<double> x(n), y(m);
1388                for(j = 0; j < n; j++)
1389                        x[j] = (j + 1) * (j + 1);
1390                y = g.Forward(0, x);
1391                // y[0] = x[0] + x[1]
1392                ok &= (y[0] == x[0] + x[1]);
1393                // y[1] = x[1] + x[2]
1394                ok &= (y[0] == x[0] + x[1]);
1395
1396                return ok;
1397        }
1398        bool not_identically_equal(void)
1399        {       bool ok = true;
1400                using CppAD::AD;
1401
1402                // independent variable vector
1403                size_t n = 5;
1404                CppAD::vector< AD<double> > ax(n);
1405                size_t j;
1406                for(j = 0; j < n; j++)
1407                        ax[j] = 1. / 3.;
1408                CppAD::Independent(ax);
1409       
1410                // dependent variable vector
1411                size_t m = 1;
1412                CppAD::vector< AD<double> > ay(m);
1413                ay[0]       = 0.;
1414                for(j = 0; j < n; j++)
1415                {       if( j % 2 == 0 )
1416                                ay[0] += ax[j];
1417                        else    ay[0] -= ax[j];
1418                }
1419                CppAD::ADFun<double> f(ax, ay);
1420       
1421                // Used to fail assert in optimize that forward mode results
1422                // are identically equal
1423                f.optimize();
1424       
1425                return ok;
1426        }
1427        // -----------------------------------------------------------------------
1428        double floor(const double& x)
1429        {       return std::floor(x); }
1430        CPPAD_DISCRETE_FUNCTION(double, floor)
1431        bool discrete_function(void)
1432        {       bool ok = true;
1433                using CppAD::vector;
1434       
1435                vector< CppAD::AD<double> > ax(1), ay(1);
1436                ax[0] = 0.0; 
1437                CppAD::Independent(ax);
1438                ay[0] =  floor(ax[0]) + floor(ax[0]); 
1439                CppAD::ADFun<double> f(ax, ay);
1440       
1441                size_t size_before = f.size_var();
1442                f.optimize(); 
1443                size_t size_after = f.size_var();
1444                ok &= size_after + 1 == size_before;
1445       
1446                vector<double> x(1), y(1);
1447                x[0] = -2.2;
1448                y    = f.Forward(0, x);
1449                ok &= y[0] == -6.0;
1450       
1451                return ok;
1452        }
1453        // ----------------------------------------------------------------
1454        void i_algo( 
1455                const CppAD::vector< CppAD::AD<double> >& ax ,
1456                      CppAD::vector< CppAD::AD<double> >& ay )
1457        {       ay[0] = 1.0 / ax[0]; }
1458        //
1459        // Test bug where atomic functions were not properly conditionally skipped.
1460        bool cond_exp_skip_atomic(void)
1461        {       bool ok = true;
1462                using CppAD::AD;
1463                using CppAD::vector;
1464       
1465                // Create a checkpoint version of the function i_algo
1466                vector< AD<double> > au(1), av(1), aw(1);
1467                au[0] = 1.0;
1468                CppAD::checkpoint<double> i_check("i_check", i_algo, au, av);
1469       
1470                // independent variable vector
1471                vector< AD<double> > ax(2), ay(1);
1472                ax[0] = 1.0;
1473                ax[1] = 2.0;
1474                Independent(ax);
1475       
1476                // call atomic function that does not get used
1477                au[0] = ax[0];
1478                i_check(au, av);
1479                au[0] = ax[1];
1480                i_check(au, aw);
1481                AD<double> zero = 0.0;
1482                ay[0] = CondExpGt(av[0], zero, av[0], aw[0]);
1483
1484                // create function object f : ax -> ay
1485                CppAD::ADFun<double> f(ax, ay);
1486
1487                // run case that skips the second call to afun
1488                // (can use trace in forward0sweep.hpp to see this).
1489                vector<double> x(2), y_before(1), y_after(1);
1490                x[0]      = 1.0;
1491                x[1]      = 2.0;
1492                y_before  = f.Forward(0, x);
1493                f.optimize();
1494                y_after   = f.Forward(0, x);
1495       
1496                ok &= y_before[0] == y_after[0];
1497               
1498                return ok;
1499        }
1500        //
1501        // Test bug where conditional dependence did not pass through
1502        // atomic functions
1503        bool cond_exp_atomic_dependence(void)
1504        {       bool ok = true;
1505                using CppAD::AD;
1506                using CppAD::vector;
1507       
1508                // Create a checkpoint version of the function i_algo
1509                vector< AD<double> > au(1), av(1), aw(1);
1510                au[0] = 1.0;
1511                CppAD::checkpoint<double> i_check("i_check", i_algo, au, av);
1512       
1513                vector< AD<double> > ax(2), ay(1);
1514                AD<double> zero = 0.0; 
1515                ax[0] = 1.0;
1516                ax[1] = 1.0;
1517                Independent(ax);
1518                av[0] = ax[0] + ax[1];
1519                i_check(av, aw);
1520                ay[0] = CondExpGt(aw[0], zero, zero, aw[0]);
1521                CppAD::ADFun<double> f;
1522                f.Dependent(ax, ay);
1523
1524                // run case that skips the second call to afun
1525                // (but not for order zero)
1526                vector<double> x(2), y_before(1), y_after(1);
1527                vector<double> dx(2), dy_before(1), dy_after(1);
1528                x[0]      = 1.0;
1529                x[1]      = 1.0;
1530                y_before  = f.Forward(0, x);
1531                dx[0]     = 2.0;
1532                dx[1]     = 2.0;
1533                dy_before = f.Forward(1, dx);
1534                f.optimize();
1535                y_after   = f.Forward(0, x);
1536                dy_after  = f.Forward(1, dx);
1537
1538                ok &= y_before[0]  == y_after[0];
1539                ok &= dy_before[0] == dy_after[0];
1540
1541                return ok;
1542        }
1543        // -----------------------------------------------------------------------
1544        // Test reverse mode conditionalay skipping commands.
1545        template <class Type>
1546        Type my_max(const CppAD::vector<Type>& arg)
1547        {       Type res = arg[0];
1548                for(size_t j = 0;j < arg.size(); j++)
1549                res = CondExpGt(res, arg[j], res, arg[j]);
1550                return res;
1551        }
1552        bool cond_exp_reverse(void)
1553        {       bool ok = true;
1554                size_t n = 3;
1555                using CppAD::vector;
1556                using CppAD::AD;
1557       
1558                vector< AD<double> > ax(n), ay(1);
1559                for(size_t j = 0; j < n; j++)
1560                        ax[j] = 1.0;
1561                Independent(ax);
1562                ay[0] = my_max(ax) + my_max(ax);
1563                CppAD::ADFun<double> f(ax, ay);
1564       
1565                f.optimize();
1566       
1567                vector<double> x(n), w(1), dx(n);
1568                for(size_t j = 0;j < n; j++)
1569                        x[j] = double(j);
1570                f.Forward(0, x);
1571                w[0] = 1.0;
1572                dx = f.Reverse(1, w);
1573                for(size_t j = 0; j < n; j++)
1574                {       if( j == n-1 )
1575                                ok &= dx[j] == 2.0;
1576                        else
1577                                ok &= dx[j] == 0.0;
1578                }
1579                return ok;
1580        }
1581        // Test case where an expression depends on both the true
1582        // and false cases (bug fixed 2014-12-22)
1583        bool cond_exp_both_true_and_false(void)
1584        {       bool ok = true;
1585                using CppAD::vector;
1586                using CppAD::AD;
1587
1588                // f(x) = x[0] + x[0] if x[0] >= 3
1589                //      = x[0] + x[1] otherwise
1590                vector< AD<double> > ax(2), ay(3);
1591                ax[0] = 1.0;
1592                ax[1] = 2.0;
1593                Independent(ax);
1594                AD<double> three(3);
1595                AD<double> value = ax[0] + ax[1];
1596                // a simple value
1597                ay[0]  = CppAD::CondExpGe(ax[0], three, value, value);
1598                // a  binary exprpression
1599                ay[1]  = CppAD::CondExpGe(ax[0], three, ax[0]-ax[1], ax[0]-ax[1]);
1600                // a unary expression
1601                ay[2]  = CppAD::CondExpGe(ax[0], three, exp(ax[0]), exp(ax[0]) );
1602                CppAD::ADFun<double> f(ax, ay);
1603                f.optimize();
1604
1605                // check case where x[0] >= 3
1606                vector<double> x(2), y(3);
1607                x[0] = 4.0;
1608                x[1] = 2.0;
1609                y    = f.Forward(0, x);
1610                ok  &= y[0] == x[0] + x[1];
1611                ok  &= y[1] == x[0] - x[1];
1612                ok  &= y[2] == exp(x[0]);
1613
1614                // check case where x[0] < 3
1615                x[0] = 1.0;
1616                x[1] = 2.0;
1617                y    = f.Forward(0, x);
1618                ok  &= y[0] == x[0] + x[1];
1619                ok  &= y[1] == x[0] - x[1];
1620                ok  &= y[2] == exp(x[0]);
1621
1622                return ok;
1623        }
1624}
1625
1626bool optimize(void)
1627{       bool ok = true;
1628
1629        atomic_sparsity_option = CppAD::atomic_base<double>::bool_sparsity_enum;
1630        for(size_t i = 0; i < 2; i++)
1631        {       // check conditional expression sparsity pattern
1632                // (used to optimize calls to atomic functions).
1633                ok     &= atomic_cond_exp_sparsity();
1634                // check optimizing out entire atomic function
1635                ok     &= atomic_cond_exp();
1636                // check optimizing out atomic arguments
1637                ok     &= atomic_no_used();
1638                ok     &= atomic_arguments();
1639                atomic_sparsity_option = 
1640                        CppAD::atomic_base<double>::set_sparsity_enum;
1641        }
1642        // check nested conditional expressions
1643        ok     &= nested_cond_exp();
1644        // check reverse dependency analysis optimization
1645        ok     &= depend_one();
1646        ok     &= depend_two();
1647        ok     &= depend_three();
1648        // check removal of duplicate expressions
1649        ok     &= duplicate_one();
1650        ok     &= duplicate_two();
1651        ok     &= duplicate_three();
1652        ok     &= duplicate_four();
1653        // convert sequence of additions to cummulative summation
1654        ok     &= cummulative_sum();
1655        ok     &= forward_csum();
1656        ok     &= reverse_csum();
1657        // sparsity patterns
1658        ok     &= forward_sparse_jacobian();
1659        ok     &= reverse_sparse_jacobian();
1660        ok     &= reverse_sparse_hessian();
1661        // check that CondExp properly detects dependencies
1662        ok     &= cond_exp_depend();
1663        // check old_atomic functions
1664        ok     &= old_atomic_test();
1665        // case where results are not identically equal
1666        ok     &= not_identically_equal();
1667        // case where a discrete function is used
1668        ok     &= discrete_function();
1669        // check conditional skip of an atomic function
1670        ok     &= cond_exp_skip_atomic();
1671        // check conditional dependence through atomic function
1672        ok     &= cond_exp_atomic_dependence();
1673        // check reverse mode conditional skipping
1674        ok     &= cond_exp_reverse();
1675        // check case where an expresion needed by both true and false case
1676        ok     &=  cond_exp_both_true_and_false();
1677        //
1678        CppAD::user_atomic<double>::clear();
1679        return ok;
1680}
Note: See TracBrowser for help on using the repository browser.