一段自动求导的代码，C++er 来看下

lost_templar

自己写了做自动数值求导的；

一次导数正常，二次就会出错，有些莫名其妙。

一次求导正确：

1. double ds( double x )
   
2. {
   
3.     return std::sin(x);
   
4. }

复制代码

求导使用：

1. derivative<0, double(double)>( ds )( 0.0 );
   

复制代码

输出 1

二次求导则会出错：

1.     auto const& dds = derivative<0, double(double)>(ds);
   
2. 
   
3.     derivative<0, double(double)>(dds)( 0.0 );
   

复制代码

依旧输出 1，不知为何


derivitative 实现的大致代码如下：

1.     template<std::size_t M, typename Dummy_Type> struct derivative;
   
2. 
   
3.     template< std::size_t M, typename R, typename... Types >
   
4.     struct derivative< M, R(Types...) >
   
5.     {
   
6.         typedef R                                               return_type;
   
7.         typedef std::function< R(Types...) >                    function_type;
   
8.         typedef typename type_at< M, Types... >::result_type    result_type;
   
9.         typedef std::size_t                                     size_type;
   
10. 
    
11.         function_type ff;
    
12. 
    
13.         template< typename Function >
    
14.         derivative( const Function& ff_ ) : ff( ff_ ) { }
    
15. 
    
16.         result_type operator()( Types... vts ) const
    
17.         {
    
18.             typedef typename stepper_value_type<result_type>::value_type value_type; //fix for complex
    
19. 
    
20.             const value_type decrease_step     = 1.6180339887498948482;
    
21.             const value_type decrease_step_2   = 2.6180339887498948482;
    
22.             const value_type safe_boundary     = 2.6180339887498948482;
    
23.             const value_type start_step        = 1.0e-8;
    
24.             const result_type x                = value_at<M, Types...>()( vts... );
    
25.             const size_type iter_depth         = 64;
    
26. 
    
27.             result_type error   = std::numeric_limits<value_type>::max();
    
28.             result_type ans     = std::numeric_limits<value_type>::quiet_NaN();
    
29.             value_type step     = start_step;
    
30. 
    
31.             auto const& lhs_f = [&step]( result_type x ) { return x - step; };
    
32.             auto const& rhs_f = [&step]( result_type x ) { return x + step; };
    
33. 
    
34.             oscillate_function< M, return_type, Types... > const lhs_of( ff, lhs_f );
    
35.             oscillate_function< M, return_type, Types... > const rhs_of( ff, rhs_f );
    
36. 
    
37.             //matrix<result_type> a( iter_depth, iter_depth );
    
38.             result_type a[iter_depth][iter_depth];
    
39. 
    
40.             a[0][0] = ( rhs_of( vts... ) - lhs_of( vts... ) ) / ( step + step );
    
41. 
    
42.             //for ( size_type i = 1; i != a.row(); ++i )
    
43.             for ( size_type i = 1; i != iter_depth; ++i )
    
44.             {
    
45.                 step /= decrease_step;
    
46.                 a[i][0] = ( rhs_of( vts... ) - lhs_of( vts... ) ) / ( step + step );
    
47. 
    
48.                 result_type factor = decrease_step_2;
    
49. 
    
50.                 for ( size_type j = 1; j <= i; ++j )
    
51.                 {
    
52.                     const result_type factor_1 = factor - result_type(1);
    
53. 
    
54.                     a[i][j]  = ( a[i][j-1] * factor - a[i-1][j-1] ) / factor_1;
    
55. 
    
56.                     factor *= decrease_step_2;
    
57. 
    
58.                     const result_type error_so_far = std::max( std::abs(a[i][j]-a[i][j-1]), std::abs(a[i][j]-a[i-1][j-1]) );
    
59. 
    
60.                     if ( error > error_so_far )
    
61.                     {
    
62.                         error = error_so_far;
    
63.                         ans = a[i][j];
    
64.                     }
    
65.                 }
    
66. 
    
67.                 if ( std::abs( a[i][i] - a[i-1][i-1] ) >=  safe_boundary * error )
    
68.                     return ans;
    
69. 
    
70.             }
    
71. 
    
72.             return ans;
    
73.         }
    
74. 
    
75.     };
    

