00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #ifndef __NUMERICALDIFFERENTIATION_H__
00017 #define __NUMERICALDIFFERENTIATION_H__
00018
00019
00020
00021 #include "vector.h"
00022 #include "matrix.h"
00023
00024
00025
00026 #include "../../parsers/tinyxml/tinyxml.h"
00027
00028
00029 namespace OpenNN
00030 {
00031
00034
00035 class NumericalDifferentiation
00036 {
00037
00038 public:
00039
00040
00041
00042 explicit NumericalDifferentiation(void);
00043
00044
00045
00046 NumericalDifferentiation(const NumericalDifferentiation&);
00047
00048
00049
00050 virtual ~NumericalDifferentiation(void);
00051
00052
00053
00054 NumericalDifferentiation& operator = (const NumericalDifferentiation&);
00055
00056
00057
00058 bool operator == (const NumericalDifferentiation&) const;
00059
00061
00062 enum NumericalDifferentiationMethod{ForwardDifferences, CentralDifferences};
00063
00064
00065
00066 const NumericalDifferentiationMethod& get_numerical_differentiation_method(void) const;
00067 std::string write_numerical_differentiation_method(void) const;
00068
00069 const unsigned int& get_precision_digits(void) const;
00070
00071 const bool& get_display(void) const;
00072
00073 void set(const NumericalDifferentiation&);
00074
00075 void set_numerical_differentiation_method(const NumericalDifferentiationMethod&);
00076 void set_numerical_differentiation_method(const std::string&);
00077
00078 void set_precision_digits(const unsigned int&);
00079
00080 void set_display(const bool&);
00081
00082 void set_default(void);
00083
00084 double calculate_h(const double&) const;
00085
00086 Vector<double> calculate_h(const Vector<double>&) const;
00087
00088
00089
00090 TiXmlElement* to_XML(void) const;
00091 void from_XML(TiXmlElement*);
00092
00093
00094
00095
00096
00097
00098
00103
00104 template<class T>
00105 double calculate_forward_differences_derivative(const T& t, double (T::*f)(const double&) const, const double& x) const
00106 {
00107 const double y = (t.*f)(x);
00108
00109 const double h = calculate_h(x);
00110
00111 const double y_forward = (t.*f)(x + h);
00112
00113 const double d = (y_forward - y)/h;
00114
00115 return(d);
00116 }
00117
00118
00119
00120
00125
00126 template<class T>
00127 double calculate_central_differences_derivative(const T& t, double (T::*f)(const double&) const , const double& x) const
00128 {
00129 const double h = calculate_h(x);
00130
00131 const double y_forward = (t.*f)(x+h);
00132
00133 const double y_backward = (t.*f)(x-h);
00134
00135 const double d = (y_forward - y_backward)/(2.0*h);
00136
00137 return(d);
00138 }
00139
00140
00141
00142
00147
00148 template<class T>
00149 double calculate_derivative(const T& t, double (T::*f)(const double&) const , const double& x) const
00150 {
00151 switch(numerical_differentiation_method)
00152 {
00153 case ForwardDifferences:
00154 {
00155 return(calculate_forward_differences_derivative(t, f, x));
00156 }
00157 break;
00158
00159 case CentralDifferences:
00160 {
00161 return(calculate_central_differences_derivative(t, f, x));
00162 }
00163 break;
00164
00165 default:
00166 {
00167 std::ostringstream buffer;
00168
00169 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
00170 << "double calculate_derivative(const T&, double (T::*f)(const double&) const , const double&) const method.\n"
00171 << "Unknown numerical differentiation method.\n";
00172
00173 throw std::logic_error(buffer.str());
00174 }
00175 break;
00176 }
00177 }
00178
00179
00180
00181
00182
00187
00188 template<class T>
00189 Vector<double> calculate_forward_differences_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00190 {
00191 const Vector<double> h = calculate_h(x);
00192
00193 const Vector<double> y = (t.*f)(x);
00194
00195 const Vector<double> x_forward = x + h;
00196 const Vector<double> y_forward = (t.*f)(x_forward);
00197
00198 const Vector<double> d = (y_forward - y)/h;
00199
00200 return(d);
00201 }
00202
00203
00204
00205
00206
00211
00212 template<class T>
00213 Vector<double> calculate_central_differences_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00214 {
00215 const Vector<double> h = calculate_h(x);
00216
00217 const Vector<double> x_forward = x + h;
00218 const Vector<double> x_backward = x - h;
00219
00220 const Vector<double> y_forward = (t.*f)(x_forward);
00221 const Vector<double> y_backward = (t.*f)(x_backward);
00222
00223 const Vector<double> y = (t.*f)(x);
00224
00225 const Vector<double> d = (y_forward - y_backward)/(h*2.0);
00226
00227 return(d);
00228 }
00229
00230
00231
00232
00233
00238
00239 template<class T>
00240 Vector<double> calculate_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00241 {
00242 switch(numerical_differentiation_method)
00243 {
00244 case ForwardDifferences:
00245 {
00246 return(calculate_forward_differences_derivative(t, f, x));
00247 }
00248 break;
00249
00250 case CentralDifferences:
00251 {
00252 return(calculate_central_differences_derivative(t, f, x));
00253 }
00254 break;
00255
00256 default:
00257 {
00258 std::ostringstream buffer;
00259
00260 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
00261 << "Vector<double> calculate_derivative(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
00262 << "Unknown numerical differentiation method.\n";
00263
00264 throw std::logic_error(buffer.str());
00265 }
00266 break;
00267 }
00268 }
00269
00270
00271
00272
00273
00280
00281 template<class T>
00282 Vector<double> calculate_forward_differences_derivative(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
00283 {
00284 const Vector<double> y = (t.*f)(dummy, x);
00285
00286 const Vector<double> h = calculate_h(x);
00287 const Vector<double> x_forward = x + h;
00288
00289 const Vector<double> y_forward = (t.*f)(dummy, x_forward);
00290
00291 const Vector<double> d = (y_forward - y)/h;
00292
00293 return(d);
00294 }
00295
00296
00297
00298
00299
00306
00307 template<class T>
00308 Vector<double> calculate_central_differences_derivative(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
00309 {
00310 const Vector<double> h = calculate_h(x);
00311
00312 const Vector<double> x_forward = x + h;
00313 const Vector<double> x_backward = x - h;
00314
00315 const Vector<double> y_forward = (t.*f)(dummy, x_forward);
00316 const Vector<double> y_backward = (t.*f)(dummy, x_backward);
00317
00318 const Vector<double> d = (y_forward - y_backward)/(h*2.0);
00319
00320 return(d);
00321 }
00322
00323
00324
00325
00326
00333
00334 template<class T>
00335 Vector<double> calculate_derivative(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
00336 {
00337 switch(numerical_differentiation_method)
00338 {
00339 case ForwardDifferences:
00340 {
00341 return(calculate_forward_differences_derivative(t, f, dummy, x));
00342 }
00343 break;
00344
00345 case CentralDifferences:
00346 {
00347 return(calculate_central_differences_derivative(t, f, dummy, x));
00348 }
00349 break;
00350
00351 default:
00352 {
00353 std::ostringstream buffer;
00354
00355 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
00356 << "Vector<double> calculate_derivative(const T&, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int&, const Vector<double>&) const method.\n"
00357 << "Unknown numerical differentiation method.\n";
00358
00359 throw std::logic_error(buffer.str());
00360 }
00361 break;
00362 }
00363 }
00364
00365
00366
00367
00368
00369
00370
00375
00376 template<class T>
00377 double calculate_forward_differences_second_derivative(const T& t, double (T::*f)(const double&) const, const double& x) const
00378 {
00379 const double h = calculate_h(x);
00380
00381 const double x_forward_2 = x + 2.0*h;
00382
00383 const double y_forward_2 = (t.*f)(x_forward_2);
00384
00385 const double x_forward = x + h;
00386
00387 const double y_forward = (t.*f)(x_forward);
00388
00389 const double y = (t.*f)(x);
00390
00391 return((y_forward_2 - 2*y_forward + y)/pow(h, 2));
00392 }
00393
00394
00395
00396
00401
00402 template<class T>
00403 double calculate_central_differences_second_derivative(const T& t, double (T::*f)(const double&) const , const double& x) const
00404 {
00405 const double h = calculate_h(x);
00406
00407 const double x_forward_2 = x + 2.0*h;
00408
00409 const double y_forward_2 = (t.*f)(x_forward_2);
00410
00411 const double x_forward = x + h;
00412
00413 const double y_forward = (t.*f)(x_forward);
00414
00415 const double y = (t.*f)(x);
00416
00417 const double x_backward = x - h;
00418
00419 const double y_backward = (t.*f)(x_backward);
00420
00421 const double x_backward_2 = x - 2.0*h;
00422
00423 const double y_backward_2 = (t.*f)(x_backward_2);
00424
00425 const double d2 = (-y_forward_2 + 16.0*y_forward -30.0*y + 16.0*y_backward - y_backward_2)/(12.0*pow(h, 2));
00426
00427 return(d2);
00428 }
00429
00430
00431
00432
00437
00438 template<class T>
00439 double calculate_second_derivative(const T& t, double (T::*f)(const double&) const , const double& x) const
00440 {
00441 switch(numerical_differentiation_method)
00442 {
00443 case ForwardDifferences:
00444 {
00445 return(calculate_forward_differences_second_derivative(t, f, x));
00446 }
00447 break;
00448
00449 case CentralDifferences:
00450 {
00451 return(calculate_central_differences_second_derivative(t, f, x));
00452 }
00453 break;
00454
00455 default:
00456 {
00457 std::ostringstream buffer;
00458
00459 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
00460 << "double calculate_second_derivative(const T& t, double (T::*f)(const double&) const, const double&) const method.\n"
00461 << "Unknown numerical differentiation method.\n";
00462
00463 throw std::logic_error(buffer.str());
00464 }
00465 break;
00466 }
00467 }
00468
00469
00470
00471
00476
00477 template<class T>
00478 Vector<double> calculate_forward_differences_second_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00479 {
00480 const Vector<double> y = (t.*f)(x);
00481
00482 const Vector<double> h = calculate_h(x);
00483
00484 const Vector<double> x_forward = x + h;
00485 const Vector<double> x_forward_2 = x + h*2.0;
00486
00487 const Vector<double> y_forward = (t.*f)(x_forward);
00488 const Vector<double> y_forward_2 = (t.*f)(x_forward_2);
00489
00490 return((y_forward_2 - y_forward*2.0 + y)/(h*h));
00491 }
00492
00493
00494
00495
00500
00501 template<class T>
00502 Vector<double> calculate_central_differences_second_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00503 {
00504 const Vector<double> h = calculate_h(x);
00505
00506 const Vector<double> x_forward = x + h;
00507 const Vector<double> x_forward_2 = x + h*2.0;
00508
00509 const Vector<double> x_backward = x - h;
00510 const Vector<double> x_backward_2 = x - h*2.0;
00511
00512 const Vector<double> y = (t.*f)(x);
00513
00514 const Vector<double> y_forward = (t.*f)(x_forward);
00515 const Vector<double> y_forward_2 = (t.*f)(x_forward_2);
00516
00517 const Vector<double> y_backward = (t.*f)(x_backward);
00518 const Vector<double> y_backward_2 = (t.*f)(x_backward_2);
00519
00520 return((y_forward_2*-1.0 + y_forward*16.0 + y*-30.0 + y_backward*16.0 + y_backward_2*-1.0)/(h*h*12.0));
00521 }
00522
00523
00524
00525
00530
00531 template<class T>
00532 Vector<double> calculate_second_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00533 {
00534 switch(numerical_differentiation_method)
00535 {
00536 case ForwardDifferences:
00537 {
00538 return(calculate_forward_differences_second_derivative(t, f, x));
00539 }
00540 break;
00541
00542 case CentralDifferences:
00543 {
00544 return(calculate_central_differences_second_derivative(t, f, x));
00545 }
00546 break;
00547
00548 default:
00549 {
00550 std::ostringstream buffer;
00551
00552 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
00553 << "Vector<double> calculate_second_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
00554 << "Unknown numerical differentiation method.\n";
00555
00556 throw std::logic_error(buffer.str());
00557 }
00558 break;
00559 }
00560 }
00561
00562
00563
00564
00571
00572 template<class T>
00573 Vector<double> calculate_forward_differences_second_derivative(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
00574 {
00575 const Vector<double> y = (t.*f)(dummy, x);
00576
00577 const Vector<double> h = calculate_h(x);
00578
00579 const Vector<double> x_forward = x + h;
00580 const Vector<double> x_forward_2 = x + h*2.0;
00581
00582 const Vector<double> y_forward = (t.*f)(dummy, x_forward);
00583 const Vector<double> y_forward_2 = (t.*f)(dummy, x_forward_2);
00584
00585 return((y_forward_2 - y_forward*2.0 + y)/(h*h));
00586 }
00587
00588
00589
00590
00597
00598 template<class T>
00599 Vector<double> calculate_central_differences_second_derivative(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
00600 {
00601 const Vector<double> h = calculate_h(x);
00602
00603 const Vector<double> x_forward = x + h;
00604 const Vector<double> x_forward_2 = x + h*2.0;
00605
00606 const Vector<double> x_backward = x - h;
00607 const Vector<double> x_backward_2 = x - h*2.0;
00608
00609 const Vector<double> y = (t.*f)(dummy, x);
00610
00611 const Vector<double> y_forward = (t.*f)(dummy, x_forward);
00612 const Vector<double> y_forward_2 = (t.*f)(dummy, x_forward_2);
00613
00614 const Vector<double> y_backward = (t.*f)(dummy, x_backward);
00615 const Vector<double> y_backward_2 = (t.*f)(dummy, x_backward_2);
00616
00617 return((y_forward_2*-1.0 + y_forward*16.0 + y*-30.0 + y_backward*16.0 + y_backward_2*-1.0)/(h*h*12.0));
00618 }
00619
00620
00621
00622
00629
00630 template<class T>
00631 Vector<double> calculate_second_derivative(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
00632 {
00633 switch(numerical_differentiation_method)
00634 {
00635 case ForwardDifferences:
00636 {
00637 return(calculate_forward_differences_second_derivative(t, f, dummy, x));
00638 }
00639 break;
00640
00641 case CentralDifferences:
00642 {
00643 return(calculate_central_differences_second_derivative(t, f, dummy, x));
00644 }
00645 break;
00646
00647 default:
00648 {
00649 std::ostringstream buffer;
00650
00651 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
00652 << "Vector<double> calculate_second_derivative(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int&, const Vector<double>&) const method.\n"
00653 << "Unknown numerical differentiation method.\n";
00654
00655 throw std::logic_error(buffer.str());
00656 }
00657 break;
00658 }
00659 }
00660
00661
00662
00663
00664
00665
00671
00672 template<class T>
00673 Vector<double> calculate_forward_differences_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00674 {
00675 const unsigned int n = x.size();
00676
00677 double h;
00678
00679 double y = (t.*f)(x);
00680
00681 Vector<double> x_forward(x);
00682
00683 double y_forward;
00684
00685 Vector<double> g(n);
00686
00687 for(unsigned int i = 0; i < n; i++)
00688 {
00689 h = calculate_h(x[i]);
00690
00691 x_forward[i] += h;
00692 y_forward = (t.*f)(x_forward);
00693 x_forward[i] -= h;
00694
00695 g[i] = (y_forward - y)/h;
00696 }
00697
00698 return(g);
00699 }
00700
00701
00702
00703
00709
00710 template<class T>
00711 Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00712 {
00713 const unsigned int n = x.size();
00714
00715 double h;
00716
00717 Vector<double> x_forward(x);
00718 Vector<double> x_backward(x);
00719
00720 double y_forward;
00721 double y_backward;
00722
00723 Vector<double> g(n);
00724
00725 for(unsigned int i = 0; i < n; i++)
00726 {
00727 h = calculate_h(x[i]);
00728
00729 x_forward[i] += h;
00730 y_forward = (t.*f)(x_forward);
00731 x_forward[i] -= h;
00732
00733 x_backward[i] -= h;
00734 y_backward = (t.*f)(x_backward);
00735 x_backward[i] += h;
00736
00737 g[i] = (y_forward - y_backward)/(2.0*h);
00738 }
00739
00740 return(g);
00741 }
00742
00743
00744
00745
00751
00752 template<class T>
00753 Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
00754 {
00755 switch(numerical_differentiation_method)
00756 {
00757 case ForwardDifferences:
00758 {
00759 return(calculate_forward_differences_gradient(t, f, x));
00760 }
00761 break;
00762
00763 case CentralDifferences:
00764 {
00765 return(calculate_central_differences_gradient(t, f, x));
00766 }
00767 break;
00768
00769 default:
00770 {
00771 std::ostringstream buffer;
00772
00773 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
00774 << "Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
00775 << "Unknown numerical differentiation method.\n";
00776
00777 throw std::logic_error(buffer.str());
00778 }
00779 break;
00780 }
00781 }
00782
00783
00784
00785
00791
00792 template<class T>
00793 Vector<double> calculate_forward_differences_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>& x) const
00794 {
00795 const unsigned int n = x.size();
00796
00797 double h;
00798
00799 double y = (t.*f)(x);
00800
00801 Vector<double> x_forward(x);
00802
00803 double y_forward;
00804
00805 Vector<double> g(n);
00806
00807 for(unsigned int i = 0; i < n; i++)
00808 {
00809 h = calculate_h(x[i]);
00810
00811 x_forward[i] += h;
00812 y_forward = (t.*f)(x_forward);
00813 x_forward[i] -= h;
00814
00815 g[i] = (y_forward - y)/h;
00816 }
00817
00818 return(g);
00819 }
00820
00821
00822
00823
00829
00830 template<class T>
00831 Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>& x) const
00832 {
00833 const unsigned int n = x.size();
00834
00835 double h;
00836
00837 Vector<double> x_forward(x);
00838 Vector<double> x_backward(x);
00839
00840 double y_forward;
00841 double y_backward;
00842
00843 Vector<double> g(n);
00844
00845 for(unsigned int i = 0; i < n; i++)
00846 {
00847 h = calculate_h(x[i]);
00848
00849 x_forward[i] += h;
00850 y_forward = (t.*f)(x_forward);
00851 x_forward[i] -= h;
00852
00853 x_backward[i] -= h;
00854 y_backward = (t.*f)(x_backward);
00855 x_backward[i] += h;
00856
00857 g[i] = (y_forward - y_backward)/(2.0*h);
00858 }
00859
00860 return(g);
00861 }
00862
00863
00864
00865
00871
00872 template<class T>
00873 Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>& x) const
00874 {
00875 switch(numerical_differentiation_method)
00876 {
00877 case ForwardDifferences:
00878 {
00879 return(calculate_forward_differences_gradient(t, f, x));
00880 }
00881 break;
00882
00883 case CentralDifferences:
00884 {
00885 return(calculate_central_differences_gradient(t, f, x));
00886 }
00887 break;
00888
00889 default:
00890 {
00891 std::ostringstream buffer;
00892
00893 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
00894 << "Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>&) const method.\n"
00895 << "Unknown numerical differentiation method.\n";
00896
00897 throw std::logic_error(buffer.str());
00898 }
00899 break;
00900 }
00901 }
00902
00903
00904
00905
00913
00914 template<class T>
00915 Vector<double> calculate_forward_differences_gradient(const T& t, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
00916 {
00917 const unsigned int n = x.size();
00918
00919 double h;
00920
00921 const double y = (t.*f)(dummy, x);
00922
00923 Vector<double> x_forward(x);
00924
00925 double y_forward;
00926
00927 Vector<double> g(n);
00928
00929 for(unsigned int i = 0; i < n; i++)
00930 {
00931 h = calculate_h(x[i]);
00932
00933 x_forward[i] += h;
00934 y_forward = (t.*f)(dummy, x_forward);
00935 x_forward[i] -= h;
00936
00937 g[i] = (y_forward - y)/h;
00938 }
00939
00940 return(g);
00941 }
00942
00943
00944
00945
00953
00954 template<class T>
00955 Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
00956 {
00957 const unsigned int n = x.size();
00958
00959 double h;
00960
00961 Vector<double> x_forward(x);
00962 Vector<double> x_backward(x);
00963
00964 double y_forward;
00965 double y_backward;
00966
00967 Vector<double> g(n);
00968
00969 for(unsigned int i = 0; i < n; i++)
00970 {
00971 h = calculate_h(x[i]);
00972
00973 x_forward[i] += h;
00974 y_forward = (t.*f)(dummy, x_forward);
00975 x_forward[i] -= h;
00976
00977 x_backward[i] -= h;
00978 y_backward = (t.*f)(dummy, x_backward);
00979 x_backward[i] += h;
00980
00981 g[i] = (y_forward - y_backward)/(2.0*h);
00982 }
00983
00984 return(g);
00985 }
00986
00987
00988
00989
00997
00998 template<class T>
00999 Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
01000 {
01001 switch(numerical_differentiation_method)
01002 {
01003 case ForwardDifferences:
01004 {
01005 return(calculate_forward_differences_gradient(t, f, dummy, x));
01006 }
01007 break;
01008
01009 case CentralDifferences:
01010 {
01011 return(calculate_central_differences_gradient(t, f, dummy, x));
01012 }
01013 break;
01014
01015 default:
01016 {
01017 std::ostringstream buffer;
01018
01019 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
01020 << "Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method.\n"
01021 << "Unknown numerical differentiation method.\n";
01022
01023 throw std::logic_error(buffer.str());
01024 }
01025 break;
01026 }
01027 }
01028
01029
01030
01031
01039
01040 template<class T>
01041 Vector<double> calculate_forward_differences_gradient(const T& t, double (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
01042 {
01043 unsigned int n = x.size();
01044
01045 double h;
01046
01047 double y = (t.*f)(dummy, x);
01048
01049 Vector<double> x_forward(x);
01050
01051 double y_forward;
01052
01053 Vector<double> g(n);
01054
01055 for(unsigned int i = 0; i < n; i++)
01056 {
01057 h = calculate_h(x[i]);
01058
01059 x_forward[i] += h;
01060 y_forward = (t.*f)(dummy, x_forward);
01061 x_forward[i] -= h;
01062
01063 g[i] = (y_forward - y)/h;
01064 }
01065
01066 return(g);
01067 }
01068
01069
01070
01071
01079
01080 template<class T>
01081 Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
01082 {
01083 unsigned int n = x.size();
01084
01085 double h;
01086
01087 Vector<double> x_forward(x);
01088 Vector<double> x_backward(x);
01089
01090 double y_forward;
01091 double y_backward;
01092
01093 Vector<double> g(n);
01094
01095 for(unsigned int i = 0; i < n; i++)
01096 {
01097 h = calculate_h(x[i]);
01098
01099 x_forward[i] += h;
01100 y_forward = (t.*f)(dummy, x_forward);
01101 x_forward[i] -= h;
01102
01103 x_backward[i] -= h;
01104 y_backward = (t.*f)(dummy, x_backward);
01105 x_backward[i] += h;
01106
01107 g[i] = (y_forward - y_backward)/(2.0*h);
01108 }
01109
01110 return(g);
01111 }
01112
01113
01114
01115
01123
01124 template<class T>
01125 Vector<double> calculate_gradient(const T& t, double (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
01126 {
01127 switch(numerical_differentiation_method)
01128 {
01129 case ForwardDifferences:
01130 {
01131 return(calculate_forward_differences_gradient(t, f, dummy, x));
01132 }
01133 break;
01134
01135 case CentralDifferences:
01136 {
01137 return(calculate_central_differences_gradient(t, f, dummy, x));
01138 }
01139 break;
01140
01141 default:
01142 {
01143 std::ostringstream buffer;
01144
01145 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
01146 << "Vector<double> calculate_gradient(const T& t, double (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int&, const Vector<double>&) const method.\n"
01147 << "Unknown numerical differentiation method.\n";
01148
01149 throw std::logic_error(buffer.str());
01150 }
01151 break;
01152 }
01153 }
01154
01155
01156
01157
01158
01159
01160
01166
01167 template<class T>
01168 Matrix<double> calculate_forward_differences_Hessian(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
01169 {
01170 unsigned int n = x.size();
01171
01172 Matrix<double> H(n, n);
01173
01174 double h_i;
01175 double h_j;
01176
01177 double y = (t.*f)(x);
01178
01179 Vector<double> x_forward_2i(x);
01180 Vector<double> x_forward_ij(x);
01181 Vector<double> x_forward_i(x);
01182 Vector<double> x_forward_j(x);
01183
01184 double y_forward_2i;
01185 double y_forward_ij;
01186 double y_forward_i;
01187 double y_forward_j;
01188
01189 for(unsigned int i = 0; i < n; i++)
01190 {
01191 h_i = calculate_h(x[i]);
01192
01193 x_forward_i[i] += h_i;
01194 y_forward_i = (t.*f)(x_forward_i);
01195 x_forward_i[i] -= h_i;
01196
01197 x_forward_2i[i] += 2.0*h_i;
01198 y_forward_2i = (t.*f)(x_forward_2i);
01199 x_forward_2i[i] -= 2.0*h_i;
01200
01201 H[i][i] = (y_forward_2i - 2*y_forward_i + y)/pow(h_i, 2);
01202
01203 for(unsigned int j = i; j < n; j++)
01204 {
01205 h_j = calculate_h(x[j]);
01206
01207 x_forward_j[j] += h_j;
01208 y_forward_j = (t.*f)(x_forward_j);
01209 x_forward_j[j] -= h_j;
01210
01211 x_forward_ij[i] += h_i;
01212 x_forward_ij[j] += h_j;
01213 y_forward_ij = (t.*f)(x_forward_ij);
01214 x_forward_ij[i] -= h_i;
01215 x_forward_ij[j] -= h_j;
01216
01217 H[i][j] = (y_forward_ij - y_forward_i - y_forward_j + y)/(h_i*h_j);
01218 }
01219 }
01220
01221 for(unsigned int i = 0; i < n; i++)
01222 {
01223 for(unsigned int j = 0; j < i; j++)
01224 {
01225 H[i][j] = H[j][i];
01226 }
01227 }
01228
01229 return(H);
01230 }
01231
01232
01233
01234
01240
01241 template<class T>
01242 Matrix<double> calculate_central_differences_Hessian(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
01243 {
01244 unsigned int n = x.size();
01245
01246 double y = (t.*f)(x);
01247
01248 Matrix<double> H(n, n);
01249
01250 double h_i;
01251 double h_j;
01252
01253 Vector<double> x_backward_2i(x);
01254 Vector<double> x_backward_i(x);
01255
01256 Vector<double> x_forward_i(x);
01257 Vector<double> x_forward_2i(x);
01258
01259 Vector<double> x_backward_ij(x);
01260 Vector<double> x_forward_ij(x);
01261
01262 Vector<double> x_backward_i_forward_j(x);
01263 Vector<double> x_forward_i_backward_j(x);
01264
01265 double y_backward_2i;
01266 double y_backward_i;
01267
01268 double y_forward_i;
01269 double y_forward_2i;
01270
01271 double y_backward_ij;
01272 double y_forward_ij;
01273
01274 double y_backward_i_forward_j;
01275 double y_forward_i_backward_j;
01276
01277 for(unsigned int i = 0; i < n; i++)
01278 {
01279 h_i = calculate_h(x[i]);
01280
01281 x_backward_2i[i] -= 2.0*h_i;
01282 y_backward_2i = (t.*f)(x_backward_2i);
01283 x_backward_2i[i] += 2.0*h_i;
01284
01285 x_backward_i[i] -= h_i;
01286 y_backward_i = (t.*f)(x_backward_i);
01287 x_backward_i[i] += h_i;
01288
01289 x_forward_i[i] += h_i;
01290 y_forward_i = (t.*f)(x_forward_i);
01291 x_forward_i[i] -= h_i;
01292
01293 x_forward_2i[i] += 2.0*h_i;
01294 y_forward_2i = (t.*f)(x_forward_2i);
01295 x_forward_2i[i] -= 2.0*h_i;
01296
01297 H[i][i] = (-y_forward_2i + 16.0*y_forward_i -30.0*y + 16.0*y_backward_i - y_backward_2i)/(12.0*pow(h_i, 2));
01298
01299 for(unsigned int j = i; j < n; j++)
01300 {
01301 h_j = calculate_h(x[j]);
01302
01303 x_backward_ij[i] -= h_i;
01304 x_backward_ij[j] -= h_j;
01305 y_backward_ij = (t.*f)(x_backward_ij);
01306 x_backward_ij[i] += h_i;
01307 x_backward_ij[j] += h_j;
01308
01309 x_forward_ij[i] += h_i;
01310 x_forward_ij[j] += h_j;
01311 y_forward_ij = (t.*f)(x_forward_ij);
01312 x_forward_ij[i] -= h_i;
01313 x_forward_ij[j] -= h_j;
01314
01315 x_backward_i_forward_j[i] -= h_i;
01316 x_backward_i_forward_j[j] += h_j;
01317 y_backward_i_forward_j = (t.*f)(x_backward_i_forward_j);
01318 x_backward_i_forward_j[i] += h_i;
01319 x_backward_i_forward_j[j] -= h_j;
01320
01321 x_forward_i_backward_j[i] += h_i;
01322 x_forward_i_backward_j[j] -= h_j;
01323 y_forward_i_backward_j = (t.*f)(x_forward_i_backward_j);
01324 x_forward_i_backward_j[i] -= h_i;
01325 x_forward_i_backward_j[j] += h_j;
01326
01327 H[i][j] = (y_forward_ij - y_forward_i_backward_j - y_backward_i_forward_j + y_backward_ij)/(4.0*h_i*h_j);
01328 }
01329 }
01330
01331 for(unsigned int i = 0; i < n; i++)
01332 {
01333 for(unsigned int j = 0; j < i; j++)
01334 {
01335 H[i][j] = H[j][i];
01336 }
01337 }
01338
01339 return(H);
01340 }
01341
01342
01343
01344
01350
01351 template<class T>
01352 Matrix<double> calculate_Hessian(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
01353 {
01354 switch(numerical_differentiation_method)
01355 {
01356 case ForwardDifferences:
01357 {
01358 return(calculate_forward_differences_Hessian(t, f, x));
01359 }
01360 break;
01361
01362 case CentralDifferences:
01363 {
01364 return(calculate_central_differences_Hessian(t, f, x));
01365 }
01366 break;
01367
01368 default:
01369 {
01370 std::ostringstream buffer;
01371
01372 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
01373 << "double calculate_Hessian(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
01374 << "Unknown numerical differentiation method.\n";
01375
01376 throw std::logic_error(buffer.str());
01377 }
01378 break;
01379 }
01380 }
01381
01382
01383
01384
01392
01393 template<class T>
01394 Matrix<double> calculate_forward_differences_Hessian(const T& t, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
01395 {
01396 unsigned int n = x.size();
01397
01398 Matrix<double> H(n, n);
01399
01400 double h_i;
01401 double h_j;
01402
01403 double y = (t.*f)(dummy, x);
01404
01405 Vector<double> x_forward_2i(x);
01406 Vector<double> x_forward_ij(x);
01407 Vector<double> x_forward_i(x);
01408 Vector<double> x_forward_j(x);
01409
01410 double y_forward_2i;
01411 double y_forward_ij;
01412 double y_forward_i;
01413 double y_forward_j;
01414
01415 for(unsigned int i = 0; i < n; i++)
01416 {
01417 h_i = calculate_h(x[i]);
01418
01419 x_forward_i[i] += h_i;
01420 y_forward_i = (t.*f)(dummy, x_forward_i);
01421 x_forward_i[i] -= h_i;
01422
01423 x_forward_2i[i] += 2.0*h_i;
01424 y_forward_2i = (t.*f)(dummy, x_forward_2i);
01425 x_forward_2i[i] -= 2.0*h_i;
01426
01427 H[i][i] = (y_forward_2i - 2*y_forward_i + y)/pow(h_i, 2);
01428
01429 for(unsigned int j = i; j < n; j++)
01430 {
01431 h_j = calculate_h(x[j]);
01432
01433 x_forward_j[j] += h_j;
01434 y_forward_j = (t.*f)(dummy, x_forward_j);
01435 x_forward_j[j] -= h_j;
01436
01437 x_forward_ij[i] += h_i;
01438 x_forward_ij[j] += h_j;
01439 y_forward_ij = (t.*f)(dummy, x_forward_ij);
01440 x_forward_ij[i] -= h_i;
01441 x_forward_ij[j] -= h_j;
01442
01443 H[i][j] = (y_forward_ij - y_forward_i - y_forward_j + y)/(h_i*h_j);
01444 }
01445 }
01446
01447 for(unsigned int i = 0; i < n; i++)
01448 {
01449 for(unsigned int j = 0; j < i; j++)
01450 {
01451 H[i][j] = H[j][i];
01452 }
01453 }
01454
01455 return(H);
01456 }
01457
01458
01459
01460
01468
01469 template<class T>
01470 Matrix<double> calculate_central_differences_Hessian(const T& t, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
01471 {
01472 unsigned int n = x.size();
01473
01474 double y = (t.*f)(dummy, x);
01475
01476 Matrix<double> H(n, n);
01477
01478 double h_i;
01479 double h_j;
01480
01481 Vector<double> x_backward_2i(x);
01482 Vector<double> x_backward_i(x);
01483
01484 Vector<double> x_forward_i(x);
01485 Vector<double> x_forward_2i(x);
01486
01487 Vector<double> x_backward_ij(x);
01488 Vector<double> x_forward_ij(x);
01489
01490 Vector<double> x_backward_i_forward_j(x);
01491 Vector<double> x_forward_i_backward_j(x);
01492
01493 double y_backward_2i;
01494 double y_backward_i;
01495
01496 double y_forward_i;
01497 double y_forward_2i;
01498
01499 double y_backward_ij;
01500 double y_forward_ij;
01501
01502 double y_backward_i_forward_j;
01503 double y_forward_i_backward_j;
01504
01505 for(unsigned int i = 0; i < n; i++)
01506 {
01507 h_i = calculate_h(x[i]);
01508
01509 x_backward_2i[i] -= 2.0*h_i;
01510 y_backward_2i = (t.*f)(dummy, x_backward_2i);
01511 x_backward_2i[i] += 2.0*h_i;
01512
01513 x_backward_i[i] -= h_i;
01514 y_backward_i = (t.*f)(dummy, x_backward_i);
01515 x_backward_i[i] += h_i;
01516
01517 x_forward_i[i] += h_i;
01518 y_forward_i = (t.*f)(dummy, x_forward_i);
01519 x_forward_i[i] -= h_i;
01520
01521 x_forward_2i[i] += 2.0*h_i;
01522 y_forward_2i = (t.*f)(dummy, x_forward_2i);
01523 x_forward_2i[i] -= 2.0*h_i;
01524
01525 H[i][i] = (-y_forward_2i + 16.0*y_forward_i -30.0*y + 16.0*y_backward_i - y_backward_2i)/(12.0*pow(h_i, 2));
01526
01527 for(unsigned int j = i; j < n; j++)
01528 {
01529 h_j = calculate_h(x[j]);
01530
01531 x_backward_ij[i] -= h_i;
01532 x_backward_ij[j] -= h_j;
01533 y_backward_ij = (t.*f)(dummy, x_backward_ij);
01534 x_backward_ij[i] += h_i;
01535 x_backward_ij[j] += h_j;
01536
01537 x_forward_ij[i] += h_i;
01538 x_forward_ij[j] += h_j;
01539 y_forward_ij = (t.*f)(dummy, x_forward_ij);
01540 x_forward_ij[i] -= h_i;
01541 x_forward_ij[j] -= h_j;
01542
01543 x_backward_i_forward_j[i] -= h_i;
01544 x_backward_i_forward_j[j] += h_j;
01545 y_backward_i_forward_j = (t.*f)(dummy, x_backward_i_forward_j);
01546 x_backward_i_forward_j[i] += h_i;
01547 x_backward_i_forward_j[j] -= h_j;
01548
01549 x_forward_i_backward_j[i] += h_i;
01550 x_forward_i_backward_j[j] -= h_j;
01551 y_forward_i_backward_j = (t.*f)(dummy, x_forward_i_backward_j);
01552 x_forward_i_backward_j[i] -= h_i;
01553 x_forward_i_backward_j[j] += h_j;
01554
01555 H[i][j] = (y_forward_ij - y_forward_i_backward_j - y_backward_i_forward_j + y_backward_ij)/(4.0*h_i*h_j);
01556 }
01557 }
01558
01559 for(unsigned int i = 0; i < n; i++)
01560 {
01561 for(unsigned int j = 0; j < i; j++)
01562 {
01563 H[i][j] = H[j][i];
01564 }
01565 }
01566
01567 return(H);
01568 }
01569
01570
01571
01572
01580
01581 template<class T>
01582 Matrix<double> calculate_Hessian(const T& t, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
01583 {
01584 switch(numerical_differentiation_method)
01585 {
01586 case ForwardDifferences:
01587 {
01588 return(calculate_forward_differences_Hessian(t, f, dummy, x));
01589 }
01590 break;
01591
01592 case CentralDifferences:
01593 {
01594 return(calculate_central_differences_Hessian(t, f, dummy, x));
01595 }
01596 break;
01597
01598 default:
01599 {
01600 std::ostringstream buffer;
01601
01602 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
01603 << "double calculate_Hessian(const T& t, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method.\n"
01604 << "Unknown numerical differentiation method.\n";
01605
01606 throw std::logic_error(buffer.str());
01607 }
01608 break;
01609 }
01610 }
01611
01612
01613
01614
01622
01623 template<class T>
01624 Matrix<double> calculate_forward_differences_Hessian(const T& t, double (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
01625 {
01626 unsigned int n = x.size();
01627
01628 Matrix<double> H(n, n);
01629
01630 double h_i;
01631 double h_j;
01632
01633 double y = (t.*f)(dummy, x);
01634
01635 Vector<double> x_forward_2i(x);
01636 Vector<double> x_forward_ij(x);
01637 Vector<double> x_forward_i(x);
01638 Vector<double> x_forward_j(x);
01639
01640 double y_forward_2i;
01641 double y_forward_ij;
01642 double y_forward_i;
01643 double y_forward_j;
01644
01645 for(unsigned int i = 0; i < n; i++)
01646 {
01647 h_i = calculate_h(x[i]);
01648
01649 x_forward_i[i] += h_i;
01650 y_forward_i = (t.*f)(dummy, x_forward_i);
01651 x_forward_i[i] -= h_i;
01652
01653 x_forward_2i[i] += 2.0*h_i;
01654 y_forward_2i = (t.*f)(dummy, x_forward_2i);
01655 x_forward_2i[i] -= 2.0*h_i;
01656
01657 H[i][i] = (y_forward_2i - 2*y_forward_i + y)/pow(h_i, 2);
01658
01659 for(unsigned int j = i; j < n; j++)
01660 {
01661 h_j = calculate_h(x[j]);
01662
01663 x_forward_j[j] += h_j;
01664 y_forward_j = (t.*f)(dummy, x_forward_j);
01665 x_forward_j[j] -= h_j;
01666
01667 x_forward_ij[i] += h_i;
01668 x_forward_ij[j] += h_j;
01669 y_forward_ij = (t.*f)(dummy, x_forward_ij);
01670 x_forward_ij[i] -= h_i;
01671 x_forward_ij[j] -= h_j;
01672
01673 H[i][j] = (y_forward_ij - y_forward_i - y_forward_j + y)/(h_i*h_j);
01674 }
01675 }
01676
01677 for(unsigned int i = 0; i < n; i++)
01678 {
01679 for(unsigned int j = 0; j < i; j++)
01680 {
01681 H[i][j] = H[j][i];
01682 }
01683 }
01684
01685 return(H);
01686 }
01687
01688
01689
01690
01698
01699 template<class T>
01700 Matrix<double> calculate_central_differences_Hessian(const T& t, double (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
01701 {
01702 unsigned int n = x.size();
01703
01704 double y = (t.*f)(dummy, x);
01705
01706 Matrix<double> H(n, n);
01707
01708 double h_i;
01709 double h_j;
01710
01711 Vector<double> x_backward_2i(x);
01712 Vector<double> x_backward_i(x);
01713
01714 Vector<double> x_forward_i(x);
01715 Vector<double> x_forward_2i(x);
01716
01717 Vector<double> x_backward_ij(x);
01718 Vector<double> x_forward_ij(x);
01719
01720 Vector<double> x_backward_i_forward_j(x);
01721 Vector<double> x_forward_i_backward_j(x);
01722
01723 double y_backward_2i;
01724 double y_backward_i;
01725
01726 double y_forward_i;
01727 double y_forward_2i;
01728
01729 double y_backward_ij;
01730 double y_forward_ij;
01731
01732 double y_backward_i_forward_j;
01733 double y_forward_i_backward_j;
01734
01735 for(unsigned int i = 0; i < n; i++)
01736 {
01737 h_i = calculate_h(x[i]);
01738
01739 x_backward_2i[i] -= 2.0*h_i;
01740 y_backward_2i = (t.*f)(dummy, x_backward_2i);
01741 x_backward_2i[i] += 2.0*h_i;
01742
01743 x_backward_i[i] -= h_i;
01744 y_backward_i = (t.*f)(dummy, x_backward_i);
01745 x_backward_i[i] += h_i;
01746
01747 x_forward_i[i] += h_i;
01748 y_forward_i = (t.*f)(dummy, x_forward_i);
01749 x_forward_i[i] -= h_i;
01750
01751 x_forward_2i[i] += 2.0*h_i;
01752 y_forward_2i = (t.*f)(dummy, x_forward_2i);
01753 x_forward_2i[i] -= 2.0*h_i;
01754
01755 H[i][i] = (-y_forward_2i + 16.0*y_forward_i -30.0*y + 16.0*y_backward_i - y_backward_2i)/(12.0*pow(h_i, 2));
01756
01757 for(unsigned int j = i; j < n; j++)
01758 {
01759 h_j = calculate_h(x[j]);
01760
01761 x_backward_ij[i] -= h_i;
01762 x_backward_ij[j] -= h_j;
01763 y_backward_ij = (t.*f)(dummy, x_backward_ij);
01764 x_backward_ij[i] += h_i;
01765 x_backward_ij[j] += h_j;
01766
01767 x_forward_ij[i] += h_i;
01768 x_forward_ij[j] += h_j;
01769 y_forward_ij = (t.*f)(dummy, x_forward_ij);
01770 x_forward_ij[i] -= h_i;
01771 x_forward_ij[j] -= h_j;
01772
01773 x_backward_i_forward_j[i] -= h_i;
01774 x_backward_i_forward_j[j] += h_j;
01775 y_backward_i_forward_j = (t.*f)(dummy, x_backward_i_forward_j);
01776 x_backward_i_forward_j[i] += h_i;
01777 x_backward_i_forward_j[j] -= h_j;
01778
01779 x_forward_i_backward_j[i] += h_i;
01780 x_forward_i_backward_j[j] -= h_j;
01781 y_forward_i_backward_j = (t.*f)(dummy, x_forward_i_backward_j);
01782 x_forward_i_backward_j[i] -= h_i;
01783 x_forward_i_backward_j[j] += h_j;
01784
01785 H[i][j] = (y_forward_ij - y_forward_i_backward_j - y_backward_i_forward_j + y_backward_ij)/(4.0*h_i*h_j);
01786 }
01787 }
01788
01789 for(unsigned int i = 0; i < n; i++)
01790 {
01791 for(unsigned int j = 0; j < i; j++)
01792 {
01793 H[i][j] = H[j][i];
01794 }
01795 }
01796
01797 return(H);
01798 }
01799
01800
01801
01802
01810
01811 template<class T>
01812 Matrix<double> calculate_Hessian(const T& t, double (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
01813 {
01814 switch(numerical_differentiation_method)
01815 {
01816 case ForwardDifferences:
01817 {
01818 return(calculate_forward_differences_Hessian(t, f, dummy, x));
01819 }
01820 break;
01821
01822 case CentralDifferences:
01823 {
01824 return(calculate_central_differences_Hessian(t, f, dummy, x));
01825 }
01826 break;
01827
01828 default:
01829 {
01830 std::ostringstream buffer;
01831
01832 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
01833 << "double calculate_Hessian(const T& t, double (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int&, const Vector<double>&) const method.\n"
01834 << "Unknown numerical differentiation method.\n";
01835
01836 throw std::logic_error(buffer.str());
01837 }
01838 break;
01839 }
01840 }
01841
01842
01843
01844
01845
01846
01852
01853 template<class T>
01854 Matrix<double> calculate_forward_differences_Jacobian(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
01855 {
01856 Vector<double> y = (t.*f)(x);
01857
01858 double h;
01859
01860 unsigned int n = x.size();
01861 unsigned int m = y.size();
01862
01863 Vector<double> x_forward(x);
01864 Vector<double> y_forward(n);
01865
01866 Matrix<double> J(m,n);
01867
01868 for(unsigned int j = 0; j < n; j++)
01869 {
01870 h = calculate_h(x[j]);
01871
01872 x_forward[j] += h;
01873 y_forward = (t.*f)(x_forward);
01874 x_forward[j] -= h;
01875
01876 for(unsigned int i = 0; i < m; i++)
01877 {
01878 J[i][j] = (y_forward[i] - y[i])/h;
01879 }
01880 }
01881
01882 return(J);
01883 }
01884
01885
01886
01887
01893
01894 template<class T>
01895 Matrix<double> calculate_central_differences_Jacobian(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
01896 {
01897 Vector<double> y = (t.*f)(x);
01898
01899 double h;
01900
01901 unsigned int n = x.size();
01902 unsigned int m = y.size();
01903
01904 Vector<double> x_forward(x);
01905 Vector<double> x_backward(x);
01906
01907 Vector<double> y_forward(n);
01908 Vector<double> y_backward(n);
01909
01910 Matrix<double> J(m,n);
01911
01912 for(unsigned int j = 0; j < n; j++)
01913 {
01914 h = calculate_h(x[j]);
01915
01916 x_backward[j] -= h;
01917 y_backward = (t.*f)(x_backward);
01918 x_backward[j] += h;
01919
01920 x_forward[j] += h;
01921 y_forward = (t.*f)(x_forward);
01922 x_forward[j] -= h;
01923
01924 for(unsigned int i = 0; i < m; i++)
01925 {
01926 J[i][j] = (y_forward[i] - y_backward[i])/(2.0*h);
01927 }
01928 }
01929
01930 return(J);
01931 }
01932
01933
01934
01935
01941
01942 template<class T>
01943 Matrix<double> calculate_Jacobian(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
01944 {
01945 switch(numerical_differentiation_method)
01946 {
01947 case ForwardDifferences:
01948 {
01949 return(calculate_forward_differences_Jacobian(t, f, x));
01950 }
01951 break;
01952
01953 case CentralDifferences:
01954 {
01955 return(calculate_central_differences_Jacobian(t, f, x));
01956 }
01957 break;
01958
01959 default:
01960 {
01961 std::ostringstream buffer;
01962
01963 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
01964 << "Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
01965 << "Unknown numerical differentiation method.\n";
01966
01967 throw std::logic_error(buffer.str());
01968 }
01969 break;
01970 }
01971 }
01972
01973
01974
01975
01982
01983 template<class T>
01984 Matrix<double> calculate_forward_differences_Jacobian(const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
01985 {
01986 Vector<double> y = (t.*f)(dummy, x);
01987
01988 double h;
01989
01990 unsigned int n = x.size();
01991 unsigned int m = y.size();
01992
01993 Vector<double> x_forward(x);
01994 Vector<double> y_forward(n);
01995
01996 Matrix<double> J(m,n);
01997
01998 for(unsigned int j = 0; j < n; j++)
01999 {
02000 h = calculate_h(x[j]);
02001
02002 x_forward[j] += h;
02003 y_forward = (t.*f)(dummy, x_forward);
02004 x_forward[j] -= h;
02005
02006 for(unsigned int i = 0; i < m; i++)
02007 {
02008 J[i][j] = (y_forward[i] - y[i])/h;
02009 }
02010 }
02011
02012 return(J);
02013 }
02014
02015
02016
02017
02024
02025 template<class T>
02026 Matrix<double> calculate_central_differences_Jacobian(const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
02027 {
02028 Vector<double> y = (t.*f)(dummy, x);
02029
02030 double h;
02031
02032 unsigned int n = x.size();
02033 unsigned int m = y.size();
02034
02035 Vector<double> x_forward(x);
02036 Vector<double> x_backward(x);
02037
02038 Vector<double> y_forward(n);
02039 Vector<double> y_backward(n);
02040
02041 Matrix<double> J(m,n);
02042
02043 for(unsigned int j = 0; j < n; j++)
02044 {
02045 h = calculate_h(x[j]);
02046
02047 x_backward[j] -= h;
02048 y_backward = (t.*f)(dummy, x_backward);
02049 x_backward[j] += h;
02050
02051 x_forward[j] += h;
02052 y_forward = (t.*f)(dummy, x_forward);
02053 x_forward[j] -= h;
02054
02055 for(unsigned int i = 0; i < m; i++)
02056 {
02057 J[i][j] = (y_forward[i] - y_backward[i])/(2.0*h);
02058 }
02059 }
02060
02061 return(J);
02062 }
02063
02064
02065
02066
02073
02074 template<class T>
02075 Matrix<double> calculate_Jacobian(const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const
02076 {
02077 switch(numerical_differentiation_method)
02078 {
02079 case ForwardDifferences:
02080 {
02081 return(calculate_forward_differences_Jacobian(t, f, dummy, x));
02082 }
02083 break;
02084
02085 case CentralDifferences:
02086 {
02087 return(calculate_central_differences_Jacobian(t, f, dummy, x));
02088 }
02089 break;
02090
02091 default:
02092 {
02093 std::ostringstream buffer;
02094
02095 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
02096 << "Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method.\n"
02097 << "Unknown numerical differentiation method.\n";
02098
02099 throw std::logic_error(buffer.str());
02100 }
02101 break;
02102 }
02103 }
02104
02105
02106
02107
02108
02116
02117 template<class T>
02118 Matrix<double> calculate_forward_differences_Jacobian(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
02119 {
02120 Vector<double> y = (t.*f)(dummy, x);
02121
02122 double h;
02123
02124 unsigned int n = x.size();
02125 unsigned int m = y.size();
02126
02127 Vector<double> x_forward(x);
02128 Vector<double> y_forward(n);
02129
02130 Matrix<double> J(m,n);
02131
02132 for(unsigned int j = 0; j < n; j++)
02133 {
02134 h = calculate_h(x[j]);
02135
02136 x_forward[j] += h;
02137 y_forward = (t.*f)(dummy, x_forward);
02138 x_forward[j] -= h;
02139
02140 for(unsigned int i = 0; i < m; i++)
02141 {
02142 J[i][j] = (y_forward[i] - y[i])/h;
02143 }
02144 }
02145
02146 return(J);
02147 }
02148
02149
02150
02151
02159
02160 template<class T>
02161 Matrix<double> calculate_central_differences_Jacobian(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
02162 {
02163 Vector<double> y = (t.*f)(dummy, x);
02164
02165 double h;
02166
02167 unsigned int n = x.size();
02168 unsigned int m = y.size();
02169
02170 Vector<double> x_forward(x);
02171 Vector<double> x_backward(x);
02172
02173 Vector<double> y_forward(n);
02174 Vector<double> y_backward(n);
02175
02176 Matrix<double> J(m,n);
02177
02178 for(unsigned int j = 0; j < n; j++)
02179 {
02180 h = calculate_h(x[j]);
02181
02182 x_backward[j] -= h;
02183 y_backward = (t.*f)(dummy, x_backward);
02184 x_backward[j] += h;
02185
02186 x_forward[j] += h;
02187 y_forward = (t.*f)(dummy, x_forward);
02188 x_forward[j] -= h;
02189
02190 for(unsigned int i = 0; i < m; i++)
02191 {
02192 J[i][j] = (y_forward[i] - y_backward[i])/(2.0*h);
02193 }
02194 }
02195
02196 return(J);
02197 }
02198
02199
02200
02201
02209
02210 template<class T>
02211 Matrix<double> calculate_Jacobian(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
02212 {
02213 switch(numerical_differentiation_method)
02214 {
02215 case ForwardDifferences:
02216 {
02217 return(calculate_forward_differences_Jacobian(t, f, dummy, x));
02218 }
02219 break;
02220
02221 case CentralDifferences:
02222 {
02223 return(calculate_central_differences_Jacobian(t, f, dummy, x));
02224 }
02225 break;
02226
02227 default:
02228 {
02229 std::ostringstream buffer;
02230
02231 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
02232 << "Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int&, const Vector<double>&) const method.\n"
02233 << "Unknown numerical differentiation method.\n";
02234
02235 throw std::logic_error(buffer.str());
02236 }
02237 break;
02238 }
02239 }
02240
02241
02242
02243
02252
02253 template<class T>
02254 Matrix<double> calculate_forward_differences_Jacobian
02255 (const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&, const Vector<double>&) const, const unsigned int& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
02256 {
02257 Vector<double> y = (t.*f)(dummy_int, dummy_vector, x);
02258
02259 double h;
02260
02261 unsigned int n = x.size();
02262 unsigned int m = y.size();
02263
02264 Vector<double> x_forward(x);
02265 Vector<double> y_forward(n);
02266
02267 Matrix<double> J(m,n);
02268
02269 for(unsigned int j = 0; j < n; j++)
02270 {
02271 h = calculate_h(x[j]);
02272
02273 x_forward[j] += h;
02274 y_forward = (t.*f)(dummy_int, dummy_vector, x_forward);
02275 x_forward[j] -= h;
02276
02277 for(unsigned int i = 0; i < m; i++)
02278 {
02279 J[i][j] = (y_forward[i] - y[i])/h;
02280 }
02281 }
02282
02283 return(J);
02284 }
02285
02286
02287
02288
02297
02298 template<class T>
02299 Matrix<double> calculate_central_differences_Jacobian
02300 (const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&, const Vector<double>&) const, const unsigned int& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
02301 {
02302 const unsigned int n = x.size();
02303
02304 const Vector<double> y = (t.*f)(dummy_int, dummy_vector, x);
02305 const unsigned int m = y.size();
02306
02307 double h;
02308
02309 Vector<double> x_forward(x);
02310 Vector<double> x_backward(x);
02311
02312 Vector<double> y_forward(n);
02313 Vector<double> y_backward(n);
02314
02315 Matrix<double> J(m,n);
02316
02317 for(unsigned int j = 0; j < n; j++)
02318 {
02319 h = calculate_h(x[j]);
02320
02321 x_backward[j] -= h;
02322 y_backward = (t.*f)(dummy_int, dummy_vector, x_backward);
02323 x_backward[j] += h;
02324
02325 x_forward[j] += h;
02326 y_forward = (t.*f)(dummy_int, dummy_vector, x_forward);
02327 x_forward[j] -= h;
02328
02329 for(unsigned int i = 0; i < m; i++)
02330 {
02331 J[i][j] = (y_forward[i] - y_backward[i])/(2.0*h);
02332 }
02333 }
02334
02335 return(J);
02336 }
02337
02338
02339
02340
02349
02350 template<class T>
02351 Matrix<double> calculate_Jacobian
02352 (const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&, const Vector<double>&) const, const unsigned int& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
02353 {
02354 switch(numerical_differentiation_method)
02355 {
02356 case ForwardDifferences:
02357 {
02358 return(calculate_forward_differences_Jacobian(t, f, dummy_int, dummy_vector, x));
02359 }
02360 break;
02361
02362 case CentralDifferences:
02363 {
02364 return(calculate_central_differences_Jacobian(t, f, dummy_int, dummy_vector, x));
02365 }
02366 break;
02367
02368 default:
02369 {
02370 std::ostringstream buffer;
02371
02372 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
02373 << "Matrix<double> calculate_Jacobian\n"
02374 << "(const T&, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&, const Vector<double>&) const, const unsigned int&, const Vector<double>&, const Vector<double>&) const method.\n"
02375 << "Unknown numerical differentiation method.\n";
02376
02377 throw std::logic_error(buffer.str());
02378 }
02379 break;
02380 }
02381 }
02382
02383
02384
02385
02394
02395 template<class T>
02396 Matrix<double> calculate_forward_differences_Jacobian
02397 (const T& t, Vector<double> (T::*f)(const unsigned int&, const unsigned int&, const Vector<double>&) const, const unsigned int& dummy_int_1, const unsigned int& dummy_int_2, const Vector<double>& x) const
02398 {
02399 const Vector<double> y = (t.*f)(dummy_int_1, dummy_int_2, x);
02400
02401 const unsigned int n = x.size();
02402 const unsigned int m = y.size();
02403
02404 double h;
02405
02406 Vector<double> x_forward(x);
02407 Vector<double> y_forward(n);
02408
02409 Matrix<double> J(m,n);
02410
02411 for(unsigned int j = 0; j < n; j++)
02412 {
02413 h = calculate_h(x[j]);
02414
02415 x_forward[j] += h;
02416 y_forward = (t.*f)(dummy_int_1, dummy_int_2, x_forward);
02417 x_forward[j] -= h;
02418
02419 for(unsigned int i = 0; i < m; i++)
02420 {
02421 J[i][j] = (y_forward[i] - y[i])/h;
02422 }
02423 }
02424
02425 return(J);
02426 }
02427
02428
02429
02430
02439
02440 template<class T>
02441 Matrix<double> calculate_central_differences_Jacobian
02442 (const T& t, Vector<double> (T::*f)(const unsigned int&, const unsigned int&, const Vector<double>&) const, const unsigned int& dummy_int_1, const unsigned int& dummy_int_2, const Vector<double>& x) const
02443 {
02444 const Vector<double> y = (t.*f)(dummy_int_1, dummy_int_2, x);
02445
02446 const unsigned int n = x.size();
02447 const unsigned int m = y.size();
02448
02449 double h;
02450
02451 Vector<double> x_forward(x);
02452 Vector<double> x_backward(x);
02453
02454 Vector<double> y_forward(n);
02455 Vector<double> y_backward(n);
02456
02457 Matrix<double> J(m,n);
02458
02459 for(unsigned int j = 0; j < n; j++)
02460 {
02461 h = calculate_h(x[j]);
02462
02463 x_backward[j] -= h;
02464 y_backward = (t.*f)(dummy_int_1, dummy_int_2, x_backward);
02465 x_backward[j] += h;
02466
02467 x_forward[j] += h;
02468 y_forward = (t.*f)(dummy_int_1, dummy_int_2, x_forward);
02469 x_forward[j] -= h;
02470
02471 for(unsigned int i = 0; i < m; i++)
02472 {
02473 J[i][j] = (y_forward[i] - y_backward[i])/(2.0*h);
02474 }
02475 }
02476
02477 return(J);
02478 }
02479
02480
02481
02482
02491
02492 template<class T>
02493 Matrix<double> calculate_Jacobian
02494 (const T& t, Vector<double> (T::*f)(const unsigned int&, const unsigned int&, const Vector<double>&) const, const unsigned int& dummy_int_1, const unsigned int& dummy_int_2, const Vector<double>& x) const
02495 {
02496 switch(numerical_differentiation_method)
02497 {
02498 case ForwardDifferences:
02499 {
02500 return(calculate_forward_differences_Jacobian(t, f, dummy_int_1, dummy_int_2, x));
02501 }
02502 break;
02503
02504 case CentralDifferences:
02505 {
02506 return(calculate_central_differences_Jacobian(t, f, dummy_int_1, dummy_int_2, x));
02507 }
02508 break;
02509
02510 default:
02511 {
02512 std::ostringstream buffer;
02513
02514 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
02515 << "Matrix<double> calculate_Jacobian\n"
02516 << "(const T&, Vector<double> (T::*f)(const unsigned int&, const unsigned int&, const Vector<double>&) const, const unsigned int&, const unsigned int&, const Vector<double>&) const method.\n"
02517 << "Unknown numerical differentiation method.\n";
02518
02519 throw std::logic_error(buffer.str());
02520 }
02521 break;
02522 }
02523 }
02524
02525
02526
02527
02528
02529
02530
02536
02537 template<class T>
02538 Vector< Matrix <double> > calculate_forward_differences_Hessian_form(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
02539 {
02540 Vector<double> y = (t.*f)(x);
02541
02542 unsigned int s = y.size();
02543 unsigned int n = x.size();
02544
02545 double h_j;
02546 double h_k;
02547
02548 Vector<double> x_forward_j(x);
02549 Vector<double> x_forward_2j(x);
02550
02551 Vector<double> x_forward_k(x);
02552 Vector<double> x_forward_jk(x);
02553
02554 Vector<double> y_forward_j(s);
02555 Vector<double> y_forward_2j(s);
02556
02557 Vector<double> y_forward_k(s);
02558 Vector<double> y_forward_jk(s);
02559
02560 Vector< Matrix<double> > H(s);
02561
02562 for(unsigned int i = 0; i < s; i++)
02563 {
02564 H[i].set(n,n);
02565
02566 for(unsigned int j = 0; j < n; j++)
02567 {
02568 h_j = calculate_h(x[j]);
02569
02570 x_forward_j[j] += h_j;
02571 y_forward_j = (t.*f)(x_forward_j);
02572 x_forward_j[j] -= h_j;
02573
02574 x_forward_2j[j] += 2.0*h_j;
02575 y_forward_2j = (t.*f)(x_forward_2j);
02576 x_forward_2j[j] -= 2.0*h_j;
02577
02578 H[i][j][j] = (y_forward_2j[i] - 2.0*y_forward_j[i] + y[i])/pow(h_j, 2);
02579
02580 for(unsigned int k = j; k < n; k++)
02581 {
02582 h_k = calculate_h(x[k]);
02583
02584 x_forward_k[k] += h_k;
02585 y_forward_k = (t.*f)(x_forward_k);
02586 x_forward_k[k] -= h_k;
02587
02588 x_forward_jk[j] += h_j;
02589 x_forward_jk[k] += h_k;
02590 y_forward_jk = (t.*f)(x_forward_jk);
02591 x_forward_jk[j] -= h_j;
02592 x_forward_jk[k] -= h_k;
02593
02594 H[i][j][k] = (y_forward_jk[i] - y_forward_j[i] - y_forward_k[i] + y[i])/(h_j*h_k);
02595 }
02596 }
02597
02598 for(unsigned int j = 0; j < n; j++)
02599 {
02600 for(unsigned int k = 0; k < j; k++)
02601 {
02602 H[i][j][k] = H[i][k][j];
02603 }
02604 }
02605 }
02606
02607 return(H);
02608 }
02609
02610
02611
02612
02618
02619 template<class T>
02620 Vector< Matrix <double> > calculate_central_differences_Hessian_form(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
02621 {
02622 Vector<double> y = (t.*f)(x);
02623
02624 unsigned int s = y.size();
02625 unsigned int n = x.size();
02626
02627 double h_j;
02628 double h_k;
02629
02630 Vector<double> x_backward_2j(x);
02631 Vector<double> x_backward_j(x);
02632
02633 Vector<double> x_forward_j(x);
02634 Vector<double> x_forward_2j(x);
02635
02636 Vector<double> x_backward_jk(x);
02637 Vector<double> x_forward_jk(x);
02638
02639 Vector<double> x_backward_j_forward_k(x);
02640 Vector<double> x_forward_j_backward_k(x);
02641
02642 Vector<double> y_backward_2j;
02643 Vector<double> y_backward_j;
02644
02645 Vector<double> y_forward_j;
02646 Vector<double> y_forward_2j;
02647
02648 Vector<double> y_backward_jk;
02649 Vector<double> y_forward_jk;
02650
02651 Vector<double> y_backward_j_forward_k;
02652 Vector<double> y_forward_j_backward_k;
02653
02654 Vector< Matrix<double> > H(s);
02655
02656 for(unsigned int i = 0; i < s; i++)
02657 {
02658 H[i].set(n,n);
02659
02660 for(unsigned int j = 0; j < n; j++)
02661 {
02662 h_j = calculate_h(x[j]);
02663
02664 x_backward_2j[j] -= 2.0*h_j;
02665 y_backward_2j = (t.*f)(x_backward_2j);
02666 x_backward_2j[j] += 2.0*h_j;
02667
02668 x_backward_j[j] -= h_j;
02669 y_backward_j = (t.*f)(x_backward_j);
02670 x_backward_j[j] += h_j;
02671
02672 x_forward_j[j] += h_j;
02673 y_forward_j = (t.*f)(x_forward_j);
02674 x_forward_j[j] -= h_j;
02675
02676 x_forward_2j[j] += 2.0*h_j;
02677 y_forward_2j = (t.*f)(x_forward_2j);
02678 x_forward_2j[j] -= 2.0*h_j;
02679
02680 H[i][j][j] = (-y_forward_2j[i] + 16.0*y_forward_j[i] -30.0*y[i] + 16.0*y_backward_j[i] - y_backward_2j[i])/(12.0*pow(h_j, 2));
02681
02682 for(unsigned int k = j; k < n; k++)
02683 {
02684 h_k = calculate_h(x[k]);
02685
02686 x_backward_jk[j] -= h_j;
02687 x_backward_jk[k] -= h_k;
02688 y_backward_jk = (t.*f)(x_backward_jk);
02689 x_backward_jk[j] += h_j;
02690 x_backward_jk[k] += h_k;
02691
02692 x_forward_jk[j] += h_j;
02693 x_forward_jk[k] += h_k;
02694 y_forward_jk = (t.*f)(x_forward_jk);
02695 x_forward_jk[j] -= h_j;
02696 x_forward_jk[k] -= h_k;
02697
02698 x_backward_j_forward_k[j] -= h_j;
02699 x_backward_j_forward_k[k] += h_k;
02700 y_backward_j_forward_k = (t.*f)(x_backward_j_forward_k);
02701 x_backward_j_forward_k[j] += h_j;
02702 x_backward_j_forward_k[k] -= h_k;
02703
02704 x_forward_j_backward_k[j] += h_j;
02705 x_forward_j_backward_k[k] -= h_k;
02706 y_forward_j_backward_k = (t.*f)(x_forward_j_backward_k);
02707 x_forward_j_backward_k[j] -= h_j;
02708 x_forward_j_backward_k[k] += h_k;
02709
02710 H[i][j][k] = (y_forward_jk[i] - y_forward_j_backward_k[i] - y_backward_j_forward_k[i] + y_backward_jk[i])/(4.0*h_j*h_k);
02711 }
02712 }
02713
02714 for(unsigned int j = 0; j < n; j++)
02715 {
02716 for(unsigned int k = 0; k < j; k++)
02717 {
02718 H[i][j][k] = H[i][k][j];
02719 }
02720 }
02721 }
02722
02723 return(H);
02724 }
02725
02726
02727
02728
02734
02735 template<class T>
02736 Vector< Matrix<double> > calculate_Hessian_form(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
02737 {
02738 switch(numerical_differentiation_method)
02739 {
02740 case ForwardDifferences:
02741 {
02742 return(calculate_forward_differences_Hessian_form(t, f, x));
02743 }
02744 break;
02745
02746 case CentralDifferences:
02747 {
02748 return(calculate_central_differences_Hessian_form(t, f, x));
02749 }
02750 break;
02751
02752 default:
02753 {
02754 std::ostringstream buffer;
02755
02756 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
02757 << "double calculate_Hessian_form(const T& t, Vector<double> (T::*f)(const Vector<double>&), const Vector<double>&) const method.\n"
02758 << "Unknown numerical differentiation method.\n";
02759
02760 throw std::logic_error(buffer.str());
02761 }
02762 break;
02763 }
02764 }
02765
02766
02767
02768
02769
02777
02778 template<class T>
02779 Vector< Matrix <double> > calculate_forward_differences_Hessian_form
02780 (const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy_vector, const Vector<double>& x) const
02781 {
02782 Vector<double> y = (t.*f)(dummy_vector, x);
02783
02784 unsigned int s = y.size();
02785 unsigned int n = x.size();
02786
02787 double h_j;
02788 double h_k;
02789
02790 Vector<double> x_forward_j(x);
02791 Vector<double> x_forward_2j(x);
02792
02793 Vector<double> x_forward_k(x);
02794 Vector<double> x_forward_jk(x);
02795
02796 Vector<double> y_forward_j(s);
02797 Vector<double> y_forward_2j(s);
02798
02799 Vector<double> y_forward_k(s);
02800 Vector<double> y_forward_jk(s);
02801
02802 Vector< Matrix<double> > H(s);
02803
02804 for(unsigned int i = 0; i < s; i++)
02805 {
02806 H[i].set(n,n);
02807
02808 for(unsigned int j = 0; j < n; j++)
02809 {
02810 h_j = calculate_h(x[j]);
02811
02812 x_forward_j[j] += h_j;
02813 y_forward_j = (t.*f)(dummy_vector, x_forward_j);
02814 x_forward_j[j] -= h_j;
02815
02816 x_forward_2j[j] += 2.0*h_j;
02817 y_forward_2j = (t.*f)(dummy_vector, x_forward_2j);
02818 x_forward_2j[j] -= 2.0*h_j;
02819
02820 H[i][j][j] = (y_forward_2j[i] - 2.0*y_forward_j[i] + y[i])/pow(h_j, 2);
02821
02822 for(unsigned int k = j; k < n; k++)
02823 {
02824 h_k = calculate_h(x[k]);
02825
02826 x_forward_k[k] += h_k;
02827 y_forward_k = (t.*f)(dummy_vector, x_forward_k);
02828 x_forward_k[k] -= h_k;
02829
02830 x_forward_jk[j] += h_j;
02831 x_forward_jk[k] += h_k;
02832 y_forward_jk = (t.*f)(dummy_vector, x_forward_jk);
02833 x_forward_jk[j] -= h_j;
02834 x_forward_jk[k] -= h_k;
02835
02836 H[i][j][k] = (y_forward_jk[i] - y_forward_j[i] - y_forward_k[i] + y[i])/(h_j*h_k);
02837 }
02838 }
02839
02840 for(unsigned int j = 0; j < n; j++)
02841 {
02842 for(unsigned int k = 0; k < j; k++)
02843 {
02844 H[i][j][k] = H[i][k][j];
02845 }
02846 }
02847 }
02848
02849 return(H);
02850 }
02851
02852
02853
02854
02855
02863
02864 template<class T>
02865 Vector< Matrix <double> > calculate_central_differences_Hessian_form
02866 (const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy_vector, const Vector<double>& x) const
02867 {
02868 Vector<double> y = (t.*f)(dummy_vector, x);
02869
02870 unsigned int s = y.size();
02871 unsigned int n = x.size();
02872
02873 double h_j;
02874 double h_k;
02875
02876 Vector<double> x_backward_2j(x);
02877 Vector<double> x_backward_j(x);
02878
02879 Vector<double> x_forward_j(x);
02880 Vector<double> x_forward_2j(x);
02881
02882 Vector<double> x_backward_jk(x);
02883 Vector<double> x_forward_jk(x);
02884
02885 Vector<double> x_backward_j_forward_k(x);
02886 Vector<double> x_forward_j_backward_k(x);
02887
02888 Vector<double> y_backward_2j;
02889 Vector<double> y_backward_j;
02890
02891 Vector<double> y_forward_j;
02892 Vector<double> y_forward_2j;
02893
02894 Vector<double> y_backward_jk;
02895 Vector<double> y_forward_jk;
02896
02897 Vector<double> y_backward_j_forward_k;
02898 Vector<double> y_forward_j_backward_k;
02899
02900 Vector< Matrix<double> > H(s);
02901
02902 for(unsigned int i = 0; i < s; i++)
02903 {
02904 H[i].set(n,n);
02905
02906 for(unsigned int j = 0; j < n; j++)
02907 {
02908 h_j = calculate_h(x[j]);
02909
02910 x_backward_2j[j] -= 2.0*h_j;
02911 y_backward_2j = (t.*f)(dummy_vector, x_backward_2j);
02912 x_backward_2j[j] += 2.0*h_j;
02913
02914 x_backward_j[j] -= h_j;
02915 y_backward_j = (t.*f)(dummy_vector, x_backward_j);
02916 x_backward_j[j] += h_j;
02917
02918 x_forward_j[j] += h_j;
02919 y_forward_j = (t.*f)(dummy_vector, x_forward_j);
02920 x_forward_j[j] -= h_j;
02921
02922 x_forward_2j[j] += 2.0*h_j;
02923 y_forward_2j = (t.*f)(dummy_vector, x_forward_2j);
02924 x_forward_2j[j] -= 2.0*h_j;
02925
02926 H[i][j][j] = (-y_forward_2j[i] + 16.0*y_forward_j[i] -30.0*y[i] + 16.0*y_backward_j[i] - y_backward_2j[i])/(12.0*pow(h_j, 2));
02927
02928 for(unsigned int k = j; k < n; k++)
02929 {
02930 h_k = calculate_h(x[k]);
02931
02932 x_backward_jk[j] -= h_j;
02933 x_backward_jk[k] -= h_k;
02934 y_backward_jk = (t.*f)(dummy_vector, x_backward_jk);
02935 x_backward_jk[j] += h_j;
02936 x_backward_jk[k] += h_k;
02937
02938 x_forward_jk[j] += h_j;
02939 x_forward_jk[k] += h_k;
02940 y_forward_jk = (t.*f)(dummy_vector, x_forward_jk);
02941 x_forward_jk[j] -= h_j;
02942 x_forward_jk[k] -= h_k;
02943
02944 x_backward_j_forward_k[j] -= h_j;
02945 x_backward_j_forward_k[k] += h_k;
02946 y_backward_j_forward_k = (t.*f)(dummy_vector, x_backward_j_forward_k);
02947 x_backward_j_forward_k[j] += h_j;
02948 x_backward_j_forward_k[k] -= h_k;
02949
02950 x_forward_j_backward_k[j] += h_j;
02951 x_forward_j_backward_k[k] -= h_k;
02952 y_forward_j_backward_k = (t.*f)(dummy_vector, x_forward_j_backward_k);
02953 x_forward_j_backward_k[j] -= h_j;
02954 x_forward_j_backward_k[k] += h_k;
02955
02956 H[i][j][k] = (y_forward_jk[i] - y_forward_j_backward_k[i] - y_backward_j_forward_k[i] + y_backward_jk[i])/(4.0*h_j*h_k);
02957 }
02958 }
02959
02960 for(unsigned int j = 0; j < n; j++)
02961 {
02962 for(unsigned int k = 0; k < j; k++)
02963 {
02964 H[i][j][k] = H[i][k][j];
02965 }
02966 }
02967 }
02968
02969 return(H);
02970 }
02971
02972
02973
02974
02975
02983
02984 template<class T>
02985 Vector< Matrix<double> > calculate_Hessian_form
02986 (const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy_vector, const Vector<double>& x) const
02987 {
02988 switch(numerical_differentiation_method)
02989 {
02990 case ForwardDifferences:
02991 {
02992 return(calculate_forward_differences_Hessian_form(t, f, dummy_vector, x));
02993 }
02994 break;
02995
02996 case CentralDifferences:
02997 {
02998 return(calculate_central_differences_Hessian_form(t, f, dummy_vector, x));
02999 }
03000 break;
03001
03002 default:
03003 {
03004 std::ostringstream buffer;
03005
03006 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
03007 << "Vector< Matrix<double> > calculate_Hessian_form\n"
03008 << "(const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method.\n"
03009 << "Unknown numerical differentiation method.\n";
03010
03011 throw std::logic_error(buffer.str());
03012 }
03013 break;
03014 }
03015 }
03016
03017
03018
03019
03027
03028 template<class T>
03029 Vector< Matrix <double> > calculate_forward_differences_Hessian_form(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
03030 {
03031 Vector<double> y = (t.*f)(dummy, x);
03032
03033 unsigned int s = y.size();
03034 unsigned int n = x.size();
03035
03036 double h_j;
03037 double h_k;
03038
03039 Vector<double> x_forward_j(x);
03040 Vector<double> x_forward_2j(x);
03041
03042 Vector<double> x_forward_k(x);
03043 Vector<double> x_forward_jk(x);
03044
03045 Vector<double> y_forward_j(s);
03046 Vector<double> y_forward_2j(s);
03047
03048 Vector<double> y_forward_k(s);
03049 Vector<double> y_forward_jk(s);
03050
03051 Vector< Matrix<double> > H(s);
03052
03053 for(unsigned int i = 0; i < s; i++)
03054 {
03055 H[i].set(n,n);
03056
03057 for(unsigned int j = 0; j < n; j++)
03058 {
03059 h_j = calculate_h(x[j]);
03060
03061 x_forward_j[j] += h_j;
03062 y_forward_j = (t.*f)(dummy, x_forward_j);
03063 x_forward_j[j] -= h_j;
03064
03065 x_forward_2j[j] += 2.0*h_j;
03066 y_forward_2j = (t.*f)(dummy, x_forward_2j);
03067 x_forward_2j[j] -= 2.0*h_j;
03068
03069 H[i][j][j] = (y_forward_2j[i] - 2.0*y_forward_j[i] + y[i])/pow(h_j, 2);
03070
03071 for(unsigned int k = j; k < n; k++)
03072 {
03073 h_k = calculate_h(x[k]);
03074
03075 x_forward_k[k] += h_k;
03076 y_forward_k = (t.*f)(dummy, x_forward_k);
03077 x_forward_k[k] -= h_k;
03078
03079 x_forward_jk[j] += h_j;
03080 x_forward_jk[k] += h_k;
03081 y_forward_jk = (t.*f)(dummy, x_forward_jk);
03082 x_forward_jk[j] -= h_j;
03083 x_forward_jk[k] -= h_k;
03084
03085 H[i][j][k] = (y_forward_jk[i] - y_forward_j[i] - y_forward_k[i] + y[i])/(h_j*h_k);
03086 }
03087 }
03088
03089 for(unsigned int j = 0; j < n; j++)
03090 {
03091 for(unsigned int k = 0; k < j; k++)
03092 {
03093 H[i][j][k] = H[i][k][j];
03094 }
03095 }
03096 }
03097
03098 return(H);
03099 }
03100
03101
03102
03103
03111
03112 template<class T>
03113 Vector< Matrix <double> > calculate_central_differences_Hessian_form(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
03114 {
03115 Vector<double> y = (t.*f)(dummy, x);
03116
03117 unsigned int s = y.size();
03118 unsigned int n = x.size();
03119
03120 double h_j;
03121 double h_k;
03122
03123 Vector<double> x_backward_2j(x);
03124 Vector<double> x_backward_j(x);
03125
03126 Vector<double> x_forward_j(x);
03127 Vector<double> x_forward_2j(x);
03128
03129 Vector<double> x_backward_jk(x);
03130 Vector<double> x_forward_jk(x);
03131
03132 Vector<double> x_backward_j_forward_k(x);
03133 Vector<double> x_forward_j_backward_k(x);
03134
03135 Vector<double> y_backward_2j;
03136 Vector<double> y_backward_j;
03137
03138 Vector<double> y_forward_j;
03139 Vector<double> y_forward_2j;
03140
03141 Vector<double> y_backward_jk;
03142 Vector<double> y_forward_jk;
03143
03144 Vector<double> y_backward_j_forward_k;
03145 Vector<double> y_forward_j_backward_k;
03146
03147 Vector< Matrix<double> > H(s);
03148
03149 for(unsigned int i = 0; i < s; i++)
03150 {
03151 H[i].set(n,n);
03152
03153 for(unsigned int j = 0; j < n; j++)
03154 {
03155 h_j = calculate_h(x[j]);
03156
03157 x_backward_2j[j] -= 2.0*h_j;
03158 y_backward_2j = (t.*f)(dummy, x_backward_2j);
03159 x_backward_2j[j] += 2.0*h_j;
03160
03161 x_backward_j[j] -= h_j;
03162 y_backward_j = (t.*f)(dummy, x_backward_j);
03163 x_backward_j[j] += h_j;
03164
03165 x_forward_j[j] += h_j;
03166 y_forward_j = (t.*f)(dummy, x_forward_j);
03167 x_forward_j[j] -= h_j;
03168
03169 x_forward_2j[j] += 2.0*h_j;
03170 y_forward_2j = (t.*f)(dummy, x_forward_2j);
03171 x_forward_2j[j] -= 2.0*h_j;
03172
03173 H[i][j][j] = (-y_forward_2j[i] + 16.0*y_forward_j[i] -30.0*y[i] + 16.0*y_backward_j[i] - y_backward_2j[i])/(12.0*pow(h_j, 2));
03174
03175 for(unsigned int k = j; k < n; k++)
03176 {
03177 h_k = calculate_h(x[k]);
03178
03179 x_backward_jk[j] -= h_j;
03180 x_backward_jk[k] -= h_k;
03181 y_backward_jk = (t.*f)(dummy, x_backward_jk);
03182 x_backward_jk[j] += h_j;
03183 x_backward_jk[k] += h_k;
03184
03185 x_forward_jk[j] += h_j;
03186 x_forward_jk[k] += h_k;
03187 y_forward_jk = (t.*f)(dummy, x_forward_jk);
03188 x_forward_jk[j] -= h_j;
03189 x_forward_jk[k] -= h_k;
03190
03191 x_backward_j_forward_k[j] -= h_j;
03192 x_backward_j_forward_k[k] += h_k;
03193 y_backward_j_forward_k = (t.*f)(dummy, x_backward_j_forward_k);
03194 x_backward_j_forward_k[j] += h_j;
03195 x_backward_j_forward_k[k] -= h_k;
03196
03197 x_forward_j_backward_k[j] += h_j;
03198 x_forward_j_backward_k[k] -= h_k;
03199 y_forward_j_backward_k = (t.*f)(dummy, x_forward_j_backward_k);
03200 x_forward_j_backward_k[j] -= h_j;
03201 x_forward_j_backward_k[k] += h_k;
03202
03203 H[i][j][k] = (y_forward_jk[i] - y_forward_j_backward_k[i] - y_backward_j_forward_k[i] + y_backward_jk[i])/(4.0*h_j*h_k);
03204 }
03205 }
03206
03207 for(unsigned int j = 0; j < n; j++)
03208 {
03209 for(unsigned int k = 0; k < j; k++)
03210 {
03211 H[i][j][k] = H[i][k][j];
03212 }
03213 }
03214 }
03215
03216 return(H);
03217 }
03218
03219
03220
03221
03229
03230 template<class T>
03231 Vector< Matrix<double> > calculate_Hessian_form(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&) const, const unsigned int& dummy, const Vector<double>& x) const
03232 {
03233 switch(numerical_differentiation_method)
03234 {
03235 case ForwardDifferences:
03236 {
03237 return(calculate_forward_differences_Hessian_form(t, f, dummy, x));
03238 }
03239 break;
03240
03241 case CentralDifferences:
03242 {
03243 return(calculate_central_differences_Hessian_form(t, f, dummy, x));
03244 }
03245 break;
03246
03247 default:
03248 {
03249 std::ostringstream buffer;
03250
03251 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
03252 << "double calculate_Hessian_form(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&), const unsigned int&, const Vector<double>&) const method.\n"
03253 << "Unknown numerical differentiation method.\n";
03254
03255 throw std::logic_error(buffer.str());
03256 }
03257 break;
03258 }
03259 }
03260
03261
03262
03263
03264
03273
03274 template<class T>
03275 Vector< Matrix <double> > calculate_forward_differences_Hessian_form
03276 (const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&, const Vector<double>&) const, const unsigned int& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
03277 {
03278 Vector<double> y = (t.*f)(dummy_int, dummy_vector, x);
03279
03280 unsigned int s = y.size();
03281 unsigned int n = x.size();
03282
03283 double h_j;
03284 double h_k;
03285
03286 Vector<double> x_forward_j(x);
03287 Vector<double> x_forward_2j(x);
03288
03289 Vector<double> x_forward_k(x);
03290 Vector<double> x_forward_jk(x);
03291
03292 Vector<double> y_forward_j(s);
03293 Vector<double> y_forward_2j(s);
03294
03295 Vector<double> y_forward_k(s);
03296 Vector<double> y_forward_jk(s);
03297
03298 Vector< Matrix<double> > H(s);
03299
03300 for(unsigned int i = 0; i < s; i++)
03301 {
03302 H[i].set(n,n);
03303
03304 for(unsigned int j = 0; j < n; j++)
03305 {
03306 h_j = calculate_h(x[j]);
03307
03308 x_forward_j[j] += h_j;
03309 y_forward_j = (t.*f)(dummy_int, dummy_vector, x_forward_j);
03310 x_forward_j[j] -= h_j;
03311
03312 x_forward_2j[j] += 2.0*h_j;
03313 y_forward_2j = (t.*f)(dummy_int, dummy_vector, x_forward_2j);
03314 x_forward_2j[j] -= 2.0*h_j;
03315
03316 H[i][j][j] = (y_forward_2j[i] - 2.0*y_forward_j[i] + y[i])/pow(h_j, 2);
03317
03318 for(unsigned int k = j; k < n; k++)
03319 {
03320 h_k = calculate_h(x[k]);
03321
03322 x_forward_k[k] += h_k;
03323 y_forward_k = (t.*f)(dummy_int, dummy_vector, x_forward_k);
03324 x_forward_k[k] -= h_k;
03325
03326 x_forward_jk[j] += h_j;
03327 x_forward_jk[k] += h_k;
03328 y_forward_jk = (t.*f)(dummy_int, dummy_vector, x_forward_jk);
03329 x_forward_jk[j] -= h_j;
03330 x_forward_jk[k] -= h_k;
03331
03332 H[i][j][k] = (y_forward_jk[i] - y_forward_j[i] - y_forward_k[i] + y[i])/(h_j*h_k);
03333 }
03334 }
03335
03336 for(unsigned int j = 0; j < n; j++)
03337 {
03338 for(unsigned int k = 0; k < j; k++)
03339 {
03340 H[i][j][k] = H[i][k][j];
03341 }
03342 }
03343 }
03344
03345 return(H);
03346 }
03347
03348
03349
03350
03351
03360
03361 template<class T>
03362 Vector< Matrix <double> > calculate_central_differences_Hessian_form
03363 (const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&, const Vector<double>&) const, const unsigned int& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
03364 {
03365 const Vector<double> y = (t.*f)(dummy_int, dummy_vector, x);
03366
03367 unsigned int s = y.size();
03368 unsigned int n = x.size();
03369
03370 double h_j;
03371 double h_k;
03372
03373 Vector<double> x_backward_2j(x);
03374 Vector<double> x_backward_j(x);
03375
03376 Vector<double> x_forward_j(x);
03377 Vector<double> x_forward_2j(x);
03378
03379 Vector<double> x_backward_jk(x);
03380 Vector<double> x_forward_jk(x);
03381
03382 Vector<double> x_backward_j_forward_k(x);
03383 Vector<double> x_forward_j_backward_k(x);
03384
03385 Vector<double> y_backward_2j;
03386 Vector<double> y_backward_j;
03387
03388 Vector<double> y_forward_j;
03389 Vector<double> y_forward_2j;
03390
03391 Vector<double> y_backward_jk;
03392 Vector<double> y_forward_jk;
03393
03394 Vector<double> y_backward_j_forward_k;
03395 Vector<double> y_forward_j_backward_k;
03396
03397 Vector< Matrix<double> > H(s);
03398
03399 for(unsigned int i = 0; i < s; i++)
03400 {
03401 H[i].set(n,n);
03402
03403 for(unsigned int j = 0; j < n; j++)
03404 {
03405 h_j = calculate_h(x[j]);
03406
03407 x_backward_2j[j] -= 2.0*h_j;
03408 y_backward_2j = (t.*f)(dummy_int, dummy_vector, x_backward_2j);
03409 x_backward_2j[j] += 2.0*h_j;
03410
03411 x_backward_j[j] -= h_j;
03412 y_backward_j = (t.*f)(dummy_int, dummy_vector, x_backward_j);
03413 x_backward_j[j] += h_j;
03414
03415 x_forward_j[j] += h_j;
03416 y_forward_j = (t.*f)(dummy_int, dummy_vector, x_forward_j);
03417 x_forward_j[j] -= h_j;
03418
03419 x_forward_2j[j] += 2.0*h_j;
03420 y_forward_2j = (t.*f)(dummy_int, dummy_vector, x_forward_2j);
03421 x_forward_2j[j] -= 2.0*h_j;
03422
03423 H[i][j][j] = (-y_forward_2j[i] + 16.0*y_forward_j[i] -30.0*y[i] + 16.0*y_backward_j[i] - y_backward_2j[i])/(12.0*pow(h_j, 2));
03424
03425 for(unsigned int k = j; k < n; k++)
03426 {
03427 h_k = calculate_h(x[k]);
03428
03429 x_backward_jk[j] -= h_j;
03430 x_backward_jk[k] -= h_k;
03431 y_backward_jk = (t.*f)(dummy_int, dummy_vector, x_backward_jk);
03432 x_backward_jk[j] += h_j;
03433 x_backward_jk[k] += h_k;
03434
03435 x_forward_jk[j] += h_j;
03436 x_forward_jk[k] += h_k;
03437 y_forward_jk = (t.*f)(dummy_int, dummy_vector, x_forward_jk);
03438 x_forward_jk[j] -= h_j;
03439 x_forward_jk[k] -= h_k;
03440
03441 x_backward_j_forward_k[j] -= h_j;
03442 x_backward_j_forward_k[k] += h_k;
03443 y_backward_j_forward_k = (t.*f)(dummy_int, dummy_vector, x_backward_j_forward_k);
03444 x_backward_j_forward_k[j] += h_j;
03445 x_backward_j_forward_k[k] -= h_k;
03446
03447 x_forward_j_backward_k[j] += h_j;
03448 x_forward_j_backward_k[k] -= h_k;
03449 y_forward_j_backward_k = (t.*f)(dummy_int, dummy_vector, x_forward_j_backward_k);
03450 x_forward_j_backward_k[j] -= h_j;
03451 x_forward_j_backward_k[k] += h_k;
03452
03453 H[i][j][k] = (y_forward_jk[i] - y_forward_j_backward_k[i] - y_backward_j_forward_k[i] + y_backward_jk[i])/(4.0*h_j*h_k);
03454 }
03455 }
03456
03457 for(unsigned int j = 0; j < n; j++)
03458 {
03459 for(unsigned int k = 0; k < j; k++)
03460 {
03461 H[i][j][k] = H[i][k][j];
03462 }
03463 }
03464 }
03465
03466 return(H);
03467 }
03468
03469
03470
03471
03472
03481
03482 template<class T>
03483 Vector< Matrix<double> > calculate_Hessian_form
03484 (const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&, const Vector<double>&) const, const unsigned int& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
03485 {
03486 switch(numerical_differentiation_method)
03487 {
03488 case ForwardDifferences:
03489 {
03490 return(calculate_forward_differences_Hessian_form(t, f, dummy_int, dummy_vector, x));
03491 }
03492 break;
03493
03494 case CentralDifferences:
03495 {
03496 return(calculate_central_differences_Hessian_form(t, f, dummy_int, dummy_vector, x));
03497 }
03498 break;
03499
03500 default:
03501 {
03502 std::ostringstream buffer;
03503
03504 buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
03505 << "Vector< Matrix<double> > calculate_Hessian_form\n"
03506 << "(const T& t, Vector<double> (T::*f)(const unsigned int&, const Vector<double>&, const Vector<double>&) const, const unsigned int&, const Vector<double>&, const Vector<double>&) const method.\n"
03507 << "Unknown numerical differentiation method.\n";
03508
03509 throw std::logic_error(buffer.str());
03510 }
03511 break;
03512 }
03513 }
03514
03515
03516 private:
03517
03519
03520 NumericalDifferentiationMethod numerical_differentiation_method;
03521
03523
03524 unsigned int precision_digits;
03525
03527
03528 bool display;
03529
03530 };
03531
03532 }
03533
03534 #endif
03535
03536
03537
03538
03539
03540
03541
03542
03543
03544
03545
03546
03547
03548
03549
03550
03551
03552
03553