source: trunk/test_more/forward_order.cpp @ 3688

Last change on this file since 3688 was 3214, checked in by bradbell, 6 years ago

merge in mul_dir branch except for multiple direcitons specific part:

check_include_omh.sh: ignore junk files.
search.sh: include compare_c directory, simplify.
whats_new_14.omh: change dates and contect from mul_dir to trunk.
svn_merge.sh: add --accept thiers-full option.

  • Property svn:keywords set to Id
File size: 5.2 KB
Line 
1/* $Id: forward_order.cpp 3214 2014-03-18 20:50:38Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 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/*
14*/
15// BEGIN C++
16# include <cppad/cppad.hpp>
17
18namespace {
19
20        double my_discrete(const double& x)
21        {       return static_cast<int> ( x ); }
22        CPPAD_DISCRETE_FUNCTION(double, my_discrete)
23
24}
25bool forward_order(void)
26{       bool ok = true;
27        using CppAD::AD;
28        using CppAD::NearEqual;
29        size_t j, k;
30        double eps = 10. * CppAD::numeric_limits<double>::epsilon();
31
32        // domain space vector
33        size_t n = 23, m = n;
34        CPPAD_TESTVECTOR(AD<double>) X(n), Y(m);
35        for(j = 0; j < n; j++)
36                X[j] = 0.0;
37
38        // declare independent variables and starting recording
39        CppAD::Independent(X);
40
41        // identity function values
42        size_t i = 0;
43        size_t identity_begin = i;
44        Y[i] = cos( acos( X[i] ) );                   i++; // AcosOp,  CosOp
45        Y[i] = sin( asin( X[i] ) );                   i++; // AsinOp,  SinOp
46        Y[i] = tan( atan( X[i] ) );                   i++; // AtanOp,  TanOp
47        Y[i] = CondExpGt(X[i], X[i-1], X[i], X[i-2]); i++; // CExpOp
48        Y[i] = X[i-1] * X[i] / X[i-1];                i++; // DivvvOp, MulvvOp
49        Y[i] = X[i] * X[i] * 1.0 / X[i];              i++; // DivpvOp
50        Y[i] = 5.0 * X[i] / 5.0;                      i++; // DivvpOp, MulpvOp
51        Y[i] = exp( log( X[i] ) );                    i++; // ExpOp,   LogOp
52        Y[i] = pow( sqrt( X[i] ), 2.0);               i++; // PowvpOp, SqrtOp
53        Y[i] = log( pow( std::exp(1.), X[i] ) );      i++; // PowpvOp
54        Y[i] = log( pow( X[i], X[i] ) ) / log( X[i]); i++; // PowvvOp
55        Y[i] = -2. - ((X[i-1] - X[i]) - 2.) + X[i-1]; i++; // Sub*Op: pv, vv, vp
56        size_t identity_end = i;
57
58        // other functions
59        Y[i] = abs( X[i] );         i++;   // AbsOp
60        Y[i] = X[i-1] + X[i] + 2.0; i++;   // AddvvOp, AddvpOp
61        Y[i] = cosh( X[i] );        i++;   // CoshOp
62        Y[i] = my_discrete( X[i] ); i++;   // DisOp
63        Y[i] = 4.0;                 i++;   // ParOp
64        Y[i] = sign( X[i] );        i++;   // SignOp
65        Y[i] = sinh( X[i] );        i++;   // SinhOp
66        Y[i] = tanh(X[i]);          i++;   // TanhOp
67
68        // VecAD operations
69        CppAD::VecAD<double> V(n);
70        AD<double> index = 1.;
71        V[index] = 3.0;
72        Y[i]     = V[index];            i++;   // StppOp, LdpOp
73        V[index] = X[0];
74        Y[i]     = V[index];            i++;   // StpvOp, LdpOp
75        index    = double(n) * X[3];
76        V[index] = X[1];
77        Y[i]     = V[index];            i++;   // StvvOp, LdvOp
78
79        // create f: X -> Y and stop tape recording
80        assert( i == m );
81        CppAD::ADFun<double> f;
82        f.Dependent(X, Y);
83
84        // initially, no values stored in f
85        ok &= f.size_order() == 0;
86
87        // Set X_j (t) = x + t
88        size_t p = 2, p1 = p+1;
89        CPPAD_TESTVECTOR(double) x(n), x_p(n * p1), y_p(m * p1);
90        for(j = 0; j < n; j++)
91        {       x[j]            = double(j) / double(n); 
92                x_p[j * p1 + 0] = x[j]; // order 0
93                x_p[j * p1 + 1] = 1.;   // order 1
94                x_p[j * p1 + 2] = 0.;   // order 2
95        }
96        // compute orders 0, 1, 2
97        y_p  = f.Forward(p, x_p);
98
99        // identity functions
100        CPPAD_TESTVECTOR(double) y(p1);
101        i = 0;
102        for(j = identity_begin; j != identity_end; j++)
103        {       y[0] = x[j];
104                y[1] = 1.0;
105                y[2] = 0.0;
106                for(k = 0; k < p1; k++)
107                        ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
108                i++;
109        }
110
111        // y_i = abs( x_i )
112        y[0] = CppAD::abs( x[i] );
113        y[1] = CppAD::sign( x[i] );
114        y[2] = 0.0;
115        for(k = 0; k < p1; k++)
116                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
117
118        // y_i = x_[i-1] + x_i + 2
119        i++;
120        y[0] = x[i-1] + x[i] + 2.0;
121        y[1] = 2.0;
122        y[2] = 0.0;
123        for(k = 0; k < p1; k++)
124                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
125
126        // y_i = cosh( x_i )
127        i++;
128        y[0] = CppAD::cosh( x[i] );
129        y[1] = CppAD::sinh( x[i] );
130        y[2] = CppAD::cosh( x[i] ) / 2.0;
131        for(k = 0; k < p1; k++)
132                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
133
134        // y_i = my_discrete( x_i )
135        i++;
136        y[0] = my_discrete( x[i] );
137        y[1] = 0.0;
138        y[2] = 0.0;
139        for(k = 0; k < p1; k++)
140                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
141
142        // y_i = 4
143        i++;
144        y[0] = 4.0;
145        y[1] = 0.0;
146        y[2] = 0.0;
147        for(k = 0; k < p1; k++)
148                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
149
150        // y_i = sign( x_i )
151        i++;
152        y[0] = CppAD::sign( x[i] );
153        y[1] = 0.0;
154        y[2] = 0.0;
155        for(k = 0; k < p1; k++)
156                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
157
158        // y_i = sinh( x_i )
159        i++;
160        y[0] = CppAD::sinh( x[i] );
161        y[1] = CppAD::cosh( x[i] );
162        y[2] = CppAD::sinh( x[i] ) / 2.0;
163        for(k = 0; k < p1; k++)
164                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
165
166        // y_i = tanh( x_i )
167        i++;
168        y[0] = CppAD::tanh( x[i] );
169        y[1] = 1.0 - y[0] * y[0];
170        y[2] = - 2.0 * y[0] * y[1] / 2.0;
171        for(k = 0; k < p1; k++)
172                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
173
174        // y_i = 3.0;
175        i++;
176        y[0] = 3.0;
177        y[1] = 0.0;
178        y[2] = 0.0;
179        for(k = 0; k < p1; k++)
180                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
181
182        // y_i = x_0
183        i++;
184        y[0] = x[0];
185        y[1] = 1.0;
186        y[2] = 0.0;
187        for(k = 0; k < p1; k++)
188                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
189
190        // y_i = x_1
191        i++;
192        y[0] = x[1];
193        y[1] = 1.0;
194        y[2] = 0.0;
195        for(k = 0; k < p1; k++)
196                ok  &= NearEqual(y[k] , y_p[i * p1 + k], eps, eps);
197
198        return ok;
199}
200// END C++
Note: See TracBrowser for help on using the repository browser.