复制代码

其中 type_at 如下实现：

1.     template< std::size_t N, typename T, typename... Types >
   
2.     struct type_at
   
3.     {
   
4.         static_assert( N < sizeof...(Types)+1, "dim size exceeds limitation!" );
   
5.         typedef typename type_at<N-1, Types...>::result_type result_type;
   
6.     };
   
7. 
   
8.     template<typename T, typename... Types>
   
9.     struct type_at< 0, T, Types...>
   
10.     {
    
11.         typedef T result_type;
    
12.     };
    

复制代码

那个 matrix 可以用个二维数组代替，应该不难理解。

lost_templar

不好意思，漏了 oscillate_function 的定义：

1.     namespace oscillate_function_private_fdspojiasldkjasasfdioj4asfd4d
   
2.     {
   
3.     template< typename R, std::size_t Bn, std::size_t N >
   
4.     struct oscillate_backward_impl
   
5.     {
   
6.         template<typename F, typename Type, typename... Types >
   
7.         R impl( F F_, Type t, Types... vts ) const
   
8.         {
   
9.             return oscillate_backward_impl<R, Bn+1, N>().impl( F_, vts..., t );
   
10.         }
    
11.     };
    
12. 
    
13.     template< typename R, std::size_t N >
    
14.     struct oscillate_backward_impl<R, N, N>
    
15.     {
    
16.         template< typename F, typename... Types >
    
17.         R impl( F F_, Types... vts ) const
    
18.         {
    
19.             return F_( vts... );
    
20.         }
    
21.     };
    
22. 
    
23.     template<typename R, std::size_t Fn, std::size_t Bn>
    
24.     struct oscillate_forward_impl
    
25.     {
    
26.         template<typename F, typename f, typename Type, typename... Types>
    
27.         R impl( F F_, f f_, Type t, Types... vts ) const
    
28.         {
    
29.             return oscillate_forward_impl< R, Fn-1, Bn >().impl( F_, f_, vts..., t );
    
30.         }
    
31.     };
    
32. 
    
33.     template<typename R, std::size_t Bn>
    
34.     struct oscillate_forward_impl< R, 0, Bn>
    
35.     {
    
36.         template<typename F, typename f, typename Type, typename... Types>
    
37.         R impl( F F_, f f_, Type t, Types... vts ) const
    
38.         {
    
39.             return oscillate_backward_impl<R, 0, Bn>().impl( F_, vts..., f_(t) );
    
40.         }
    
41.     };
    
42. 
    
43.     }//namespace oscillate_function_private_fdspojiasldkjasasfdioj4asfd4d
    
44. 
    
45.     template< std::size_t N, typename R, typename... Types >
    
46.     struct oscillate_function
    
47.     {
    
48.         typedef R return_type;
    
49.         typedef typename type_at< N, Types... >::result_type oscillate_type;
    
50.         typedef std::function<R(Types...)> function_type;
    
51.         typedef std::function<oscillate_type(oscillate_type)> oscillate_function_type;
    
52. 
    
53.         static_assert( N < sizeof...(Types), "dim size exceeds arguments limitation." );
    
54. 
    
55.         function_type F;
    
56.         oscillate_function_type f;
    
57. 
    
58.         oscillate_function( const function_type& F_, const oscillate_function_type& f_ ) : F( F_ ), f( f_ ) {}
    
59. 
    
60.         return_type operator()( Types... vts ) const
    
61.         {
    
62.             using namespace oscillate_function_private_fdspojiasldkjasasfdioj4asfd4d;
    
63.             return oscillate_forward_impl< R, N, sizeof...(vts)-N-1 >().impl( F, f, vts... );
    
64.         }
    
65. 
    
66.     };
    

复制代码

lost_templar

好吧，还有个 value_at

1.     template< std::size_t N, typename T, typename... Types >
   
2.     struct value_at
   
3.     {
   
4.         typedef typename type_at< N, T, Types...>::result_type result_type;
   
5. 
   
6.         result_type operator()( T, Types... vts ) const
   
7.         {
   
8.             return value_at<N-1, Types...>()( vts... );
   
9.         }
   
10.     };
    
11. 
    
