source: trunk/test_more/div.cpp @ 2506

Last change on this file since 2506 was 2506, checked in by bradbell, 8 years ago

Change Licenses: CPL-1.0 -> EPL-1.0, GPL-2.0->GPL-3.0

  • Property svn:keywords set to Id
File size: 5.3 KB
Line 
1/* $Id: div.cpp 2506 2012-10-24 19:36:49Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 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
13/*
14Two old Div examples now used just for valiadation testing
15*/
16
17# include <cppad/cppad.hpp>
18
19namespace { // BEGIN empty namespace
20
21bool DivTestOne(void)
22{       bool ok = true;
23
24        using namespace CppAD;
25
26        // assign some parameters
27        AD<double> zero = 0.;
28        AD<double>  one = 1.;
29
30        // independent variable vector, indices, values, and declaration
31        CPPAD_TESTVECTOR(AD<double>) U(2);
32        size_t s = 0;
33        size_t t = 1;
34        U[s]     = 2.;
35        U[t]     = 3.;
36        Independent(U);
37
38        // dependent variable vector and indices
39        CPPAD_TESTVECTOR(AD<double>) Z(6);
40        size_t x = 0;
41        size_t y = 1;
42        size_t z = 2;
43        size_t u = 3;
44        size_t v = 4;
45        size_t w = 5;
46
47        // dependent variables
48        Z[x] = U[s]   / U[t];   // AD<double> / AD<double>
49        Z[y] = Z[x]   / 4.;     // AD<double> / double
50        Z[z] = 5. / Z[y];       //     double / AD<double>
51        Z[u] =  Z[z] / one;     // division by a parameter equal to one
52        Z[v] =  Z[z] / 1.;      // division by a double equal to one
53        Z[w] =  zero / Z[z];    // division into a parameter equal to zero
54
55        // check division into a zero valued parameter results in a parameter
56        // (must do this before creating f because it erases the tape)
57        ok &= Parameter(Z[w]);
58
59        // create f : U -> Z and vectors used for derivative calculations
60        ADFun<double> f(U, Z);
61        CPPAD_TESTVECTOR(double) q( f.Domain() );
62        CPPAD_TESTVECTOR(double) r( f.Range() );
63 
64        // check parameter flag
65        ok &= f.Parameter(w);
66
67        // check values
68        ok &= NearEqual( Z[x] , 2. / 3. ,           1e-10 , 1e-10);
69        ok &= NearEqual( Z[y] , 2. / ( 3. * 4. ) ,  1e-10 , 1e-10);
70        ok &= NearEqual( Z[z] , 5. * 3. * 4. / 2. , 1e-10 , 1e-10);
71        ok &= ( Z[w] == 0. );
72        ok &= ( Z[u] == Z[z] );
73
74        // forward computation of partials w.r.t. s
75        q[s] = 1.;
76        q[t] = 0.;
77        r    = f.Forward(1, q);
78        ok &= NearEqual(r[x], 1./U[t],                1e-10 , 1e-10); // dx/ds
79        ok &= NearEqual(r[y], 1./(U[t]*4.),           1e-10 , 1e-10); // dy/ds
80        ok &= NearEqual(r[z], -5.*U[t]*4./(U[s]*U[s]),1e-10 , 1e-10); // dz/ds
81        ok &= ( r[u] == r[z] );                                       // du/ds
82        ok &= ( r[v] == r[z] );                                       // dv/ds
83        ok &= ( r[w] == 0. );                                         // dw/ds
84
85        // forward computation in the direction (1, 1)
86        q[s] = 1.;
87        q[t] = 1.;
88        r    = f.Forward(1, q);
89        ok  &= NearEqual(r[x], 1./U[t] - U[s]/(U[t] * U[t]), 1e-10, 1e-10);
90
91        // second order reverse mode computation
92        CPPAD_TESTVECTOR(double) Q( f.Domain() * 2 );
93        r[x] = 1.;
94        r[y] = r[z] = r[u] = r[v] = r[w] = 0.;
95        Q    = f.Reverse(2, r);
96        ok  &= NearEqual(
97                Q[s * f.Domain() + 1], 
98                - 1. / (U[t] * U[t]), 
99                1e-10, 
100                1e-10
101        );
102
103        return ok;
104}
105
106bool DivTestTwo(void)
107{       bool ok = true;
108        using namespace CppAD;
109
110        // independent variable vector
111        double u0 = .5;
112        CPPAD_TESTVECTOR(AD<double>) U(1);
113        U[0]      = u0;
114        Independent(U);
115
116        AD<double> a = U[0] / 1.; // AD<double> / double
117        AD<double> b = a  / 2;    // AD<double> / int
118        AD<double> c = 3. / b;    // double     / AD<double>
119        AD<double> d = 4  / c;    // int        / AD<double>
120
121        // dependent variable vector
122        CPPAD_TESTVECTOR(AD<double>) Z(1);
123        Z[0] = U[0] * U[0] / d;   // AD<double> / AD<double>
124
125        // create f: U -> Z and vectors used for derivative calculations
126        ADFun<double> f(U, Z); 
127        CPPAD_TESTVECTOR(double) v(1);
128        CPPAD_TESTVECTOR(double) w(1);
129
130        // check value
131        ok &= NearEqual(Value(Z[0]) , u0*u0/(4/(3/(u0/2))),  1e-10 , 1e-10);
132
133        // forward computation of partials w.r.t. u
134        size_t j;
135        size_t p     = 5;
136        double jfac  = 1.;
137        v[0]         = 1.;
138        double value = 6. / 4.;
139        for(j = 1; j < p; j++)
140        {
141                jfac *= j;
142                w     = f.Forward(j, v);       
143                ok &= NearEqual(jfac*w[0], value, 1e-10 , 1e-10); // d^jz/du^j
144                v[0]  = 0.;
145                value = 0.;
146        }
147
148        // reverse computation of partials of Taylor coefficients
149        CPPAD_TESTVECTOR(double) r(p); 
150        w[0]  = 1.;
151        r     = f.Reverse(p, w);
152        jfac  = 1.;
153        value = 6. / 4.;
154        for(j = 0; j < p; j++)
155        {
156                ok &= NearEqual(jfac*r[j], value, 1e-10 , 1e-10); // d^jz/du^j
157                jfac *= (j + 1);
158                value = 0.;
159        }
160
161        return ok;
162}
163
164bool DivTestThree(void)
165{       bool ok = true;
166        using namespace CppAD;
167
168        // more testing of variable / variable case
169        double x0 = 2.;
170        double x1 = 3.;
171        size_t n  = 2;
172        CPPAD_TESTVECTOR(AD<double>) X(n);
173        X[0]      = x0;
174        X[1]      = x1;
175        Independent(X);
176        size_t m  = 1;
177        CPPAD_TESTVECTOR(AD<double>) Y(m);
178        Y[0]      = X[0] / X[1];
179        ADFun<double> f(X, Y);
180
181        CPPAD_TESTVECTOR(double) dx(n), dy(m);
182        double check;
183        dx[0] = 1.;
184        dx[1] = 1.;
185        dy    = f.Forward(1, dx);
186        check = 1. / x1 - x0 / (x1 * x1); 
187        ok   &= NearEqual(dy[0], check, 1e-10 , 1e-10);
188
189        CPPAD_TESTVECTOR(double) w(m), dw(n);
190        w[0]  = 1.;
191        dw    = f.Reverse(1, w);
192        check = 1. / x1;
193        ok   &= NearEqual(dw[0], check, 1e-10 , 1e-10);
194        check = - x0 / (x1 * x1);
195        ok   &= NearEqual(dw[1], check, 1e-10 , 1e-10);
196
197        return ok;
198}
199
200} // END empty namespace
201
202bool Div(void)
203{       bool ok = true;
204        ok &= DivTestOne();
205        ok &= DivTestTwo(); 
206        return ok;
207}
Note: See TracBrowser for help on using the repository browser.