亚洲av成人无遮挡网站在线观看,少妇性bbb搡bbb爽爽爽,亚洲av日韩精品久久久久久,兔费看少妇性l交大片免费,无码少妇一区二区三区

Chinaunix

標題: 一段自動求導(dǎo)的代碼，C++er 來看下 [打印本頁]
作者: lost_templar    時間: 2013-04-11 21:01
標題: 一段自動求導(dǎo)的代碼，C++er 來看下
本帖最后由 lost_templar 于 2013-04-12 17:11 編輯

自己寫了做自動數(shù)值求導(dǎo)的；

一次導(dǎo)數(shù)正常，二次就會出錯，有些莫名其妙。

一次求導(dǎo)正確：
  1. double ds( double x )
  2. {
  3.     return std::sin(x);
  4. }
復(fù)制代碼
求導(dǎo)使用：
  1. derivative<0, double(double)>( ds )( 0.0 );
復(fù)制代碼
輸出 1

二次求導(dǎo)則會出錯：
  1.     auto const& dds = derivative<0, double(double)>(ds);

  3.     derivative<0, double(double)>(dds)( 0.0 );
復(fù)制代碼
依舊輸出 1，不知為何


derivitative 實現(xiàn)的大致代碼如下：
  1.     template<std::size_t M, typename Dummy_Type> struct derivative;

  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;

  11.         function_type ff;

  13.         template< typename Function >
  14.         derivative( const Function& ff_ ) : ff( ff_ ) { }

  16.         result_type operator()( Types... vts ) const
  17.         {
  18.             typedef typename stepper_value_type<result_type>::value_type value_type; //fix for complex

  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;

  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;

  31.             auto const& lhs_f = [&step]( result_type x ) { return x - step; };
  32.             auto const& rhs_f = [&step]( result_type x ) { return x + step; };

  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 );

  37.             //matrix<result_type> a( iter_depth, iter_depth );
  38.             result_type a[iter_depth][iter_depth];

  40.             a[0][0] = ( rhs_of( vts... ) - lhs_of( vts... ) ) / ( step + step );

  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 );

  48.                 result_type factor = decrease_step_2;

  50.                 for ( size_type j = 1; j <= i; ++j )
  51.                 {
  52.                     const result_type factor_1 = factor - result_type(1);

  54.                     a[i][j]  = ( a[i][j-1] * factor - a[i-1][j-1] ) / factor_1;

  56.                     factor *= decrease_step_2;

  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]) );

  60.                     if ( error > error_so_far )
  61.                     {
  62.                         error = error_so_far;
  63.                         ans = a[i][j];
  64.                     }
  65.                 }

  67.                 if ( std::abs( a[i][i] - a[i-1][i-1] ) >=  safe_boundary * error )
  68.                     return ans;

  70.             }

  72.             return ans;
  73.         }

  75.     };
復(fù)制代碼
其中 type_at 如下實現(xiàn)：
  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.     };

  8.     template<typename T, typename... Types>
  9.     struct type_at< 0, T, Types...>
  10.     {
  11.         typedef T result_type;
  12.     };
復(fù)制代碼
那個 matrix 可以用個二維數(shù)組代替，應(yīng)該不難理解。

作者: lost_templar    時間: 2013-04-11 21:13
不好意思，漏了 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.     };

  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.     };

  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.     };

  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.     };

  43.     }//namespace oscillate_function_private_fdspojiasldkjasasfdioj4asfd4d

  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;

  53.         static_assert( N < sizeof...(Types), "dim size exceeds arguments limitation." );

  55.         function_type F;
  56.         oscillate_function_type f;

  58.         oscillate_function( const function_type& F_, const oscillate_function_type& f_ ) : F( F_ ), f( f_ ) {}

  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.         }

  66.     };
復(fù)制代碼
作者: lost_templar    時間: 2013-04-11 21:15
本帖最后由 lost_templar 于 2013-04-11 21:28 編輯

好吧，還有個 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;

  6.         result_type operator()( T, Types... vts ) const
  7.         {
  8.             return value_at<N-1, Types...>()( vts... );
  9.         }
  10.     };

  12.     template< typename T, typename... Types >
  13.     struct value_at< 0, T, Types... >
  14.     {
  15.         typedef T result_type;

  17.         result_type operator()( T vt, Types... ) const
  18.         {
  19.             return vt;
  20.         }
  21.     };
復(fù)制代碼
淫呢？.......

作者: 幻の上帝    時間: 2013-04-12 16:12
給matrix定義吧，用二維數(shù)組代替編譯不過……
作者: 幻の上帝    時間: 2013-04-12 16:12
還有stepper_value_type呢。
typedef typename stepper_value_type<result_type>::value_type value_type;
雖然這行用typedef double value_type;直接代替好像問題不大？

作者: lost_templar    時間: 2013-04-12 16:33
本帖最后由 lost_templar 于 2013-06-10 05:59 編輯

post deleted.........
作者: lost_templar    時間: 2013-04-12 16:35
回復(fù) 5# 幻の上帝


    steper_value_type 可以直接用 value_type 代替；