12.     template< typename T, typename... Types >
    
13.     struct value_at< 0, T, Types... >
    
14.     {
    
15.         typedef T result_type;
    
16. 
    
17.         result_type operator()( T vt, Types... ) const
    
18.         {
    
19.             return vt;
    
20.         }
    
21.     };
    

复制代码

淫呢？.......

幻の上帝

给matrix定义吧，用二维数组代替编译不过……

幻の上帝

还有stepper_value_type呢。
typedef typename stepper_value_type<result_type>::value_type value_type;
虽然这行用typedef double value_type;直接代替好像问题不大？

lost_templar

post deleted.........

lost_templar

回复 5# 幻の上帝


    steper_value_type 可以直接用 value_type 代替；

至于矩阵那个，可以修改下两行通过：

1.             oscillate_function< M, return_type, Types... > const rhs_of( ff, rhs_f );
   
2. 
   
3.             //matrix<result_type> a( iter_depth, iter_depth );
   
4.             result_type a[iter_depth][iter_depth];
   
5. 
   
6.             a[0][0] = ( rhs_of( vts... ) - lhs_of( vts... ) ) / ( step + step );
   
7. 
   
8.             //for ( size_type i = 1; i != a.row(); ++i )
   
9.             for ( size_type i = 1; i != iter_depth; ++i )
   
