numerical_differentiation.h
1// OpenNN: Open Neural Networks Library
2// www.opennn.net
3//
4// 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
5//
6// Artificial Intelligence Techniques SL
7// artelnics@artelnics.com
8
9#ifndef NUMERICALDIFFERENTIATION_H
10#define NUMERICALDIFFERENTIATION_H
11
12// System includes
13
14#include<iostream>
15#include<vector>
16#include<vector>
17#include<limits>
18#include<cstddef>
19
20// OpenNN includes
21
22#include "config.h"
23
24using namespace std;
25using namespace Eigen;
26
27namespace OpenNN
28{
29
32
33class NumericalDifferentiation
34{
35
36public:
37
38 // Constructors
39
40 explicit NumericalDifferentiation();
41
42 // Destructor
43
44 virtual ~NumericalDifferentiation();
45
47
48 const Index& get_precision_digits() const;
49
50 const bool& get_display() const;
51
52 void set(const NumericalDifferentiation&);
53
54 void set_precision_digits(const Index&);
55
56 void set_display(const bool&);
57
58 void set_default();
59
60 type calculate_eta() const;
61
62 type calculate_h(const type&) const;
63 Tensor<type, 1> calculate_h(const Tensor<type, 1>&) const;
64 Tensor<type, 2> calculate_h(const Tensor<type, 2>&) const;
65 Tensor<type, 4> calculate_h(const Tensor<type, 4>&) const;
66
67 template<class T>
68 type calculate_derivatives(const T& t, type(T::*f)(const type&) const , const type& x) const
69 {
70 const type h = calculate_h(x);
71
72 const type y_forward = (t.*f)(x+h);
73
74 const type y_backward = (t.*f)(x-h);
75
76 const type d = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
77
78 return d;
79 }
80
81
86
87 template<class T>
88 Tensor<type, 1> calculate_derivatives(const T& t, Tensor<type, 1>(T::*f)(const Tensor<type, 1>&) const, const Tensor<type, 1>& x) const
89 {
90
91 const Tensor<type, 1> h = calculate_h(x);
92
93 const Tensor<type, 1> x_forward = x + h;
94 const Tensor<type, 1> x_backward = x - h;
95
96 const Tensor<type, 1> y_forward = (t.*f)(x_forward);
97 const Tensor<type, 1> y_backward = (t.*f)(x_backward);
98
99 const Tensor<type, 1> y = (t.*f)(x);
100
101 const Tensor<type, 1> d = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
102
103 return d;
104 }
105
106
107 template<class T>
108 Tensor<type, 2> calculate_derivatives(const T& t, Tensor<type, 2>(T::*f)(const Tensor<type, 2>&) const, const Tensor<type, 2>& x) const
109 {
110 const Tensor<type, 2> h = calculate_h(x);
111
112 const Tensor<type, 2> x_forward = x + h;
113 const Tensor<type, 2> x_backward = x - h;
114
115 const Tensor<type, 2> y_forward = (t.*f)(x_forward);
116 const Tensor<type, 2> y_backward = (t.*f)(x_backward);
117
118 const Tensor<type, 2> y = (t.*f)(x);
119
120 const Tensor<type, 2> d = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
121
122 return d;
123 }
124
125
126 template<class T>
127 Tensor<type, 4> calculate_derivatives(const T& t, Tensor<type, 4>(T::*f)(const Tensor<type, 4>&) const, const Tensor<type, 4>& x) const
128 {
129 const Tensor<type, 4> h = calculate_h(x);
130
131 const Tensor<type, 4> x_forward = x + h;
132 const Tensor<type, 4> x_backward = x - h;
133
134 const Tensor<type, 4> y_forward = (t.*f)(x_forward);
135 const Tensor<type, 4> y_backward = (t.*f)(x_backward);
136
137 const Tensor<type, 4> y = (t.*f)(x);
138
139 const Tensor<type, 4> d = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
140
141 return d;
142 }
143
144
151
152 template<class T>
153 Tensor<type, 1> calculate_derivatives(const T& t, Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 1>&) const, const Index& dummy, const Tensor<type, 1>& x) const
154 {
155 const Tensor<type, 1> h = calculate_h(x);
156
157 const Tensor<type, 1> x_forward = x + h;
158 const Tensor<type, 1> x_backward = x - h;
159
160 const Tensor<type, 1> y_forward = (t.*f)(dummy, x_forward);
161 const Tensor<type, 1> y_backward = (t.*f)(dummy, x_backward);
162
163 const Tensor<type, 1> d = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
164
165 return d;
166 }
167
168
169 template<class T>
170 Tensor<type, 2> calculate_derivatives(const T& t, void(T::*f)(const Tensor<type, 2>&, Tensor<type, 2>&) const, const Index& dummy, const Tensor<type, 2>& x) const
171 {
172 const Index rn = x.dimension(0);
173 const Index cn = x.dimension(1);
174
175 const Tensor<type, 2> h = calculate_h(x);
176
177 const Tensor<type, 2> x_forward = x + h;
178 const Tensor<type, 2> x_backward = x - h;
179
180 Tensor<type, 2> y_forward(rn,cn);
181 (t.*f)(x_forward, y_forward);
182 Tensor<type, 2> y_backward(rn,cn);
183 (t.*f)(x_backward, y_backward);
184
185 const Tensor<type, 2> d = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
186
187 return d;
188 }
189
190
191 template<class T>
192 Tensor<type, 4> calculate_derivatives(const T& t, void(T::*f)(const Tensor<type, 4>&, Tensor<type, 4>&) const, const Index& dummy, const Tensor<type, 4>& x) const
193 {
194 const Index rn = x.dimension(0);
195 const Index cn = x.dimension(1);
196 const Index kn = x.dimension(2);
197 const Index in = x.dimension(3);
198
199 const Tensor<type, 4> h = calculate_h(x);
200
201 const Tensor<type, 4> x_forward = x + h;
202 const Tensor<type, 4> x_backward = x - h;
203
204 Tensor<type, 4> y_forward(rn,cn, kn, in);
205 (t.*f)(x_forward, y_forward);
206 Tensor<type, 4> y_backward(rn,cn, kn, in);
207 (t.*f)(x_backward, y_backward);
208
209 const Tensor<type, 4> d = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
210
211 return d;
212 }
213
214
219
220 template<class T>
221 type calculate_second_derivatives(const T& t, type(T::*f)(const type&) const , const type& x) const
222 {
223 const type h = calculate_h(x);
224
225 const type x_forward_2 = x + static_cast<type>(2.0)*h;
226
227 const type y_forward_2 = (t.*f)(x_forward_2);
228
229 const type x_forward = x + h;
230
231 const type y_forward = (t.*f)(x_forward);
232
233 const type y = (t.*f)(x);
234
235 const type x_backward = x - h;
236
237 const type y_backward = (t.*f)(x_backward);
238
239 const type x_backward_2 = x - static_cast<type>(2.0)*h;
240
241 const type y_backward_2 = (t.*f)(x_backward_2);
242
243 const type d2 = (-y_forward_2 + type(16.0)*y_forward - type(30.0)*y + type(16.0)*y_backward - y_backward_2)/(type(12.0)*pow(h, type(2)));
244
245 return d2;
246 }
247
248
253
254 template<class T>
255 Tensor<type, 1> calculate_second_derivatives(const T& t, Tensor<type, 1>(T::*f)(const Tensor<type, 1>&) const, const Tensor<type, 1>& x) const
256 {
257 const Tensor<type, 1> h = calculate_h(x);
258
259 const Tensor<type, 1> x_forward = x + h;
260 const Tensor<type, 1> x_forward_2 = x + h*static_cast<type>(2.0);
261
262 const Tensor<type, 1> x_backward = x - h;
263 const Tensor<type, 1> x_backward_2 = x - h*static_cast<type>(2.0);
264
265 const Tensor<type, 1> y = (t.*f)(x);
266
267 const Tensor<type, 1> y_forward = (t.*f)(x_forward);
268 const Tensor<type, 1> y_forward_2 = (t.*f)(x_forward_2);
269
270 const Tensor<type, 1> y_backward = (t.*f)(x_backward);
271 const Tensor<type, 1> y_backward_2 = (t.*f)(x_backward_2);
272
273 return ((y_forward_2*type(-1) + y_forward*type(16) + y*type(-30) + y_backward*type(16) + y_backward_2*type(-1))/(h*h*type(12)));
274 }
275
276
283
284 template<class T>
285 Tensor<type, 1> calculate_second_derivatives(const T& t,
286 Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 1>&) const,
287 const Index& dummy, const Tensor<type, 1>& x) const
288 {
289 const Tensor<type, 1> h = calculate_h(x);
290
291 const Tensor<type, 1> x_forward = x + h;
292 const Tensor<type, 1> x_forward_2 = x + h*static_cast<type>(2.0);
293
294 const Tensor<type, 1> x_backward = x - h;
295 const Tensor<type, 1> x_backward_2 = x - h*static_cast<type>(2.0);
296
297 const Tensor<type, 1> y = (t.*f)(dummy, x);
298
299 const Tensor<type, 1> y_forward = (t.*f)(dummy, x_forward);
300 const Tensor<type, 1> y_forward_2 = (t.*f)(dummy, x_forward_2);
301
302 const Tensor<type, 1> y_backward = (t.*f)(dummy, x_backward);
303 const Tensor<type, 1> y_backward_2 = (t.*f)(dummy, x_backward_2);
304
305 return (y_forward_2*type(-1) + y_forward*type(16) + y*type(-30) + y_backward*type(16) + y_backward_2*type(-1))/(h*h*type(12));
306 }
307
308
314
315 template<class T>
316 Tensor<type, 1> calculate_gradient(const T& t, type (T::*f)(const Tensor<type, 1>&) const, const Tensor<type, 1>& x) const
317 {
318
319 const Index n = x.size();
320
321 type h;
322
323 Tensor<type, 1> x_forward(x);
324 Tensor<type, 1> x_backward(x);
325
326 type y_forward;
327 type y_backward;
328
329 Tensor<type, 1> g(n);
330
331 for(Index i = 0; i < n; i++)
332 {
333 h = calculate_h(x(i));
334
335 x_forward(i) += h;
336
337 y_forward = (t.*f)(x_forward);
338
339 x_forward(i) -= h;
340 x_backward(i) -= h;
341
342 y_backward = (t.*f)(x_backward);
343 x_backward(i) += h;
344
345 g(i) = (y_forward - y_backward)/(type(2.0)*h);
346
347 }
348
349 return g;
350 }
351
357
358 template<class T>
359 Tensor<type, 1> calculate_gradient(const T& t, type(T::*f)(const Tensor<type, 1>&), const Tensor<type, 1>& x) const
360 {
361 const Index n = x.size();
362
363 type h;
364
365 Tensor<type, 1> x_forward(x);
366 Tensor<type, 1> x_backward(x);
367
368 type y_forward;
369 type y_backward;
370
371 Tensor<type, 1> g(n);
372
373 for(Index i = 0; i < n; i++)
374 {
375 h = calculate_h(x(i));
376
377 x_forward(i) += h;
378 y_forward = (t.*f)(x_forward);
379 x_forward(i) -= h;
380
381 x_backward(i) -= h;
382 y_backward = (t.*f)(x_backward);
383 x_backward(i) += h;
384
385 g(i) = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
386 }
387
388 return g;
389 }
390
398
399 template<class T>
400 Tensor<type, 1> calculate_gradient(const T& t, type(T::*f)(const Tensor<type, 1>&, const Tensor<type, 1>&) const, const Tensor<type, 1>& dummy, const Tensor<type, 1>& x) const
401 {
402 const Index n = x.size();
403
404 type h;
405
406 Tensor<type, 1> x_forward(x);
407 Tensor<type, 1> x_backward(x);
408
409 type y_forward;
410 type y_backward;
411
412 Tensor<type, 1> g(n);
413
414 for(Index i = 0; i < n; i++)
415 {
416 h = calculate_h(x(i));
417
418 x_forward(i) += h;
419 y_forward = (t.*f)(dummy, x_forward);
420 x_forward(i) -= h;
421
422 x_backward(i) -= h;
423 y_backward = (t.*f)(dummy, x_backward);
424 x_backward(i) += h;
425
426 g(i) = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
427 }
428
429 return g;
430 }
431
439
440 template<class T>
441 Tensor<type, 1> calculate_gradient(const T& t, type(T::*f)(const Index&, const Tensor<type, 1>&) const, const Index& dummy, const Tensor<type, 1>& x) const
442 {
443 const Index n = x.size();
444
445 type h;
446
447 Tensor<type, 1> x_forward(x);
448 Tensor<type, 1> x_backward(x);
449
450 type y_forward;
451 type y_backward;
452
453 Tensor<type, 1> g(n);
454
455 for(Index i = 0; i < n; i++)
456 {
457 h = calculate_h(x(i));
458
459 x_forward(i) += h;
460 y_forward = (t.*f)(dummy, x_forward);
461 x_forward(i) -= h;
462
463 x_backward(i) -= h;
464 y_backward = (t.*f)(dummy, x_backward);
465 x_backward(i) += h;
466
467 g(i) = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
468 }
469
470 return g;
471 }
472
473
481
482 template<class T>
483 Tensor<type, 1> calculate_gradient(const T& t, Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 2>&) const, const Index& dummy, const Tensor<type, 2>& x) const
484 {
485 const Index n = x.size();
486
487 type h;
488
489 Tensor<type, 1> x_forward(x);
490 Tensor<type, 1> x_backward(x);
491
492 type y_forward;
493 type y_backward;
494
495 Tensor<type, 1> g(n);
496
497 for(Index i = 0; i < n; i++)
498 {
499 h = calculate_h(x(i));
500
501 x_forward(i) += h;
502 y_forward = (t.*f)(dummy, x_forward);
503 x_forward(i) -= h;
504
505 x_backward(i) -= h;
506 y_backward = (t.*f)(dummy, x_backward);
507 x_backward(i) += h;
508
509 g(i) = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
510 }
511
512 return g;
513 }
514
515
523
524 template<class T>
525 Tensor<type, 1> calculate_gradient(const T& t, type(T::*f)(const Tensor<Index, 1>&, const Tensor<type, 1>&) const, const Tensor<Index, 1>& dummy, const Tensor<type, 1>& x) const
526 {
527 const Index n = x.size();
528
529 type h;
530 Tensor<type, 1> x_forward(x);
531 Tensor<type, 1> x_backward(x);
532
533 type y_forward;
534 type y_backward;
535
536 Tensor<type, 1> g(n);
537
538 for(Index i = 0; i < n; i++)
539 {
540 h = calculate_h(x(i));
541
542 x_forward(i) += h;
543 y_forward = (t.*f)(dummy, x_forward);
544 x_forward(i) -= h;
545
546 x_backward(i) -= h;
547 y_backward = (t.*f)(dummy, x_backward);
548 x_backward(i) += h;
549
550 g(i) = (y_forward - y_backward)/(h*type(2));
551 }
552
553 return g;
554 }
555
556
557 template<class T>
558 Tensor<type, 2> calculate_gradient_matrix(const T& t, Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 2>&) const, const Index& integer, const Tensor<type, 2>& x) const
559 {
560 const Index rows_number = x.dimension(0);
561 const Index columns_number = x.dimension(1);
562
563 Tensor<type, 2> gradient(rows_number, columns_number);
564
565 type h;
566 Tensor<type, 2> x_forward(x);
567 Tensor<type, 2> x_backward(x);
568
569 type y_forward;
570 type y_backward;
571
572 for(Index i = 0; i < rows_number; i++)
573 {
574 for(Index j = 0; j < columns_number; j++)
575 {
576 h = calculate_h(x(i,j));
577
578 x_forward(i,j) += h;
579 y_forward = (t.*f)(integer, x_forward)(i);
580 x_forward(i,j) -= h;
581
582 x_backward(i,j) -= h;
583 y_backward = (t.*f)(integer, x_backward)(i);
584 x_backward(i,j) += h;
585
586 gradient(i,j) = (y_forward - y_backward)/(static_cast<type>(2.0)*h);
587 }
588 }
589
590 return gradient;
591 }
592
593
599
600 template<class T>
601 Tensor<type, 2> calculate_hessian(const T& t, type(T::*f)(const Tensor<type, 1>&) const, const Tensor<type, 1>& x) const
602 {
603 const Index n = x.size();
604
605 type y = (t.*f)(x);
606
607 Tensor<type, 2> H(n, n);
608
609 type h_i;
610 type h_j;
611
612 Tensor<type, 1> x_backward_2i(x);
613 Tensor<type, 1> x_backward_i(x);
614
615 Tensor<type, 1> x_forward_i(x);
616 Tensor<type, 1> x_forward_2i(x);
617
618 Tensor<type, 1> x_backward_ij(x);
619 Tensor<type, 1> x_forward_ij(x);
620
621 Tensor<type, 1> x_backward_i_forward_j(x);
622 Tensor<type, 1> x_forward_i_backward_j(x);
623
624 type y_backward_2i;
625 type y_backward_i;
626
627 type y_forward_i;
628 type y_forward_2i;
629
630 type y_backward_ij;
631 type y_forward_ij;
632
633 type y_backward_i_forward_j;
634 type y_forward_i_backward_j;
635
636 for(Index i = 0; i < n; i++)
637 {
638 h_i = calculate_h(x(i));
639
640 x_backward_2i(i) -= static_cast<type>(2.0)*h_i;
641 y_backward_2i = (t.*f)(x_backward_2i);
642 x_backward_2i(i) += static_cast<type>(2.0)*h_i;
643
644 x_backward_i(i) -= h_i;
645 y_backward_i = (t.*f)(x_backward_i);
646 x_backward_i(i) += h_i;
647
648 x_forward_i(i) += h_i;
649 y_forward_i = (t.*f)(x_forward_i);
650 x_forward_i(i) -= h_i;
651
652 x_forward_2i(i) += static_cast<type>(2.0)*h_i;
653 y_forward_2i = (t.*f)(x_forward_2i);
654 x_forward_2i(i) -= static_cast<type>(2.0)*h_i;
655
656 H(i,i) = (-y_forward_2i + type(16.0)*y_forward_i - type(30.0)*y + type(16.0)*y_backward_i - y_backward_2i)/(type(12.0)*pow(h_i, type(2)));
657
658 for(Index j = i; j < n; j++)
659 {
660 h_j = calculate_h(x(j));
661
662 x_backward_ij(i) -= h_i;
663 x_backward_ij(j) -= h_j;
664 y_backward_ij = (t.*f)(x_backward_ij);
665 x_backward_ij(i) += h_i;
666 x_backward_ij(j) += h_j;
667
668 x_forward_ij(i) += h_i;
669 x_forward_ij(j) += h_j;
670 y_forward_ij = (t.*f)(x_forward_ij);
671 x_forward_ij(i) -= h_i;
672 x_forward_ij(j) -= h_j;
673
674 x_backward_i_forward_j(i) -= h_i;
675 x_backward_i_forward_j(j) += h_j;
676 y_backward_i_forward_j = (t.*f)(x_backward_i_forward_j);
677 x_backward_i_forward_j(i) += h_i;
678 x_backward_i_forward_j(j) -= h_j;
679
680 x_forward_i_backward_j(i) += h_i;
681 x_forward_i_backward_j(j) -= h_j;
682 y_forward_i_backward_j = (t.*f)(x_forward_i_backward_j);
683 x_forward_i_backward_j(i) -= h_i;
684 x_forward_i_backward_j(j) += h_j;
685
686 H(i,j) = (y_forward_ij - y_forward_i_backward_j - y_backward_i_forward_j + y_backward_ij)/(type(4.0)*h_i*h_j);
687 }
688 }
689
690 for(Index i = 0; i < n; i++)
691 {
692 for(Index j = 0; j < i; j++)
693 {
694 H(i,j) = H(j,i);
695 }
696 }
697
698 return H;
699 }
700
701
709
710 template<class T>
711 Tensor<type, 2> calculate_hessian(const T& t, type(T::*f)(const Tensor<type, 1>&, const Tensor<type, 1>&) const, const Tensor<type, 1>& dummy, const Tensor<type, 1>& x) const
712 {
713 const Index n = x.size();
714
715 type y = (t.*f)(dummy, x);
716
717 Tensor<type, 2> H(n, n);
718
719 type h_i;
720 type h_j;
721
722 Tensor<type, 1> x_backward_2i(x);
723 Tensor<type, 1> x_backward_i(x);
724
725 Tensor<type, 1> x_forward_i(x);
726 Tensor<type, 1> x_forward_2i(x);
727
728 Tensor<type, 1> x_backward_ij(x);
729 Tensor<type, 1> x_forward_ij(x);
730
731 Tensor<type, 1> x_backward_i_forward_j(x);
732 Tensor<type, 1> x_forward_i_backward_j(x);
733
734 type y_backward_2i;
735 type y_backward_i;
736
737 type y_forward_i;
738 type y_forward_2i;
739
740 type y_backward_ij;
741 type y_forward_ij;
742
743 type y_backward_i_forward_j;
744 type y_forward_i_backward_j;
745
746 for(Index i = 0; i < n; i++)
747 {
748 h_i = calculate_h(x(i));
749
750 x_backward_2i(i) -= static_cast<type>(2.0)*h_i;
751 y_backward_2i = (t.*f)(dummy, x_backward_2i);
752 x_backward_2i(i) += static_cast<type>(2.0)*h_i;
753
754 x_backward_i(i) -= h_i;
755 y_backward_i = (t.*f)(dummy, x_backward_i);
756 x_backward_i(i) += h_i;
757
758 x_forward_i(i) += h_i;
759 y_forward_i = (t.*f)(dummy, x_forward_i);
760 x_forward_i(i) -= h_i;
761
762 x_forward_2i(i) += static_cast<type>(2.0)*h_i;
763 y_forward_2i = (t.*f)(dummy, x_forward_2i);
764 x_forward_2i(i) -= static_cast<type>(2.0)*h_i;
765
766 H(i,i) = (-y_forward_2i + type(16.0)*y_forward_i - type(30.0)*y + type(16.0)*y_backward_i - y_backward_2i)/(type(12.0)*pow(h_i, type(2)));
767
768 for(Index j = i; j < n; j++)
769 {
770 h_j = calculate_h(x(j));
771
772 x_backward_ij(i) -= h_i;
773 x_backward_ij(j) -= h_j;
774 y_backward_ij = (t.*f)(dummy, x_backward_ij);
775 x_backward_ij(i) += h_i;
776 x_backward_ij(j) += h_j;
777
778 x_forward_ij(i) += h_i;
779 x_forward_ij(j) += h_j;
780 y_forward_ij = (t.*f)(dummy, x_forward_ij);
781 x_forward_ij(i) -= h_i;
782 x_forward_ij(j) -= h_j;
783
784 x_backward_i_forward_j(i) -= h_i;
785 x_backward_i_forward_j(j) += h_j;
786 y_backward_i_forward_j = (t.*f)(dummy, x_backward_i_forward_j);
787 x_backward_i_forward_j(i) += h_i;
788 x_backward_i_forward_j(j) -= h_j;
789
790 x_forward_i_backward_j(i) += h_i;
791 x_forward_i_backward_j(j) -= h_j;
792 y_forward_i_backward_j = (t.*f)(dummy, x_forward_i_backward_j);
793 x_forward_i_backward_j(i) -= h_i;
794 x_forward_i_backward_j(j) += h_j;
795
796 H(i,j) = (y_forward_ij - y_forward_i_backward_j - y_backward_i_forward_j + y_backward_ij)/(type(4.0)*h_i*h_j);
797 }
798 }
799
800 for(Index i = 0; i < n; i++)
801 {
802 for(Index j = 0; j < i; j++)
803 {
804 H(i,j) = H(j,i);
805 }
806 }
807
808 return H;
809 }
810
818
819 template<class T>
820 Tensor<type, 2> calculate_hessian(const T& t, type(T::*f)(const Index&, const Tensor<type, 1>&) const, const Index& dummy, const Tensor<type, 1>& x) const
821 {
822 const Index n = x.size();
823
824 type y = (t.*f)(dummy, x);
825
826 Tensor<type, 2> H(n, n);
827
828 type h_i;
829 type h_j;
830
831 Tensor<type, 1> x_backward_2i(x);
832 Tensor<type, 1> x_backward_i(x);
833
834 Tensor<type, 1> x_forward_i(x);
835 Tensor<type, 1> x_forward_2i(x);
836
837 Tensor<type, 1> x_backward_ij(x);
838 Tensor<type, 1> x_forward_ij(x);
839
840 Tensor<type, 1> x_backward_i_forward_j(x);
841 Tensor<type, 1> x_forward_i_backward_j(x);
842
843 type y_backward_2i;
844 type y_backward_i;
845
846 type y_forward_i;
847 type y_forward_2i;
848
849 type y_backward_ij;
850 type y_forward_ij;
851
852 type y_backward_i_forward_j;
853 type y_forward_i_backward_j;
854
855 for(Index i = 0; i < n; i++)
856 {
857 h_i = calculate_h(x(i));
858
859 x_backward_2i(i) -= static_cast<type>(2.0)*h_i;
860 y_backward_2i = (t.*f)(dummy, x_backward_2i);
861 x_backward_2i(i) += static_cast<type>(2.0)*h_i;
862
863 x_backward_i(i) -= h_i;
864 y_backward_i = (t.*f)(dummy, x_backward_i);
865 x_backward_i(i) += h_i;
866
867 x_forward_i(i) += h_i;
868 y_forward_i = (t.*f)(dummy, x_forward_i);
869 x_forward_i(i) -= h_i;
870
871 x_forward_2i(i) += static_cast<type>(2.0)*h_i;
872 y_forward_2i = (t.*f)(dummy, x_forward_2i);
873 x_forward_2i(i) -= static_cast<type>(2.0)*h_i;
874
875 H(i,i) = (-y_forward_2i + type(16.0)*y_forward_i - type(30.0)*y + type(16.0)*y_backward_i - y_backward_2i)/(type(12.0)*pow(h_i, type(2)));
876
877 for(Index j = i; j < n; j++)
878 {
879 h_j = calculate_h(x(j));
880
881 x_backward_ij(i) -= h_i;
882 x_backward_ij(j) -= h_j;
883 y_backward_ij = (t.*f)(dummy, x_backward_ij);
884 x_backward_ij(i) += h_i;
885 x_backward_ij(j) += h_j;
886
887 x_forward_ij(i) += h_i;
888 x_forward_ij(j) += h_j;
889 y_forward_ij = (t.*f)(dummy, x_forward_ij);
890 x_forward_ij(i) -= h_i;
891 x_forward_ij(j) -= h_j;
892
893 x_backward_i_forward_j(i) -= h_i;
894 x_backward_i_forward_j(j) += h_j;
895 y_backward_i_forward_j = (t.*f)(dummy, x_backward_i_forward_j);
896 x_backward_i_forward_j(i) += h_i;
897 x_backward_i_forward_j(j) -= h_j;
898
899 x_forward_i_backward_j(i) += h_i;
900 x_forward_i_backward_j(j) -= h_j;
901 y_forward_i_backward_j = (t.*f)(dummy, x_forward_i_backward_j);
902 x_forward_i_backward_j(i) -= h_i;
903 x_forward_i_backward_j(j) += h_j;
904
905 H(i,j) = (y_forward_ij - y_forward_i_backward_j - y_backward_i_forward_j + y_backward_ij)/(type(4.0)*h_i*h_j);
906 }
907 }
908
909 for(Index i = 0; i < n; i++)
910 {
911 for(Index j = 0; j < i; j++)
912 {
913 H(i,j) = H(j,i);
914 }
915 }
916
917 return H;
918 }
919
925
926 template<class T>
927 Tensor<type, 2> calculate_Jacobian(const T& t, Tensor<type, 1>(T::*f)(const Tensor<type, 1>&) const, const Tensor<type, 1>& x) const
928 {
929 Tensor<type, 1> y = (t.*f)(x);
930
931 type h;
932
933 const Index n = x.size();
934 Index m = y.size();
935
936 Tensor<type, 1> x_forward(x);
937 Tensor<type, 1> x_backward(x);
938
939 Tensor<type, 1> y_forward(n);
940 Tensor<type, 1> y_backward(n);
941
942 Tensor<type, 2> J(m,n);
943
944 for(Index j = 0; j < n; j++)
945 {
946 h = calculate_h(x(j));
947
948 x_backward(j) -= h;
949 y_backward = (t.*f)(x_backward);
950 x_backward(j) += h;
951
952 x_forward(j) += h;
953 y_forward = (t.*f)(x_forward);
954 x_forward(j) -= h;
955
956 for(Index i = 0; i < m; i++)
957 {
958 J(i,j) = (y_forward(i) - y_backward(i))/(static_cast<type>(2.0)*h);
959 }
960 }
961
962 return J;
963 }
964
965
972
973 template<class T>
974 Tensor<type, 2> calculate_Jacobian(const T& t, Tensor<type, 1>(T::*f)(const Tensor<type, 1>&, const Tensor<type, 1>&) const, const Tensor<type, 1>& dummy, const Tensor<type, 1>& x) const
975 {
976 Tensor<type, 1> y = (t.*f)(dummy, x);
977
978 type h;
979
980 const Index n = x.size();
981 Index m = y.size();
982
983 Tensor<type, 1> x_forward(x);
984 Tensor<type, 1> x_backward(x);
985
986 Tensor<type, 1> y_forward(n);
987 Tensor<type, 1> y_backward(n);
988
989 Tensor<type, 2> J(m,n);
990
991 for(Index j = 0; j < n; j++)
992 {
993 h = calculate_h(x(j));
994
995 x_backward(j) -= h;
996 y_backward = (t.*f)(dummy, x_backward);
997 x_backward(j) += h;
998
999 x_forward(j) += h;
1000 y_forward = (t.*f)(dummy, x_forward);
1001 x_forward(j) -= h;
1002
1003 for(Index i = 0; i < m; i++)
1004 {
1005 J(i,j) = (y_forward(i) - y_backward(i))/(static_cast<type>(2.0)*h);
1006 }
1007 }
1008
1009 return J;
1010 }
1011
1019
1020 template<class T>
1021 Tensor<type, 2> calculate_Jacobian(const T& t, Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 1>&) const, const Index& dummy, const Tensor<type, 1>& x) const
1022 {
1023 Tensor<type, 1> y = (t.*f)(dummy, x);
1024
1025 type h;
1026
1027 const Index n = x.size();
1028 Index m = y.size();
1029
1030 Tensor<type, 1> x_forward(x);
1031 Tensor<type, 1> x_backward(x);
1032
1033 Tensor<type, 1> y_forward(n);
1034 Tensor<type, 1> y_backward(n);
1035
1036 Tensor<type, 2> J(m,n);
1037
1038 for(Index j = 0; j < n; j++)
1039 {
1040 h = calculate_h(x(j));
1041
1042 x_backward(j) -= h;
1043 y_backward = (t.*f)(dummy, x_backward);
1044 x_backward(j) += h;
1045
1046 x_forward(j) += h;
1047 y_forward = (t.*f)(dummy, x_forward);
1048 x_forward(j) -= h;
1049
1050 for(Index i = 0; i < m; i++)
1051 {
1052 J(i,j) = (y_forward(i) - y_backward(i))/(static_cast<type>(2.0)*h);
1053 }
1054 }
1055
1056 return J;
1057 }
1058
1067
1068 template<class T>
1069 Tensor<type, 2> calculate_Jacobian
1070 (const T& t, Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 1>&, const Tensor<type, 1>&) const, const Index& dummy_int, const Tensor<type, 1>& dummy_vector, const Tensor<type, 1>& x) const
1071 {
1072 const Index n = x.size();
1073
1074 const Tensor<type, 1> y = (t.*f)(dummy_int, dummy_vector, x);
1075 const Index m = y.size();
1076
1077 type h;
1078
1079 Tensor<type, 1> x_forward(x);
1080 Tensor<type, 1> x_backward(x);
1081
1082 Tensor<type, 1> y_forward(n);
1083 Tensor<type, 1> y_backward(n);
1084
1085 Tensor<type, 2> J(m,n);
1086
1087 for(Index j = 0; j < n; j++)
1088 {
1089 h = calculate_h(x(j));
1090
1091 x_backward(j) -= h;
1092 y_backward = (t.*f)(dummy_int, dummy_vector, x_backward);
1093 x_backward(j) += h;
1094
1095 x_forward(j) += h;
1096 y_forward = (t.*f)(dummy_int, dummy_vector, x_forward);
1097 x_forward(j) -= h;
1098
1099 for(Index i = 0; i < m; i++)
1100 {
1101 J(i,j) = (y_forward(i) - y_backward(i))/(static_cast<type>(2.0)*h);
1102 }
1103 }
1104
1105 return J;
1106 }
1107
1116
1117 template<class T>
1118 Tensor<type, 2> calculate_Jacobian
1119 (const T& t, Tensor<type, 1>(T::*f)(const Index&, const Index&, const Tensor<type, 1>&) const, const Index& dummy_int_1, const Index& dummy_int_2, const Tensor<type, 1>& x) const
1120 {
1121 const Tensor<type, 1> y = (t.*f)(dummy_int_1, dummy_int_2, x);
1122
1123 const Index n = x.size();
1124 const Index m = y.size();
1125
1126 type h;
1127
1128 Tensor<type, 1> x_forward(x);
1129 Tensor<type, 1> x_backward(x);
1130
1131 Tensor<type, 1> y_forward(n);
1132 Tensor<type, 1> y_backward(n);
1133
1134 Tensor<type, 2> J(m,n);
1135
1136 for(Index j = 0; j < n; j++)
1137 {
1138 h = calculate_h(x(j));
1139
1140 x_backward(j) -= h;
1141 y_backward = (t.*f)(dummy_int_1, dummy_int_2, x_backward);
1142 x_backward(j) += h;
1143
1144 x_forward(j) += h;
1145 y_forward = (t.*f)(dummy_int_1, dummy_int_2, x_forward);
1146 x_forward(j) -= h;
1147
1148 for(Index i = 0; i < m; i++)
1149 {
1150 J(i,j) = (y_forward(i) - y_backward(i))/(static_cast<type>(2.0)*h);
1151 }
1152 }
1153
1154 return J;
1155 }
1156
1162
1163 template<class T>
1164 Tensor<Tensor<type, 2>, 1> calculate_hessian(const T& t, Tensor<type, 1>(T::*f)(const Tensor<type, 1>&) const, const Tensor<type, 1>& x) const
1165 {
1166 Tensor<type, 1> y = (t.*f)(x);
1167
1168 Index s = y.size();
1169 const Index n = x.size();
1170
1171 type h_j;
1172 type h_k;
1173
1174 Tensor<type, 1> x_backward_2j(x);
1175 Tensor<type, 1> x_backward_j(x);
1176
1177 Tensor<type, 1> x_forward_j(x);
1178 Tensor<type, 1> x_forward_2j(x);
1179
1180 Tensor<type, 1> x_backward_jk(x);
1181 Tensor<type, 1> x_forward_jk(x);
1182
1183 Tensor<type, 1> x_backward_j_forward_k(x);
1184 Tensor<type, 1> x_forward_j_backward_k(x);
1185
1186 Tensor<type, 1> y_backward_2j;
1187 Tensor<type, 1> y_backward_j;
1188
1189 Tensor<type, 1> y_forward_j;
1190 Tensor<type, 1> y_forward_2j;
1191
1192 Tensor<type, 1> y_backward_jk;
1193 Tensor<type, 1> y_forward_jk;
1194
1195 Tensor<type, 1> y_backward_j_forward_k;
1196 Tensor<type, 1> y_forward_j_backward_k;
1197
1198 Tensor<Tensor<type, 2>, 1> H(s);
1199
1200 for(Index i = 0; i < s; i++)
1201 {
1202 H(i).resize(n,n);
1203
1204 for(Index j = 0; j < n; j++)
1205 {
1206 h_j = calculate_h(x(j));
1207
1208 x_backward_2j(j) -= static_cast<type>(2.0)*h_j;
1209 y_backward_2j = (t.*f)(x_backward_2j);
1210 x_backward_2j(j) += static_cast<type>(2.0)*h_j;
1211
1212 x_backward_j(j) -= h_j;
1213 y_backward_j = (t.*f)(x_backward_j);
1214 x_backward_j(j) += h_j;
1215
1216 x_forward_j(j) += h_j;
1217 y_forward_j = (t.*f)(x_forward_j);
1218 x_forward_j(j) -= h_j;
1219
1220 x_forward_2j(j) += static_cast<type>(2.0)*h_j;
1221 y_forward_2j = (t.*f)(x_forward_2j);
1222 x_forward_2j(j) -= static_cast<type>(2.0)*h_j;
1223
1224 H(i)(j,j) = (-y_forward_2j(i) + type(16.0)*y_forward_j(i) - type(30.0)*y(i) + type(16.0)*y_backward_j(i) - y_backward_2j(i))/(type(12.0)*pow(h_j, type(2)));
1225
1226 for(Index k = j; k < n; k++)
1227 {
1228 h_k = calculate_h(x[k]);
1229
1230 x_backward_jk(j) -= h_j;
1231 x_backward_jk[k] -= h_k;
1232 y_backward_jk = (t.*f)(x_backward_jk);
1233 x_backward_jk(j) += h_j;
1234 x_backward_jk[k] += h_k;
1235
1236 x_forward_jk(j) += h_j;
1237 x_forward_jk[k] += h_k;
1238 y_forward_jk = (t.*f)(x_forward_jk);
1239 x_forward_jk(j) -= h_j;
1240 x_forward_jk[k] -= h_k;
1241
1242 x_backward_j_forward_k(j) -= h_j;
1243 x_backward_j_forward_k[k] += h_k;
1244 y_backward_j_forward_k = (t.*f)(x_backward_j_forward_k);
1245 x_backward_j_forward_k(j) += h_j;
1246 x_backward_j_forward_k[k] -= h_k;
1247
1248 x_forward_j_backward_k(j) += h_j;
1249 x_forward_j_backward_k[k] -= h_k;
1250 y_forward_j_backward_k = (t.*f)(x_forward_j_backward_k);
1251 x_forward_j_backward_k(j) -= h_j;
1252 x_forward_j_backward_k[k] += h_k;
1253
1254 H(i)(j,k) = (y_forward_jk(i) - y_forward_j_backward_k(i) - y_backward_j_forward_k(i) + y_backward_jk(i))/(type(4.0)*h_j*h_k);
1255 }
1256 }
1257
1258 for(Index j = 0; j < n; j++)
1259 {
1260 for(Index k = 0; k < j; k++)
1261 {
1262 H(i)(j,k) = H(i)(k,j);
1263 }
1264 }
1265 }
1266
1267 return H;
1268 }
1269
1270
1278
1279 template<class T>
1280 Tensor<Tensor<type, 2>, 1> calculate_hessian
1281 (const T& t, Tensor<type, 1>(T::*f)(const Tensor<type, 1>&, const Tensor<type, 1>&) const, const Tensor<type, 1>& dummy_vector, const Tensor<type, 1>& x) const
1282 {
1283 Tensor<type, 1> y = (t.*f)(dummy_vector, x);
1284
1285 Index s = y.size();
1286 const Index n = x.size();
1287
1288 type h_j;
1289 type h_k;
1290
1291 Tensor<type, 1> x_backward_2j(x);
1292 Tensor<type, 1> x_backward_j(x);
1293
1294 Tensor<type, 1> x_forward_j(x);
1295 Tensor<type, 1> x_forward_2j(x);
1296
1297 Tensor<type, 1> x_backward_jk(x);
1298 Tensor<type, 1> x_forward_jk(x);
1299
1300 Tensor<type, 1> x_backward_j_forward_k(x);
1301 Tensor<type, 1> x_forward_j_backward_k(x);
1302
1303 Tensor<type, 1> y_backward_2j;
1304 Tensor<type, 1> y_backward_j;
1305
1306 Tensor<type, 1> y_forward_j;
1307 Tensor<type, 1> y_forward_2j;
1308
1309 Tensor<type, 1> y_backward_jk;
1310 Tensor<type, 1> y_forward_jk;
1311
1312 Tensor<type, 1> y_backward_j_forward_k;
1313 Tensor<type, 1> y_forward_j_backward_k;
1314
1315 Tensor<Tensor<type, 2>, 1> H(s);
1316
1317 for(Index i = 0; i < s; i++)
1318 {
1319 // H(i).set(n,n);
1320 H(i).resize(n,n);
1321
1322 for(Index j = 0; j < n; j++)
1323 {
1324 h_j = calculate_h(x(j));
1325
1326 x_backward_2j(j) -= static_cast<type>(2.0)*h_j;
1327 y_backward_2j = (t.*f)(dummy_vector, x_backward_2j);
1328 x_backward_2j(j) += static_cast<type>(2.0)*h_j;
1329
1330 x_backward_j(j) -= h_j;
1331 y_backward_j = (t.*f)(dummy_vector, x_backward_j);
1332 x_backward_j(j) += h_j;
1333
1334 x_forward_j(j) += h_j;
1335 y_forward_j = (t.*f)(dummy_vector, x_forward_j);
1336 x_forward_j(j) -= h_j;
1337
1338 x_forward_2j(j) += static_cast<type>(2.0)*h_j;
1339 y_forward_2j = (t.*f)(dummy_vector, x_forward_2j);
1340 x_forward_2j(j) -= static_cast<type>(2.0)*h_j;
1341
1342 H(i)(j,j) = (-y_forward_2j(i) + type(16.0)*y_forward_j(i) - type(30.0)*y(i) + type(16.0)*y_backward_j(i) - y_backward_2j(i))/(type(12.0)*pow(h_j, type(2)));
1343
1344 for(Index k = j; k < n; k++)
1345 {
1346 h_k = calculate_h(x[k]);
1347
1348 x_backward_jk(j) -= h_j;
1349 x_backward_jk[k] -= h_k;
1350 y_backward_jk = (t.*f)(dummy_vector, x_backward_jk);
1351 x_backward_jk(j) += h_j;
1352 x_backward_jk[k] += h_k;
1353
1354 x_forward_jk(j) += h_j;
1355 x_forward_jk[k] += h_k;
1356 y_forward_jk = (t.*f)(dummy_vector, x_forward_jk);
1357 x_forward_jk(j) -= h_j;
1358 x_forward_jk[k] -= h_k;
1359
1360 x_backward_j_forward_k(j) -= h_j;
1361 x_backward_j_forward_k[k] += h_k;
1362 y_backward_j_forward_k = (t.*f)(dummy_vector, x_backward_j_forward_k);
1363 x_backward_j_forward_k(j) += h_j;
1364 x_backward_j_forward_k[k] -= h_k;
1365
1366 x_forward_j_backward_k(j) += h_j;
1367 x_forward_j_backward_k[k] -= h_k;
1368 y_forward_j_backward_k = (t.*f)(dummy_vector, x_forward_j_backward_k);
1369 x_forward_j_backward_k(j) -= h_j;
1370 x_forward_j_backward_k[k] += h_k;
1371
1372 H(i)(j,k) = (y_forward_jk(i) - y_forward_j_backward_k(i) - y_backward_j_forward_k(i) + y_backward_jk(i))/(type(4.0)*h_j*h_k);
1373 }
1374 }
1375
1376 for(Index j = 0; j < n; j++)
1377 {
1378 for(Index k = 0; k < j; k++)
1379 {
1380 H(i)(j,k) = H(i)(k,j);
1381 }
1382 }
1383 }
1384
1385 return H;
1386 }
1387
1395
1396 template<class T>
1397 Tensor<Tensor<type, 2>, 1> calculate_hessian(const T& t, Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 1>&) const, const Index& dummy, const Tensor<type, 1>& x) const
1398 {
1399 Tensor<type, 1> y = (t.*f)(dummy, x);
1400
1401 Index s = y.size();
1402 const Index n = x.size();
1403
1404 type h_j;
1405 type h_k;
1406
1407 Tensor<type, 1> x_backward_2j(x);
1408 Tensor<type, 1> x_backward_j(x);
1409
1410 Tensor<type, 1> x_forward_j(x);
1411 Tensor<type, 1> x_forward_2j(x);
1412
1413 Tensor<type, 1> x_backward_jk(x);
1414 Tensor<type, 1> x_forward_jk(x);
1415
1416 Tensor<type, 1> x_backward_j_forward_k(x);
1417 Tensor<type, 1> x_forward_j_backward_k(x);
1418
1419 Tensor<type, 1> y_backward_2j;
1420 Tensor<type, 1> y_backward_j;
1421
1422 Tensor<type, 1> y_forward_j;
1423 Tensor<type, 1> y_forward_2j;
1424
1425 Tensor<type, 1> y_backward_jk;
1426 Tensor<type, 1> y_forward_jk;
1427
1428 Tensor<type, 1> y_backward_j_forward_k;
1429 Tensor<type, 1> y_forward_j_backward_k;
1430
1431 Tensor<Tensor<type, 2>, 1> H(s);
1432
1433 for(Index i = 0; i < s; i++)
1434 {
1435 H(i).resize(n,n);
1436
1437 for(Index j = 0; j < n; j++)
1438 {
1439 h_j = calculate_h(x(j));
1440
1441 x_backward_2j(j) -= static_cast<type>(2.0)*h_j;
1442 y_backward_2j = (t.*f)(dummy, x_backward_2j);
1443 x_backward_2j(j) += static_cast<type>(2.0)*h_j;
1444
1445 x_backward_j(j) -= h_j;
1446 y_backward_j = (t.*f)(dummy, x_backward_j);
1447 x_backward_j(j) += h_j;
1448
1449 x_forward_j(j) += h_j;
1450 y_forward_j = (t.*f)(dummy, x_forward_j);
1451 x_forward_j(j) -= h_j;
1452
1453 x_forward_2j(j) += static_cast<type>(2.0)*h_j;
1454 y_forward_2j = (t.*f)(dummy, x_forward_2j);
1455 x_forward_2j(j) -= static_cast<type>(2.0)*h_j;
1456
1457 H(i)(j,j) = (-y_forward_2j(i) + type(16.0)*y_forward_j(i) - type(30.0)*y(i) + type(16.0)*y_backward_j(i) - y_backward_2j(i))/(type(12.0)*pow(h_j, type(2)));
1458
1459 for(Index k = j; k < n; k++)
1460 {
1461 h_k = calculate_h(x[k]);
1462
1463 x_backward_jk(j) -= h_j;
1464 x_backward_jk[k] -= h_k;
1465 y_backward_jk = (t.*f)(dummy, x_backward_jk);
1466 x_backward_jk(j) += h_j;
1467 x_backward_jk[k] += h_k;
1468
1469 x_forward_jk(j) += h_j;
1470 x_forward_jk[k] += h_k;
1471 y_forward_jk = (t.*f)(dummy, x_forward_jk);
1472 x_forward_jk(j) -= h_j;
1473 x_forward_jk[k] -= h_k;
1474
1475 x_backward_j_forward_k(j) -= h_j;
1476 x_backward_j_forward_k[k] += h_k;
1477 y_backward_j_forward_k = (t.*f)(dummy, x_backward_j_forward_k);
1478 x_backward_j_forward_k(j) += h_j;
1479 x_backward_j_forward_k[k] -= h_k;
1480
1481 x_forward_j_backward_k(j) += h_j;
1482 x_forward_j_backward_k[k] -= h_k;
1483 y_forward_j_backward_k = (t.*f)(dummy, x_forward_j_backward_k);
1484 x_forward_j_backward_k(j) -= h_j;
1485 x_forward_j_backward_k[k] += h_k;
1486
1487 H(i)(j,k) = (y_forward_jk(i) - y_forward_j_backward_k(i) - y_backward_j_forward_k(i) + y_backward_jk(i))/(type(4.0)*h_j*h_k);
1488 }
1489 }
1490
1491 for(Index j = 0; j < n; j++)
1492 {
1493 for(Index k = 0; k < j; k++)
1494 {
1495 H(i)(j,k) = H(i)(k,j);
1496 }
1497 }
1498 }
1499
1500 return H;
1501 }
1502
1511
1512 template<class T>
1513 Tensor<Tensor<type, 2>, 1> calculate_hessian
1514 (const T& t, Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 1>&, const Tensor<type, 1>&) const, const Index& dummy_int, const Tensor<type, 1>& dummy_vector, const Tensor<type, 1>& x) const
1515 {
1516 const Tensor<type, 1> y = (t.*f)(dummy_int, dummy_vector, x);
1517
1518 Index s = y.size();
1519 const Index n = x.size();
1520
1521 type h_j;
1522 type h_k;
1523
1524 Tensor<type, 1> x_backward_2j(x);
1525 Tensor<type, 1> x_backward_j(x);
1526
1527 Tensor<type, 1> x_forward_j(x);
1528 Tensor<type, 1> x_forward_2j(x);
1529
1530 Tensor<type, 1> x_backward_jk(x);
1531 Tensor<type, 1> x_forward_jk(x);
1532
1533 Tensor<type, 1> x_backward_j_forward_k(x);
1534 Tensor<type, 1> x_forward_j_backward_k(x);
1535
1536 Tensor<type, 1> y_backward_2j;
1537 Tensor<type, 1> y_backward_j;
1538
1539 Tensor<type, 1> y_forward_j;
1540 Tensor<type, 1> y_forward_2j;
1541
1542 Tensor<type, 1> y_backward_jk;
1543 Tensor<type, 1> y_forward_jk;
1544
1545 Tensor<type, 1> y_backward_j_forward_k;
1546 Tensor<type, 1> y_forward_j_backward_k;
1547
1548 Tensor<Tensor<type, 2>, 1> H(s);
1549
1550 for(Index i = 0; i < s; i++)
1551 {
1552 H(i).resize(n,n);
1553
1554 for(Index j = 0; j < n; j++)
1555 {
1556 h_j = calculate_h(x(j));
1557
1558 x_backward_2j(j) -= static_cast<type>(2.0)*h_j;
1559 y_backward_2j = (t.*f)(dummy_int, dummy_vector, x_backward_2j);
1560 x_backward_2j(j) += static_cast<type>(2.0)*h_j;
1561
1562 x_backward_j(j) -= h_j;
1563 y_backward_j = (t.*f)(dummy_int, dummy_vector, x_backward_j);
1564 x_backward_j(j) += h_j;
1565
1566 x_forward_j(j) += h_j;
1567 y_forward_j = (t.*f)(dummy_int, dummy_vector, x_forward_j);
1568 x_forward_j(j) -= h_j;
1569
1570 x_forward_2j(j) += static_cast<type>(2.0)*h_j;
1571 y_forward_2j = (t.*f)(dummy_int, dummy_vector, x_forward_2j);
1572 x_forward_2j(j) -= static_cast<type>(2.0)*h_j;
1573
1574 H(i)(j,j) = (-y_forward_2j(i) + type(16.0)*y_forward_j(i) - type(30.0)*y(i) + type(16.0)*y_backward_j(i) - y_backward_2j(i))/(type(12.0)*pow(h_j, type(2)));
1575
1576 for(Index k = j; k < n; k++)
1577 {
1578 h_k = calculate_h(x[k]);
1579
1580 x_backward_jk(j) -= h_j;
1581 x_backward_jk[k] -= h_k;
1582 y_backward_jk = (t.*f)(dummy_int, dummy_vector, x_backward_jk);
1583 x_backward_jk(j) += h_j;
1584 x_backward_jk[k] += h_k;
1585
1586 x_forward_jk(j) += h_j;
1587 x_forward_jk[k] += h_k;
1588 y_forward_jk = (t.*f)(dummy_int, dummy_vector, x_forward_jk);
1589 x_forward_jk(j) -= h_j;
1590 x_forward_jk[k] -= h_k;
1591
1592 x_backward_j_forward_k(j) -= h_j;
1593 x_backward_j_forward_k[k] += h_k;
1594 y_backward_j_forward_k = (t.*f)(dummy_int, dummy_vector, x_backward_j_forward_k);
1595 x_backward_j_forward_k(j) += h_j;
1596 x_backward_j_forward_k[k] -= h_k;
1597
1598 x_forward_j_backward_k(j) += h_j;
1599 x_forward_j_backward_k[k] -= h_k;
1600 y_forward_j_backward_k = (t.*f)(dummy_int, dummy_vector, x_forward_j_backward_k);
1601 x_forward_j_backward_k(j) -= h_j;
1602 x_forward_j_backward_k[k] += h_k;
1603
1604 H(i)(j,k) = (y_forward_jk(i) - y_forward_j_backward_k(i) - y_backward_j_forward_k(i) + y_backward_jk(i))/(type(4.0)*h_j*h_k);
1605 }
1606 }
1607
1608 for(Index j = 0; j < n; j++)
1609 {
1610 for(Index k = 0; k < j; k++)
1611 {
1612 H(i)(j,k) = H(i)(k,j);
1613 }
1614 }
1615 }
1616
1617 return H;
1618 }
1619
1628
1629 template<class T>
1630 Tensor< Tensor<type, 2>, 2> calculate_hessian_matrices
1631 (const T& t, Tensor<type, 1>(T::*f)(const Index&, const Tensor<type, 1>&, const Tensor<type, 1>&) const, const Index& dummy_int, const Tensor<type, 1>& dummy_vector, const Tensor<type, 1>& x) const
1632 {
1633 const Tensor<type, 1> y = (t.*f)(dummy_int, dummy_vector, x);
1634
1635 Index s = y.size();
1636 const Index n = x.size();
1637
1638 type h_j;
1639 type h_k;
1640
1641 Tensor<type, 1> x_backward_2j(x);
1642 Tensor<type, 1> x_backward_j(x);
1643
1644 Tensor<type, 1> x_forward_j(x);
1645 Tensor<type, 1> x_forward_2j(x);
1646
1647 Tensor<type, 1> x_backward_jk(x);
1648 Tensor<type, 1> x_forward_jk(x);
1649
1650 Tensor<type, 1> x_backward_j_forward_k(x);
1651 Tensor<type, 1> x_forward_j_backward_k(x);
1652
1653 Tensor<type, 1> y_backward_2j;
1654 Tensor<type, 1> y_backward_j;
1655
1656 Tensor<type, 1> y_forward_j;
1657 Tensor<type, 1> y_forward_2j;
1658
1659 Tensor<type, 1> y_backward_jk;
1660 Tensor<type, 1> y_forward_jk;
1661
1662 Tensor<type, 1> y_backward_j_forward_k;
1663 Tensor<type, 1> y_forward_j_backward_k;
1664
1665 Tensor< Tensor<type, 2>, 2> H;
1666
1667 H.resize(n,s);
1668
1669 for(Index i = 0; i < s; i++)
1670 {
1671 for(Index t = 0; t < n; t++)
1672 {
1673 H(i,t).resize(n,s);
1674
1675 for(Index j = 0; j < n; j++)
1676 {
1677 h_j = calculate_h(x(j));
1678
1679 x_backward_2j(j) -= static_cast<type>(2.0)*h_j;
1680 y_backward_2j = (t.*f)(dummy_int, dummy_vector, x_backward_2j);
1681 x_backward_2j(j) += static_cast<type>(2.0)*h_j;
1682
1683 x_backward_j(j) -= h_j;
1684 y_backward_j = (t.*f)(dummy_int, dummy_vector, x_backward_j);
1685 x_backward_j(j) += h_j;
1686
1687 x_forward_j(j) += h_j;
1688 y_forward_j = (t.*f)(dummy_int, dummy_vector, x_forward_j);
1689 x_forward_j(j) -= h_j;
1690
1691 x_forward_2j(j) += static_cast<type>(2.0)*h_j;
1692 y_forward_2j = (t.*f)(dummy_int, dummy_vector, x_forward_2j);
1693 x_forward_2j(j) -= static_cast<type>(2.0)*h_j;
1694
1695 H(i)(j,j) = (-y_forward_2j(i) + type(16.0)*y_forward_j(i) - type(30.0)*y(i) + type(16.0)*y_backward_j(i) - y_backward_2j(i))/(type(12.0)*pow(h_j, type(2)));
1696
1697 for(Index k = j; k < s; k++)
1698 {
1699 h_k = calculate_h(x[k]);
1700
1701 x_backward_jk(j) -= h_j;
1702 x_backward_jk[k] -= h_k;
1703 y_backward_jk = (t.*f)(dummy_int, dummy_vector, x_backward_jk);
1704 x_backward_jk(j) += h_j;
1705 x_backward_jk[k] += h_k;
1706
1707 x_forward_jk(j) += h_j;
1708 x_forward_jk[k] += h_k;
1709 y_forward_jk = (t.*f)(dummy_int, dummy_vector, x_forward_jk);
1710 x_forward_jk(j) -= h_j;
1711 x_forward_jk[k] -= h_k;
1712
1713 x_backward_j_forward_k(j) -= h_j;
1714 x_backward_j_forward_k[k] += h_k;
1715 y_backward_j_forward_k = (t.*f)(dummy_int, dummy_vector, x_backward_j_forward_k);
1716 x_backward_j_forward_k(j) += h_j;
1717 x_backward_j_forward_k[k] -= h_k;
1718
1719 x_forward_j_backward_k(j) += h_j;
1720 x_forward_j_backward_k[k] -= h_k;
1721 y_forward_j_backward_k = (t.*f)(dummy_int, dummy_vector, x_forward_j_backward_k);
1722 x_forward_j_backward_k(j) -= h_j;
1723 x_forward_j_backward_k[k] += h_k;
1724
1725 H(i,t)(j,k) = (y_forward_jk(i) - y_forward_j_backward_k(i) - y_backward_j_forward_k(i) + y_backward_jk(i))/(type(4.0)*h_j*h_k);
1726 }
1727 }
1728
1729 for(Index j = 0; j < n; j++)
1730 {
1731 for(Index k = 0; k < j; k++)
1732 {
1733 H(i,t)(j,k) = H(i)(k,j);
1734 }
1735 }
1736 }
1737 }
1738
1739 return H;
1740 }
1741
1742private:
1743
1745
1746 Index precision_digits;
1747
1749
1750 bool display = true;
1751
1752};
1753
1754}
1755
1756#endif
1757
1758
1759// OpenNN: Open Neural Networks Library.
1760// Copyright(C) 2005-2021 Artificial Intelligence Techniques, SL.
1761//
1762// This library is free software; you can redistribute it and/or
1763// modify it under the terms of the GNU Lesser General Public
1764// License as published by the Free Software Foundation; either
1765// version 2.1 of the License, or any later version.
1766//
1767// This library is distributed in the hope that it will be useful,
1768// but WITHOUT ANY WARRANTY; without even the implied warranty of
1769// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1770// Lesser General Public License for more details.
1771
1772// You should have received a copy of the GNU Lesser General Public
1773// License along with this library; if not, write to the Free Software
1774// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
OpenNN::NumericalDifferentiation::calculate_gradient
Tensor< type, 1 > calculate_gradient(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Tensor< type, 2 > &) const, const Index &dummy, const Tensor< type, 2 > &x) const
Definition: numerical_differentiation.h:483
OpenNN::NumericalDifferentiation::calculate_hessian
Tensor< type, 2 > calculate_hessian(const T &t, type(T::*f)(const Index &, const Tensor< type, 1 > &) const, const Index &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:820
OpenNN::NumericalDifferentiation::get_precision_digits
const Index & get_precision_digits() const
Enumeration of available methods for numerical differentiation.
Definition: numerical_differentiation.cpp:32
OpenNN::NumericalDifferentiation::calculate_hessian
Tensor< type, 2 > calculate_hessian(const T &t, type(T::*f)(const Tensor< type, 1 > &) const, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:601
OpenNN::NumericalDifferentiation::calculate_Jacobian
Tensor< type, 2 > calculate_Jacobian(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Index &, const Tensor< type, 1 > &) const, const Index &dummy_int_1, const Index &dummy_int_2, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:1119
OpenNN::NumericalDifferentiation::get_display
const bool & get_display() const
Returns the flag used by this class for displaying or not displaying warnings.
Definition: numerical_differentiation.cpp:40
OpenNN::NumericalDifferentiation::calculate_second_derivatives
type calculate_second_derivatives(const T &t, type(T::*f)(const type &) const, const type &x) const
Definition: numerical_differentiation.h:221
OpenNN::NumericalDifferentiation::calculate_gradient
Tensor< type, 1 > calculate_gradient(const T &t, type(T::*f)(const Tensor< type, 1 > &, const Tensor< type, 1 > &) const, const Tensor< type, 1 > &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:400
OpenNN::NumericalDifferentiation::calculate_gradient
Tensor< type, 1 > calculate_gradient(const T &t, type(T::*f)(const Index &, const Tensor< type, 1 > &) const, const Index &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:441
OpenNN::NumericalDifferentiation::calculate_derivatives
Tensor< type, 1 > calculate_derivatives(const T &t, Tensor< type, 1 >(T::*f)(const Tensor< type, 1 > &) const, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:88
OpenNN::NumericalDifferentiation::calculate_Jacobian
Tensor< type, 2 > calculate_Jacobian(const T &t, Tensor< type, 1 >(T::*f)(const Tensor< type, 1 > &) const, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:927
OpenNN::NumericalDifferentiation::display
bool display
Flag for displaying warning messages from this class.
Definition: numerical_differentiation.h:1750
OpenNN::NumericalDifferentiation::calculate_gradient
Tensor< type, 1 > calculate_gradient(const T &t, type(T::*f)(const Tensor< type, 1 > &), const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:359
OpenNN::NumericalDifferentiation::calculate_hessian
Tensor< Tensor< type, 2 >, 1 > calculate_hessian(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Tensor< type, 1 > &, const Tensor< type, 1 > &) const, const Index &dummy_int, const Tensor< type, 1 > &dummy_vector, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:1514
OpenNN::NumericalDifferentiation::calculate_Jacobian
Tensor< type, 2 > calculate_Jacobian(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Tensor< type, 1 > &) const, const Index &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:1021
OpenNN::NumericalDifferentiation::calculate_derivatives
Tensor< type, 1 > calculate_derivatives(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Tensor< type, 1 > &) const, const Index &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:153
OpenNN::NumericalDifferentiation::calculate_gradient
Tensor< type, 1 > calculate_gradient(const T &t, type(T::*f)(const Tensor< Index, 1 > &, const Tensor< type, 1 > &) const, const Tensor< Index, 1 > &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:525
OpenNN::NumericalDifferentiation::calculate_second_derivatives
Tensor< type, 1 > calculate_second_derivatives(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Tensor< type, 1 > &) const, const Index &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:285
OpenNN::NumericalDifferentiation::calculate_second_derivatives
Tensor< type, 1 > calculate_second_derivatives(const T &t, Tensor< type, 1 >(T::*f)(const Tensor< type, 1 > &) const, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:255
OpenNN::NumericalDifferentiation::precision_digits
Index precision_digits
Number of precision digits.
Definition: numerical_differentiation.h:1746
OpenNN::NumericalDifferentiation::calculate_hessian
Tensor< Tensor< type, 2 >, 1 > calculate_hessian(const T &t, Tensor< type, 1 >(T::*f)(const Tensor< type, 1 > &, const Tensor< type, 1 > &) const, const Tensor< type, 1 > &dummy_vector, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:1281
OpenNN::NumericalDifferentiation::calculate_hessian
Tensor< Tensor< type, 2 >, 1 > calculate_hessian(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Tensor< type, 1 > &) const, const Index &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:1397
OpenNN::NumericalDifferentiation::calculate_hessian
Tensor< Tensor< type, 2 >, 1 > calculate_hessian(const T &t, Tensor< type, 1 >(T::*f)(const Tensor< type, 1 > &) const, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:1164
OpenNN::NumericalDifferentiation::calculate_hessian_matrices
Tensor< Tensor< type, 2 >, 2 > calculate_hessian_matrices(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Tensor< type, 1 > &, const Tensor< type, 1 > &) const, const Index &dummy_int, const Tensor< type, 1 > &dummy_vector, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:1631
OpenNN::NumericalDifferentiation::calculate_Jacobian
Tensor< type, 2 > calculate_Jacobian(const T &t, Tensor< type, 1 >(T::*f)(const Index &, const Tensor< type, 1 > &, const Tensor< type, 1 > &) const, const Index &dummy_int, const Tensor< type, 1 > &dummy_vector, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:1070
OpenNN::NumericalDifferentiation::calculate_gradient
Tensor< type, 1 > calculate_gradient(const T &t, type(T::*f)(const Tensor< type, 1 > &) const, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:316
OpenNN::NumericalDifferentiation::calculate_hessian
Tensor< type, 2 > calculate_hessian(const T &t, type(T::*f)(const Tensor< type, 1 > &, const Tensor< type, 1 > &) const, const Tensor< type, 1 > &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:711
OpenNN::NumericalDifferentiation::calculate_Jacobian
Tensor< type, 2 > calculate_Jacobian(const T &t, Tensor< type, 1 >(T::*f)(const Tensor< type, 1 > &, const Tensor< type, 1 > &) const, const Tensor< type, 1 > &dummy, const Tensor< type, 1 > &x) const
Definition: numerical_differentiation.h:974
half_float::pow
half pow(half x, half y)
Definition: half.hpp:3427
std
Extensions to the C++ standard library.
Definition: half.hpp:2325