至于矩陣那個，可以修改下兩行通過：
  1.             oscillate_function< M, return_type, Types... > const rhs_of( ff, rhs_f );

  3.             //matrix<result_type> a( iter_depth, iter_depth );
  4.             result_type a[iter_depth][iter_depth];

  6.             a[0][0] = ( rhs_of( vts... ) - lhs_of( vts... ) ) / ( step + step );

  8.             //for ( size_type i = 1; i != a.row(); ++i )
  9.             for ( size_type i = 1; i != iter_depth; ++i )
  10.             {
復(fù)制代碼
作者: 幻の上帝    時間: 2013-04-13 08:44
似乎迭代一次就在std::abs( a[i][i] - a[i-1][i-1] ) >=  safe_boundary * error后面停了，是預(yù)期嗎？（不過就算不停結(jié)果也一樣不對……）
另外matrix\functional下min.hpp和max.hpp在debug時有問題， assert( m.size() ) > 0 ; 看來得是assert( m.size() > 0 );。

作者: lost_templar    時間: 2013-04-13 20:09
本帖最后由 lost_templar 于 2013-04-13 20:10 編輯
幻の上帝 發(fā)表于 2013-04-13 08:44
似乎迭代一次就在std::abs( a - a ) >=  safe_boundary * error后面停了，是預(yù)期嗎？（不過就算不停結(jié)果也一 ...


停頓是預(yù)期的行為；

后邊那兩個 assert 是我的失誤，那兩個函數(shù)自從寫出來就從來沒有用過，也沒有測試過{:3_190:}

另外我檢查過 d(sin x)/dx 的輸出，
  1. double ds( double x )
  2. {
  3.     return std::sin(x);
  4. }
復(fù)制代碼
  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";
復(fù)制代碼
是一個預(yù)期的 cos 函數(shù)：
  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
復(fù)制代碼
作者: lost_templar    時間: 2013-04-13 20:38
回復(fù) 8# 幻の上帝


    我覺得可能是問題出在用 derivative<0, double(double)> 這樣一個類型的參數(shù)來初始化另外一個 derivative<0, double(double)> 這樣一個類型的時候：
  1. auto const& dds = derivative<0, double(double)>(ds);
  2. derivative<0, double(double)>(dds)( 0.0 )
復(fù)制代碼
上邊代碼的第二句本來預(yù)期的行為是調(diào)用這個帶 template 的 constructor：
  1. struct derivative
  2. {
  3.        function_type ff;

  5.         template< typename Function >
  6.         derivative( const Function& ff_ ) : ff( ff_ ) { }
  7. };
復(fù)制代碼
但是在實際調(diào)用的是一個缺省的隱藏構(gòu)造函數(shù)：
  1. struct derivative
  2. {
  3.     derivative( const derivative& other );
  4. };
復(fù)制代碼
因此得到錯誤的結(jié)果

作者: lost_templar    時間: 2013-04-13 20:59
本帖最后由 lost_templar 于 2013-04-13 21:04 編輯

回復(fù) 8# 幻の上帝


    確實是 copy constructor 被觸發(fā)了的原因，我加入這段代碼：
  1.         derivative& operator = ( const derivative& other )
  2.         {
  3.             ff = [&](Types... vts) -> return_type { return other( vts... ); };
  4.         }

  6.         derivative( const derivative& other )
  7.         {
  8.             ff = [&](Types... vts) -> return_type { return other( vts... ); };
  9.         }
復(fù)制代碼
到類中之后，
  1.     auto const& dds = derivative<0, double(double)>(ds);
  2.     std::cout << "d( d sin(x) / dx) / dx at (0.0) is " << derivative<0, double(double)>(dds)( 0.0 ) << "\n";
復(fù)制代碼
輸出正確的結(jié)果 0， 但是觸發(fā)了 move constructor 的這段代碼依舊輸出 1：

  2. std::cout << "d( d sin(x) / dx) / dx at (0.0) is " << derivative<0, double(double)>(derivative<0, double(double)>(ds) )( 0.0 ) << "\n";
復(fù)制代碼
盡管我已經(jīng)把 move constructor 和 move assignment operator 寫為：

  2.         derivative( derivative&& other )
  3.         {
  4.             ff = [&](Types... vts) -> return_type { return other( vts... ); };
  5.         }
  6.         derivative& operator = ( derivative&& other )
  7.         {
  8.             ff = [&](Types... vts) -> return_type { return other( vts... ); };
  9.         }
復(fù)制代碼
而如果將 move constructor 和 move assignment operator 定義為 delete，則上邊的代碼不能通過編譯 {:3_193:}
作者: lost_templar    時間: 2013-04-13 21:32
回復(fù) 5# 幻の上帝


    我已經(jīng)在 http://134.60.10.69/gitweb/ 更新了代碼，現(xiàn)在只有
  1. derivative<0, double(double)>(derivative<0, double(double)>(ds) )( 0.0 )
復(fù)制代碼
這個搞不定了，它非要去調(diào)用 move constructor
作者: lost_templar    時間: 2013-04-13 23:42
回復(fù) 12# lost_templar

覺得自己進入誤區(qū)了，對于這樣的代碼
  1. derivative a(something);
  2. A b = a;
復(fù)制代碼
讓 a 和 b 具有不同的行為，本身就是很邪惡的{:3_182:}

   
作者: 幻の上帝    時間: 2013-04-14 08:54
回復(fù) 13# lost_templar

……我也試出來了，是copy ctor的問題。再包裝一層std::function或lambda傳進去就正常了。
確實，如果自己重寫copy ctor，可以解決一部分問題，不過會讓代碼的意義不太明確。
可以學(xué)標準庫的做法，考慮寫一個模版（比如叫make_derivative）轉(zhuǎn)發(fā)參數(shù)，并對derivative提供重載，這樣就區(qū)分開來了。而且可以讓用戶不用照抄模版實參。
或者提供顯式的tag放在構(gòu)造函數(shù)里，然后提供tag的重載，不過可能寫起來比較麻煩一點。
另外建議考慮一下使用std::result_of或decltype推導(dǎo)返回類型。
作者: lost_templar    時間: 2013-04-14 16:35
回復(fù) 14# 幻の上帝

非常感謝！
非常有用的建議，我會仔細地研究一下標準庫，再看下怎么設(shè)計這里的接口。


   



歡迎光臨 Chinaunix (http://www.72891.cn/) Powered by Discuz! X3.2