10.             {
    

复制代码

幻の上帝

似乎迭代一次就在std::abs( a[i][i] - a[i-1][i-1] ) >=  safe_boundary * error后面停了，是预期吗？（不过就算不停结果也一样不对……）
另外matrix\functional下min.hpp和max.hpp在debug时有问题， assert( m.size() ) > 0 ; 看来得是assert( m.size() > 0 );。

lost_templar

幻の上帝 发表于 2013-04-13 08:44
似乎迭代一次就在std::abs( a - a ) >=  safe_boundary * error后面停了，是预期吗？（不过就算不停结果也一 ...

停顿是预期的行为；

后边那两个 assert 是我的失误，那两个函数自从写出来就从来没有用过，也没有测试过{:3_190:}

另外我检查过 d(sin x)/dx 的输出，

1. double ds( double x )
   
2. {
   
3.     return std::sin(x);
   
4. }

复制代码

1.       for ( double d_start = -0.1; d_start <= 0.1; d_start += 0.001 )
   
2.           std::cout << derivative<0, double(double)>( ds )( d_start ) << "\n";

复制代码

是一个预期的 cos 函数：

1. 0.995004
   
2. 0.995104
   
3. 0.995202
   
4. 0.995299
   
5. 0.995396
   
6. 0.995491
   
7. 0.995585
   
8. 0.995679
   
9. 0.995771
   
10. 0.995862
    
11. 0.995953
    
12. 0.996042
    
13. 0.99613
    
14. 0.996218
    
15. 0.996304
    
16. 0.99639
    
17. 0.996474
    
18. 0.996557
    
19. 0.99664
    
20. 0.996721
    
21. 0.996802
    
22. 0.996881
    
23. 0.99696
    
24. 0.997037
    
25. 0.997113
    
26. 0.997189
    
27. 0.997263
    
28. 0.997337
    
29. 0.997409
    
30. 0.997481
    
31. 0.997551
    
32. 0.99762
    
33. 0.997689
    
34. 0.997756
    
35. 0.997823
    
36. 0.997888
    
37. 0.997953
    
38. 0.998016
    
39. 0.998079
    
40. 0.99814
    
41. 0.998201
    
42. 0.99826
    
43. 0.998318
    
44. 0.998376
    
45. 0.998432
    
46. 0.998488
    
47. 0.998542
    
48. 0.998596
    
49. 0.998648
    
50. 0.9987
    
51. 0.99875
    
52. 0.9988
    
53. 0.998848
    
54. 0.998896
    
55. 0.998942
    
56. 0.998988
    
57. 0.999032
    
58. 0.999076
    
59. 0.999118
    
60. 0.99916
    
61. 0.9992
    
62. 0.99924
    
63. 0.999278
    
64. 0.999316
    
65. 0.999352
    
66. 0.999388
    
67. 0.999422
    
68. 0.999456
    
69. 0.999488
    
70. 0.99952
    
71. 0.99955
    
72. 0.99958
    
73. 0.999608
    
74. 0.999636
    
75. 0.999662
    
76. 0.999688
    
77. 0.999712
    
78. 0.999736
    
79. 0.999758
    
80. 0.99978
    
81. 0.9998
    
82. 0.99982
    
83. 0.999838
    
84. 0.999856
    
85. 0.999872
    
86. 0.999888
    
87. 0.999902
    
88. 0.999916
    
89. 0.999928
    
90. 0.99994
    
91. 0.99995
    
92. 0.99996
    
93. 0.999968
    
94. 0.999976
    
95. 0.999982
    
96. 0.999987
    
97. 0.999992
    
98. 0.999995
    
99. 0.999998
    
100. 0.999999
     
101. 1
     
102. 0.999999
     
103. 0.999998
     
104. 0.999995
     
105. 0.999992
     
106. 0.999987
     
107. 0.999982
     
108. 0.999976
     
109. 0.999968
     
110. 0.99996
     
111. 0.99995
     
112. 0.99994
     
113. 0.999928
     
114. 0.999916
     
115. 0.999902
     
116. 0.999888
     
117. 0.999872
     
118. 0.999856
     
119. 0.999838
     
120. 0.99982
     
121. 0.9998
     
122. 0.99978
     
123. 0.999758
     
124. 0.999736
     
125. 0.999712
     
126. 0.999688
     
127. 0.999662
     
128. 0.999636
     
129. 0.999608
     
130. 0.99958
     
131. 0.99955
     
132. 0.99952
     
133. 0.999488
     
134. 0.999456
     
135. 0.999422
     
136. 0.999388
     
137. 0.999352
     
138. 0.999316
     
139. 0.999278
     
140. 0.99924
     
141. 0.9992
     
142. 0.99916
     
143. 0.999118
     
144. 0.999076
     
145. 0.999032
     
146. 0.998988
     
147. 0.998942
     
148. 0.998896
     
149. 0.998848
     
150. 0.9988
     
151. 0.99875
     
152. 0.9987
     
153. 0.998648
     
154. 0.998596
     
155. 0.998542
     
156. 0.998488
     
157. 0.998432
     
158. 0.998376
     
159. 0.998318
     
160. 0.99826
     
161. 0.998201
     
162. 0.99814
     
163. 0.998079
     
164. 0.998016
     
165. 0.997953
     
166. 0.997888
     
167. 0.997823
     
168. 0.997756
     
169. 0.997689
     
170. 0.99762
     
171. 0.997551
     
172. 0.997481
     
173. 0.997409
     
174. 0.997337
     
175. 0.997263
     
176. 0.997189
     
177. 0.997113
     
178. 0.997037
     
179. 0.99696
     
180. 0.996881
     
181. 0.996802
     
182. 0.996721
     
183. 0.99664
     
184. 0.996557
     
185. 0.996474
     
186. 0.99639
     
187. 0.996304
     
188. 0.996218
     
189. 0.99613
     
190. 0.996042
     
191. 0.995953
     
192. 0.995862
     
193. 0.995771
     
194. 0.995679
     
195. 0.995585
     
196. 0.995491
     
197. 0.995396
     
198. 0.995299
     
199. 0.995202
     
200. 0.995104

复制代码

lost_templar

回复 8# 幻の上帝


    我觉得可能是问题出在用 derivative<0, double(double)> 这样一个类型的参数来初始化另外一个 derivative<0, double(double)> 这样一个类型的时候：

1. auto const& dds = derivative<0, double(double)>(ds);
   
2. derivative<0, double(double)>(dds)( 0.0 )

复制代码

上边代码的第二句本来预期的行为是调用这个带 template 的 constructor：

1. struct derivative
   
2. {
   
3.        function_type ff;
   
4. 
   
5.         template< typename Function >
   
6.         derivative( const Function& ff_ ) : ff( ff_ ) { }
   
7. };

复制代码

但是在实际调用的是一个缺省的隐藏构造函数：

1. struct derivative
   
2. {
   
3.     derivative( const derivative& other );
   
4. };

复制代码

因此得到错误的结果