Merge pull request #6352 from fraserharris:patch-1
[platform/upstream/opencv.git] / modules / ml / src / ann_mlp.cpp
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 //  By downloading, copying, installing or using the software you agree to this license.
6 //  If you do not agree to this license, do not download, install,
7 //  copy or use the software.
8 //
9 //
10 //                        Intel License Agreement
11 //
12 // Copyright (C) 2000, Intel Corporation, all rights reserved.
13 // Third party copyrights are property of their respective owners.
14 //
15 // Redistribution and use in source and binary forms, with or without modification,
16 // are permitted provided that the following conditions are met:
17 //
18 //   * Redistribution's of source code must retain the above copyright notice,
19 //     this list of conditions and the following disclaimer.
20 //
21 //   * Redistribution's in binary form must reproduce the above copyright notice,
22 //     this list of conditions and the following disclaimer in the documentation
23 //     and/or other materials provided with the distribution.
24 //
25 //   * The name of Intel Corporation may not be used to endorse or promote products
26 //     derived from this software without specific prior written permission.
27 //
28 // This software is provided by the copyright holders and contributors "as is" and
29 // any express or implied warranties, including, but not limited to, the implied
30 // warranties of merchantability and fitness for a particular purpose are disclaimed.
31 // In no event shall the Intel Corporation or contributors be liable for any direct,
32 // indirect, incidental, special, exemplary, or consequential damages
33 // (including, but not limited to, procurement of substitute goods or services;
34 // loss of use, data, or profits; or business interruption) however caused
35 // and on any theory of liability, whether in contract, strict liability,
36 // or tort (including negligence or otherwise) arising in any way out of
37 // the use of this software, even if advised of the possibility of such damage.
38 //
39 //M*/
40
41 #include "precomp.hpp"
42
43 namespace cv { namespace ml {
44
45 struct AnnParams
46 {
47     AnnParams()
48     {
49         termCrit = TermCriteria( TermCriteria::COUNT + TermCriteria::EPS, 1000, 0.01 );
50         trainMethod = ANN_MLP::RPROP;
51         bpDWScale = bpMomentScale = 0.1;
52         rpDW0 = 0.1; rpDWPlus = 1.2; rpDWMinus = 0.5;
53         rpDWMin = FLT_EPSILON; rpDWMax = 50.;
54     }
55
56     TermCriteria termCrit;
57     int trainMethod;
58
59     double bpDWScale;
60     double bpMomentScale;
61
62     double rpDW0;
63     double rpDWPlus;
64     double rpDWMinus;
65     double rpDWMin;
66     double rpDWMax;
67 };
68
69 template <typename T>
70 inline T inBounds(T val, T min_val, T max_val)
71 {
72     return std::min(std::max(val, min_val), max_val);
73 }
74
75 class ANN_MLPImpl : public ANN_MLP
76 {
77 public:
78     ANN_MLPImpl()
79     {
80         clear();
81         setActivationFunction( SIGMOID_SYM, 0, 0 );
82         setLayerSizes(Mat());
83         setTrainMethod(ANN_MLP::RPROP, 0.1, FLT_EPSILON);
84     }
85
86     virtual ~ANN_MLPImpl() {}
87
88     CV_IMPL_PROPERTY(TermCriteria, TermCriteria, params.termCrit)
89     CV_IMPL_PROPERTY(double, BackpropWeightScale, params.bpDWScale)
90     CV_IMPL_PROPERTY(double, BackpropMomentumScale, params.bpMomentScale)
91     CV_IMPL_PROPERTY(double, RpropDW0, params.rpDW0)
92     CV_IMPL_PROPERTY(double, RpropDWPlus, params.rpDWPlus)
93     CV_IMPL_PROPERTY(double, RpropDWMinus, params.rpDWMinus)
94     CV_IMPL_PROPERTY(double, RpropDWMin, params.rpDWMin)
95     CV_IMPL_PROPERTY(double, RpropDWMax, params.rpDWMax)
96
97     void clear()
98     {
99         min_val = max_val = min_val1 = max_val1 = 0.;
100         rng = RNG((uint64)-1);
101         weights.clear();
102         trained = false;
103         max_buf_sz = 1 << 12;
104     }
105
106     int layer_count() const { return (int)layer_sizes.size(); }
107
108     void setTrainMethod(int method, double param1, double param2)
109     {
110         if (method != ANN_MLP::RPROP && method != ANN_MLP::BACKPROP)
111             method = ANN_MLP::RPROP;
112         params.trainMethod = method;
113         if(method == ANN_MLP::RPROP )
114         {
115             if( param1 < FLT_EPSILON )
116                 param1 = 1.;
117             params.rpDW0 = param1;
118             params.rpDWMin = std::max( param2, 0. );
119         }
120         else if(method == ANN_MLP::BACKPROP )
121         {
122             if( param1 <= 0 )
123                 param1 = 0.1;
124             params.bpDWScale = inBounds<double>(param1, 1e-3, 1.);
125             if( param2 < 0 )
126                 param2 = 0.1;
127             params.bpMomentScale = std::min( param2, 1. );
128         }
129     }
130
131     int getTrainMethod() const
132     {
133         return params.trainMethod;
134     }
135
136     void setActivationFunction(int _activ_func, double _f_param1, double _f_param2 )
137     {
138         if( _activ_func < 0 || _activ_func > GAUSSIAN )
139             CV_Error( CV_StsOutOfRange, "Unknown activation function" );
140
141         activ_func = _activ_func;
142
143         switch( activ_func )
144         {
145         case SIGMOID_SYM:
146             max_val = 0.95; min_val = -max_val;
147             max_val1 = 0.98; min_val1 = -max_val1;
148             if( fabs(_f_param1) < FLT_EPSILON )
149                 _f_param1 = 2./3;
150             if( fabs(_f_param2) < FLT_EPSILON )
151                 _f_param2 = 1.7159;
152             break;
153         case GAUSSIAN:
154             max_val = 1.; min_val = 0.05;
155             max_val1 = 1.; min_val1 = 0.02;
156             if( fabs(_f_param1) < FLT_EPSILON )
157                 _f_param1 = 1.;
158             if( fabs(_f_param2) < FLT_EPSILON )
159                 _f_param2 = 1.;
160             break;
161         default:
162             min_val = max_val = min_val1 = max_val1 = 0.;
163             _f_param1 = 1.;
164             _f_param2 = 0.;
165         }
166
167         f_param1 = _f_param1;
168         f_param2 = _f_param2;
169     }
170
171
172     void init_weights()
173     {
174         int i, j, k, l_count = layer_count();
175
176         for( i = 1; i < l_count; i++ )
177         {
178             int n1 = layer_sizes[i-1];
179             int n2 = layer_sizes[i];
180             double val = 0, G = n2 > 2 ? 0.7*pow((double)n1,1./(n2-1)) : 1.;
181             double* w = weights[i].ptr<double>();
182
183             // initialize weights using Nguyen-Widrow algorithm
184             for( j = 0; j < n2; j++ )
185             {
186                 double s = 0;
187                 for( k = 0; k <= n1; k++ )
188                 {
189                     val = rng.uniform(0., 1.)*2-1.;
190                     w[k*n2 + j] = val;
191                     s += fabs(val);
192                 }
193
194                 if( i < l_count - 1 )
195                 {
196                     s = 1./(s - fabs(val));
197                     for( k = 0; k <= n1; k++ )
198                         w[k*n2 + j] *= s;
199                     w[n1*n2 + j] *= G*(-1+j*2./n2);
200                 }
201             }
202         }
203     }
204
205     Mat getLayerSizes() const
206     {
207         return Mat_<int>(layer_sizes, true);
208     }
209
210     void setLayerSizes( InputArray _layer_sizes )
211     {
212         clear();
213
214         _layer_sizes.copyTo(layer_sizes);
215         int l_count = layer_count();
216
217         weights.resize(l_count + 2);
218         max_lsize = 0;
219
220         if( l_count > 0 )
221         {
222             for( int i = 0; i < l_count; i++ )
223             {
224                 int n = layer_sizes[i];
225                 if( n < 1 + (0 < i && i < l_count-1))
226                     CV_Error( CV_StsOutOfRange,
227                              "there should be at least one input and one output "
228                              "and every hidden layer must have more than 1 neuron" );
229                 max_lsize = std::max( max_lsize, n );
230                 if( i > 0 )
231                     weights[i].create(layer_sizes[i-1]+1, n, CV_64F);
232             }
233
234             int ninputs = layer_sizes.front();
235             int noutputs = layer_sizes.back();
236             weights[0].create(1, ninputs*2, CV_64F);
237             weights[l_count].create(1, noutputs*2, CV_64F);
238             weights[l_count+1].create(1, noutputs*2, CV_64F);
239         }
240     }
241
242     float predict( InputArray _inputs, OutputArray _outputs, int ) const
243     {
244         if( !trained )
245             CV_Error( CV_StsError, "The network has not been trained or loaded" );
246
247         Mat inputs = _inputs.getMat();
248         int type = inputs.type(), l_count = layer_count();
249         int n = inputs.rows, dn0 = n;
250
251         CV_Assert( (type == CV_32F || type == CV_64F) && inputs.cols == layer_sizes[0] );
252         int noutputs = layer_sizes[l_count-1];
253         Mat outputs;
254
255         int min_buf_sz = 2*max_lsize;
256         int buf_sz = n*min_buf_sz;
257
258         if( buf_sz > max_buf_sz )
259         {
260             dn0 = max_buf_sz/min_buf_sz;
261             dn0 = std::max( dn0, 1 );
262             buf_sz = dn0*min_buf_sz;
263         }
264
265         cv::AutoBuffer<double> _buf(buf_sz+noutputs);
266         double* buf = _buf;
267
268         if( !_outputs.needed() )
269         {
270             CV_Assert( n == 1 );
271             outputs = Mat(n, noutputs, type, buf + buf_sz);
272         }
273         else
274         {
275             _outputs.create(n, noutputs, type);
276             outputs = _outputs.getMat();
277         }
278
279         int dn = 0;
280         for( int i = 0; i < n; i += dn )
281         {
282             dn = std::min( dn0, n - i );
283
284             Mat layer_in = inputs.rowRange(i, i + dn);
285             Mat layer_out( dn, layer_in.cols, CV_64F, buf);
286
287             scale_input( layer_in, layer_out );
288             layer_in = layer_out;
289
290             for( int j = 1; j < l_count; j++ )
291             {
292                 double* data = buf + ((j&1) ? max_lsize*dn0 : 0);
293                 int cols = layer_sizes[j];
294
295                 layer_out = Mat(dn, cols, CV_64F, data);
296                 Mat w = weights[j].rowRange(0, layer_in.cols);
297                 gemm(layer_in, w, 1, noArray(), 0, layer_out);
298                 calc_activ_func( layer_out, weights[j] );
299
300                 layer_in = layer_out;
301             }
302
303             layer_out = outputs.rowRange(i, i + dn);
304             scale_output( layer_in, layer_out );
305         }
306
307         if( n == 1 )
308         {
309             int maxIdx[] = {0, 0};
310             minMaxIdx(outputs, 0, 0, 0, maxIdx);
311             return (float)(maxIdx[0] + maxIdx[1]);
312         }
313
314         return 0.f;
315     }
316
317     void scale_input( const Mat& _src, Mat& _dst ) const
318     {
319         int cols = _src.cols;
320         const double* w = weights[0].ptr<double>();
321
322         if( _src.type() == CV_32F )
323         {
324             for( int i = 0; i < _src.rows; i++ )
325             {
326                 const float* src = _src.ptr<float>(i);
327                 double* dst = _dst.ptr<double>(i);
328                 for( int j = 0; j < cols; j++ )
329                     dst[j] = src[j]*w[j*2] + w[j*2+1];
330             }
331         }
332         else
333         {
334             for( int i = 0; i < _src.rows; i++ )
335             {
336                 const float* src = _src.ptr<float>(i);
337                 double* dst = _dst.ptr<double>(i);
338                 for( int j = 0; j < cols; j++ )
339                     dst[j] = src[j]*w[j*2] + w[j*2+1];
340             }
341         }
342     }
343
344     void scale_output( const Mat& _src, Mat& _dst ) const
345     {
346         int cols = _src.cols;
347         const double* w = weights[layer_count()].ptr<double>();
348
349         if( _dst.type() == CV_32F )
350         {
351             for( int i = 0; i < _src.rows; i++ )
352             {
353                 const double* src = _src.ptr<double>(i);
354                 float* dst = _dst.ptr<float>(i);
355                 for( int j = 0; j < cols; j++ )
356                     dst[j] = (float)(src[j]*w[j*2] + w[j*2+1]);
357             }
358         }
359         else
360         {
361             for( int i = 0; i < _src.rows; i++ )
362             {
363                 const double* src = _src.ptr<double>(i);
364                 double* dst = _dst.ptr<double>(i);
365                 for( int j = 0; j < cols; j++ )
366                     dst[j] = src[j]*w[j*2] + w[j*2+1];
367             }
368         }
369     }
370
371     void calc_activ_func( Mat& sums, const Mat& w ) const
372     {
373         const double* bias = w.ptr<double>(w.rows-1);
374         int i, j, n = sums.rows, cols = sums.cols;
375         double scale = 0, scale2 = f_param2;
376
377         switch( activ_func )
378         {
379             case IDENTITY:
380                 scale = 1.;
381                 break;
382             case SIGMOID_SYM:
383                 scale = -f_param1;
384                 break;
385             case GAUSSIAN:
386                 scale = -f_param1*f_param1;
387                 break;
388             default:
389                 ;
390         }
391
392         CV_Assert( sums.isContinuous() );
393
394         if( activ_func != GAUSSIAN )
395         {
396             for( i = 0; i < n; i++ )
397             {
398                 double* data = sums.ptr<double>(i);
399                 for( j = 0; j < cols; j++ )
400                     data[j] = (data[j] + bias[j])*scale;
401             }
402
403             if( activ_func == IDENTITY )
404                 return;
405         }
406         else
407         {
408             for( i = 0; i < n; i++ )
409             {
410                 double* data = sums.ptr<double>(i);
411                 for( j = 0; j < cols; j++ )
412                 {
413                     double t = data[j] + bias[j];
414                     data[j] = t*t*scale;
415                 }
416             }
417         }
418
419         exp( sums, sums );
420
421         if( sums.isContinuous() )
422         {
423             cols *= n;
424             n = 1;
425         }
426
427         switch( activ_func )
428         {
429             case SIGMOID_SYM:
430                 for( i = 0; i < n; i++ )
431                 {
432                     double* data = sums.ptr<double>(i);
433                     for( j = 0; j < cols; j++ )
434                     {
435                         double t = scale2*(1. - data[j])/(1. + data[j]);
436                         data[j] = t;
437                     }
438                 }
439                 break;
440
441             case GAUSSIAN:
442                 for( i = 0; i < n; i++ )
443                 {
444                     double* data = sums.ptr<double>(i);
445                     for( j = 0; j < cols; j++ )
446                         data[j] = scale2*data[j];
447                 }
448                 break;
449
450             default:
451                 ;
452         }
453     }
454
455     void calc_activ_func_deriv( Mat& _xf, Mat& _df, const Mat& w ) const
456     {
457         const double* bias = w.ptr<double>(w.rows-1);
458         int i, j, n = _xf.rows, cols = _xf.cols;
459
460         if( activ_func == IDENTITY )
461         {
462             for( i = 0; i < n; i++ )
463             {
464                 double* xf = _xf.ptr<double>(i);
465                 double* df = _df.ptr<double>(i);
466
467                 for( j = 0; j < cols; j++ )
468                 {
469                     xf[j] += bias[j];
470                     df[j] = 1;
471                 }
472             }
473         }
474         else if( activ_func == GAUSSIAN )
475         {
476             double scale = -f_param1*f_param1;
477             double scale2 = scale*f_param2;
478             for( i = 0; i < n; i++ )
479             {
480                 double* xf = _xf.ptr<double>(i);
481                 double* df = _df.ptr<double>(i);
482
483                 for( j = 0; j < cols; j++ )
484                 {
485                     double t = xf[j] + bias[j];
486                     df[j] = t*2*scale2;
487                     xf[j] = t*t*scale;
488                 }
489             }
490             exp( _xf, _xf );
491
492             for( i = 0; i < n; i++ )
493             {
494                 double* xf = _xf.ptr<double>(i);
495                 double* df = _df.ptr<double>(i);
496
497                 for( j = 0; j < cols; j++ )
498                     df[j] *= xf[j];
499             }
500         }
501         else
502         {
503             double scale = f_param1;
504             double scale2 = f_param2;
505
506             for( i = 0; i < n; i++ )
507             {
508                 double* xf = _xf.ptr<double>(i);
509                 double* df = _df.ptr<double>(i);
510
511                 for( j = 0; j < cols; j++ )
512                 {
513                     xf[j] = (xf[j] + bias[j])*scale;
514                     df[j] = -fabs(xf[j]);
515                 }
516             }
517
518             exp( _df, _df );
519
520             // ((1+exp(-ax))^-1)'=a*((1+exp(-ax))^-2)*exp(-ax);
521             // ((1-exp(-ax))/(1+exp(-ax)))'=(a*exp(-ax)*(1+exp(-ax)) + a*exp(-ax)*(1-exp(-ax)))/(1+exp(-ax))^2=
522             // 2*a*exp(-ax)/(1+exp(-ax))^2
523             scale *= 2*f_param2;
524             for( i = 0; i < n; i++ )
525             {
526                 double* xf = _xf.ptr<double>(i);
527                 double* df = _df.ptr<double>(i);
528
529                 for( j = 0; j < cols; j++ )
530                 {
531                     int s0 = xf[j] > 0 ? 1 : -1;
532                     double t0 = 1./(1. + df[j]);
533                     double t1 = scale*df[j]*t0*t0;
534                     t0 *= scale2*(1. - df[j])*s0;
535                     df[j] = t1;
536                     xf[j] = t0;
537                 }
538             }
539         }
540     }
541
542     void calc_input_scale( const Mat& inputs, int flags )
543     {
544         bool reset_weights = (flags & UPDATE_WEIGHTS) == 0;
545         bool no_scale = (flags & NO_INPUT_SCALE) != 0;
546         double* scale = weights[0].ptr<double>();
547         int count = inputs.rows;
548
549         if( reset_weights )
550         {
551             int i, j, vcount = layer_sizes[0];
552             int type = inputs.type();
553             double a = no_scale ? 1. : 0.;
554
555             for( j = 0; j < vcount; j++ )
556                 scale[2*j] = a, scale[j*2+1] = 0.;
557
558             if( no_scale )
559                 return;
560
561             for( i = 0; i < count; i++ )
562             {
563                 const uchar* p = inputs.ptr(i);
564                 const float* f = (const float*)p;
565                 const double* d = (const double*)p;
566                 for( j = 0; j < vcount; j++ )
567                 {
568                     double t = type == CV_32F ? (double)f[j] : d[j];
569                     scale[j*2] += t;
570                     scale[j*2+1] += t*t;
571                 }
572             }
573
574             for( j = 0; j < vcount; j++ )
575             {
576                 double s = scale[j*2], s2 = scale[j*2+1];
577                 double m = s/count, sigma2 = s2/count - m*m;
578                 scale[j*2] = sigma2 < DBL_EPSILON ? 1 : 1./sqrt(sigma2);
579                 scale[j*2+1] = -m*scale[j*2];
580             }
581         }
582     }
583
584     void calc_output_scale( const Mat& outputs, int flags )
585     {
586         int i, j, vcount = layer_sizes.back();
587         int type = outputs.type();
588         double m = min_val, M = max_val, m1 = min_val1, M1 = max_val1;
589         bool reset_weights = (flags & UPDATE_WEIGHTS) == 0;
590         bool no_scale = (flags & NO_OUTPUT_SCALE) != 0;
591         int l_count = layer_count();
592         double* scale = weights[l_count].ptr<double>();
593         double* inv_scale = weights[l_count+1].ptr<double>();
594         int count = outputs.rows;
595
596         if( reset_weights )
597         {
598             double a0 = no_scale ? 1 : DBL_MAX, b0 = no_scale ? 0 : -DBL_MAX;
599
600             for( j = 0; j < vcount; j++ )
601             {
602                 scale[2*j] = inv_scale[2*j] = a0;
603                 scale[j*2+1] = inv_scale[2*j+1] = b0;
604             }
605
606             if( no_scale )
607                 return;
608         }
609
610         for( i = 0; i < count; i++ )
611         {
612             const uchar* p = outputs.ptr(i);
613             const float* f = (const float*)p;
614             const double* d = (const double*)p;
615
616             for( j = 0; j < vcount; j++ )
617             {
618                 double t = type == CV_32F ? (double)f[j] : d[j];
619
620                 if( reset_weights )
621                 {
622                     double mj = scale[j*2], Mj = scale[j*2+1];
623                     if( mj > t ) mj = t;
624                     if( Mj < t ) Mj = t;
625
626                     scale[j*2] = mj;
627                     scale[j*2+1] = Mj;
628                 }
629                 else if( !no_scale )
630                 {
631                     t = t*inv_scale[j*2] + inv_scale[2*j+1];
632                     if( t < m1 || t > M1 )
633                         CV_Error( CV_StsOutOfRange,
634                                  "Some of new output training vector components run exceed the original range too much" );
635                 }
636             }
637         }
638
639         if( reset_weights )
640             for( j = 0; j < vcount; j++ )
641             {
642                 // map mj..Mj to m..M
643                 double mj = scale[j*2], Mj = scale[j*2+1];
644                 double a, b;
645                 double delta = Mj - mj;
646                 if( delta < DBL_EPSILON )
647                     a = 1, b = (M + m - Mj - mj)*0.5;
648                 else
649                     a = (M - m)/delta, b = m - mj*a;
650                 inv_scale[j*2] = a; inv_scale[j*2+1] = b;
651                 a = 1./a; b = -b*a;
652                 scale[j*2] = a; scale[j*2+1] = b;
653             }
654     }
655
656     void prepare_to_train( const Mat& inputs, const Mat& outputs,
657                            Mat& sample_weights, int flags )
658     {
659         if( layer_sizes.empty() )
660             CV_Error( CV_StsError,
661                      "The network has not been created. Use method create or the appropriate constructor" );
662
663         if( (inputs.type() != CV_32F && inputs.type() != CV_64F) ||
664             inputs.cols != layer_sizes[0] )
665             CV_Error( CV_StsBadArg,
666                      "input training data should be a floating-point matrix with "
667                      "the number of rows equal to the number of training samples and "
668                      "the number of columns equal to the size of 0-th (input) layer" );
669
670         if( (outputs.type() != CV_32F && outputs.type() != CV_64F) ||
671             outputs.cols != layer_sizes.back() )
672             CV_Error( CV_StsBadArg,
673                      "output training data should be a floating-point matrix with "
674                      "the number of rows equal to the number of training samples and "
675                      "the number of columns equal to the size of last (output) layer" );
676
677         if( inputs.rows != outputs.rows )
678             CV_Error( CV_StsUnmatchedSizes, "The numbers of input and output samples do not match" );
679
680         Mat temp;
681         double s = sum(sample_weights)[0];
682         sample_weights.convertTo(temp, CV_64F, 1./s);
683         sample_weights = temp;
684
685         calc_input_scale( inputs, flags );
686         calc_output_scale( outputs, flags );
687     }
688
689     bool train( const Ptr<TrainData>& trainData, int flags )
690     {
691         const int MAX_ITER = 1000;
692         const double DEFAULT_EPSILON = FLT_EPSILON;
693
694         // initialize training data
695         Mat inputs = trainData->getTrainSamples();
696         Mat outputs = trainData->getTrainResponses();
697         Mat sw = trainData->getTrainSampleWeights();
698         prepare_to_train( inputs, outputs, sw, flags );
699
700         // ... and link weights
701         if( !(flags & UPDATE_WEIGHTS) )
702             init_weights();
703
704         TermCriteria termcrit;
705         termcrit.type = TermCriteria::COUNT + TermCriteria::EPS;
706         termcrit.maxCount = std::max((params.termCrit.type & CV_TERMCRIT_ITER ? params.termCrit.maxCount : MAX_ITER), 1);
707         termcrit.epsilon = std::max((params.termCrit.type & CV_TERMCRIT_EPS ? params.termCrit.epsilon : DEFAULT_EPSILON), DBL_EPSILON);
708
709         int iter = params.trainMethod == ANN_MLP::BACKPROP ?
710             train_backprop( inputs, outputs, sw, termcrit ) :
711             train_rprop( inputs, outputs, sw, termcrit );
712
713         trained = iter > 0;
714         return trained;
715     }
716
717     int train_backprop( const Mat& inputs, const Mat& outputs, const Mat& _sw, TermCriteria termCrit )
718     {
719         int i, j, k;
720         double prev_E = DBL_MAX*0.5, E = 0;
721         int itype = inputs.type(), otype = outputs.type();
722
723         int count = inputs.rows;
724
725         int iter = -1, max_iter = termCrit.maxCount*count;
726         double epsilon = termCrit.epsilon*count;
727
728         int l_count = layer_count();
729         int ivcount = layer_sizes[0];
730         int ovcount = layer_sizes.back();
731
732         // allocate buffers
733         vector<vector<double> > x(l_count);
734         vector<vector<double> > df(l_count);
735         vector<Mat> dw(l_count);
736
737         for( i = 0; i < l_count; i++ )
738         {
739             int n = layer_sizes[i];
740             x[i].resize(n+1);
741             df[i].resize(n);
742             dw[i] = Mat::zeros(weights[i].size(), CV_64F);
743         }
744
745         Mat _idx_m(1, count, CV_32S);
746         int* _idx = _idx_m.ptr<int>();
747         for( i = 0; i < count; i++ )
748             _idx[i] = i;
749
750         AutoBuffer<double> _buf(max_lsize*2);
751         double* buf[] = { _buf, (double*)_buf + max_lsize };
752
753         const double* sw = _sw.empty() ? 0 : _sw.ptr<double>();
754
755         // run back-propagation loop
756         /*
757          y_i = w_i*x_{i-1}
758          x_i = f(y_i)
759          E = 1/2*||u - x_N||^2
760          grad_N = (x_N - u)*f'(y_i)
761          dw_i(t) = momentum*dw_i(t-1) + dw_scale*x_{i-1}*grad_i
762          w_i(t+1) = w_i(t) + dw_i(t)
763          grad_{i-1} = w_i^t*grad_i
764         */
765         for( iter = 0; iter < max_iter; iter++ )
766         {
767             int idx = iter % count;
768             double sweight = sw ? count*sw[idx] : 1.;
769
770             if( idx == 0 )
771             {
772                 //printf("%d. E = %g\n", iter/count, E);
773                 if( fabs(prev_E - E) < epsilon )
774                     break;
775                 prev_E = E;
776                 E = 0;
777
778                 // shuffle indices
779                 for( i = 0; i < count; i++ )
780                 {
781                     j = rng.uniform(0, count);
782                     k = rng.uniform(0, count);
783                     std::swap(_idx[j], _idx[k]);
784                 }
785             }
786
787             idx = _idx[idx];
788
789             const uchar* x0data_p = inputs.ptr(idx);
790             const float* x0data_f = (const float*)x0data_p;
791             const double* x0data_d = (const double*)x0data_p;
792
793             double* w = weights[0].ptr<double>();
794             for( j = 0; j < ivcount; j++ )
795                 x[0][j] = (itype == CV_32F ? (double)x0data_f[j] : x0data_d[j])*w[j*2] + w[j*2 + 1];
796
797             Mat x1( 1, ivcount, CV_64F, &x[0][0] );
798
799             // forward pass, compute y[i]=w*x[i-1], x[i]=f(y[i]), df[i]=f'(y[i])
800             for( i = 1; i < l_count; i++ )
801             {
802                 int n = layer_sizes[i];
803                 Mat x2(1, n, CV_64F, &x[i][0] );
804                 Mat _w = weights[i].rowRange(0, x1.cols);
805                 gemm(x1, _w, 1, noArray(), 0, x2);
806                 Mat _df(1, n, CV_64F, &df[i][0] );
807                 calc_activ_func_deriv( x2, _df, weights[i] );
808                 x1 = x2;
809             }
810
811             Mat grad1( 1, ovcount, CV_64F, buf[l_count&1] );
812             w = weights[l_count+1].ptr<double>();
813
814             // calculate error
815             const uchar* udata_p = outputs.ptr(idx);
816             const float* udata_f = (const float*)udata_p;
817             const double* udata_d = (const double*)udata_p;
818
819             double* gdata = grad1.ptr<double>();
820             for( k = 0; k < ovcount; k++ )
821             {
822                 double t = (otype == CV_32F ? (double)udata_f[k] : udata_d[k])*w[k*2] + w[k*2+1] - x[l_count-1][k];
823                 gdata[k] = t*sweight;
824                 E += t*t;
825             }
826             E *= sweight;
827
828             // backward pass, update weights
829             for( i = l_count-1; i > 0; i-- )
830             {
831                 int n1 = layer_sizes[i-1], n2 = layer_sizes[i];
832                 Mat _df(1, n2, CV_64F, &df[i][0]);
833                 multiply( grad1, _df, grad1 );
834                 Mat _x(n1+1, 1, CV_64F, &x[i-1][0]);
835                 x[i-1][n1] = 1.;
836                 gemm( _x, grad1, params.bpDWScale, dw[i], params.bpMomentScale, dw[i] );
837                 add( weights[i], dw[i], weights[i] );
838                 if( i > 1 )
839                 {
840                     Mat grad2(1, n1, CV_64F, buf[i&1]);
841                     Mat _w = weights[i].rowRange(0, n1);
842                     gemm( grad1, _w, 1, noArray(), 0, grad2, GEMM_2_T );
843                     grad1 = grad2;
844                 }
845             }
846         }
847
848         iter /= count;
849         return iter;
850     }
851
852     struct RPropLoop : public ParallelLoopBody
853     {
854         RPropLoop(ANN_MLPImpl* _ann,
855                   const Mat& _inputs, const Mat& _outputs, const Mat& _sw,
856                   int _dcount0, vector<Mat>& _dEdw, double* _E)
857         {
858             ann = _ann;
859             inputs = _inputs;
860             outputs = _outputs;
861             sw = _sw.ptr<double>();
862             dcount0 = _dcount0;
863             dEdw = &_dEdw;
864             pE = _E;
865         }
866
867         ANN_MLPImpl* ann;
868         vector<Mat>* dEdw;
869         Mat inputs, outputs;
870         const double* sw;
871         int dcount0;
872         double* pE;
873
874         void operator()( const Range& range ) const
875         {
876             double inv_count = 1./inputs.rows;
877             int ivcount = ann->layer_sizes.front();
878             int ovcount = ann->layer_sizes.back();
879             int itype = inputs.type(), otype = outputs.type();
880             int count = inputs.rows;
881             int i, j, k, l_count = ann->layer_count();
882             vector<vector<double> > x(l_count);
883             vector<vector<double> > df(l_count);
884             vector<double> _buf(ann->max_lsize*dcount0*2);
885             double* buf[] = { &_buf[0], &_buf[ann->max_lsize*dcount0] };
886             double E = 0;
887
888             for( i = 0; i < l_count; i++ )
889             {
890                 x[i].resize(ann->layer_sizes[i]*dcount0);
891                 df[i].resize(ann->layer_sizes[i]*dcount0);
892             }
893
894             for( int si = range.start; si < range.end; si++ )
895             {
896                 int i0 = si*dcount0, i1 = std::min((si + 1)*dcount0, count);
897                 int dcount = i1 - i0;
898                 const double* w = ann->weights[0].ptr<double>();
899
900                 // grab and preprocess input data
901                 for( i = 0; i < dcount; i++ )
902                 {
903                     const uchar* x0data_p = inputs.ptr(i0 + i);
904                     const float* x0data_f = (const float*)x0data_p;
905                     const double* x0data_d = (const double*)x0data_p;
906
907                     double* xdata = &x[0][i*ivcount];
908                     for( j = 0; j < ivcount; j++ )
909                         xdata[j] = (itype == CV_32F ? (double)x0data_f[j] : x0data_d[j])*w[j*2] + w[j*2+1];
910                 }
911                 Mat x1(dcount, ivcount, CV_64F, &x[0][0]);
912
913                 // forward pass, compute y[i]=w*x[i-1], x[i]=f(y[i]), df[i]=f'(y[i])
914                 for( i = 1; i < l_count; i++ )
915                 {
916                     Mat x2( dcount, ann->layer_sizes[i], CV_64F, &x[i][0] );
917                     Mat _w = ann->weights[i].rowRange(0, x1.cols);
918                     gemm( x1, _w, 1, noArray(), 0, x2 );
919                     Mat _df( x2.size(), CV_64F, &df[i][0] );
920                     ann->calc_activ_func_deriv( x2, _df, ann->weights[i] );
921                     x1 = x2;
922                 }
923
924                 Mat grad1(dcount, ovcount, CV_64F, buf[l_count & 1]);
925
926                 w = ann->weights[l_count+1].ptr<double>();
927
928                 // calculate error
929                 for( i = 0; i < dcount; i++ )
930                 {
931                     const uchar* udata_p = outputs.ptr(i0+i);
932                     const float* udata_f = (const float*)udata_p;
933                     const double* udata_d = (const double*)udata_p;
934
935                     const double* xdata = &x[l_count-1][i*ovcount];
936                     double* gdata = grad1.ptr<double>(i);
937                     double sweight = sw ? sw[si+i] : inv_count, E1 = 0;
938
939                     for( j = 0; j < ovcount; j++ )
940                     {
941                         double t = (otype == CV_32F ? (double)udata_f[j] : udata_d[j])*w[j*2] + w[j*2+1] - xdata[j];
942                         gdata[j] = t*sweight;
943                         E1 += t*t;
944                     }
945                     E += sweight*E1;
946                 }
947
948                 for( i = l_count-1; i > 0; i-- )
949                 {
950                     int n1 = ann->layer_sizes[i-1], n2 = ann->layer_sizes[i];
951                     Mat _df(dcount, n2, CV_64F, &df[i][0]);
952                     multiply(grad1, _df, grad1);
953
954                     {
955                         AutoLock lock(ann->mtx);
956                         Mat _dEdw = dEdw->at(i).rowRange(0, n1);
957                         x1 = Mat(dcount, n1, CV_64F, &x[i-1][0]);
958                         gemm(x1, grad1, 1, _dEdw, 1, _dEdw, GEMM_1_T);
959
960                         // update bias part of dEdw
961                         double* dst = dEdw->at(i).ptr<double>(n1);
962                         for( k = 0; k < dcount; k++ )
963                         {
964                             const double* src = grad1.ptr<double>(k);
965                             for( j = 0; j < n2; j++ )
966                                 dst[j] += src[j];
967                         }
968                     }
969
970                     Mat grad2( dcount, n1, CV_64F, buf[i&1] );
971                     if( i > 1 )
972                     {
973                         Mat _w = ann->weights[i].rowRange(0, n1);
974                         gemm(grad1, _w, 1, noArray(), 0, grad2, GEMM_2_T);
975                     }
976                     grad1 = grad2;
977                 }
978             }
979             {
980                 AutoLock lock(ann->mtx);
981                 *pE += E;
982             }
983         }
984     };
985
986     int train_rprop( const Mat& inputs, const Mat& outputs, const Mat& _sw, TermCriteria termCrit )
987     {
988         const int max_buf_size = 1 << 16;
989         int i, iter = -1, count = inputs.rows;
990
991         double prev_E = DBL_MAX*0.5;
992
993         int max_iter = termCrit.maxCount;
994         double epsilon = termCrit.epsilon;
995         double dw_plus = params.rpDWPlus;
996         double dw_minus = params.rpDWMinus;
997         double dw_min = params.rpDWMin;
998         double dw_max = params.rpDWMax;
999
1000         int l_count = layer_count();
1001
1002         // allocate buffers
1003         vector<Mat> dw(l_count), dEdw(l_count), prev_dEdw_sign(l_count);
1004
1005         int total = 0;
1006         for( i = 0; i < l_count; i++ )
1007         {
1008             total += layer_sizes[i];
1009             dw[i].create(weights[i].size(), CV_64F);
1010             dw[i].setTo(Scalar::all(params.rpDW0));
1011             prev_dEdw_sign[i] = Mat::zeros(weights[i].size(), CV_8S);
1012             dEdw[i] = Mat::zeros(weights[i].size(), CV_64F);
1013         }
1014
1015         int dcount0 = max_buf_size/(2*total);
1016         dcount0 = std::max( dcount0, 1 );
1017         dcount0 = std::min( dcount0, count );
1018         int chunk_count = (count + dcount0 - 1)/dcount0;
1019
1020         // run rprop loop
1021         /*
1022          y_i(t) = w_i(t)*x_{i-1}(t)
1023          x_i(t) = f(y_i(t))
1024          E = sum_over_all_samples(1/2*||u - x_N||^2)
1025          grad_N = (x_N - u)*f'(y_i)
1026
1027          std::min(dw_i{jk}(t)*dw_plus, dw_max), if dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) > 0
1028          dw_i{jk}(t) = std::max(dw_i{jk}(t)*dw_minus, dw_min), if dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) < 0
1029          dw_i{jk}(t-1) else
1030
1031          if (dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) < 0)
1032          dE/dw_i{jk}(t)<-0
1033          else
1034          w_i{jk}(t+1) = w_i{jk}(t) + dw_i{jk}(t)
1035          grad_{i-1}(t) = w_i^t(t)*grad_i(t)
1036          */
1037         for( iter = 0; iter < max_iter; iter++ )
1038         {
1039             double E = 0;
1040
1041             for( i = 0; i < l_count; i++ )
1042                 dEdw[i].setTo(Scalar::all(0));
1043
1044             // first, iterate through all the samples and compute dEdw
1045             RPropLoop invoker(this, inputs, outputs, _sw, dcount0, dEdw, &E);
1046             parallel_for_(Range(0, chunk_count), invoker);
1047             //invoker(Range(0, chunk_count));
1048
1049             // now update weights
1050             for( i = 1; i < l_count; i++ )
1051             {
1052                 int n1 = layer_sizes[i-1], n2 = layer_sizes[i];
1053                 for( int k = 0; k <= n1; k++ )
1054                 {
1055                     CV_Assert(weights[i].size() == Size(n2, n1+1));
1056                     double* wk = weights[i].ptr<double>(k);
1057                     double* dwk = dw[i].ptr<double>(k);
1058                     double* dEdwk = dEdw[i].ptr<double>(k);
1059                     schar* prevEk = prev_dEdw_sign[i].ptr<schar>(k);
1060
1061                     for( int j = 0; j < n2; j++ )
1062                     {
1063                         double Eval = dEdwk[j];
1064                         double dval = dwk[j];
1065                         double wval = wk[j];
1066                         int s = CV_SIGN(Eval);
1067                         int ss = prevEk[j]*s;
1068                         if( ss > 0 )
1069                         {
1070                             dval *= dw_plus;
1071                             dval = std::min( dval, dw_max );
1072                             dwk[j] = dval;
1073                             wk[j] = wval + dval*s;
1074                         }
1075                         else if( ss < 0 )
1076                         {
1077                             dval *= dw_minus;
1078                             dval = std::max( dval, dw_min );
1079                             prevEk[j] = 0;
1080                             dwk[j] = dval;
1081                             wk[j] = wval + dval*s;
1082                         }
1083                         else
1084                         {
1085                             prevEk[j] = (schar)s;
1086                             wk[j] = wval + dval*s;
1087                         }
1088                         dEdwk[j] = 0.;
1089                     }
1090                 }
1091             }
1092
1093             //printf("%d. E = %g\n", iter, E);
1094             if( fabs(prev_E - E) < epsilon )
1095                 break;
1096             prev_E = E;
1097         }
1098
1099         return iter;
1100     }
1101
1102     void write_params( FileStorage& fs ) const
1103     {
1104         const char* activ_func_name = activ_func == IDENTITY ? "IDENTITY" :
1105                                       activ_func == SIGMOID_SYM ? "SIGMOID_SYM" :
1106                                       activ_func == GAUSSIAN ? "GAUSSIAN" : 0;
1107
1108         if( activ_func_name )
1109             fs << "activation_function" << activ_func_name;
1110         else
1111             fs << "activation_function_id" << activ_func;
1112
1113         if( activ_func != IDENTITY )
1114         {
1115             fs << "f_param1" << f_param1;
1116             fs << "f_param2" << f_param2;
1117         }
1118
1119         fs << "min_val" << min_val << "max_val" << max_val << "min_val1" << min_val1 << "max_val1" << max_val1;
1120
1121         fs << "training_params" << "{";
1122         if( params.trainMethod == ANN_MLP::BACKPROP )
1123         {
1124             fs << "train_method" << "BACKPROP";
1125             fs << "dw_scale" << params.bpDWScale;
1126             fs << "moment_scale" << params.bpMomentScale;
1127         }
1128         else if( params.trainMethod == ANN_MLP::RPROP )
1129         {
1130             fs << "train_method" << "RPROP";
1131             fs << "dw0" << params.rpDW0;
1132             fs << "dw_plus" << params.rpDWPlus;
1133             fs << "dw_minus" << params.rpDWMinus;
1134             fs << "dw_min" << params.rpDWMin;
1135             fs << "dw_max" << params.rpDWMax;
1136         }
1137         else
1138             CV_Error(CV_StsError, "Unknown training method");
1139
1140         fs << "term_criteria" << "{";
1141         if( params.termCrit.type & TermCriteria::EPS )
1142             fs << "epsilon" << params.termCrit.epsilon;
1143         if( params.termCrit.type & TermCriteria::COUNT )
1144             fs << "iterations" << params.termCrit.maxCount;
1145         fs << "}" << "}";
1146     }
1147
1148     void write( FileStorage& fs ) const
1149     {
1150         if( layer_sizes.empty() )
1151             return;
1152         int i, l_count = layer_count();
1153
1154         writeFormat(fs);
1155         fs << "layer_sizes" << layer_sizes;
1156
1157         write_params( fs );
1158
1159         size_t esz = weights[0].elemSize();
1160
1161         fs << "input_scale" << "[";
1162         fs.writeRaw("d", weights[0].ptr(), weights[0].total()*esz);
1163
1164         fs << "]" << "output_scale" << "[";
1165         fs.writeRaw("d", weights[l_count].ptr(), weights[l_count].total()*esz);
1166
1167         fs << "]" << "inv_output_scale" << "[";
1168         fs.writeRaw("d", weights[l_count+1].ptr(), weights[l_count+1].total()*esz);
1169
1170         fs << "]" << "weights" << "[";
1171         for( i = 1; i < l_count; i++ )
1172         {
1173             fs << "[";
1174             fs.writeRaw("d", weights[i].ptr(), weights[i].total()*esz);
1175             fs << "]";
1176         }
1177         fs << "]";
1178     }
1179
1180     void read_params( const FileNode& fn )
1181     {
1182         String activ_func_name = (String)fn["activation_function"];
1183         if( !activ_func_name.empty() )
1184         {
1185             activ_func = activ_func_name == "SIGMOID_SYM" ? SIGMOID_SYM :
1186                          activ_func_name == "IDENTITY" ? IDENTITY :
1187                          activ_func_name == "GAUSSIAN" ? GAUSSIAN : -1;
1188             CV_Assert( activ_func >= 0 );
1189         }
1190         else
1191             activ_func = (int)fn["activation_function_id"];
1192
1193         f_param1 = (double)fn["f_param1"];
1194         f_param2 = (double)fn["f_param2"];
1195
1196         setActivationFunction( activ_func, f_param1, f_param2 );
1197
1198         min_val = (double)fn["min_val"];
1199         max_val = (double)fn["max_val"];
1200         min_val1 = (double)fn["min_val1"];
1201         max_val1 = (double)fn["max_val1"];
1202
1203         FileNode tpn = fn["training_params"];
1204         params = AnnParams();
1205
1206         if( !tpn.empty() )
1207         {
1208             String tmethod_name = (String)tpn["train_method"];
1209
1210             if( tmethod_name == "BACKPROP" )
1211             {
1212                 params.trainMethod = ANN_MLP::BACKPROP;
1213                 params.bpDWScale = (double)tpn["dw_scale"];
1214                 params.bpMomentScale = (double)tpn["moment_scale"];
1215             }
1216             else if( tmethod_name == "RPROP" )
1217             {
1218                 params.trainMethod = ANN_MLP::RPROP;
1219                 params.rpDW0 = (double)tpn["dw0"];
1220                 params.rpDWPlus = (double)tpn["dw_plus"];
1221                 params.rpDWMinus = (double)tpn["dw_minus"];
1222                 params.rpDWMin = (double)tpn["dw_min"];
1223                 params.rpDWMax = (double)tpn["dw_max"];
1224             }
1225             else
1226                 CV_Error(CV_StsParseError, "Unknown training method (should be BACKPROP or RPROP)");
1227
1228             FileNode tcn = tpn["term_criteria"];
1229             if( !tcn.empty() )
1230             {
1231                 FileNode tcn_e = tcn["epsilon"];
1232                 FileNode tcn_i = tcn["iterations"];
1233                 params.termCrit.type = 0;
1234                 if( !tcn_e.empty() )
1235                 {
1236                     params.termCrit.type |= TermCriteria::EPS;
1237                     params.termCrit.epsilon = (double)tcn_e;
1238                 }
1239                 if( !tcn_i.empty() )
1240                 {
1241                     params.termCrit.type |= TermCriteria::COUNT;
1242                     params.termCrit.maxCount = (int)tcn_i;
1243                 }
1244             }
1245         }
1246     }
1247
1248     void read( const FileNode& fn )
1249     {
1250         clear();
1251
1252         vector<int> _layer_sizes;
1253         readVectorOrMat(fn["layer_sizes"], _layer_sizes);
1254         setLayerSizes( _layer_sizes );
1255
1256         int i, l_count = layer_count();
1257         read_params(fn);
1258
1259         size_t esz = weights[0].elemSize();
1260
1261         FileNode w = fn["input_scale"];
1262         w.readRaw("d", weights[0].ptr(), weights[0].total()*esz);
1263
1264         w = fn["output_scale"];
1265         w.readRaw("d", weights[l_count].ptr(), weights[l_count].total()*esz);
1266
1267         w = fn["inv_output_scale"];
1268         w.readRaw("d", weights[l_count+1].ptr(), weights[l_count+1].total()*esz);
1269
1270         FileNodeIterator w_it = fn["weights"].begin();
1271
1272         for( i = 1; i < l_count; i++, ++w_it )
1273             (*w_it).readRaw("d", weights[i].ptr(), weights[i].total()*esz);
1274         trained = true;
1275     }
1276
1277     Mat getWeights(int layerIdx) const
1278     {
1279         CV_Assert( 0 <= layerIdx && layerIdx < (int)weights.size() );
1280         return weights[layerIdx];
1281     }
1282
1283     bool isTrained() const
1284     {
1285         return trained;
1286     }
1287
1288     bool isClassifier() const
1289     {
1290         return false;
1291     }
1292
1293     int getVarCount() const
1294     {
1295         return layer_sizes.empty() ? 0 : layer_sizes[0];
1296     }
1297
1298     String getDefaultName() const
1299     {
1300         return "opencv_ml_ann_mlp";
1301     }
1302
1303     vector<int> layer_sizes;
1304     vector<Mat> weights;
1305     double f_param1, f_param2;
1306     double min_val, max_val, min_val1, max_val1;
1307     int activ_func;
1308     int max_lsize, max_buf_sz;
1309     AnnParams params;
1310     RNG rng;
1311     Mutex mtx;
1312     bool trained;
1313 };
1314
1315
1316 Ptr<ANN_MLP> ANN_MLP::create()
1317 {
1318     return makePtr<ANN_MLPImpl>();
1319 }
1320
1321 Ptr<ANN_MLP> ANN_MLP::load(const String& filepath)
1322 {
1323     FileStorage fs;
1324     fs.open(filepath, FileStorage::READ);
1325
1326     Ptr<ANN_MLP> ann = makePtr<ANN_MLPImpl>();
1327
1328     ((ANN_MLPImpl*)ann.get())->read(fs.getFirstTopLevelNode());
1329     return ann;
1330 }
1331
1332
1333     }}
1334
1335 /* End of file. */