Codebase list cppad / ed313bd
master: rosen_34.cpp: use dynamic parameters and ADFun swap. fun_construct.hpp: fix a bug in swap. Brad Bell 3 years ago
3 changed file(s) with 107 addition(s) and 83 deletion(s). Raw diff Collapse all Expand all
+82
-72
example/general/rosen_34.cpp less more
4444
4545 namespace {
4646 class Fun {
47 private:
48 const bool use_x_;
49 CppAD::ADFun<double> ode_ind_;
50 CppAD::ADFun<double> ode_dep_;
4751 public:
4852 // constructor
49 Fun(bool use_x_) : use_x(use_x_)
53 Fun(bool use_x) : use_x_(use_x)
5054 { }
5155
5256 // compute f(t, x) both for double and AD<double>
5357 template <class Scalar>
5458 void Ode(
55 const Scalar &t,
59 const Scalar &t,
5660 const CPPAD_TESTVECTOR(Scalar) &x,
5761 CPPAD_TESTVECTOR(Scalar) &f)
5862 { size_t n = x.size();
5963 Scalar ti(1);
6064 f[0] = Scalar(1);
61 size_t i;
62 for(i = 1; i < n; i++)
65 for(size_t i = 1; i < n; i++)
6366 { ti *= t;
64 // convert int(size_t) to avoid warning
65 // on _MSC_VER systems
66 if( use_x )
67 f[i] = int(i+1) * x[i-1];
67 if( use_x_ )
68 f[i] = Scalar(i+1) * x[i-1];
6869 else
69 f[i] = int(i+1) * ti;
70 f[i] = Scalar(i+1) * ti;
7071 }
7172 }
7273
7374 // compute partial of f(t, x) w.r.t. t using AD
7475 void Ode_ind(
75 const double &t,
76 const double &t,
7677 const CPPAD_TESTVECTOR(double) &x,
7778 CPPAD_TESTVECTOR(double) &f_t)
7879 { using namespace CppAD;
7980
8081 size_t n = x.size();
81 CPPAD_TESTVECTOR(AD<double>) T(1);
82 CPPAD_TESTVECTOR(AD<double>) X(n);
83 CPPAD_TESTVECTOR(AD<double>) F(n);
84
85 // set argument values
86 T[0] = t;
87 size_t i;
88 for(i = 0; i < n; i++)
89 X[i] = x[i];
90
91 // declare independent variables
92 Independent(T);
93
94 // compute f(t, x)
95 this->Ode(T[0], X, F);
96
97 // define AD function object
98 ADFun<double> fun(T, F);
99
82 if( ode_ind_.size_var() == 0 )
83 { // record function that evaluates f(t, x)
84 // with t as independent variable and x as dynamcic parameter
85 CPPAD_TESTVECTOR(AD<double>) at(1);
86 CPPAD_TESTVECTOR(AD<double>) ax(n);
87 CPPAD_TESTVECTOR(AD<double>) af(n);
88
89 // set argument values
90 at[0] = t;
91 size_t i;
92 for(i = 0; i < n; i++)
93 ax[i] = x[i];
94
95 // declare independent variables and dynamic parameters
96 size_t abort_op_index = 0;
97 bool record_compare = false;
98 Independent(at, abort_op_index, record_compare, ax);
99
100 // compute f(t, x)
101 this->Ode(at[0], ax, af);
102
103 // define AD function object
104 ADFun<double> fun(at, af);
105
106 // store result in ode_ind_ so can be re-used
107 ode_ind_.swap(fun);
108 assert( ode_ind_.size_var() != 0 );
109 assert( fun.size_var() == 0 );
110 }
100111 // compute partial of f w.r.t t
101 CPPAD_TESTVECTOR(double) dt(1);
102 dt[0] = 1.;
103 f_t = fun.Forward(1, dt);
112 CPPAD_TESTVECTOR(double) t_vec(1);
113 t_vec[0] = t;
114 ode_ind_.new_dynamic(x);
115 f_t = ode_ind_.Jacobian(t_vec); // partial f(t, x) w.r.t. t
104116 }
105117
106118 // compute partial of f(t, x) w.r.t. x using AD
107119 void Ode_dep(
108 const double &t,
120 const double &t,
109121 const CPPAD_TESTVECTOR(double) &x,
110122 CPPAD_TESTVECTOR(double) &f_x)
111123 { using namespace CppAD;
112124
113125 size_t n = x.size();
114 CPPAD_TESTVECTOR(AD<double>) T(1);
115 CPPAD_TESTVECTOR(AD<double>) X(n);
116 CPPAD_TESTVECTOR(AD<double>) F(n);
117
118 // set argument values
119 T[0] = t;
120 size_t i, j;
121 for(i = 0; i < n; i++)
122 X[i] = x[i];
123
124 // declare independent variables
125 Independent(X);
126
127 // compute f(t, x)
128 this->Ode(T[0], X, F);
129
130 // define AD function object
131 ADFun<double> fun(X, F);
132
126 if( ode_dep_.size_var() == 0 )
127 { // record function that evaluates f(t, x)
128 // with x as independent variable and t as dynamcic parameter
129 CPPAD_TESTVECTOR(AD<double>) at(1);
130 CPPAD_TESTVECTOR(AD<double>) ax(n);
131 CPPAD_TESTVECTOR(AD<double>) af(n);
132
133 // set argument values
134 at[0] = t;
135 for(size_t i = 0; i < n; i++)
136 ax[i] = x[i];
137
138 // declare independent variables
139 size_t abort_op_index = 0;
140 bool record_compare = false;
141 Independent(ax, abort_op_index, record_compare, at);
142
143 // compute f(t, x)
144 this->Ode(at[0], ax, af);
145
146 // define AD function object
147 ADFun<double> fun(ax, af);
148
149 // store result in ode_dep_ so can be re-used
150 ode_dep_.swap(fun);
151 assert( ode_ind_.size_var() != 0 );
152 assert( fun.size_var() == 0 );
153 }
133154 // compute partial of f w.r.t x
134 CPPAD_TESTVECTOR(double) dx(n);
135 CPPAD_TESTVECTOR(double) df(n);
136 for(j = 0; j < n; j++)
137 dx[j] = 0.;
138 for(j = 0; j < n; j++)
139 { dx[j] = 1.;
140 df = fun.Forward(1, dx);
141 for(i = 0; i < n; i++)
142 f_x [i * n + j] = df[i];
143 dx[j] = 0.;
144 }
145 }
146
147 private:
148 const bool use_x;
155 CPPAD_TESTVECTOR(double) t_vec(1), dx(n), df(n);
156 t_vec[0] = t;
157 ode_dep_.new_dynamic(t_vec);
158 f_x = ode_dep_.Jacobian(x); // partial f(t, x) w.r.t. x
159 }
160
149161
150162 };
151163 }
152164
153165 bool rosen_34(void)
154166 { bool ok = true; // initial return value
155 size_t i; // temporary indices
156167
157168 using CppAD::NearEqual;
158169 double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
164175
165176 // xi = X(0)
166177 CPPAD_TESTVECTOR(double) xi(n);
167 for(i = 0; i <n; i++)
178 for(size_t i = 0; i <n; i++)
168179 xi[i] = 0.;
169180
170 size_t use_x;
171 for( use_x = 0; use_x < 2; use_x++)
181 for(size_t use_x = 0; use_x < 2; use_x++)
172182 { // function object depends on value of use_x
173183 Fun F(use_x > 0);
174184
177187 xf = CppAD::Rosen34(F, M, ti, tf, xi, e);
178188
179189 double check = tf;
180 for(i = 0; i < n; i++)
190 for(size_t i = 0; i < n; i++)
181191 { // check that error is always positive
182192 ok &= (e[i] >= 0.);
183193 // 4th order method is exact for i < 4
+11
-11
include/cppad/core/fun_construct.hpp less more
367367 function_name_.swap( f.function_name_ );
368368 //
369369 // bool objects
370 exceed_collision_limit_ = f.exceed_collision_limit_;
371 base2ad_return_value_ = false; // f might be, but this is not
372 has_been_optimized_ = f.has_been_optimized_;
373 check_for_nan_ = f.check_for_nan_;
370 std::swap( exceed_collision_limit_ , f.exceed_collision_limit_);
371 std::swap( base2ad_return_value_ , f.base2ad_return_value_);
372 std::swap( has_been_optimized_ , f.has_been_optimized_);
373 std::swap( check_for_nan_ , f.check_for_nan_);
374374 //
375375 // size_t objects
376 compare_change_count_ = f.compare_change_count_;
377 compare_change_number_ = f.compare_change_number_;
378 compare_change_op_index_ = f.compare_change_op_index_;
379 num_order_taylor_ = f.num_order_taylor_;
380 cap_order_taylor_ = f.cap_order_taylor_;
381 num_direction_taylor_ = f.num_direction_taylor_;
382 num_var_tape_ = f.num_var_tape_;
376 std::swap( compare_change_count_ , f.compare_change_count_);
377 std::swap( compare_change_number_ , f.compare_change_number_);
378 std::swap( compare_change_op_index_ , f.compare_change_op_index_);
379 std::swap( num_order_taylor_ , f.num_order_taylor_);
380 std::swap( cap_order_taylor_ , f.cap_order_taylor_);
381 std::swap( num_direction_taylor_ , f.num_direction_taylor_);
382 std::swap( num_var_tape_ , f.num_var_tape_);
383383 //
384384 // pod_vector objects
385385 ind_taddr_.swap( f.ind_taddr_);
+14
-0
omh/appendix/whats_new/20.omh less more
4747
4848 $section Changes and Additions to CppAD During 2020$$
4949
50 $head 12-20$$
51 $list number$$
52 The $cref rosen_34.cpp$$ example was changed to use
53 $cref/dynamic/independent/dynamic/$$ parameters so it does not
54 need to record a new function for every argument value.
55 The new version of this example also uses the
56 ADFun $cref/swap/FunConstruct/swap/$$ operation to avoid
57 making a separate copy of an ADFun object.
58 $lnext
59 Fix a bug in the
60 ADFun $cref/swap/FunConstruct/swap/$$, which was added on 12-19.
61 $lend
62
63
5064 $head 12-19$$
5165 Merge in $code cg$$ branch which had the following changes:
5266