1 /*M///////////////////////////////////////////////////////////////////////////////////////
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
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.
10 // Intel License Agreement
12 // Copyright (C) 2000, Intel Corporation, all rights reserved.
13 // Third party copyrights are property of their respective owners.
15 // Redistribution and use in source and binary forms, with or without modification,
16 // are permitted provided that the following conditions are met:
18 // * Redistribution's of source code must retain the above copyright notice,
19 // this list of conditions and the following disclaimer.
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.
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.
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.
41 #include "precomp.hpp"
43 namespace cv { namespace ml {
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.;
56 TermCriteria termCrit;
70 inline T inBounds(T val, T min_val, T max_val)
72 return std::min(std::max(val, min_val), max_val);
75 class ANN_MLPImpl : public ANN_MLP
81 setActivationFunction( SIGMOID_SYM, 0, 0 );
83 setTrainMethod(ANN_MLP::RPROP, 0.1, FLT_EPSILON);
86 virtual ~ANN_MLPImpl() {}
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)
99 min_val = max_val = min_val1 = max_val1 = 0.;
100 rng = RNG((uint64)-1);
103 max_buf_sz = 1 << 12;
106 int layer_count() const { return (int)layer_sizes.size(); }
108 void setTrainMethod(int method, double param1, double param2)
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 )
115 if( param1 < FLT_EPSILON )
117 params.rpDW0 = param1;
118 params.rpDWMin = std::max( param2, 0. );
120 else if(method == ANN_MLP::BACKPROP )
124 params.bpDWScale = inBounds<double>(param1, 1e-3, 1.);
127 params.bpMomentScale = std::min( param2, 1. );
131 int getTrainMethod() const
133 return params.trainMethod;
136 void setActivationFunction(int _activ_func, double _f_param1, double _f_param2 )
138 if( _activ_func < 0 || _activ_func > GAUSSIAN )
139 CV_Error( CV_StsOutOfRange, "Unknown activation function" );
141 activ_func = _activ_func;
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 )
150 if( fabs(_f_param2) < FLT_EPSILON )
154 max_val = 1.; min_val = 0.05;
155 max_val1 = 1.; min_val1 = 0.02;
156 if( fabs(_f_param1) < FLT_EPSILON )
158 if( fabs(_f_param2) < FLT_EPSILON )
162 min_val = max_val = min_val1 = max_val1 = 0.;
167 f_param1 = _f_param1;
168 f_param2 = _f_param2;
174 int i, j, k, l_count = layer_count();
176 for( i = 1; i < l_count; i++ )
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>();
183 // initialize weights using Nguyen-Widrow algorithm
184 for( j = 0; j < n2; j++ )
187 for( k = 0; k <= n1; k++ )
189 val = rng.uniform(0., 1.)*2-1.;
194 if( i < l_count - 1 )
196 s = 1./(s - fabs(val));
197 for( k = 0; k <= n1; k++ )
199 w[n1*n2 + j] *= G*(-1+j*2./n2);
205 Mat getLayerSizes() const
207 return Mat_<int>(layer_sizes, true);
210 void setLayerSizes( InputArray _layer_sizes )
214 _layer_sizes.copyTo(layer_sizes);
215 int l_count = layer_count();
217 weights.resize(l_count + 2);
222 for( int i = 0; i < l_count; i++ )
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 );
231 weights[i].create(layer_sizes[i-1]+1, n, CV_64F);
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);
242 float predict( InputArray _inputs, OutputArray _outputs, int ) const
245 CV_Error( CV_StsError, "The network has not been trained or loaded" );
247 Mat inputs = _inputs.getMat();
248 int type = inputs.type(), l_count = layer_count();
249 int n = inputs.rows, dn0 = n;
251 CV_Assert( (type == CV_32F || type == CV_64F) && inputs.cols == layer_sizes[0] );
252 int noutputs = layer_sizes[l_count-1];
255 int min_buf_sz = 2*max_lsize;
256 int buf_sz = n*min_buf_sz;
258 if( buf_sz > max_buf_sz )
260 dn0 = max_buf_sz/min_buf_sz;
261 dn0 = std::max( dn0, 1 );
262 buf_sz = dn0*min_buf_sz;
265 cv::AutoBuffer<double> _buf(buf_sz+noutputs);
268 if( !_outputs.needed() )
271 outputs = Mat(n, noutputs, type, buf + buf_sz);
275 _outputs.create(n, noutputs, type);
276 outputs = _outputs.getMat();
280 for( int i = 0; i < n; i += dn )
282 dn = std::min( dn0, n - i );
284 Mat layer_in = inputs.rowRange(i, i + dn);
285 Mat layer_out( dn, layer_in.cols, CV_64F, buf);
287 scale_input( layer_in, layer_out );
288 layer_in = layer_out;
290 for( int j = 1; j < l_count; j++ )
292 double* data = buf + ((j&1) ? max_lsize*dn0 : 0);
293 int cols = layer_sizes[j];
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] );
300 layer_in = layer_out;
303 layer_out = outputs.rowRange(i, i + dn);
304 scale_output( layer_in, layer_out );
309 int maxIdx[] = {0, 0};
310 minMaxIdx(outputs, 0, 0, 0, maxIdx);
311 return (float)(maxIdx[0] + maxIdx[1]);
317 void scale_input( const Mat& _src, Mat& _dst ) const
319 int cols = _src.cols;
320 const double* w = weights[0].ptr<double>();
322 if( _src.type() == CV_32F )
324 for( int i = 0; i < _src.rows; i++ )
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];
334 for( int i = 0; i < _src.rows; i++ )
336 const double* src = _src.ptr<double>(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];
344 void scale_output( const Mat& _src, Mat& _dst ) const
346 int cols = _src.cols;
347 const double* w = weights[layer_count()].ptr<double>();
349 if( _dst.type() == CV_32F )
351 for( int i = 0; i < _src.rows; i++ )
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]);
361 for( int i = 0; i < _src.rows; i++ )
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];
371 void calc_activ_func( Mat& sums, const Mat& w ) const
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;
386 scale = -f_param1*f_param1;
392 CV_Assert( sums.isContinuous() );
394 if( activ_func != GAUSSIAN )
396 for( i = 0; i < n; i++ )
398 double* data = sums.ptr<double>(i);
399 for( j = 0; j < cols; j++ )
400 data[j] = (data[j] + bias[j])*scale;
403 if( activ_func == IDENTITY )
408 for( i = 0; i < n; i++ )
410 double* data = sums.ptr<double>(i);
411 for( j = 0; j < cols; j++ )
413 double t = data[j] + bias[j];
421 if( sums.isContinuous() )
430 for( i = 0; i < n; i++ )
432 double* data = sums.ptr<double>(i);
433 for( j = 0; j < cols; j++ )
435 if(!cvIsInf(data[j]))
437 double t = scale2*(1. - data[j])/(1. + data[j]);
449 for( i = 0; i < n; i++ )
451 double* data = sums.ptr<double>(i);
452 for( j = 0; j < cols; j++ )
453 data[j] = scale2*data[j];
462 void calc_activ_func_deriv( Mat& _xf, Mat& _df, const Mat& w ) const
464 const double* bias = w.ptr<double>(w.rows-1);
465 int i, j, n = _xf.rows, cols = _xf.cols;
467 if( activ_func == IDENTITY )
469 for( i = 0; i < n; i++ )
471 double* xf = _xf.ptr<double>(i);
472 double* df = _df.ptr<double>(i);
474 for( j = 0; j < cols; j++ )
481 else if( activ_func == GAUSSIAN )
483 double scale = -f_param1*f_param1;
484 double scale2 = scale*f_param2;
485 for( i = 0; i < n; i++ )
487 double* xf = _xf.ptr<double>(i);
488 double* df = _df.ptr<double>(i);
490 for( j = 0; j < cols; j++ )
492 double t = xf[j] + bias[j];
499 for( i = 0; i < n; i++ )
501 double* xf = _xf.ptr<double>(i);
502 double* df = _df.ptr<double>(i);
504 for( j = 0; j < cols; j++ )
510 double scale = f_param1;
511 double scale2 = f_param2;
513 for( i = 0; i < n; i++ )
515 double* xf = _xf.ptr<double>(i);
516 double* df = _df.ptr<double>(i);
518 for( j = 0; j < cols; j++ )
520 xf[j] = (xf[j] + bias[j])*scale;
521 df[j] = -fabs(xf[j]);
527 // ((1+exp(-ax))^-1)'=a*((1+exp(-ax))^-2)*exp(-ax);
528 // ((1-exp(-ax))/(1+exp(-ax)))'=(a*exp(-ax)*(1+exp(-ax)) + a*exp(-ax)*(1-exp(-ax)))/(1+exp(-ax))^2=
529 // 2*a*exp(-ax)/(1+exp(-ax))^2
531 for( i = 0; i < n; i++ )
533 double* xf = _xf.ptr<double>(i);
534 double* df = _df.ptr<double>(i);
536 for( j = 0; j < cols; j++ )
538 int s0 = xf[j] > 0 ? 1 : -1;
539 double t0 = 1./(1. + df[j]);
540 double t1 = scale*df[j]*t0*t0;
541 t0 *= scale2*(1. - df[j])*s0;
549 void calc_input_scale( const Mat& inputs, int flags )
551 bool reset_weights = (flags & UPDATE_WEIGHTS) == 0;
552 bool no_scale = (flags & NO_INPUT_SCALE) != 0;
553 double* scale = weights[0].ptr<double>();
554 int count = inputs.rows;
558 int i, j, vcount = layer_sizes[0];
559 int type = inputs.type();
560 double a = no_scale ? 1. : 0.;
562 for( j = 0; j < vcount; j++ )
563 scale[2*j] = a, scale[j*2+1] = 0.;
568 for( i = 0; i < count; i++ )
570 const uchar* p = inputs.ptr(i);
571 const float* f = (const float*)p;
572 const double* d = (const double*)p;
573 for( j = 0; j < vcount; j++ )
575 double t = type == CV_32F ? (double)f[j] : d[j];
581 for( j = 0; j < vcount; j++ )
583 double s = scale[j*2], s2 = scale[j*2+1];
584 double m = s/count, sigma2 = s2/count - m*m;
585 scale[j*2] = sigma2 < DBL_EPSILON ? 1 : 1./sqrt(sigma2);
586 scale[j*2+1] = -m*scale[j*2];
591 void calc_output_scale( const Mat& outputs, int flags )
593 int i, j, vcount = layer_sizes.back();
594 int type = outputs.type();
595 double m = min_val, M = max_val, m1 = min_val1, M1 = max_val1;
596 bool reset_weights = (flags & UPDATE_WEIGHTS) == 0;
597 bool no_scale = (flags & NO_OUTPUT_SCALE) != 0;
598 int l_count = layer_count();
599 double* scale = weights[l_count].ptr<double>();
600 double* inv_scale = weights[l_count+1].ptr<double>();
601 int count = outputs.rows;
605 double a0 = no_scale ? 1 : DBL_MAX, b0 = no_scale ? 0 : -DBL_MAX;
607 for( j = 0; j < vcount; j++ )
609 scale[2*j] = inv_scale[2*j] = a0;
610 scale[j*2+1] = inv_scale[2*j+1] = b0;
617 for( i = 0; i < count; i++ )
619 const uchar* p = outputs.ptr(i);
620 const float* f = (const float*)p;
621 const double* d = (const double*)p;
623 for( j = 0; j < vcount; j++ )
625 double t = type == CV_32F ? (double)f[j] : d[j];
629 double mj = scale[j*2], Mj = scale[j*2+1];
638 t = t*inv_scale[j*2] + inv_scale[2*j+1];
639 if( t < m1 || t > M1 )
640 CV_Error( CV_StsOutOfRange,
641 "Some of new output training vector components run exceed the original range too much" );
647 for( j = 0; j < vcount; j++ )
649 // map mj..Mj to m..M
650 double mj = scale[j*2], Mj = scale[j*2+1];
652 double delta = Mj - mj;
653 if( delta < DBL_EPSILON )
654 a = 1, b = (M + m - Mj - mj)*0.5;
656 a = (M - m)/delta, b = m - mj*a;
657 inv_scale[j*2] = a; inv_scale[j*2+1] = b;
659 scale[j*2] = a; scale[j*2+1] = b;
663 void prepare_to_train( const Mat& inputs, const Mat& outputs,
664 Mat& sample_weights, int flags )
666 if( layer_sizes.empty() )
667 CV_Error( CV_StsError,
668 "The network has not been created. Use method create or the appropriate constructor" );
670 if( (inputs.type() != CV_32F && inputs.type() != CV_64F) ||
671 inputs.cols != layer_sizes[0] )
672 CV_Error( CV_StsBadArg,
673 "input 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 0-th (input) layer" );
677 if( (outputs.type() != CV_32F && outputs.type() != CV_64F) ||
678 outputs.cols != layer_sizes.back() )
679 CV_Error( CV_StsBadArg,
680 "output training data should be a floating-point matrix with "
681 "the number of rows equal to the number of training samples and "
682 "the number of columns equal to the size of last (output) layer" );
684 if( inputs.rows != outputs.rows )
685 CV_Error( CV_StsUnmatchedSizes, "The numbers of input and output samples do not match" );
688 double s = sum(sample_weights)[0];
689 sample_weights.convertTo(temp, CV_64F, 1./s);
690 sample_weights = temp;
692 calc_input_scale( inputs, flags );
693 calc_output_scale( outputs, flags );
696 bool train( const Ptr<TrainData>& trainData, int flags )
698 const int MAX_ITER = 1000;
699 const double DEFAULT_EPSILON = FLT_EPSILON;
701 // initialize training data
702 Mat inputs = trainData->getTrainSamples();
703 Mat outputs = trainData->getTrainResponses();
704 Mat sw = trainData->getTrainSampleWeights();
705 prepare_to_train( inputs, outputs, sw, flags );
707 // ... and link weights
708 if( !(flags & UPDATE_WEIGHTS) )
711 TermCriteria termcrit;
712 termcrit.type = TermCriteria::COUNT + TermCriteria::EPS;
713 termcrit.maxCount = std::max((params.termCrit.type & CV_TERMCRIT_ITER ? params.termCrit.maxCount : MAX_ITER), 1);
714 termcrit.epsilon = std::max((params.termCrit.type & CV_TERMCRIT_EPS ? params.termCrit.epsilon : DEFAULT_EPSILON), DBL_EPSILON);
716 int iter = params.trainMethod == ANN_MLP::BACKPROP ?
717 train_backprop( inputs, outputs, sw, termcrit ) :
718 train_rprop( inputs, outputs, sw, termcrit );
724 int train_backprop( const Mat& inputs, const Mat& outputs, const Mat& _sw, TermCriteria termCrit )
727 double prev_E = DBL_MAX*0.5, E = 0;
728 int itype = inputs.type(), otype = outputs.type();
730 int count = inputs.rows;
732 int iter = -1, max_iter = termCrit.maxCount*count;
733 double epsilon = termCrit.epsilon*count;
735 int l_count = layer_count();
736 int ivcount = layer_sizes[0];
737 int ovcount = layer_sizes.back();
740 vector<vector<double> > x(l_count);
741 vector<vector<double> > df(l_count);
742 vector<Mat> dw(l_count);
744 for( i = 0; i < l_count; i++ )
746 int n = layer_sizes[i];
749 dw[i] = Mat::zeros(weights[i].size(), CV_64F);
752 Mat _idx_m(1, count, CV_32S);
753 int* _idx = _idx_m.ptr<int>();
754 for( i = 0; i < count; i++ )
757 AutoBuffer<double> _buf(max_lsize*2);
758 double* buf[] = { _buf, (double*)_buf + max_lsize };
760 const double* sw = _sw.empty() ? 0 : _sw.ptr<double>();
762 // run back-propagation loop
766 E = 1/2*||u - x_N||^2
767 grad_N = (x_N - u)*f'(y_i)
768 dw_i(t) = momentum*dw_i(t-1) + dw_scale*x_{i-1}*grad_i
769 w_i(t+1) = w_i(t) + dw_i(t)
770 grad_{i-1} = w_i^t*grad_i
772 for( iter = 0; iter < max_iter; iter++ )
774 int idx = iter % count;
775 double sweight = sw ? count*sw[idx] : 1.;
779 //printf("%d. E = %g\n", iter/count, E);
780 if( fabs(prev_E - E) < epsilon )
786 for( i = 0; i < count; i++ )
788 j = rng.uniform(0, count);
789 k = rng.uniform(0, count);
790 std::swap(_idx[j], _idx[k]);
796 const uchar* x0data_p = inputs.ptr(idx);
797 const float* x0data_f = (const float*)x0data_p;
798 const double* x0data_d = (const double*)x0data_p;
800 double* w = weights[0].ptr<double>();
801 for( j = 0; j < ivcount; j++ )
802 x[0][j] = (itype == CV_32F ? (double)x0data_f[j] : x0data_d[j])*w[j*2] + w[j*2 + 1];
804 Mat x1( 1, ivcount, CV_64F, &x[0][0] );
806 // forward pass, compute y[i]=w*x[i-1], x[i]=f(y[i]), df[i]=f'(y[i])
807 for( i = 1; i < l_count; i++ )
809 int n = layer_sizes[i];
810 Mat x2(1, n, CV_64F, &x[i][0] );
811 Mat _w = weights[i].rowRange(0, x1.cols);
812 gemm(x1, _w, 1, noArray(), 0, x2);
813 Mat _df(1, n, CV_64F, &df[i][0] );
814 calc_activ_func_deriv( x2, _df, weights[i] );
818 Mat grad1( 1, ovcount, CV_64F, buf[l_count&1] );
819 w = weights[l_count+1].ptr<double>();
822 const uchar* udata_p = outputs.ptr(idx);
823 const float* udata_f = (const float*)udata_p;
824 const double* udata_d = (const double*)udata_p;
826 double* gdata = grad1.ptr<double>();
827 for( k = 0; k < ovcount; k++ )
829 double t = (otype == CV_32F ? (double)udata_f[k] : udata_d[k])*w[k*2] + w[k*2+1] - x[l_count-1][k];
830 gdata[k] = t*sweight;
835 // backward pass, update weights
836 for( i = l_count-1; i > 0; i-- )
838 int n1 = layer_sizes[i-1], n2 = layer_sizes[i];
839 Mat _df(1, n2, CV_64F, &df[i][0]);
840 multiply( grad1, _df, grad1 );
841 Mat _x(n1+1, 1, CV_64F, &x[i-1][0]);
843 gemm( _x, grad1, params.bpDWScale, dw[i], params.bpMomentScale, dw[i] );
844 add( weights[i], dw[i], weights[i] );
847 Mat grad2(1, n1, CV_64F, buf[i&1]);
848 Mat _w = weights[i].rowRange(0, n1);
849 gemm( grad1, _w, 1, noArray(), 0, grad2, GEMM_2_T );
859 struct RPropLoop : public ParallelLoopBody
861 RPropLoop(ANN_MLPImpl* _ann,
862 const Mat& _inputs, const Mat& _outputs, const Mat& _sw,
863 int _dcount0, vector<Mat>& _dEdw, double* _E)
868 sw = _sw.ptr<double>();
881 void operator()( const Range& range ) const
883 double inv_count = 1./inputs.rows;
884 int ivcount = ann->layer_sizes.front();
885 int ovcount = ann->layer_sizes.back();
886 int itype = inputs.type(), otype = outputs.type();
887 int count = inputs.rows;
888 int i, j, k, l_count = ann->layer_count();
889 vector<vector<double> > x(l_count);
890 vector<vector<double> > df(l_count);
891 vector<double> _buf(ann->max_lsize*dcount0*2);
892 double* buf[] = { &_buf[0], &_buf[ann->max_lsize*dcount0] };
895 for( i = 0; i < l_count; i++ )
897 x[i].resize(ann->layer_sizes[i]*dcount0);
898 df[i].resize(ann->layer_sizes[i]*dcount0);
901 for( int si = range.start; si < range.end; si++ )
903 int i0 = si*dcount0, i1 = std::min((si + 1)*dcount0, count);
904 int dcount = i1 - i0;
905 const double* w = ann->weights[0].ptr<double>();
907 // grab and preprocess input data
908 for( i = 0; i < dcount; i++ )
910 const uchar* x0data_p = inputs.ptr(i0 + i);
911 const float* x0data_f = (const float*)x0data_p;
912 const double* x0data_d = (const double*)x0data_p;
914 double* xdata = &x[0][i*ivcount];
915 for( j = 0; j < ivcount; j++ )
916 xdata[j] = (itype == CV_32F ? (double)x0data_f[j] : x0data_d[j])*w[j*2] + w[j*2+1];
918 Mat x1(dcount, ivcount, CV_64F, &x[0][0]);
920 // forward pass, compute y[i]=w*x[i-1], x[i]=f(y[i]), df[i]=f'(y[i])
921 for( i = 1; i < l_count; i++ )
923 Mat x2( dcount, ann->layer_sizes[i], CV_64F, &x[i][0] );
924 Mat _w = ann->weights[i].rowRange(0, x1.cols);
925 gemm( x1, _w, 1, noArray(), 0, x2 );
926 Mat _df( x2.size(), CV_64F, &df[i][0] );
927 ann->calc_activ_func_deriv( x2, _df, ann->weights[i] );
931 Mat grad1(dcount, ovcount, CV_64F, buf[l_count & 1]);
933 w = ann->weights[l_count+1].ptr<double>();
936 for( i = 0; i < dcount; i++ )
938 const uchar* udata_p = outputs.ptr(i0+i);
939 const float* udata_f = (const float*)udata_p;
940 const double* udata_d = (const double*)udata_p;
942 const double* xdata = &x[l_count-1][i*ovcount];
943 double* gdata = grad1.ptr<double>(i);
944 double sweight = sw ? sw[si+i] : inv_count, E1 = 0;
946 for( j = 0; j < ovcount; j++ )
948 double t = (otype == CV_32F ? (double)udata_f[j] : udata_d[j])*w[j*2] + w[j*2+1] - xdata[j];
949 gdata[j] = t*sweight;
955 for( i = l_count-1; i > 0; i-- )
957 int n1 = ann->layer_sizes[i-1], n2 = ann->layer_sizes[i];
958 Mat _df(dcount, n2, CV_64F, &df[i][0]);
959 multiply(grad1, _df, grad1);
962 AutoLock lock(ann->mtx);
963 Mat _dEdw = dEdw->at(i).rowRange(0, n1);
964 x1 = Mat(dcount, n1, CV_64F, &x[i-1][0]);
965 gemm(x1, grad1, 1, _dEdw, 1, _dEdw, GEMM_1_T);
967 // update bias part of dEdw
968 double* dst = dEdw->at(i).ptr<double>(n1);
969 for( k = 0; k < dcount; k++ )
971 const double* src = grad1.ptr<double>(k);
972 for( j = 0; j < n2; j++ )
977 Mat grad2( dcount, n1, CV_64F, buf[i&1] );
980 Mat _w = ann->weights[i].rowRange(0, n1);
981 gemm(grad1, _w, 1, noArray(), 0, grad2, GEMM_2_T);
987 AutoLock lock(ann->mtx);
993 int train_rprop( const Mat& inputs, const Mat& outputs, const Mat& _sw, TermCriteria termCrit )
995 const int max_buf_size = 1 << 16;
996 int i, iter = -1, count = inputs.rows;
998 double prev_E = DBL_MAX*0.5;
1000 int max_iter = termCrit.maxCount;
1001 double epsilon = termCrit.epsilon;
1002 double dw_plus = params.rpDWPlus;
1003 double dw_minus = params.rpDWMinus;
1004 double dw_min = params.rpDWMin;
1005 double dw_max = params.rpDWMax;
1007 int l_count = layer_count();
1010 vector<Mat> dw(l_count), dEdw(l_count), prev_dEdw_sign(l_count);
1013 for( i = 0; i < l_count; i++ )
1015 total += layer_sizes[i];
1016 dw[i].create(weights[i].size(), CV_64F);
1017 dw[i].setTo(Scalar::all(params.rpDW0));
1018 prev_dEdw_sign[i] = Mat::zeros(weights[i].size(), CV_8S);
1019 dEdw[i] = Mat::zeros(weights[i].size(), CV_64F);
1022 int dcount0 = max_buf_size/(2*total);
1023 dcount0 = std::max( dcount0, 1 );
1024 dcount0 = std::min( dcount0, count );
1025 int chunk_count = (count + dcount0 - 1)/dcount0;
1029 y_i(t) = w_i(t)*x_{i-1}(t)
1031 E = sum_over_all_samples(1/2*||u - x_N||^2)
1032 grad_N = (x_N - u)*f'(y_i)
1034 std::min(dw_i{jk}(t)*dw_plus, dw_max), if dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) > 0
1035 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
1038 if (dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) < 0)
1041 w_i{jk}(t+1) = w_i{jk}(t) + dw_i{jk}(t)
1042 grad_{i-1}(t) = w_i^t(t)*grad_i(t)
1044 for( iter = 0; iter < max_iter; iter++ )
1048 for( i = 0; i < l_count; i++ )
1049 dEdw[i].setTo(Scalar::all(0));
1051 // first, iterate through all the samples and compute dEdw
1052 RPropLoop invoker(this, inputs, outputs, _sw, dcount0, dEdw, &E);
1053 parallel_for_(Range(0, chunk_count), invoker);
1054 //invoker(Range(0, chunk_count));
1056 // now update weights
1057 for( i = 1; i < l_count; i++ )
1059 int n1 = layer_sizes[i-1], n2 = layer_sizes[i];
1060 for( int k = 0; k <= n1; k++ )
1062 CV_Assert(weights[i].size() == Size(n2, n1+1));
1063 double* wk = weights[i].ptr<double>(k);
1064 double* dwk = dw[i].ptr<double>(k);
1065 double* dEdwk = dEdw[i].ptr<double>(k);
1066 schar* prevEk = prev_dEdw_sign[i].ptr<schar>(k);
1068 for( int j = 0; j < n2; j++ )
1070 double Eval = dEdwk[j];
1071 double dval = dwk[j];
1072 double wval = wk[j];
1073 int s = CV_SIGN(Eval);
1074 int ss = prevEk[j]*s;
1078 dval = std::min( dval, dw_max );
1080 wk[j] = wval + dval*s;
1085 dval = std::max( dval, dw_min );
1088 wk[j] = wval + dval*s;
1092 prevEk[j] = (schar)s;
1093 wk[j] = wval + dval*s;
1100 //printf("%d. E = %g\n", iter, E);
1101 if( fabs(prev_E - E) < epsilon )
1109 void write_params( FileStorage& fs ) const
1111 const char* activ_func_name = activ_func == IDENTITY ? "IDENTITY" :
1112 activ_func == SIGMOID_SYM ? "SIGMOID_SYM" :
1113 activ_func == GAUSSIAN ? "GAUSSIAN" : 0;
1115 if( activ_func_name )
1116 fs << "activation_function" << activ_func_name;
1118 fs << "activation_function_id" << activ_func;
1120 if( activ_func != IDENTITY )
1122 fs << "f_param1" << f_param1;
1123 fs << "f_param2" << f_param2;
1126 fs << "min_val" << min_val << "max_val" << max_val << "min_val1" << min_val1 << "max_val1" << max_val1;
1128 fs << "training_params" << "{";
1129 if( params.trainMethod == ANN_MLP::BACKPROP )
1131 fs << "train_method" << "BACKPROP";
1132 fs << "dw_scale" << params.bpDWScale;
1133 fs << "moment_scale" << params.bpMomentScale;
1135 else if( params.trainMethod == ANN_MLP::RPROP )
1137 fs << "train_method" << "RPROP";
1138 fs << "dw0" << params.rpDW0;
1139 fs << "dw_plus" << params.rpDWPlus;
1140 fs << "dw_minus" << params.rpDWMinus;
1141 fs << "dw_min" << params.rpDWMin;
1142 fs << "dw_max" << params.rpDWMax;
1145 CV_Error(CV_StsError, "Unknown training method");
1147 fs << "term_criteria" << "{";
1148 if( params.termCrit.type & TermCriteria::EPS )
1149 fs << "epsilon" << params.termCrit.epsilon;
1150 if( params.termCrit.type & TermCriteria::COUNT )
1151 fs << "iterations" << params.termCrit.maxCount;
1155 void write( FileStorage& fs ) const
1157 if( layer_sizes.empty() )
1159 int i, l_count = layer_count();
1162 fs << "layer_sizes" << layer_sizes;
1166 size_t esz = weights[0].elemSize();
1168 fs << "input_scale" << "[";
1169 fs.writeRaw("d", weights[0].ptr(), weights[0].total()*esz);
1171 fs << "]" << "output_scale" << "[";
1172 fs.writeRaw("d", weights[l_count].ptr(), weights[l_count].total()*esz);
1174 fs << "]" << "inv_output_scale" << "[";
1175 fs.writeRaw("d", weights[l_count+1].ptr(), weights[l_count+1].total()*esz);
1177 fs << "]" << "weights" << "[";
1178 for( i = 1; i < l_count; i++ )
1181 fs.writeRaw("d", weights[i].ptr(), weights[i].total()*esz);
1187 void read_params( const FileNode& fn )
1189 String activ_func_name = (String)fn["activation_function"];
1190 if( !activ_func_name.empty() )
1192 activ_func = activ_func_name == "SIGMOID_SYM" ? SIGMOID_SYM :
1193 activ_func_name == "IDENTITY" ? IDENTITY :
1194 activ_func_name == "GAUSSIAN" ? GAUSSIAN : -1;
1195 CV_Assert( activ_func >= 0 );
1198 activ_func = (int)fn["activation_function_id"];
1200 f_param1 = (double)fn["f_param1"];
1201 f_param2 = (double)fn["f_param2"];
1203 setActivationFunction( activ_func, f_param1, f_param2 );
1205 min_val = (double)fn["min_val"];
1206 max_val = (double)fn["max_val"];
1207 min_val1 = (double)fn["min_val1"];
1208 max_val1 = (double)fn["max_val1"];
1210 FileNode tpn = fn["training_params"];
1211 params = AnnParams();
1215 String tmethod_name = (String)tpn["train_method"];
1217 if( tmethod_name == "BACKPROP" )
1219 params.trainMethod = ANN_MLP::BACKPROP;
1220 params.bpDWScale = (double)tpn["dw_scale"];
1221 params.bpMomentScale = (double)tpn["moment_scale"];
1223 else if( tmethod_name == "RPROP" )
1225 params.trainMethod = ANN_MLP::RPROP;
1226 params.rpDW0 = (double)tpn["dw0"];
1227 params.rpDWPlus = (double)tpn["dw_plus"];
1228 params.rpDWMinus = (double)tpn["dw_minus"];
1229 params.rpDWMin = (double)tpn["dw_min"];
1230 params.rpDWMax = (double)tpn["dw_max"];
1233 CV_Error(CV_StsParseError, "Unknown training method (should be BACKPROP or RPROP)");
1235 FileNode tcn = tpn["term_criteria"];
1238 FileNode tcn_e = tcn["epsilon"];
1239 FileNode tcn_i = tcn["iterations"];
1240 params.termCrit.type = 0;
1241 if( !tcn_e.empty() )
1243 params.termCrit.type |= TermCriteria::EPS;
1244 params.termCrit.epsilon = (double)tcn_e;
1246 if( !tcn_i.empty() )
1248 params.termCrit.type |= TermCriteria::COUNT;
1249 params.termCrit.maxCount = (int)tcn_i;
1255 void read( const FileNode& fn )
1259 vector<int> _layer_sizes;
1260 readVectorOrMat(fn["layer_sizes"], _layer_sizes);
1261 setLayerSizes( _layer_sizes );
1263 int i, l_count = layer_count();
1266 size_t esz = weights[0].elemSize();
1268 FileNode w = fn["input_scale"];
1269 w.readRaw("d", weights[0].ptr(), weights[0].total()*esz);
1271 w = fn["output_scale"];
1272 w.readRaw("d", weights[l_count].ptr(), weights[l_count].total()*esz);
1274 w = fn["inv_output_scale"];
1275 w.readRaw("d", weights[l_count+1].ptr(), weights[l_count+1].total()*esz);
1277 FileNodeIterator w_it = fn["weights"].begin();
1279 for( i = 1; i < l_count; i++, ++w_it )
1280 (*w_it).readRaw("d", weights[i].ptr(), weights[i].total()*esz);
1284 Mat getWeights(int layerIdx) const
1286 CV_Assert( 0 <= layerIdx && layerIdx < (int)weights.size() );
1287 return weights[layerIdx];
1290 bool isTrained() const
1295 bool isClassifier() const
1300 int getVarCount() const
1302 return layer_sizes.empty() ? 0 : layer_sizes[0];
1305 String getDefaultName() const
1307 return "opencv_ml_ann_mlp";
1310 vector<int> layer_sizes;
1311 vector<Mat> weights;
1312 double f_param1, f_param2;
1313 double min_val, max_val, min_val1, max_val1;
1315 int max_lsize, max_buf_sz;
1323 Ptr<ANN_MLP> ANN_MLP::create()
1325 return makePtr<ANN_MLPImpl>();
1328 Ptr<ANN_MLP> ANN_MLP::load(const String& filepath)
1331 fs.open(filepath, FileStorage::READ);
1332 CV_Assert(fs.isOpened());
1333 Ptr<ANN_MLP> ann = makePtr<ANN_MLPImpl>();
1335 ((ANN_MLPImpl*)ann.get())->read(fs.getFirstTopLevelNode());