OpenNN  2.2
Open Neural Networks Library
numerical_differentiation.h
1 /****************************************************************************************************************/
2 /* */
3 /* OpenNN: Open Neural Networks Library */
4 /* www.artelnics.com/opennn */
5 /* */
6 /* N U M E R I C A L D I F F E R E N T I A T I O N C L A S S H E A D E R */
7 /* */
8 /* Roberto Lopez */
9 /* Artelnics - Making intelligent use of data */
11 /* */
12 /****************************************************************************************************************/
13 
14 #ifndef __NUMERICALDIFFERENTIATION_H__
15 #define __NUMERICALDIFFERENTIATION_H__
16 
17 // System includes
18 
19 #include<iostream>
20 
21 // OpenNN includes
22 
23 #include "vector.h"
24 #include "matrix.h"
25 
26 // TinyXml includes
27 
28 #include "../tinyxml2/tinyxml2.h"
29 
30 
31 namespace OpenNN
32 {
33 
36 
38 {
39 
40 public:
41 
42  // DEFAULT CONSTRUCTOR
43 
44  explicit NumericalDifferentiation(void);
45 
46  // COPY CONSTRUCTOR
47 
49 
50  // DESTRUCTOR
51 
52  virtual ~NumericalDifferentiation(void);
53 
54  // ASSIGNMENT OPERATOR
55 
57 
58  // EQUAL TO OPERATOR
59 
60  bool operator == (const NumericalDifferentiation&) const;
61 
63 
64  enum NumericalDifferentiationMethod{ForwardDifferences, CentralDifferences};
65 
66  // METHODS
67 
69  std::string write_numerical_differentiation_method(void) const;
70 
71  const size_t& get_precision_digits(void) const;
72 
73  const bool& get_display(void) const;
74 
75  void set(const NumericalDifferentiation&);
76 
78  void set_numerical_differentiation_method(const std::string&);
79 
80  void set_precision_digits(const size_t&);
81 
82  void set_display(const bool&);
83 
84  void set_default(void);
85 
86  double calculate_h(const double&) const;
87 
89 
90  // Serialization methods
91 
92  tinyxml2::XMLDocument* to_XML(void) const;
93  void from_XML(const tinyxml2::XMLDocument&);
94 
95 
96  // DERIVATIVE METHODS
97 
98  // double calculate_forward_differences_derivative(const T&, double (T::*f)(const double&) const , double) const method
99 
104 
105  template<class T>
106  double calculate_forward_differences_derivative(const T& t, double (T::*f)(const double&) const, const double& x) const
107  {
108  const double y = (t.*f)(x);
109 
110  const double h = calculate_h(x);
111 
112  const double y_forward = (t.*f)(x + h);
113 
114  const double d = (y_forward - y)/h;
115 
116  return(d);
117  }
118 
119 
120  // double calculate_central_differences_derivative(const T&, double (T::*f)(const double&) const , double) const method
121 
126 
127  template<class T>
128  double calculate_central_differences_derivative(const T& t, double (T::*f)(const double&) const , const double& x) const
129  {
130  const double h = calculate_h(x);
131 
132  const double y_forward = (t.*f)(x+h);
133 
134  const double y_backward = (t.*f)(x-h);
135 
136  const double d = (y_forward - y_backward)/(2.0*h);
137 
138  return(d);
139  }
140 
141 
142  // double calculate_derivative(const T&, double (T::*f)(const double&) const , double) const method
143 
148 
149  template<class T>
150  double calculate_derivative(const T& t, double (T::*f)(const double&) const , const double& x) const
151  {
153  {
154  case ForwardDifferences:
155  {
157  }
158  break;
159 
160  case CentralDifferences:
161  {
163  }
164  break;
165 
166  default:
167  {
168  std::ostringstream buffer;
169 
170  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
171  << "double calculate_derivative(const T&, double (T::*f)(const double&) const , const double&) const method.\n"
172  << "Unknown numerical differentiation method.\n";
173 
174  throw std::logic_error(buffer.str());
175  }
176  break;
177  }
178  }
179 
180 
181  // Vector<double> calculate_forward_differences_derivative
182  //(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
183 
188 
189  template<class T>
191  {
192  const Vector<double> h = calculate_h(x);
193 
194  const Vector<double> y = (t.*f)(x);
195 
196  const Vector<double> x_forward = x + h;
197  const Vector<double> y_forward = (t.*f)(x_forward);
198 
199  const Vector<double> d = (y_forward - y)/h;
200 
201  return(d);
202  }
203 
204 
205  // Vector<double> calculate_central_differences_derivative
206  //(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
207 
212 
213  template<class T>
215  {
216  const Vector<double> h = calculate_h(x);
217 
218  const Vector<double> x_forward = x + h;
219  const Vector<double> x_backward = x - h;
220 
221  const Vector<double> y_forward = (t.*f)(x_forward);
222  const Vector<double> y_backward = (t.*f)(x_backward);
223 
224  const Vector<double> y = (t.*f)(x);
225 
226  const Vector<double> d = (y_forward - y_backward)/(h*2.0);
227 
228  return(d);
229  }
230 
231 
232  // Vector<double> calculate_derivative
233  //(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
234 
239 
240  template<class T>
241  Vector<double> calculate_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
242  {
244  {
245  case ForwardDifferences:
246  {
248  }
249  break;
250 
251  case CentralDifferences:
252  {
254  }
255  break;
256 
257  default:
258  {
259  std::ostringstream buffer;
260 
261  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
262  << "Vector<double> calculate_derivative(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
263  << "Unknown numerical differentiation method.\n";
264 
265  throw std::logic_error(buffer.str());
266  }
267  break;
268  }
269  }
270 
271 
272  // Vector<double> calculate_forward_differences_derivative
273  //(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
274 
281 
282  template<class T>
283  Vector<double> calculate_forward_differences_derivative(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
284  {
285  const Vector<double> y = (t.*f)(dummy, x);
286 
287  const Vector<double> h = calculate_h(x);
288  const Vector<double> x_forward = x + h;
289 
290  const Vector<double> y_forward = (t.*f)(dummy, x_forward);
291 
292  const Vector<double> d = (y_forward - y)/h;
293 
294  return(d);
295  }
296 
297 
298  // Vector<double> calculate_central_differences_derivative
299  //(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
300 
307 
308  template<class T>
309  Vector<double> calculate_central_differences_derivative(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
310  {
311  const Vector<double> h = calculate_h(x);
312 
313  const Vector<double> x_forward = x + h;
314  const Vector<double> x_backward = x - h;
315 
316  const Vector<double> y_forward = (t.*f)(dummy, x_forward);
317  const Vector<double> y_backward = (t.*f)(dummy, x_backward);
318 
319  const Vector<double> d = (y_forward - y_backward)/(h*2.0);
320 
321  return(d);
322  }
323 
324 
325  // Vector<double> calculate_derivative
326  //(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
327 
334 
335  template<class T>
336  Vector<double> calculate_derivative(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
337  {
339  {
340  case ForwardDifferences:
341  {
342  return(calculate_forward_differences_derivative(t, f, dummy, x));
343  }
344  break;
345 
346  case CentralDifferences:
347  {
348  return(calculate_central_differences_derivative(t, f, dummy, x));
349  }
350  break;
351 
352  default:
353  {
354  std::ostringstream buffer;
355 
356  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
357  << "Vector<double> calculate_derivative(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method.\n"
358  << "Unknown numerical differentiation method.\n";
359 
360  throw std::logic_error(buffer.str());
361  }
362  break;
363  }
364  }
365 
366 
367  // SECOND DERIVATIVE METHODS
368 
369 
370  // double calculate_forward_differences_second_derivative(const T&, double (T::*f)(const double&) const , const double&) const method
371 
376 
377  template<class T>
378  double calculate_forward_differences_second_derivative(const T& t, double (T::*f)(const double&) const, const double& x) const
379  {
380  const double h = calculate_h(x);
381 
382  const double x_forward_2 = x + 2.0*h;
383 
384  const double y_forward_2 = (t.*f)(x_forward_2);
385 
386  const double x_forward = x + h;
387 
388  const double y_forward = (t.*f)(x_forward);
389 
390  const double y = (t.*f)(x);
391 
392  return((y_forward_2 - 2*y_forward + y)/pow(h, 2));
393  }
394 
395 
396  // double calculate_central_differences_second_derivative(const T&, double (T::*f)(const double&) const , const double&) const method
397 
402 
403  template<class T>
404  double calculate_central_differences_second_derivative(const T& t, double (T::*f)(const double&) const , const double& x) const
405  {
406  const double h = calculate_h(x);
407 
408  const double x_forward_2 = x + 2.0*h;
409 
410  const double y_forward_2 = (t.*f)(x_forward_2);
411 
412  const double x_forward = x + h;
413 
414  const double y_forward = (t.*f)(x_forward);
415 
416  const double y = (t.*f)(x);
417 
418  const double x_backward = x - h;
419 
420  const double y_backward = (t.*f)(x_backward);
421 
422  const double x_backward_2 = x - 2.0*h;
423 
424  const double y_backward_2 = (t.*f)(x_backward_2);
425 
426  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));
427 
428  return(d2);
429  }
430 
431 
432  // double calculate_second_derivative(const T&, double (T::*f)(const double&) const , const double&) const method
433 
438 
439  template<class T>
440  double calculate_second_derivative(const T& t, double (T::*f)(const double&) const , const double& x) const
441  {
443  {
444  case ForwardDifferences:
445  {
447  }
448  break;
449 
450  case CentralDifferences:
451  {
453  }
454  break;
455 
456  default:
457  {
458  std::ostringstream buffer;
459 
460  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
461  << "double calculate_second_derivative(const T& t, double (T::*f)(const double&) const, const double&) const method.\n"
462  << "Unknown numerical differentiation method.\n";
463 
464  throw std::logic_error(buffer.str());
465  }
466  break;
467  }
468  }
469 
470 
471  // Vector<double> calculate_forward_differences_second_derivative(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
472 
477 
478  template<class T>
480  {
481  const Vector<double> y = (t.*f)(x);
482 
483  const Vector<double> h = calculate_h(x);
484 
485  const Vector<double> x_forward = x + h;
486  const Vector<double> x_forward_2 = x + h*2.0;
487 
488  const Vector<double> y_forward = (t.*f)(x_forward);
489  const Vector<double> y_forward_2 = (t.*f)(x_forward_2);
490 
491  return((y_forward_2 - y_forward*2.0 + y)/(h*h));
492  }
493 
494 
495  // Vector<double> calculate_central_differences_second_derivative(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
496 
501 
502  template<class T>
504  {
505  const Vector<double> h = calculate_h(x);
506 
507  const Vector<double> x_forward = x + h;
508  const Vector<double> x_forward_2 = x + h*2.0;
509 
510  const Vector<double> x_backward = x - h;
511  const Vector<double> x_backward_2 = x - h*2.0;
512 
513  const Vector<double> y = (t.*f)(x);
514 
515  const Vector<double> y_forward = (t.*f)(x_forward);
516  const Vector<double> y_forward_2 = (t.*f)(x_forward_2);
517 
518  const Vector<double> y_backward = (t.*f)(x_backward);
519  const Vector<double> y_backward_2 = (t.*f)(x_backward_2);
520 
521  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));
522  }
523 
524 
525  // Vector<double> calculate_second_derivative(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
526 
531 
532  template<class T>
533  Vector<double> calculate_second_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
534  {
536  {
537  case ForwardDifferences:
538  {
540  }
541  break;
542 
543  case CentralDifferences:
544  {
546  }
547  break;
548 
549  default:
550  {
551  std::ostringstream buffer;
552 
553  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
554  << "Vector<double> calculate_second_derivative(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
555  << "Unknown numerical differentiation method.\n";
556 
557  throw std::logic_error(buffer.str());
558  }
559  break;
560  }
561  }
562 
563 
564  // Vector<double> calculate_forward_differences_second_derivative(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
565 
572 
573  template<class T>
574  Vector<double> calculate_forward_differences_second_derivative(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
575  {
576  const Vector<double> y = (t.*f)(dummy, x);
577 
578  const Vector<double> h = calculate_h(x);
579 
580  const Vector<double> x_forward = x + h;
581  const Vector<double> x_forward_2 = x + h*2.0;
582 
583  const Vector<double> y_forward = (t.*f)(dummy, x_forward);
584  const Vector<double> y_forward_2 = (t.*f)(dummy, x_forward_2);
585 
586  return((y_forward_2 - y_forward*2.0 + y)/(h*h));
587  }
588 
589 
590  // Vector<double> calculate_central_differences_second_derivative(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
591 
598 
599  template<class T>
600  Vector<double> calculate_central_differences_second_derivative(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
601  {
602  const Vector<double> h = calculate_h(x);
603 
604  const Vector<double> x_forward = x + h;
605  const Vector<double> x_forward_2 = x + h*2.0;
606 
607  const Vector<double> x_backward = x - h;
608  const Vector<double> x_backward_2 = x - h*2.0;
609 
610  const Vector<double> y = (t.*f)(dummy, x);
611 
612  const Vector<double> y_forward = (t.*f)(dummy, x_forward);
613  const Vector<double> y_forward_2 = (t.*f)(dummy, x_forward_2);
614 
615  const Vector<double> y_backward = (t.*f)(dummy, x_backward);
616  const Vector<double> y_backward_2 = (t.*f)(dummy, x_backward_2);
617 
618  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));
619  }
620 
621 
622  // Vector<double> calculate_second_derivative(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
623 
630 
631  template<class T>
632  Vector<double> calculate_second_derivative(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
633  {
635  {
636  case ForwardDifferences:
637  {
638  return(calculate_forward_differences_second_derivative(t, f, dummy, x));
639  }
640  break;
641 
642  case CentralDifferences:
643  {
644  return(calculate_central_differences_second_derivative(t, f, dummy, x));
645  }
646  break;
647 
648  default:
649  {
650  std::ostringstream buffer;
651 
652  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
653  << "Vector<double> calculate_second_derivative(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method.\n"
654  << "Unknown numerical differentiation method.\n";
655 
656  throw std::logic_error(buffer.str());
657  }
658  break;
659  }
660  }
661 
662 
663  // GRADIENT METHODS
664 
665  // Vector<double> calculate_forward_differences_gradient(const T&, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
666 
672 
673  template<class T>
674  Vector<double> calculate_forward_differences_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
675  {
676  const size_t n = x.size();
677 
678  double h;
679 
680  double y = (t.*f)(x);
681 
682  Vector<double> x_forward(x);
683 
684  double y_forward;
685 
686  Vector<double> g(n);
687 
688  for(size_t i = 0; i < n; i++)
689  {
690  h = calculate_h(x[i]);
691 
692  x_forward[i] += h;
693  y_forward = (t.*f)(x_forward);
694  x_forward[i] -= h;
695 
696  g[i] = (y_forward - y)/h;
697  }
698 
699  return(g);
700  }
701 
702 
703  // Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const method
704 
710 
711  template<class T>
712  Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
713  {
714  const size_t n = x.size();
715 
716  double h;
717 
718  Vector<double> x_forward(x);
719  Vector<double> x_backward(x);
720 
721  double y_forward;
722  double y_backward;
723 
724  Vector<double> g(n);
725 
726  for(size_t i = 0; i < n; i++)
727  {
728  h = calculate_h(x[i]);
729 
730  x_forward[i] += h;
731  y_forward = (t.*f)(x_forward);
732  x_forward[i] -= h;
733 
734  x_backward[i] -= h;
735  y_backward = (t.*f)(x_backward);
736  x_backward[i] += h;
737 
738  g[i] = (y_forward - y_backward)/(2.0*h);
739  }
740 
741  return(g);
742  }
743 
744 
745  // Vector<double> calculate_gradient(const T&, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
746 
752 
753  template<class T>
754  Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
755  {
757  {
758  case ForwardDifferences:
759  {
761  }
762  break;
763 
764  case CentralDifferences:
765  {
767  }
768  break;
769 
770  default:
771  {
772  std::ostringstream buffer;
773 
774  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
775  << "Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
776  << "Unknown numerical differentiation method.\n";
777 
778  throw std::logic_error(buffer.str());
779  }
780  break;
781  }
782  }
783 
784 
785  // Vector<double> calculate_forward_differences_gradient(const T&, double (T::*f)(const Vector<double>&), const Vector<double>&) const method
786 
792 
793  template<class T>
794  Vector<double> calculate_forward_differences_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>& x) const
795  {
796  const size_t n = x.size();
797 
798  double h;
799 
800  double y = (t.*f)(x);
801 
802  Vector<double> x_forward(x);
803 
804  double y_forward;
805 
806  Vector<double> g(n);
807 
808  for(size_t i = 0; i < n; i++)
809  {
810  h = calculate_h(x[i]);
811 
812  x_forward[i] += h;
813  y_forward = (t.*f)(x_forward);
814  x_forward[i] -= h;
815 
816  g[i] = (y_forward - y)/h;
817  }
818 
819  return(g);
820  }
821 
822 
823  // Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>& x) const method
824 
830 
831  template<class T>
832  Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>& x) const
833  {
834  const size_t n = x.size();
835 
836  double h;
837 
838  Vector<double> x_forward(x);
839  Vector<double> x_backward(x);
840 
841  double y_forward;
842  double y_backward;
843 
844  Vector<double> g(n);
845 
846  for(size_t i = 0; i < n; i++)
847  {
848  h = calculate_h(x[i]);
849 
850  x_forward[i] += h;
851  y_forward = (t.*f)(x_forward);
852  x_forward[i] -= h;
853 
854  x_backward[i] -= h;
855  y_backward = (t.*f)(x_backward);
856  x_backward[i] += h;
857 
858  g[i] = (y_forward - y_backward)/(2.0*h);
859  }
860 
861  return(g);
862  }
863 
864 
865  // Vector<double> calculate_gradient(const T&, double (T::*f)(const Vector<double>&), const Vector<double>&) const method
866 
872 
873  template<class T>
874  Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>& x) const
875  {
877  {
878  case ForwardDifferences:
879  {
881  }
882  break;
883 
884  case CentralDifferences:
885  {
887  }
888  break;
889 
890  default:
891  {
892  std::ostringstream buffer;
893 
894  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
895  << "Vector<double> calculate_gradient(const T& t, double (T::*f)(const Vector<double>&), const Vector<double>&) const method.\n"
896  << "Unknown numerical differentiation method.\n";
897 
898  throw std::logic_error(buffer.str());
899  }
900  break;
901  }
902  }
903 
904 
905  // Vector<double> calculate_forward_differences_gradient(const T&, double (T::*f)(const Vector<double>&, const Vector<double>&), const Vector<double>&, const Vector<double>&) const
906 
914 
915  template<class T>
916  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
917  {
918  const size_t n = x.size();
919 
920  double h;
921 
922  const double y = (t.*f)(dummy, x);
923 
924  Vector<double> x_forward(x);
925 
926  double y_forward;
927 
928  Vector<double> g(n);
929 
930  for(size_t i = 0; i < n; i++)
931  {
932  h = calculate_h(x[i]);
933 
934  x_forward[i] += h;
935  y_forward = (t.*f)(dummy, x_forward);
936  x_forward[i] -= h;
937 
938  g[i] = (y_forward - y)/h;
939  }
940 
941  return(g);
942  }
943 
944 
945  // Vector<double> calculate_central_differences_gradient(const T&, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method
946 
954 
955  template<class T>
956  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
957  {
958  const size_t n = x.size();
959 
960  double h;
961 
962  Vector<double> x_forward(x);
963  Vector<double> x_backward(x);
964 
965  double y_forward;
966  double y_backward;
967 
968  Vector<double> g(n);
969 
970  for(size_t i = 0; i < n; i++)
971  {
972  h = calculate_h(x[i]);
973 
974  x_forward[i] += h;
975  y_forward = (t.*f)(dummy, x_forward);
976  x_forward[i] -= h;
977 
978  x_backward[i] -= h;
979  y_backward = (t.*f)(dummy, x_backward);
980  x_backward[i] += h;
981 
982  g[i] = (y_forward - y_backward)/(2.0*h);
983  }
984 
985  return(g);
986  }
987 
988 
989  // Vector<double> calculate_gradient(const T&, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method
990 
998 
999  template<class T>
1000  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
1001  {
1003  {
1004  case ForwardDifferences:
1005  {
1006  return(calculate_forward_differences_gradient(t, f, dummy, x));
1007  }
1008  break;
1009 
1010  case CentralDifferences:
1011  {
1012  return(calculate_central_differences_gradient(t, f, dummy, x));
1013  }
1014  break;
1015 
1016  default:
1017  {
1018  std::ostringstream buffer;
1019 
1020  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
1021  << "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"
1022  << "Unknown numerical differentiation method.\n";
1023 
1024  throw std::logic_error(buffer.str());
1025  }
1026  break;
1027  }
1028  }
1029 
1030 
1031  // Vector<double> calculate_forward_differences_gradient(const T&, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
1032 
1040 
1041  template<class T>
1042  Vector<double> calculate_forward_differences_gradient(const T& t, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
1043  {
1044  size_t n = x.size();
1045 
1046  double h;
1047 
1048  double y = (t.*f)(dummy, x);
1049 
1050  Vector<double> x_forward(x);
1051 
1052  double y_forward;
1053 
1054  Vector<double> g(n);
1055 
1056  for(size_t i = 0; i < n; i++)
1057  {
1058  h = calculate_h(x[i]);
1059 
1060  x_forward[i] += h;
1061  y_forward = (t.*f)(dummy, x_forward);
1062  x_forward[i] -= h;
1063 
1064  g[i] = (y_forward - y)/h;
1065  }
1066 
1067  return(g);
1068  }
1069 
1070 
1071  // Vector<double> calculate_central_differences_gradient(const T&, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
1072 
1080 
1081  template<class T>
1082  Vector<double> calculate_central_differences_gradient(const T& t, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
1083  {
1084  size_t n = x.size();
1085 
1086  double h;
1087 
1088  Vector<double> x_forward(x);
1089  Vector<double> x_backward(x);
1090 
1091  double y_forward;
1092  double y_backward;
1093 
1094  Vector<double> g(n);
1095 
1096  for(size_t i = 0; i < n; i++)
1097  {
1098  h = calculate_h(x[i]);
1099 
1100  x_forward[i] += h;
1101  y_forward = (t.*f)(dummy, x_forward);
1102  x_forward[i] -= h;
1103 
1104  x_backward[i] -= h;
1105  y_backward = (t.*f)(dummy, x_backward);
1106  x_backward[i] += h;
1107 
1108  g[i] = (y_forward - y_backward)/(2.0*h);
1109  }
1110 
1111  return(g);
1112  }
1113 
1114 
1115  // Vector<double> calculate_gradient(const T&, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
1116 
1124 
1125  template<class T>
1126  Vector<double> calculate_gradient(const T& t, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
1127  {
1129  {
1130  case ForwardDifferences:
1131  {
1132  return(calculate_forward_differences_gradient(t, f, dummy, x));
1133  }
1134  break;
1135 
1136  case CentralDifferences:
1137  {
1138  return(calculate_central_differences_gradient(t, f, dummy, x));
1139  }
1140  break;
1141 
1142  default:
1143  {
1144  std::ostringstream buffer;
1145 
1146  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
1147  << "Vector<double> calculate_gradient(const T& t, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method.\n"
1148  << "Unknown numerical differentiation method.\n";
1149 
1150  throw std::logic_error(buffer.str());
1151  }
1152  break;
1153  }
1154  }
1155 
1156 
1157  // HESSIAN METHODS
1158 
1159 
1160  // Matrix<double> calculate_forward_differences_Hessian(const T&, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
1161 
1167 
1168  template<class T>
1169  Matrix<double> calculate_forward_differences_Hessian(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
1170  {
1171  size_t n = x.size();
1172 
1173  Matrix<double> H(n, n);
1174 
1175  double h_i;
1176  double h_j;
1177 
1178  double y = (t.*f)(x);
1179 
1180  Vector<double> x_forward_2i(x);
1181  Vector<double> x_forward_ij(x);
1182  Vector<double> x_forward_i(x);
1183  Vector<double> x_forward_j(x);
1184 
1185  double y_forward_2i;
1186  double y_forward_ij;
1187  double y_forward_i;
1188  double y_forward_j;
1189 
1190  for(size_t i = 0; i < n; i++)
1191  {
1192  h_i = calculate_h(x[i]);
1193 
1194  x_forward_i[i] += h_i;
1195  y_forward_i = (t.*f)(x_forward_i);
1196  x_forward_i[i] -= h_i;
1197 
1198  x_forward_2i[i] += 2.0*h_i;
1199  y_forward_2i = (t.*f)(x_forward_2i);
1200  x_forward_2i[i] -= 2.0*h_i;
1201 
1202  H(i,i) = (y_forward_2i - 2*y_forward_i + y)/pow(h_i, 2);
1203 
1204  for(size_t j = i; j < n; j++)
1205  {
1206  h_j = calculate_h(x[j]);
1207 
1208  x_forward_j[j] += h_j;
1209  y_forward_j = (t.*f)(x_forward_j);
1210  x_forward_j[j] -= h_j;
1211 
1212  x_forward_ij[i] += h_i;
1213  x_forward_ij[j] += h_j;
1214  y_forward_ij = (t.*f)(x_forward_ij);
1215  x_forward_ij[i] -= h_i;
1216  x_forward_ij[j] -= h_j;
1217 
1218  H(i,j) = (y_forward_ij - y_forward_i - y_forward_j + y)/(h_i*h_j);
1219  }
1220  }
1221 
1222  for(size_t i = 0; i < n; i++)
1223  {
1224  for(size_t j = 0; j < i; j++)
1225  {
1226  H(i,j) = H(j,i);
1227  }
1228  }
1229 
1230  return(H);
1231  }
1232 
1233 
1234  // Matrix<double> calculate_central_differences_Hessian(const T&, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
1235 
1241 
1242  template<class T>
1243  Matrix<double> calculate_central_differences_Hessian(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
1244  {
1245  size_t n = x.size();
1246 
1247  double y = (t.*f)(x);
1248 
1249  Matrix<double> H(n, n);
1250 
1251  double h_i;
1252  double h_j;
1253 
1254  Vector<double> x_backward_2i(x);
1255  Vector<double> x_backward_i(x);
1256 
1257  Vector<double> x_forward_i(x);
1258  Vector<double> x_forward_2i(x);
1259 
1260  Vector<double> x_backward_ij(x);
1261  Vector<double> x_forward_ij(x);
1262 
1263  Vector<double> x_backward_i_forward_j(x);
1264  Vector<double> x_forward_i_backward_j(x);
1265 
1266  double y_backward_2i;
1267  double y_backward_i;
1268 
1269  double y_forward_i;
1270  double y_forward_2i;
1271 
1272  double y_backward_ij;
1273  double y_forward_ij;
1274 
1275  double y_backward_i_forward_j;
1276  double y_forward_i_backward_j;
1277 
1278  for(size_t i = 0; i < n; i++)
1279  {
1280  h_i = calculate_h(x[i]);
1281 
1282  x_backward_2i[i] -= 2.0*h_i;
1283  y_backward_2i = (t.*f)(x_backward_2i);
1284  x_backward_2i[i] += 2.0*h_i;
1285 
1286  x_backward_i[i] -= h_i;
1287  y_backward_i = (t.*f)(x_backward_i);
1288  x_backward_i[i] += h_i;
1289 
1290  x_forward_i[i] += h_i;
1291  y_forward_i = (t.*f)(x_forward_i);
1292  x_forward_i[i] -= h_i;
1293 
1294  x_forward_2i[i] += 2.0*h_i;
1295  y_forward_2i = (t.*f)(x_forward_2i);
1296  x_forward_2i[i] -= 2.0*h_i;
1297 
1298  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));
1299 
1300  for(size_t j = i; j < n; j++)
1301  {
1302  h_j = calculate_h(x[j]);
1303 
1304  x_backward_ij[i] -= h_i;
1305  x_backward_ij[j] -= h_j;
1306  y_backward_ij = (t.*f)(x_backward_ij);
1307  x_backward_ij[i] += h_i;
1308  x_backward_ij[j] += h_j;
1309 
1310  x_forward_ij[i] += h_i;
1311  x_forward_ij[j] += h_j;
1312  y_forward_ij = (t.*f)(x_forward_ij);
1313  x_forward_ij[i] -= h_i;
1314  x_forward_ij[j] -= h_j;
1315 
1316  x_backward_i_forward_j[i] -= h_i;
1317  x_backward_i_forward_j[j] += h_j;
1318  y_backward_i_forward_j = (t.*f)(x_backward_i_forward_j);
1319  x_backward_i_forward_j[i] += h_i;
1320  x_backward_i_forward_j[j] -= h_j;
1321 
1322  x_forward_i_backward_j[i] += h_i;
1323  x_forward_i_backward_j[j] -= h_j;
1324  y_forward_i_backward_j = (t.*f)(x_forward_i_backward_j);
1325  x_forward_i_backward_j[i] -= h_i;
1326  x_forward_i_backward_j[j] += h_j;
1327 
1328  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);
1329  }
1330  }
1331 
1332  for(size_t i = 0; i < n; i++)
1333  {
1334  for(size_t j = 0; j < i; j++)
1335  {
1336  H(i,j) = H(j,i);
1337  }
1338  }
1339 
1340  return(H);
1341  }
1342 
1343 
1344  // Matrix<double> calculate_Hessian(const T&, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
1345 
1351 
1352  template<class T>
1353  Matrix<double> calculate_Hessian(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
1354  {
1356  {
1357  case ForwardDifferences:
1358  {
1359  return(calculate_forward_differences_Hessian(t, f, x));
1360  }
1361  break;
1362 
1363  case CentralDifferences:
1364  {
1365  return(calculate_central_differences_Hessian(t, f, x));
1366  }
1367  break;
1368 
1369  default:
1370  {
1371  std::ostringstream buffer;
1372 
1373  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
1374  << "double calculate_Hessian(const T& t, double (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
1375  << "Unknown numerical differentiation method.\n";
1376 
1377  throw std::logic_error(buffer.str());
1378  }
1379  break;
1380  }
1381  }
1382 
1383 
1384  // Matrix<double> calculate_forward_differences_Hessian(const T&, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy, const Vector<double>& x) const method
1385 
1393 
1394  template<class T>
1395  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
1396  {
1397  size_t n = x.size();
1398 
1399  Matrix<double> H(n, n);
1400 
1401  double h_i;
1402  double h_j;
1403 
1404  double y = (t.*f)(dummy, x);
1405 
1406  Vector<double> x_forward_2i(x);
1407  Vector<double> x_forward_ij(x);
1408  Vector<double> x_forward_i(x);
1409  Vector<double> x_forward_j(x);
1410 
1411  double y_forward_2i;
1412  double y_forward_ij;
1413  double y_forward_i;
1414  double y_forward_j;
1415 
1416  for(size_t i = 0; i < n; i++)
1417  {
1418  h_i = calculate_h(x[i]);
1419 
1420  x_forward_i[i] += h_i;
1421  y_forward_i = (t.*f)(dummy, x_forward_i);
1422  x_forward_i[i] -= h_i;
1423 
1424  x_forward_2i[i] += 2.0*h_i;
1425  y_forward_2i = (t.*f)(dummy, x_forward_2i);
1426  x_forward_2i[i] -= 2.0*h_i;
1427 
1428  H(i,i) = (y_forward_2i - 2*y_forward_i + y)/pow(h_i, 2);
1429 
1430  for(size_t j = i; j < n; j++)
1431  {
1432  h_j = calculate_h(x[j]);
1433 
1434  x_forward_j[j] += h_j;
1435  y_forward_j = (t.*f)(dummy, x_forward_j);
1436  x_forward_j[j] -= h_j;
1437 
1438  x_forward_ij[i] += h_i;
1439  x_forward_ij[j] += h_j;
1440  y_forward_ij = (t.*f)(dummy, x_forward_ij);
1441  x_forward_ij[i] -= h_i;
1442  x_forward_ij[j] -= h_j;
1443 
1444  H(i,j) = (y_forward_ij - y_forward_i - y_forward_j + y)/(h_i*h_j);
1445  }
1446  }
1447 
1448  for(size_t i = 0; i < n; i++)
1449  {
1450  for(size_t j = 0; j < i; j++)
1451  {
1452  H(i,j) = H(j,i);
1453  }
1454  }
1455 
1456  return(H);
1457  }
1458 
1459 
1460  // Matrix<double> calculate_central_differences_Hessian(const T&, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method
1461 
1469 
1470  template<class T>
1471  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
1472  {
1473  size_t n = x.size();
1474 
1475  double y = (t.*f)(dummy, x);
1476 
1477  Matrix<double> H(n, n);
1478 
1479  double h_i;
1480  double h_j;
1481 
1482  Vector<double> x_backward_2i(x);
1483  Vector<double> x_backward_i(x);
1484 
1485  Vector<double> x_forward_i(x);
1486  Vector<double> x_forward_2i(x);
1487 
1488  Vector<double> x_backward_ij(x);
1489  Vector<double> x_forward_ij(x);
1490 
1491  Vector<double> x_backward_i_forward_j(x);
1492  Vector<double> x_forward_i_backward_j(x);
1493 
1494  double y_backward_2i;
1495  double y_backward_i;
1496 
1497  double y_forward_i;
1498  double y_forward_2i;
1499 
1500  double y_backward_ij;
1501  double y_forward_ij;
1502 
1503  double y_backward_i_forward_j;
1504  double y_forward_i_backward_j;
1505 
1506  for(size_t i = 0; i < n; i++)
1507  {
1508  h_i = calculate_h(x[i]);
1509 
1510  x_backward_2i[i] -= 2.0*h_i;
1511  y_backward_2i = (t.*f)(dummy, x_backward_2i);
1512  x_backward_2i[i] += 2.0*h_i;
1513 
1514  x_backward_i[i] -= h_i;
1515  y_backward_i = (t.*f)(dummy, x_backward_i);
1516  x_backward_i[i] += h_i;
1517 
1518  x_forward_i[i] += h_i;
1519  y_forward_i = (t.*f)(dummy, x_forward_i);
1520  x_forward_i[i] -= h_i;
1521 
1522  x_forward_2i[i] += 2.0*h_i;
1523  y_forward_2i = (t.*f)(dummy, x_forward_2i);
1524  x_forward_2i[i] -= 2.0*h_i;
1525 
1526  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));
1527 
1528  for(size_t j = i; j < n; j++)
1529  {
1530  h_j = calculate_h(x[j]);
1531 
1532  x_backward_ij[i] -= h_i;
1533  x_backward_ij[j] -= h_j;
1534  y_backward_ij = (t.*f)(dummy, x_backward_ij);
1535  x_backward_ij[i] += h_i;
1536  x_backward_ij[j] += h_j;
1537 
1538  x_forward_ij[i] += h_i;
1539  x_forward_ij[j] += h_j;
1540  y_forward_ij = (t.*f)(dummy, x_forward_ij);
1541  x_forward_ij[i] -= h_i;
1542  x_forward_ij[j] -= h_j;
1543 
1544  x_backward_i_forward_j[i] -= h_i;
1545  x_backward_i_forward_j[j] += h_j;
1546  y_backward_i_forward_j = (t.*f)(dummy, x_backward_i_forward_j);
1547  x_backward_i_forward_j[i] += h_i;
1548  x_backward_i_forward_j[j] -= h_j;
1549 
1550  x_forward_i_backward_j[i] += h_i;
1551  x_forward_i_backward_j[j] -= h_j;
1552  y_forward_i_backward_j = (t.*f)(dummy, x_forward_i_backward_j);
1553  x_forward_i_backward_j[i] -= h_i;
1554  x_forward_i_backward_j[j] += h_j;
1555 
1556  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);
1557  }
1558  }
1559 
1560  for(size_t i = 0; i < n; i++)
1561  {
1562  for(size_t j = 0; j < i; j++)
1563  {
1564  H(i,j) = H(j,i);
1565  }
1566  }
1567 
1568  return(H);
1569  }
1570 
1571 
1572  // Matrix<double> calculate_Hessian(const T&, double (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method
1573 
1581 
1582  template<class T>
1583  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
1584  {
1586  {
1587  case ForwardDifferences:
1588  {
1589  return(calculate_forward_differences_Hessian(t, f, dummy, x));
1590  }
1591  break;
1592 
1593  case CentralDifferences:
1594  {
1595  return(calculate_central_differences_Hessian(t, f, dummy, x));
1596  }
1597  break;
1598 
1599  default:
1600  {
1601  std::ostringstream buffer;
1602 
1603  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
1604  << "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"
1605  << "Unknown numerical differentiation method.\n";
1606 
1607  throw std::logic_error(buffer.str());
1608  }
1609  break;
1610  }
1611  }
1612 
1613 
1614  // Matrix<double> calculate_forward_differences_Hessian(const T&, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
1615 
1623 
1624  template<class T>
1625  Matrix<double> calculate_forward_differences_Hessian(const T& t, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
1626  {
1627  size_t n = x.size();
1628 
1629  Matrix<double> H(n, n);
1630 
1631  double h_i;
1632  double h_j;
1633 
1634  double y = (t.*f)(dummy, x);
1635 
1636  Vector<double> x_forward_2i(x);
1637  Vector<double> x_forward_ij(x);
1638  Vector<double> x_forward_i(x);
1639  Vector<double> x_forward_j(x);
1640 
1641  double y_forward_2i;
1642  double y_forward_ij;
1643  double y_forward_i;
1644  double y_forward_j;
1645 
1646  for(size_t i = 0; i < n; i++)
1647  {
1648  h_i = calculate_h(x[i]);
1649 
1650  x_forward_i[i] += h_i;
1651  y_forward_i = (t.*f)(dummy, x_forward_i);
1652  x_forward_i[i] -= h_i;
1653 
1654  x_forward_2i[i] += 2.0*h_i;
1655  y_forward_2i = (t.*f)(dummy, x_forward_2i);
1656  x_forward_2i[i] -= 2.0*h_i;
1657 
1658  H(i,i) = (y_forward_2i - 2*y_forward_i + y)/pow(h_i, 2);
1659 
1660  for(size_t j = i; j < n; j++)
1661  {
1662  h_j = calculate_h(x[j]);
1663 
1664  x_forward_j[j] += h_j;
1665  y_forward_j = (t.*f)(dummy, x_forward_j);
1666  x_forward_j[j] -= h_j;
1667 
1668  x_forward_ij[i] += h_i;
1669  x_forward_ij[j] += h_j;
1670  y_forward_ij = (t.*f)(dummy, x_forward_ij);
1671  x_forward_ij[i] -= h_i;
1672  x_forward_ij[j] -= h_j;
1673 
1674  H(i,j) = (y_forward_ij - y_forward_i - y_forward_j + y)/(h_i*h_j);
1675  }
1676  }
1677 
1678  for(size_t i = 0; i < n; i++)
1679  {
1680  for(size_t j = 0; j < i; j++)
1681  {
1682  H(i,j) = H(j,i);
1683  }
1684  }
1685 
1686  return(H);
1687  }
1688 
1689 
1690  // Matrix<double> calculate_central_differences_Hessian(const T&, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
1691 
1699 
1700  template<class T>
1701  Matrix<double> calculate_central_differences_Hessian(const T& t, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
1702  {
1703  size_t n = x.size();
1704 
1705  double y = (t.*f)(dummy, x);
1706 
1707  Matrix<double> H(n, n);
1708 
1709  double h_i;
1710  double h_j;
1711 
1712  Vector<double> x_backward_2i(x);
1713  Vector<double> x_backward_i(x);
1714 
1715  Vector<double> x_forward_i(x);
1716  Vector<double> x_forward_2i(x);
1717 
1718  Vector<double> x_backward_ij(x);
1719  Vector<double> x_forward_ij(x);
1720 
1721  Vector<double> x_backward_i_forward_j(x);
1722  Vector<double> x_forward_i_backward_j(x);
1723 
1724  double y_backward_2i;
1725  double y_backward_i;
1726 
1727  double y_forward_i;
1728  double y_forward_2i;
1729 
1730  double y_backward_ij;
1731  double y_forward_ij;
1732 
1733  double y_backward_i_forward_j;
1734  double y_forward_i_backward_j;
1735 
1736  for(size_t i = 0; i < n; i++)
1737  {
1738  h_i = calculate_h(x[i]);
1739 
1740  x_backward_2i[i] -= 2.0*h_i;
1741  y_backward_2i = (t.*f)(dummy, x_backward_2i);
1742  x_backward_2i[i] += 2.0*h_i;
1743 
1744  x_backward_i[i] -= h_i;
1745  y_backward_i = (t.*f)(dummy, x_backward_i);
1746  x_backward_i[i] += h_i;
1747 
1748  x_forward_i[i] += h_i;
1749  y_forward_i = (t.*f)(dummy, x_forward_i);
1750  x_forward_i[i] -= h_i;
1751 
1752  x_forward_2i[i] += 2.0*h_i;
1753  y_forward_2i = (t.*f)(dummy, x_forward_2i);
1754  x_forward_2i[i] -= 2.0*h_i;
1755 
1756  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));
1757 
1758  for(size_t j = i; j < n; j++)
1759  {
1760  h_j = calculate_h(x[j]);
1761 
1762  x_backward_ij[i] -= h_i;
1763  x_backward_ij[j] -= h_j;
1764  y_backward_ij = (t.*f)(dummy, x_backward_ij);
1765  x_backward_ij[i] += h_i;
1766  x_backward_ij[j] += h_j;
1767 
1768  x_forward_ij[i] += h_i;
1769  x_forward_ij[j] += h_j;
1770  y_forward_ij = (t.*f)(dummy, x_forward_ij);
1771  x_forward_ij[i] -= h_i;
1772  x_forward_ij[j] -= h_j;
1773 
1774  x_backward_i_forward_j[i] -= h_i;
1775  x_backward_i_forward_j[j] += h_j;
1776  y_backward_i_forward_j = (t.*f)(dummy, x_backward_i_forward_j);
1777  x_backward_i_forward_j[i] += h_i;
1778  x_backward_i_forward_j[j] -= h_j;
1779 
1780  x_forward_i_backward_j[i] += h_i;
1781  x_forward_i_backward_j[j] -= h_j;
1782  y_forward_i_backward_j = (t.*f)(dummy, x_forward_i_backward_j);
1783  x_forward_i_backward_j[i] -= h_i;
1784  x_forward_i_backward_j[j] += h_j;
1785 
1786  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);
1787  }
1788  }
1789 
1790  for(size_t i = 0; i < n; i++)
1791  {
1792  for(size_t j = 0; j < i; j++)
1793  {
1794  H(i,j) = H(j,i);
1795  }
1796  }
1797 
1798  return(H);
1799  }
1800 
1801 
1802  // Matrix<double> calculate_Hessian(const T&, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
1803 
1811 
1812  template<class T>
1813  Matrix<double> calculate_Hessian(const T& t, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
1814  {
1816  {
1817  case ForwardDifferences:
1818  {
1819  return(calculate_forward_differences_Hessian(t, f, dummy, x));
1820  }
1821  break;
1822 
1823  case CentralDifferences:
1824  {
1825  return(calculate_central_differences_Hessian(t, f, dummy, x));
1826  }
1827  break;
1828 
1829  default:
1830  {
1831  std::ostringstream buffer;
1832 
1833  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
1834  << "double calculate_Hessian(const T& t, double (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method.\n"
1835  << "Unknown numerical differentiation method.\n";
1836 
1837  throw std::logic_error(buffer.str());
1838  }
1839  break;
1840  }
1841  }
1842 
1843 
1844  // JACOBIAN METHODS
1845 
1846  // Matrix<double> calculate_forward_differences_Jacobian(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const
1847 
1853 
1854  template<class T>
1856  {
1857  Vector<double> y = (t.*f)(x);
1858 
1859  double h;
1860 
1861  size_t n = x.size();
1862  size_t m = y.size();
1863 
1864  Vector<double> x_forward(x);
1865  Vector<double> y_forward(n);
1866 
1867  Matrix<double> J(m,n);
1868 
1869  for(size_t j = 0; j < n; j++)
1870  {
1871  h = calculate_h(x[j]);
1872 
1873  x_forward[j] += h;
1874  y_forward = (t.*f)(x_forward);
1875  x_forward[j] -= h;
1876 
1877  for(size_t i = 0; i < m; i++)
1878  {
1879  J(i,j) = (y_forward[i] - y[i])/h;
1880  }
1881  }
1882 
1883  return(J);
1884  }
1885 
1886 
1887  // Matrix<double> calculate_central_differences_Jacobian(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
1888 
1894 
1895  template<class T>
1897  {
1898  Vector<double> y = (t.*f)(x);
1899 
1900  double h;
1901 
1902  size_t n = x.size();
1903  size_t m = y.size();
1904 
1905  Vector<double> x_forward(x);
1906  Vector<double> x_backward(x);
1907 
1908  Vector<double> y_forward(n);
1909  Vector<double> y_backward(n);
1910 
1911  Matrix<double> J(m,n);
1912 
1913  for(size_t j = 0; j < n; j++)
1914  {
1915  h = calculate_h(x[j]);
1916 
1917  x_backward[j] -= h;
1918  y_backward = (t.*f)(x_backward);
1919  x_backward[j] += h;
1920 
1921  x_forward[j] += h;
1922  y_forward = (t.*f)(x_forward);
1923  x_forward[j] -= h;
1924 
1925  for(size_t i = 0; i < m; i++)
1926  {
1927  J(i,j) = (y_forward[i] - y_backward[i])/(2.0*h);
1928  }
1929  }
1930 
1931  return(J);
1932  }
1933 
1934 
1935  // Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
1936 
1942 
1943  template<class T>
1944  Matrix<double> calculate_Jacobian(const T& t, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>& x) const
1945  {
1947  {
1948  case ForwardDifferences:
1949  {
1950  return(calculate_forward_differences_Jacobian(t, f, x));
1951  }
1952  break;
1953 
1954  case CentralDifferences:
1955  {
1956  return(calculate_central_differences_Jacobian(t, f, x));
1957  }
1958  break;
1959 
1960  default:
1961  {
1962  std::ostringstream buffer;
1963 
1964  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
1965  << "Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method.\n"
1966  << "Unknown numerical differentiation method.\n";
1967 
1968  throw std::logic_error(buffer.str());
1969  }
1970  break;
1971  }
1972  }
1973 
1974 
1975  // Matrix<double> calculate_forward_differences_Jacobian(const T&, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const
1976 
1983 
1984  template<class T>
1986  {
1987  Vector<double> y = (t.*f)(dummy, x);
1988 
1989  double h;
1990 
1991  size_t n = x.size();
1992  size_t m = y.size();
1993 
1994  Vector<double> x_forward(x);
1995  Vector<double> y_forward(n);
1996 
1997  Matrix<double> J(m,n);
1998 
1999  for(size_t j = 0; j < n; j++)
2000  {
2001  h = calculate_h(x[j]);
2002 
2003  x_forward[j] += h;
2004  y_forward = (t.*f)(dummy, x_forward);
2005  x_forward[j] -= h;
2006 
2007  for(size_t i = 0; i < m; i++)
2008  {
2009  J(i,j) = (y_forward[i] - y[i])/h;
2010  }
2011  }
2012 
2013  return(J);
2014  }
2015 
2016 
2017  // Matrix<double> calculate_central_differences_Jacobian(const T&, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method
2018 
2025 
2026  template<class T>
2028  {
2029  Vector<double> y = (t.*f)(dummy, x);
2030 
2031  double h;
2032 
2033  size_t n = x.size();
2034  size_t m = y.size();
2035 
2036  Vector<double> x_forward(x);
2037  Vector<double> x_backward(x);
2038 
2039  Vector<double> y_forward(n);
2040  Vector<double> y_backward(n);
2041 
2042  Matrix<double> J(m,n);
2043 
2044  for(size_t j = 0; j < n; j++)
2045  {
2046  h = calculate_h(x[j]);
2047 
2048  x_backward[j] -= h;
2049  y_backward = (t.*f)(dummy, x_backward);
2050  x_backward[j] += h;
2051 
2052  x_forward[j] += h;
2053  y_forward = (t.*f)(dummy, x_forward);
2054  x_forward[j] -= h;
2055 
2056  for(size_t i = 0; i < m; i++)
2057  {
2058  J(i,j) = (y_forward[i] - y_backward[i])/(2.0*h);
2059  }
2060  }
2061 
2062  return(J);
2063  }
2064 
2065 
2066  // 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
2067 
2074 
2075  template<class T>
2076  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
2077  {
2079  {
2080  case ForwardDifferences:
2081  {
2082  return(calculate_forward_differences_Jacobian(t, f, dummy, x));
2083  }
2084  break;
2085 
2086  case CentralDifferences:
2087  {
2088  return(calculate_central_differences_Jacobian(t, f, dummy, x));
2089  }
2090  break;
2091 
2092  default:
2093  {
2094  std::ostringstream buffer;
2095 
2096  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
2097  << "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"
2098  << "Unknown numerical differentiation method.\n";
2099 
2100  throw std::logic_error(buffer.str());
2101  }
2102  break;
2103  }
2104  }
2105 
2106 
2107 
2108  // Matrix<double> calculate_forward_differences_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const
2109 
2117 
2118  template<class T>
2119  Matrix<double> calculate_forward_differences_Jacobian(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
2120  {
2121  Vector<double> y = (t.*f)(dummy, x);
2122 
2123  double h;
2124 
2125  size_t n = x.size();
2126  size_t m = y.size();
2127 
2128  Vector<double> x_forward(x);
2129  Vector<double> y_forward(n);
2130 
2131  Matrix<double> J(m,n);
2132 
2133  for(size_t j = 0; j < n; j++)
2134  {
2135  h = calculate_h(x[j]);
2136 
2137  x_forward[j] += h;
2138  y_forward = (t.*f)(dummy, x_forward);
2139  x_forward[j] -= h;
2140 
2141  for(size_t i = 0; i < m; i++)
2142  {
2143  J(i,j) = (y_forward[i] - y[i])/h;
2144  }
2145  }
2146 
2147  return(J);
2148  }
2149 
2150 
2151  // Matrix<double> calculate_central_differences_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const
2152 
2160 
2161  template<class T>
2162  Matrix<double> calculate_central_differences_Jacobian(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
2163  {
2164  Vector<double> y = (t.*f)(dummy, x);
2165 
2166  double h;
2167 
2168  size_t n = x.size();
2169  size_t m = y.size();
2170 
2171  Vector<double> x_forward(x);
2172  Vector<double> x_backward(x);
2173 
2174  Vector<double> y_forward(n);
2175  Vector<double> y_backward(n);
2176 
2177  Matrix<double> J(m,n);
2178 
2179  for(size_t j = 0; j < n; j++)
2180  {
2181  h = calculate_h(x[j]);
2182 
2183  x_backward[j] -= h;
2184  y_backward = (t.*f)(dummy, x_backward);
2185  x_backward[j] += h;
2186 
2187  x_forward[j] += h;
2188  y_forward = (t.*f)(dummy, x_forward);
2189  x_forward[j] -= h;
2190 
2191  for(size_t i = 0; i < m; i++)
2192  {
2193  J(i,j) = (y_forward[i] - y_backward[i])/(2.0*h);
2194  }
2195  }
2196 
2197  return(J);
2198  }
2199 
2200 
2201  // Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
2202 
2210 
2211  template<class T>
2212  Matrix<double> calculate_Jacobian(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
2213  {
2215  {
2216  case ForwardDifferences:
2217  {
2218  return(calculate_forward_differences_Jacobian(t, f, dummy, x));
2219  }
2220  break;
2221 
2222  case CentralDifferences:
2223  {
2224  return(calculate_central_differences_Jacobian(t, f, dummy, x));
2225  }
2226  break;
2227 
2228  default:
2229  {
2230  std::ostringstream buffer;
2231 
2232  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
2233  << "Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method.\n"
2234  << "Unknown numerical differentiation method.\n";
2235 
2236  throw std::logic_error(buffer.str());
2237  }
2238  break;
2239  }
2240  }
2241 
2242 
2243  // Matrix<double> calculate_forward_differences_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method
2244 
2253 
2254  template<class T>
2256  (const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
2257  {
2258  Vector<double> y = (t.*f)(dummy_int, dummy_vector, x);
2259 
2260  double h;
2261 
2262  size_t n = x.size();
2263  size_t m = y.size();
2264 
2265  Vector<double> x_forward(x);
2266  Vector<double> y_forward(n);
2267 
2268  Matrix<double> J(m,n);
2269 
2270  for(size_t j = 0; j < n; j++)
2271  {
2272  h = calculate_h(x[j]);
2273 
2274  x_forward[j] += h;
2275  y_forward = (t.*f)(dummy_int, dummy_vector, x_forward);
2276  x_forward[j] -= h;
2277 
2278  for(size_t i = 0; i < m; i++)
2279  {
2280  J(i,j) = (y_forward[i] - y[i])/h;
2281  }
2282  }
2283 
2284  return(J);
2285  }
2286 
2287 
2288  // Matrix<double> calculate_central_differences_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method
2289 
2298 
2299  template<class T>
2301  (const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
2302  {
2303  const size_t n = x.size();
2304 
2305  const Vector<double> y = (t.*f)(dummy_int, dummy_vector, x);
2306  const size_t m = y.size();
2307 
2308  double h;
2309 
2310  Vector<double> x_forward(x);
2311  Vector<double> x_backward(x);
2312 
2313  Vector<double> y_forward(n);
2314  Vector<double> y_backward(n);
2315 
2316  Matrix<double> J(m,n);
2317 
2318  for(size_t j = 0; j < n; j++)
2319  {
2320  h = calculate_h(x[j]);
2321 
2322  x_backward[j] -= h;
2323  y_backward = (t.*f)(dummy_int, dummy_vector, x_backward);
2324  x_backward[j] += h;
2325 
2326  x_forward[j] += h;
2327  y_forward = (t.*f)(dummy_int, dummy_vector, x_forward);
2328  x_forward[j] -= h;
2329 
2330  for(size_t i = 0; i < m; i++)
2331  {
2332  J(i,j) = (y_forward[i] - y_backward[i])/(2.0*h);
2333  }
2334  }
2335 
2336  return(J);
2337  }
2338 
2339 
2340  // Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method
2341 
2350 
2351  template<class T>
2353  (const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
2354  {
2356  {
2357  case ForwardDifferences:
2358  {
2359  return(calculate_forward_differences_Jacobian(t, f, dummy_int, dummy_vector, x));
2360  }
2361  break;
2362 
2363  case CentralDifferences:
2364  {
2365  return(calculate_central_differences_Jacobian(t, f, dummy_int, dummy_vector, x));
2366  }
2367  break;
2368 
2369  default:
2370  {
2371  std::ostringstream buffer;
2372 
2373  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
2374  << "Matrix<double> calculate_Jacobian\n"
2375  << "(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method.\n"
2376  << "Unknown numerical differentiation method.\n";
2377 
2378  throw std::logic_error(buffer.str());
2379  }
2380  break;
2381  }
2382  }
2383 
2384 
2385  // Matrix<double> calculate_forward_differences_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method
2386 
2395 
2396  template<class T>
2398  (const T& t, Vector<double> (T::*f)(const size_t&, const size_t&, const Vector<double>&) const, const size_t& dummy_int_1, const size_t& dummy_int_2, const Vector<double>& x) const
2399  {
2400  const Vector<double> y = (t.*f)(dummy_int_1, dummy_int_2, x);
2401 
2402  const size_t n = x.size();
2403  const size_t m = y.size();
2404 
2405  double h;
2406 
2407  Vector<double> x_forward(x);
2408  Vector<double> y_forward(n);
2409 
2410  Matrix<double> J(m,n);
2411 
2412  for(size_t j = 0; j < n; j++)
2413  {
2414  h = calculate_h(x[j]);
2415 
2416  x_forward[j] += h;
2417  y_forward = (t.*f)(dummy_int_1, dummy_int_2, x_forward);
2418  x_forward[j] -= h;
2419 
2420  for(size_t i = 0; i < m; i++)
2421  {
2422  J(i,j) = (y_forward[i] - y[i])/h;
2423  }
2424  }
2425 
2426  return(J);
2427  }
2428 
2429 
2430  // Matrix<double> calculate_central_differences_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method
2431 
2440 
2441  template<class T>
2443  (const T& t, Vector<double> (T::*f)(const size_t&, const size_t&, const Vector<double>&) const, const size_t& dummy_int_1, const size_t& dummy_int_2, const Vector<double>& x) const
2444  {
2445  const Vector<double> y = (t.*f)(dummy_int_1, dummy_int_2, x);
2446 
2447  const size_t n = x.size();
2448  const size_t m = y.size();
2449 
2450  double h;
2451 
2452  Vector<double> x_forward(x);
2453  Vector<double> x_backward(x);
2454 
2455  Vector<double> y_forward(n);
2456  Vector<double> y_backward(n);
2457 
2458  Matrix<double> J(m,n);
2459 
2460  for(size_t j = 0; j < n; j++)
2461  {
2462  h = calculate_h(x[j]);
2463 
2464  x_backward[j] -= h;
2465  y_backward = (t.*f)(dummy_int_1, dummy_int_2, x_backward);
2466  x_backward[j] += h;
2467 
2468  x_forward[j] += h;
2469  y_forward = (t.*f)(dummy_int_1, dummy_int_2, x_forward);
2470  x_forward[j] -= h;
2471 
2472  for(size_t i = 0; i < m; i++)
2473  {
2474  J(i,j) = (y_forward[i] - y_backward[i])/(2.0*h);
2475  }
2476  }
2477 
2478  return(J);
2479  }
2480 
2481 
2482  // Matrix<double> calculate_Jacobian(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method
2483 
2492 
2493  template<class T>
2495  (const T& t, Vector<double> (T::*f)(const size_t&, const size_t&, const Vector<double>&) const, const size_t& dummy_int_1, const size_t& dummy_int_2, const Vector<double>& x) const
2496  {
2498  {
2499  case ForwardDifferences:
2500  {
2501  return(calculate_forward_differences_Jacobian(t, f, dummy_int_1, dummy_int_2, x));
2502  }
2503  break;
2504 
2505  case CentralDifferences:
2506  {
2507  return(calculate_central_differences_Jacobian(t, f, dummy_int_1, dummy_int_2, x));
2508  }
2509  break;
2510 
2511  default:
2512  {
2513  std::ostringstream buffer;
2514 
2515  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
2516  << "Matrix<double> calculate_Jacobian\n"
2517  << "(const T&, Vector<double> (T::*f)(const size_t&, const size_t&, const Vector<double>&) const, const size_t&, const size_t&, const Vector<double>&) const method.\n"
2518  << "Unknown numerical differentiation method.\n";
2519 
2520  throw std::logic_error(buffer.str());
2521  }
2522  break;
2523  }
2524  }
2525 
2526 
2527  // HESSIAN FORM METHODS
2528 
2529 
2530  // Vector< Matrix <double> > calculate_forward_differences_Hessian_form(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
2531 
2537 
2538  template<class T>
2540  {
2541  Vector<double> y = (t.*f)(x);
2542 
2543  size_t s = y.size();
2544  size_t n = x.size();
2545 
2546  double h_j;
2547  double h_k;
2548 
2549  Vector<double> x_forward_j(x);
2550  Vector<double> x_forward_2j(x);
2551 
2552  Vector<double> x_forward_k(x);
2553  Vector<double> x_forward_jk(x);
2554 
2555  Vector<double> y_forward_j(s);
2556  Vector<double> y_forward_2j(s);
2557 
2558  Vector<double> y_forward_k(s);
2559  Vector<double> y_forward_jk(s);
2560 
2561  Vector< Matrix<double> > H(s);
2562 
2563  for(size_t i = 0; i < s; i++)
2564  {
2565  H[i].set(n,n);
2566 
2567  for(size_t j = 0; j < n; j++)
2568  {
2569  h_j = calculate_h(x[j]);
2570 
2571  x_forward_j[j] += h_j;
2572  y_forward_j = (t.*f)(x_forward_j);
2573  x_forward_j[j] -= h_j;
2574 
2575  x_forward_2j[j] += 2.0*h_j;
2576  y_forward_2j = (t.*f)(x_forward_2j);
2577  x_forward_2j[j] -= 2.0*h_j;
2578 
2579  H[i](j,j) = (y_forward_2j[i] - 2.0*y_forward_j[i] + y[i])/pow(h_j, 2);
2580 
2581  for(size_t k = j; k < n; k++)
2582  {
2583  h_k = calculate_h(x[k]);
2584 
2585  x_forward_k[k] += h_k;
2586  y_forward_k = (t.*f)(x_forward_k);
2587  x_forward_k[k] -= h_k;
2588 
2589  x_forward_jk[j] += h_j;
2590  x_forward_jk[k] += h_k;
2591  y_forward_jk = (t.*f)(x_forward_jk);
2592  x_forward_jk[j] -= h_j;
2593  x_forward_jk[k] -= h_k;
2594 
2595  H[i](j,k) = (y_forward_jk[i] - y_forward_j[i] - y_forward_k[i] + y[i])/(h_j*h_k);
2596  }
2597  }
2598 
2599  for(size_t j = 0; j < n; j++)
2600  {
2601  for(size_t k = 0; k < j; k++)
2602  {
2603  H[i](j,k) = H[i](k,j);
2604  }
2605  }
2606  }
2607 
2608  return(H);
2609  }
2610 
2611 
2612  // Vector< Matrix <double> > calculate_central_differences_Hessian_form(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
2613 
2619 
2620  template<class T>
2622  {
2623  Vector<double> y = (t.*f)(x);
2624 
2625  size_t s = y.size();
2626  size_t n = x.size();
2627 
2628  double h_j;
2629  double h_k;
2630 
2631  Vector<double> x_backward_2j(x);
2632  Vector<double> x_backward_j(x);
2633 
2634  Vector<double> x_forward_j(x);
2635  Vector<double> x_forward_2j(x);
2636 
2637  Vector<double> x_backward_jk(x);
2638  Vector<double> x_forward_jk(x);
2639 
2640  Vector<double> x_backward_j_forward_k(x);
2641  Vector<double> x_forward_j_backward_k(x);
2642 
2643  Vector<double> y_backward_2j;
2644  Vector<double> y_backward_j;
2645 
2646  Vector<double> y_forward_j;
2647  Vector<double> y_forward_2j;
2648 
2649  Vector<double> y_backward_jk;
2650  Vector<double> y_forward_jk;
2651 
2652  Vector<double> y_backward_j_forward_k;
2653  Vector<double> y_forward_j_backward_k;
2654 
2655  Vector< Matrix<double> > H(s);
2656 
2657  for(size_t i = 0; i < s; i++)
2658  {
2659  H[i].set(n,n);
2660 
2661  for(size_t j = 0; j < n; j++)
2662  {
2663  h_j = calculate_h(x[j]);
2664 
2665  x_backward_2j[j] -= 2.0*h_j;
2666  y_backward_2j = (t.*f)(x_backward_2j);
2667  x_backward_2j[j] += 2.0*h_j;
2668 
2669  x_backward_j[j] -= h_j;
2670  y_backward_j = (t.*f)(x_backward_j);
2671  x_backward_j[j] += h_j;
2672 
2673  x_forward_j[j] += h_j;
2674  y_forward_j = (t.*f)(x_forward_j);
2675  x_forward_j[j] -= h_j;
2676 
2677  x_forward_2j[j] += 2.0*h_j;
2678  y_forward_2j = (t.*f)(x_forward_2j);
2679  x_forward_2j[j] -= 2.0*h_j;
2680 
2681  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));
2682 
2683  for(size_t k = j; k < n; k++)
2684  {
2685  h_k = calculate_h(x[k]);
2686 
2687  x_backward_jk[j] -= h_j;
2688  x_backward_jk[k] -= h_k;
2689  y_backward_jk = (t.*f)(x_backward_jk);
2690  x_backward_jk[j] += h_j;
2691  x_backward_jk[k] += h_k;
2692 
2693  x_forward_jk[j] += h_j;
2694  x_forward_jk[k] += h_k;
2695  y_forward_jk = (t.*f)(x_forward_jk);
2696  x_forward_jk[j] -= h_j;
2697  x_forward_jk[k] -= h_k;
2698 
2699  x_backward_j_forward_k[j] -= h_j;
2700  x_backward_j_forward_k[k] += h_k;
2701  y_backward_j_forward_k = (t.*f)(x_backward_j_forward_k);
2702  x_backward_j_forward_k[j] += h_j;
2703  x_backward_j_forward_k[k] -= h_k;
2704 
2705  x_forward_j_backward_k[j] += h_j;
2706  x_forward_j_backward_k[k] -= h_k;
2707  y_forward_j_backward_k = (t.*f)(x_forward_j_backward_k);
2708  x_forward_j_backward_k[j] -= h_j;
2709  x_forward_j_backward_k[k] += h_k;
2710 
2711  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);
2712  }
2713  }
2714 
2715  for(size_t j = 0; j < n; j++)
2716  {
2717  for(size_t k = 0; k < j; k++)
2718  {
2719  H[i](j,k) = H[i](k,j);
2720  }
2721  }
2722  }
2723 
2724  return(H);
2725  }
2726 
2727 
2728  // Vector< Matrix<double> > calculate_Hessian_form(const T&, Vector<double> (T::*f)(const Vector<double>&) const, const Vector<double>&) const method
2729 
2735 
2736  template<class T>
2738  {
2740  {
2741  case ForwardDifferences:
2742  {
2744  }
2745  break;
2746 
2747  case CentralDifferences:
2748  {
2750  }
2751  break;
2752 
2753  default:
2754  {
2755  std::ostringstream buffer;
2756 
2757  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
2758  << "double calculate_Hessian_form(const T& t, Vector<double> (T::*f)(const Vector<double>&), const Vector<double>&) const method.\n"
2759  << "Unknown numerical differentiation method.\n";
2760 
2761  throw std::logic_error(buffer.str());
2762  }
2763  break;
2764  }
2765  }
2766 
2767 
2768  // Vector< Matrix <double> > calculate_forward_differences_Hessian_form
2769  // (const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const
2770 
2778 
2779  template<class T>
2781  (const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy_vector, const Vector<double>& x) const
2782  {
2783  Vector<double> y = (t.*f)(dummy_vector, x);
2784 
2785  size_t s = y.size();
2786  size_t n = x.size();
2787 
2788  double h_j;
2789  double h_k;
2790 
2791  Vector<double> x_forward_j(x);
2792  Vector<double> x_forward_2j(x);
2793 
2794  Vector<double> x_forward_k(x);
2795  Vector<double> x_forward_jk(x);
2796 
2797  Vector<double> y_forward_j(s);
2798  Vector<double> y_forward_2j(s);
2799 
2800  Vector<double> y_forward_k(s);
2801  Vector<double> y_forward_jk(s);
2802 
2803  Vector< Matrix<double> > H(s);
2804 
2805  for(size_t i = 0; i < s; i++)
2806  {
2807  H[i].set(n,n);
2808 
2809  for(size_t j = 0; j < n; j++)
2810  {
2811  h_j = calculate_h(x[j]);
2812 
2813  x_forward_j[j] += h_j;
2814  y_forward_j = (t.*f)(dummy_vector, x_forward_j);
2815  x_forward_j[j] -= h_j;
2816 
2817  x_forward_2j[j] += 2.0*h_j;
2818  y_forward_2j = (t.*f)(dummy_vector, x_forward_2j);
2819  x_forward_2j[j] -= 2.0*h_j;
2820 
2821  H[i](j,j) = (y_forward_2j[i] - 2.0*y_forward_j[i] + y[i])/pow(h_j, 2);
2822 
2823  for(size_t k = j; k < n; k++)
2824  {
2825  h_k = calculate_h(x[k]);
2826 
2827  x_forward_k[k] += h_k;
2828  y_forward_k = (t.*f)(dummy_vector, x_forward_k);
2829  x_forward_k[k] -= h_k;
2830 
2831  x_forward_jk[j] += h_j;
2832  x_forward_jk[k] += h_k;
2833  y_forward_jk = (t.*f)(dummy_vector, x_forward_jk);
2834  x_forward_jk[j] -= h_j;
2835  x_forward_jk[k] -= h_k;
2836 
2837  H[i](j,k) = (y_forward_jk[i] - y_forward_j[i] - y_forward_k[i] + y[i])/(h_j*h_k);
2838  }
2839  }
2840 
2841  for(size_t j = 0; j < n; j++)
2842  {
2843  for(size_t k = 0; k < j; k++)
2844  {
2845  H[i](j,k) = H[i](k,j);
2846  }
2847  }
2848  }
2849 
2850  return(H);
2851  }
2852 
2853 
2854  // Vector< Matrix <double> > calculate_central_differences_Hessian_form
2855  // (const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method
2856 
2864 
2865  template<class T>
2867  (const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy_vector, const Vector<double>& x) const
2868  {
2869  Vector<double> y = (t.*f)(dummy_vector, x);
2870 
2871  size_t s = y.size();
2872  size_t n = x.size();
2873 
2874  double h_j;
2875  double h_k;
2876 
2877  Vector<double> x_backward_2j(x);
2878  Vector<double> x_backward_j(x);
2879 
2880  Vector<double> x_forward_j(x);
2881  Vector<double> x_forward_2j(x);
2882 
2883  Vector<double> x_backward_jk(x);
2884  Vector<double> x_forward_jk(x);
2885 
2886  Vector<double> x_backward_j_forward_k(x);
2887  Vector<double> x_forward_j_backward_k(x);
2888 
2889  Vector<double> y_backward_2j;
2890  Vector<double> y_backward_j;
2891 
2892  Vector<double> y_forward_j;
2893  Vector<double> y_forward_2j;
2894 
2895  Vector<double> y_backward_jk;
2896  Vector<double> y_forward_jk;
2897 
2898  Vector<double> y_backward_j_forward_k;
2899  Vector<double> y_forward_j_backward_k;
2900 
2901  Vector< Matrix<double> > H(s);
2902 
2903  for(size_t i = 0; i < s; i++)
2904  {
2905  H[i].set(n,n);
2906 
2907  for(size_t j = 0; j < n; j++)
2908  {
2909  h_j = calculate_h(x[j]);
2910 
2911  x_backward_2j[j] -= 2.0*h_j;
2912  y_backward_2j = (t.*f)(dummy_vector, x_backward_2j);
2913  x_backward_2j[j] += 2.0*h_j;
2914 
2915  x_backward_j[j] -= h_j;
2916  y_backward_j = (t.*f)(dummy_vector, x_backward_j);
2917  x_backward_j[j] += h_j;
2918 
2919  x_forward_j[j] += h_j;
2920  y_forward_j = (t.*f)(dummy_vector, x_forward_j);
2921  x_forward_j[j] -= h_j;
2922 
2923  x_forward_2j[j] += 2.0*h_j;
2924  y_forward_2j = (t.*f)(dummy_vector, x_forward_2j);
2925  x_forward_2j[j] -= 2.0*h_j;
2926 
2927  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));
2928 
2929  for(size_t k = j; k < n; k++)
2930  {
2931  h_k = calculate_h(x[k]);
2932 
2933  x_backward_jk[j] -= h_j;
2934  x_backward_jk[k] -= h_k;
2935  y_backward_jk = (t.*f)(dummy_vector, x_backward_jk);
2936  x_backward_jk[j] += h_j;
2937  x_backward_jk[k] += h_k;
2938 
2939  x_forward_jk[j] += h_j;
2940  x_forward_jk[k] += h_k;
2941  y_forward_jk = (t.*f)(dummy_vector, x_forward_jk);
2942  x_forward_jk[j] -= h_j;
2943  x_forward_jk[k] -= h_k;
2944 
2945  x_backward_j_forward_k[j] -= h_j;
2946  x_backward_j_forward_k[k] += h_k;
2947  y_backward_j_forward_k = (t.*f)(dummy_vector, x_backward_j_forward_k);
2948  x_backward_j_forward_k[j] += h_j;
2949  x_backward_j_forward_k[k] -= h_k;
2950 
2951  x_forward_j_backward_k[j] += h_j;
2952  x_forward_j_backward_k[k] -= h_k;
2953  y_forward_j_backward_k = (t.*f)(dummy_vector, x_forward_j_backward_k);
2954  x_forward_j_backward_k[j] -= h_j;
2955  x_forward_j_backward_k[k] += h_k;
2956 
2957  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);
2958  }
2959  }
2960 
2961  for(size_t j = 0; j < n; j++)
2962  {
2963  for(size_t k = 0; k < j; k++)
2964  {
2965  H[i](j,k) = H[i](k,j);
2966  }
2967  }
2968  }
2969 
2970  return(H);
2971  }
2972 
2973 
2974  // Vector< Matrix<double> > calculate_Hessian_form
2975  // (const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method
2976 
2984 
2985  template<class T>
2987  (const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>& dummy_vector, const Vector<double>& x) const
2988  {
2990  {
2991  case ForwardDifferences:
2992  {
2993  return(calculate_forward_differences_Hessian_form(t, f, dummy_vector, x));
2994  }
2995  break;
2996 
2997  case CentralDifferences:
2998  {
2999  return(calculate_central_differences_Hessian_form(t, f, dummy_vector, x));
3000  }
3001  break;
3002 
3003  default:
3004  {
3005  std::ostringstream buffer;
3006 
3007  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
3008  << "Vector< Matrix<double> > calculate_Hessian_form\n"
3009  << "(const T& t, Vector<double> (T::*f)(const Vector<double>&, const Vector<double>&) const, const Vector<double>&, const Vector<double>&) const method.\n"
3010  << "Unknown numerical differentiation method.\n";
3011 
3012  throw std::logic_error(buffer.str());
3013  }
3014  break;
3015  }
3016  }
3017 
3018 
3019  // Vector< Matrix <double> > calculate_forward_differences_Hessian_form(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
3020 
3028 
3029  template<class T>
3030  Vector< Matrix <double> > calculate_forward_differences_Hessian_form(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
3031  {
3032  Vector<double> y = (t.*f)(dummy, x);
3033 
3034  size_t s = y.size();
3035  size_t n = x.size();
3036 
3037  double h_j;
3038  double h_k;
3039 
3040  Vector<double> x_forward_j(x);
3041  Vector<double> x_forward_2j(x);
3042 
3043  Vector<double> x_forward_k(x);
3044  Vector<double> x_forward_jk(x);
3045 
3046  Vector<double> y_forward_j(s);
3047  Vector<double> y_forward_2j(s);
3048 
3049  Vector<double> y_forward_k(s);
3050  Vector<double> y_forward_jk(s);
3051 
3052  Vector< Matrix<double> > H(s);
3053 
3054  for(size_t i = 0; i < s; i++)
3055  {
3056  H[i].set(n,n);
3057 
3058  for(size_t j = 0; j < n; j++)
3059  {
3060  h_j = calculate_h(x[j]);
3061 
3062  x_forward_j[j] += h_j;
3063  y_forward_j = (t.*f)(dummy, x_forward_j);
3064  x_forward_j[j] -= h_j;
3065 
3066  x_forward_2j[j] += 2.0*h_j;
3067  y_forward_2j = (t.*f)(dummy, x_forward_2j);
3068  x_forward_2j[j] -= 2.0*h_j;
3069 
3070  H[i](j,j) = (y_forward_2j[i] - 2.0*y_forward_j[i] + y[i])/pow(h_j, 2);
3071 
3072  for(size_t k = j; k < n; k++)
3073  {
3074  h_k = calculate_h(x[k]);
3075 
3076  x_forward_k[k] += h_k;
3077  y_forward_k = (t.*f)(dummy, x_forward_k);
3078  x_forward_k[k] -= h_k;
3079 
3080  x_forward_jk[j] += h_j;
3081  x_forward_jk[k] += h_k;
3082  y_forward_jk = (t.*f)(dummy, x_forward_jk);
3083  x_forward_jk[j] -= h_j;
3084  x_forward_jk[k] -= h_k;
3085 
3086  H[i](j,k) = (y_forward_jk[i] - y_forward_j[i] - y_forward_k[i] + y[i])/(h_j*h_k);
3087  }
3088  }
3089 
3090  for(size_t j = 0; j < n; j++)
3091  {
3092  for(size_t k = 0; k < j; k++)
3093  {
3094  H[i](j,k) = H[i](k,j);
3095  }
3096  }
3097  }
3098 
3099  return(H);
3100  }
3101 
3102 
3103  // Vector< Matrix <double> > calculate_central_differences_Hessian_form(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
3104 
3112 
3113  template<class T>
3114  Vector< Matrix <double> > calculate_central_differences_Hessian_form(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
3115  {
3116  Vector<double> y = (t.*f)(dummy, x);
3117 
3118  size_t s = y.size();
3119  size_t n = x.size();
3120 
3121  double h_j;
3122  double h_k;
3123 
3124  Vector<double> x_backward_2j(x);
3125  Vector<double> x_backward_j(x);
3126 
3127  Vector<double> x_forward_j(x);
3128  Vector<double> x_forward_2j(x);
3129 
3130  Vector<double> x_backward_jk(x);
3131  Vector<double> x_forward_jk(x);
3132 
3133  Vector<double> x_backward_j_forward_k(x);
3134  Vector<double> x_forward_j_backward_k(x);
3135 
3136  Vector<double> y_backward_2j;
3137  Vector<double> y_backward_j;
3138 
3139  Vector<double> y_forward_j;
3140  Vector<double> y_forward_2j;
3141 
3142  Vector<double> y_backward_jk;
3143  Vector<double> y_forward_jk;
3144 
3145  Vector<double> y_backward_j_forward_k;
3146  Vector<double> y_forward_j_backward_k;
3147 
3148  Vector< Matrix<double> > H(s);
3149 
3150  for(size_t i = 0; i < s; i++)
3151  {
3152  H[i].set(n,n);
3153 
3154  for(size_t j = 0; j < n; j++)
3155  {
3156  h_j = calculate_h(x[j]);
3157 
3158  x_backward_2j[j] -= 2.0*h_j;
3159  y_backward_2j = (t.*f)(dummy, x_backward_2j);
3160  x_backward_2j[j] += 2.0*h_j;
3161 
3162  x_backward_j[j] -= h_j;
3163  y_backward_j = (t.*f)(dummy, x_backward_j);
3164  x_backward_j[j] += h_j;
3165 
3166  x_forward_j[j] += h_j;
3167  y_forward_j = (t.*f)(dummy, x_forward_j);
3168  x_forward_j[j] -= h_j;
3169 
3170  x_forward_2j[j] += 2.0*h_j;
3171  y_forward_2j = (t.*f)(dummy, x_forward_2j);
3172  x_forward_2j[j] -= 2.0*h_j;
3173 
3174  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));
3175 
3176  for(size_t k = j; k < n; k++)
3177  {
3178  h_k = calculate_h(x[k]);
3179 
3180  x_backward_jk[j] -= h_j;
3181  x_backward_jk[k] -= h_k;
3182  y_backward_jk = (t.*f)(dummy, x_backward_jk);
3183  x_backward_jk[j] += h_j;
3184  x_backward_jk[k] += h_k;
3185 
3186  x_forward_jk[j] += h_j;
3187  x_forward_jk[k] += h_k;
3188  y_forward_jk = (t.*f)(dummy, x_forward_jk);
3189  x_forward_jk[j] -= h_j;
3190  x_forward_jk[k] -= h_k;
3191 
3192  x_backward_j_forward_k[j] -= h_j;
3193  x_backward_j_forward_k[k] += h_k;
3194  y_backward_j_forward_k = (t.*f)(dummy, x_backward_j_forward_k);
3195  x_backward_j_forward_k[j] += h_j;
3196  x_backward_j_forward_k[k] -= h_k;
3197 
3198  x_forward_j_backward_k[j] += h_j;
3199  x_forward_j_backward_k[k] -= h_k;
3200  y_forward_j_backward_k = (t.*f)(dummy, x_forward_j_backward_k);
3201  x_forward_j_backward_k[j] -= h_j;
3202  x_forward_j_backward_k[k] += h_k;
3203 
3204  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);
3205  }
3206  }
3207 
3208  for(size_t j = 0; j < n; j++)
3209  {
3210  for(size_t k = 0; k < j; k++)
3211  {
3212  H[i](j,k) = H[i](k,j);
3213  }
3214  }
3215  }
3216 
3217  return(H);
3218  }
3219 
3220 
3221  // Vector< Matrix<double> > calculate_Hessian_form(const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t&, const Vector<double>&) const method
3222 
3230 
3231  template<class T>
3232  Vector< Matrix<double> > calculate_Hessian_form(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&) const, const size_t& dummy, const Vector<double>& x) const
3233  {
3235  {
3236  case ForwardDifferences:
3237  {
3238  return(calculate_forward_differences_Hessian_form(t, f, dummy, x));
3239  }
3240  break;
3241 
3242  case CentralDifferences:
3243  {
3244  return(calculate_central_differences_Hessian_form(t, f, dummy, x));
3245  }
3246  break;
3247 
3248  default:
3249  {
3250  std::ostringstream buffer;
3251 
3252  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
3253  << "double calculate_Hessian_form(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&), const size_t&, const Vector<double>&) const method.\n"
3254  << "Unknown numerical differentiation method.\n";
3255 
3256  throw std::logic_error(buffer.str());
3257  }
3258  break;
3259  }
3260  }
3261 
3262 
3263  // Vector< Matrix <double> > calculate_forward_differences_Hessian_form
3264  // (const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const
3265 
3274 
3275  template<class T>
3277  (const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
3278  {
3279  Vector<double> y = (t.*f)(dummy_int, dummy_vector, x);
3280 
3281  size_t s = y.size();
3282  size_t n = x.size();
3283 
3284  double h_j;
3285  double h_k;
3286 
3287  Vector<double> x_forward_j(x);
3288  Vector<double> x_forward_2j(x);
3289 
3290  Vector<double> x_forward_k(x);
3291  Vector<double> x_forward_jk(x);
3292 
3293  Vector<double> y_forward_j(s);
3294  Vector<double> y_forward_2j(s);
3295 
3296  Vector<double> y_forward_k(s);
3297  Vector<double> y_forward_jk(s);
3298 
3299  Vector< Matrix<double> > H(s);
3300 
3301  for(size_t i = 0; i < s; i++)
3302  {
3303  H[i].set(n,n);
3304 
3305  for(size_t j = 0; j < n; j++)
3306  {
3307  h_j = calculate_h(x[j]);
3308 
3309  x_forward_j[j] += h_j;
3310  y_forward_j = (t.*f)(dummy_int, dummy_vector, x_forward_j);
3311  x_forward_j[j] -= h_j;
3312 
3313  x_forward_2j[j] += 2.0*h_j;
3314  y_forward_2j = (t.*f)(dummy_int, dummy_vector, x_forward_2j);
3315  x_forward_2j[j] -= 2.0*h_j;
3316 
3317  H[i](j,j) = (y_forward_2j[i] - 2.0*y_forward_j[i] + y[i])/pow(h_j, 2);
3318 
3319  for(size_t k = j; k < n; k++)
3320  {
3321  h_k = calculate_h(x[k]);
3322 
3323  x_forward_k[k] += h_k;
3324  y_forward_k = (t.*f)(dummy_int, dummy_vector, x_forward_k);
3325  x_forward_k[k] -= h_k;
3326 
3327  x_forward_jk[j] += h_j;
3328  x_forward_jk[k] += h_k;
3329  y_forward_jk = (t.*f)(dummy_int, dummy_vector, x_forward_jk);
3330  x_forward_jk[j] -= h_j;
3331  x_forward_jk[k] -= h_k;
3332 
3333  H[i](j,k) = (y_forward_jk[i] - y_forward_j[i] - y_forward_k[i] + y[i])/(h_j*h_k);
3334  }
3335  }
3336 
3337  for(size_t j = 0; j < n; j++)
3338  {
3339  for(size_t k = 0; k < j; k++)
3340  {
3341  H[i](j,k) = H[i](k,j);
3342  }
3343  }
3344  }
3345 
3346  return(H);
3347  }
3348 
3349 
3350  // Vector< Matrix <double> > calculate_central_differences_Hessian_form
3351  // (const T&, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method
3352 
3361 
3362  template<class T>
3364  (const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
3365  {
3366  const Vector<double> y = (t.*f)(dummy_int, dummy_vector, x);
3367 
3368  size_t s = y.size();
3369  size_t n = x.size();
3370 
3371  double h_j;
3372  double h_k;
3373 
3374  Vector<double> x_backward_2j(x);
3375  Vector<double> x_backward_j(x);
3376 
3377  Vector<double> x_forward_j(x);
3378  Vector<double> x_forward_2j(x);
3379 
3380  Vector<double> x_backward_jk(x);
3381  Vector<double> x_forward_jk(x);
3382 
3383  Vector<double> x_backward_j_forward_k(x);
3384  Vector<double> x_forward_j_backward_k(x);
3385 
3386  Vector<double> y_backward_2j;
3387  Vector<double> y_backward_j;
3388 
3389  Vector<double> y_forward_j;
3390  Vector<double> y_forward_2j;
3391 
3392  Vector<double> y_backward_jk;
3393  Vector<double> y_forward_jk;
3394 
3395  Vector<double> y_backward_j_forward_k;
3396  Vector<double> y_forward_j_backward_k;
3397 
3398  Vector< Matrix<double> > H(s);
3399 
3400  for(size_t i = 0; i < s; i++)
3401  {
3402  H[i].set(n,n);
3403 
3404  for(size_t j = 0; j < n; j++)
3405  {
3406  h_j = calculate_h(x[j]);
3407 
3408  x_backward_2j[j] -= 2.0*h_j;
3409  y_backward_2j = (t.*f)(dummy_int, dummy_vector, x_backward_2j);
3410  x_backward_2j[j] += 2.0*h_j;
3411 
3412  x_backward_j[j] -= h_j;
3413  y_backward_j = (t.*f)(dummy_int, dummy_vector, x_backward_j);
3414  x_backward_j[j] += h_j;
3415 
3416  x_forward_j[j] += h_j;
3417  y_forward_j = (t.*f)(dummy_int, dummy_vector, x_forward_j);
3418  x_forward_j[j] -= h_j;
3419 
3420  x_forward_2j[j] += 2.0*h_j;
3421  y_forward_2j = (t.*f)(dummy_int, dummy_vector, x_forward_2j);
3422  x_forward_2j[j] -= 2.0*h_j;
3423 
3424  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));
3425 
3426  for(size_t k = j; k < n; k++)
3427  {
3428  h_k = calculate_h(x[k]);
3429 
3430  x_backward_jk[j] -= h_j;
3431  x_backward_jk[k] -= h_k;
3432  y_backward_jk = (t.*f)(dummy_int, dummy_vector, x_backward_jk);
3433  x_backward_jk[j] += h_j;
3434  x_backward_jk[k] += h_k;
3435 
3436  x_forward_jk[j] += h_j;
3437  x_forward_jk[k] += h_k;
3438  y_forward_jk = (t.*f)(dummy_int, dummy_vector, x_forward_jk);
3439  x_forward_jk[j] -= h_j;
3440  x_forward_jk[k] -= h_k;
3441 
3442  x_backward_j_forward_k[j] -= h_j;
3443  x_backward_j_forward_k[k] += h_k;
3444  y_backward_j_forward_k = (t.*f)(dummy_int, dummy_vector, x_backward_j_forward_k);
3445  x_backward_j_forward_k[j] += h_j;
3446  x_backward_j_forward_k[k] -= h_k;
3447 
3448  x_forward_j_backward_k[j] += h_j;
3449  x_forward_j_backward_k[k] -= h_k;
3450  y_forward_j_backward_k = (t.*f)(dummy_int, dummy_vector, x_forward_j_backward_k);
3451  x_forward_j_backward_k[j] -= h_j;
3452  x_forward_j_backward_k[k] += h_k;
3453 
3454  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);
3455  }
3456  }
3457 
3458  for(size_t j = 0; j < n; j++)
3459  {
3460  for(size_t k = 0; k < j; k++)
3461  {
3462  H[i](j,k) = H[i](k,j);
3463  }
3464  }
3465  }
3466 
3467  return(H);
3468  }
3469 
3470 
3471  // Vector< Matrix<double> > calculate_Hessian_form
3472  // (const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method
3473 
3482 
3483  template<class T>
3485  (const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t& dummy_int, const Vector<double>& dummy_vector, const Vector<double>& x) const
3486  {
3488  {
3489  case ForwardDifferences:
3490  {
3491  return(calculate_forward_differences_Hessian_form(t, f, dummy_int, dummy_vector, x));
3492  }
3493  break;
3494 
3495  case CentralDifferences:
3496  {
3497  return(calculate_central_differences_Hessian_form(t, f, dummy_int, dummy_vector, x));
3498  }
3499  break;
3500 
3501  default:
3502  {
3503  std::ostringstream buffer;
3504 
3505  buffer << "OpenNN Exception: NumericalDifferentiation class.\n"
3506  << "Vector< Matrix<double> > calculate_Hessian_form\n"
3507  << "(const T& t, Vector<double> (T::*f)(const size_t&, const Vector<double>&, const Vector<double>&) const, const size_t&, const Vector<double>&, const Vector<double>&) const method.\n"
3508  << "Unknown numerical differentiation method.\n";
3509 
3510  throw std::logic_error(buffer.str());
3511  }
3512  break;
3513  }
3514  }
3515 
3516 
3517 private:
3518 
3520 
3522 
3524 
3526 
3528 
3529  bool display;
3530 
3531 };
3532 
3533 }
3534 
3535 #endif
3536 
3537 
3538 // OpenNN: Open Neural Networks Library.
3539 // Copyright (c) 2005-2015 Roberto Lopez.
3540 //
3541 // This library is free software; you can redistribute it and/or
3542 // modify it under the terms of the GNU Lesser General Public
3543 // License as published by the Free Software Foundation; either
3544 // version 2.1 of the License, or any later version.
3545 //
3546 // This library is distributed in the hope that it will be useful,
3547 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3548 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3549 // Lesser General Public License for more details.
3550 
3551 // You should have received a copy of the GNU Lesser General Public
3552 // License along with this library; if not, write to the Free Software
3553 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
3554 
Matrix< double > calculate_forward_differences_Hessian(const T &t, double(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Matrix< double > calculate_Hessian(const T &t, double(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
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
Matrix< double > calculate_central_differences_Jacobian(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Vector< double > calculate_central_differences_gradient(const T &t, double(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
size_t precision_digits
Number of precision digits.
bool operator==(const NumericalDifferentiation &) const
Vector< double > calculate_central_differences_derivative(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Vector< Matrix< double > > calculate_Hessian_form(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Vector< double > calculate_central_differences_second_derivative(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Matrix< double > calculate_central_differences_Jacobian(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Matrix< double > calculate_Jacobian(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
double calculate_second_derivative(const T &t, double(T::*f)(const double &) const , const double &x) const
Vector< Matrix< double > > calculate_central_differences_Hessian_form(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
double calculate_h(const double &) const
Vector< double > calculate_gradient(const T &t, double(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Vector< Matrix< double > > calculate_forward_differences_Hessian_form(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Vector< Matrix< double > > calculate_forward_differences_Hessian_form(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
double calculate_forward_differences_second_derivative(const T &t, double(T::*f)(const double &) const, const double &x) const
Matrix< double > calculate_forward_differences_Jacobian(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Vector< double > calculate_forward_differences_second_derivative(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
void set(void)
Sets the size of a vector to zero.
Definition: vector.h:656
void set_numerical_differentiation_method(const NumericalDifferentiationMethod &)
Matrix< double > calculate_forward_differences_Hessian(const T &t, double(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Matrix< double > calculate_forward_differences_Jacobian(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Matrix< double > calculate_Jacobian(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Vector< double > calculate_central_differences_second_derivative(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
double calculate_forward_differences_derivative(const T &t, double(T::*f)(const double &) const, const double &x) const
Vector< Matrix< double > > calculate_Hessian_form(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
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
const bool & get_display(void) const
Returns the flag used by this class for displaying or not displaying warnings.
virtual ~NumericalDifferentiation(void)
Destructor.
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
NumericalDifferentiationMethod numerical_differentiation_method
Numerical differentiation method variable.
Vector< double > calculate_gradient(const T &t, double(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Vector< double > calculate_derivative(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
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
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
Vector< double > calculate_second_derivative(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Vector< double > calculate_forward_differences_second_derivative(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Vector< double > calculate_forward_differences_gradient(const T &t, double(T::*f)(const Vector< double > &), const Vector< double > &x) const
Matrix< double > calculate_central_differences_Hessian(const T &t, double(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
void set(const NumericalDifferentiation &)
Matrix< double > calculate_Hessian(const T &t, double(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Vector< double > calculate_central_differences_gradient(const T &t, double(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Vector< Matrix< double > > calculate_central_differences_Hessian_form(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
Matrix< double > calculate_central_differences_Hessian(const T &t, double(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
const size_t & get_precision_digits(void) const
Returns the number of precision digits required for the derivatives.
Vector< double > calculate_central_differences_gradient(const T &t, double(T::*f)(const Vector< double > &), const Vector< double > &x) const
Vector< double > calculate_forward_differences_derivative(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Vector< double > calculate_derivative(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
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
Vector< double > calculate_gradient(const T &t, double(T::*f)(const Vector< double > &), const Vector< double > &x) const
double calculate_central_differences_derivative(const T &t, double(T::*f)(const double &) const , const double &x) const
double calculate_derivative(const T &t, double(T::*f)(const double &) const , const double &x) const
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
std::string write_numerical_differentiation_method(void) const
Returns a string with the name of the method to be used for numerical differentiation.
Vector< double > calculate_forward_differences_gradient(const T &t, double(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
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
bool display
Flag for displaying warning messages from this class.
Vector< double > calculate_forward_differences_derivative(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
NumericalDifferentiationMethod
Enumeration of available methods for numerical differentiation.
Vector< double > calculate_forward_differences_gradient(const T &t, double(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
Vector< double > calculate_second_derivative(const T &t, Vector< double >(T::*f)(const size_t &, const Vector< double > &) const, const size_t &dummy, const Vector< double > &x) const
const NumericalDifferentiationMethod & get_numerical_differentiation_method(void) const
Returns the method used for numerical differentiation (forward differences or central differences)...
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
tinyxml2::XMLDocument * to_XML(void) const
Serializes this numerical differentiation object into a XML document->
double calculate_central_differences_second_derivative(const T &t, double(T::*f)(const double &) const , const double &x) const
void from_XML(const tinyxml2::XMLDocument &)
Vector< double > calculate_central_differences_derivative(const T &t, Vector< double >(T::*f)(const Vector< double > &) const, const Vector< double > &x) const
NumericalDifferentiation & operator=(const NumericalDifferentiation &)