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 {
AnnParamscv::ml::AnnParams47     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         initialT=10;finalT=0.1,coolingRatio=0.95;itePerStep=10;
55         rEnergy = cv::RNG(12345);
56     }
57 
58     TermCriteria termCrit;
59     int trainMethod;
60 
61     double bpDWScale;
62     double bpMomentScale;
63 
64     double rpDW0;
65     double rpDWPlus;
66     double rpDWMinus;
67     double rpDWMin;
68     double rpDWMax;
69 
70     double initialT;
71     double finalT;
72     double coolingRatio;
73     int itePerStep;
74     RNG rEnergy;
75 };
76 
77 template <typename T>
inBounds(T val,T min_val,T max_val)78 inline T inBounds(T val, T min_val, T max_val)
79 {
80     return std::min(std::max(val, min_val), max_val);
81 }
82 
83 class SimulatedAnnealingANN_MLP
84 {
85 protected:
86     ml::ANN_MLP& nn;
87     Ptr<ml::TrainData> data;
88     int nbVariables;
89     vector<double*> adrVariables;
90     RNG rVar;
91     RNG rIndex;
92     double varTmp;
93     int index;
94 public:
SimulatedAnnealingANN_MLP(ml::ANN_MLP & x,const Ptr<ml::TrainData> & d)95     SimulatedAnnealingANN_MLP(ml::ANN_MLP& x, const Ptr<ml::TrainData>& d) : nn(x), data(d), varTmp(0.0), index(0)
96     {
97         initVarMap();
98     }
~SimulatedAnnealingANN_MLP()99     ~SimulatedAnnealingANN_MLP() {}
100 
changeState()101     void changeState()
102     {
103         index = rIndex.uniform(0, nbVariables);
104         double dv = rVar.uniform(-1.0, 1.0);
105         varTmp = *adrVariables[index];
106         *adrVariables[index] = dv;
107     }
108 
reverseState()109     void reverseState()
110     {
111         *adrVariables[index] = varTmp;
112     }
113 
energy() const114     double energy() const { return nn.calcError(data, false, noArray()); }
115 
116 protected:
initVarMap()117     void initVarMap()
118     {
119         Mat l = nn.getLayerSizes();
120         nbVariables = 0;
121         adrVariables.clear();
122         for (int i = 1; i < l.rows-1; i++)
123         {
124             Mat w = nn.getWeights(i);
125             for (int j = 0; j < w.rows; j++)
126             {
127                 for (int k = 0; k < w.cols; k++, nbVariables++)
128                 {
129                     if (j == w.rows - 1)
130                     {
131                         adrVariables.push_back(&w.at<double>(w.rows - 1, k));
132                     }
133                     else
134                     {
135                         adrVariables.push_back(&w.at<double>(j, k));
136                     }
137                 }
138             }
139         }
140     }
141 
142 };
143 
144 class ANN_MLPImpl CV_FINAL : public ANN_MLP
145 {
146 public:
ANN_MLPImpl()147     ANN_MLPImpl()
148     {
149         clear();
150         setActivationFunction( SIGMOID_SYM, 0, 0);
151         setLayerSizes(Mat());
152         setTrainMethod(ANN_MLP::RPROP, 0.1, FLT_EPSILON);
153     }
154 
~ANN_MLPImpl()155     virtual ~ANN_MLPImpl() CV_OVERRIDE {}
156 
getTermCriteria() const157     inline TermCriteria getTermCriteria() const CV_OVERRIDE { return params.termCrit; }
setTermCriteria(TermCriteria val)158     inline void setTermCriteria(TermCriteria val) CV_OVERRIDE { params.termCrit = val; }
getBackpropWeightScale() const159     inline double getBackpropWeightScale() const CV_OVERRIDE { return params.bpDWScale; }
setBackpropWeightScale(double val)160     inline void setBackpropWeightScale(double val) CV_OVERRIDE { params.bpDWScale = val; }
getBackpropMomentumScale() const161     inline double getBackpropMomentumScale() const CV_OVERRIDE { return params.bpMomentScale; }
setBackpropMomentumScale(double val)162     inline void setBackpropMomentumScale(double val) CV_OVERRIDE { params.bpMomentScale = val; }
getRpropDW0() const163     inline double getRpropDW0() const CV_OVERRIDE { return params.rpDW0; }
setRpropDW0(double val)164     inline void setRpropDW0(double val) CV_OVERRIDE { params.rpDW0 = val; }
getRpropDWPlus() const165     inline double getRpropDWPlus() const CV_OVERRIDE { return params.rpDWPlus; }
setRpropDWPlus(double val)166     inline void setRpropDWPlus(double val) CV_OVERRIDE { params.rpDWPlus = val; }
getRpropDWMinus() const167     inline double getRpropDWMinus() const CV_OVERRIDE { return params.rpDWMinus; }
setRpropDWMinus(double val)168     inline void setRpropDWMinus(double val) CV_OVERRIDE { params.rpDWMinus = val; }
getRpropDWMin() const169     inline double getRpropDWMin() const CV_OVERRIDE { return params.rpDWMin; }
setRpropDWMin(double val)170     inline void setRpropDWMin(double val) CV_OVERRIDE { params.rpDWMin = val; }
getRpropDWMax() const171     inline double getRpropDWMax() const CV_OVERRIDE { return params.rpDWMax; }
setRpropDWMax(double val)172     inline void setRpropDWMax(double val) CV_OVERRIDE { params.rpDWMax = val; }
getAnnealInitialT() const173     inline double getAnnealInitialT() const CV_OVERRIDE { return params.initialT; }
setAnnealInitialT(double val)174     inline void setAnnealInitialT(double val) CV_OVERRIDE { params.initialT = val; }
getAnnealFinalT() const175     inline double getAnnealFinalT() const CV_OVERRIDE { return params.finalT; }
setAnnealFinalT(double val)176     inline void setAnnealFinalT(double val) CV_OVERRIDE { params.finalT = val; }
getAnnealCoolingRatio() const177     inline double getAnnealCoolingRatio() const CV_OVERRIDE { return params.coolingRatio; }
setAnnealCoolingRatio(double val)178     inline void setAnnealCoolingRatio(double val) CV_OVERRIDE { params.coolingRatio = val; }
getAnnealItePerStep() const179     inline int getAnnealItePerStep() const CV_OVERRIDE { return params.itePerStep; }
setAnnealItePerStep(int val)180     inline void setAnnealItePerStep(int val) CV_OVERRIDE { params.itePerStep = val; }
181     // disabled getAnnealEnergyRNG()
setAnnealEnergyRNG(const RNG & val)182     inline void setAnnealEnergyRNG(const RNG& val) CV_OVERRIDE { params.rEnergy = val; }
183 
clear()184     void clear() CV_OVERRIDE
185     {
186         min_val = max_val = min_val1 = max_val1 = 0.;
187         rng = RNG((uint64)-1);
188         weights.clear();
189         trained = false;
190         max_buf_sz = 1 << 12;
191     }
192 
layer_count() const193     int layer_count() const { return (int)layer_sizes.size(); }
194 
setTrainMethod(int method,double param1,double param2)195     void setTrainMethod(int method, double param1, double param2) CV_OVERRIDE
196     {
197         if (method != ANN_MLP::RPROP && method != ANN_MLP::BACKPROP && method != ANN_MLP::ANNEAL)
198             method = ANN_MLP::RPROP;
199         params.trainMethod = method;
200         if(method == ANN_MLP::RPROP )
201         {
202             if( param1 < FLT_EPSILON )
203                 param1 = 1.;
204             params.rpDW0 = param1;
205             params.rpDWMin = std::max( param2, 0. );
206         }
207         else if (method == ANN_MLP::BACKPROP)
208         {
209             if (param1 <= 0)
210                 param1 = 0.1;
211             params.bpDWScale = inBounds<double>(param1, 1e-3, 1.);
212             if (param2 < 0)
213                 param2 = 0.1;
214             params.bpMomentScale = std::min(param2, 1.);
215         }
216     }
217 
getTrainMethod() const218     int getTrainMethod() const CV_OVERRIDE
219     {
220         return params.trainMethod;
221     }
222 
setActivationFunction(int _activ_func,double _f_param1,double _f_param2)223     void setActivationFunction(int _activ_func, double _f_param1, double _f_param2) CV_OVERRIDE
224     {
225         if( _activ_func < 0 || _activ_func > LEAKYRELU)
226             CV_Error( CV_StsOutOfRange, "Unknown activation function" );
227 
228         activ_func = _activ_func;
229 
230         switch( activ_func )
231         {
232         case SIGMOID_SYM:
233             max_val = 0.95; min_val = -max_val;
234             max_val1 = 0.98; min_val1 = -max_val1;
235             if( fabs(_f_param1) < FLT_EPSILON )
236                 _f_param1 = 2./3;
237             if( fabs(_f_param2) < FLT_EPSILON )
238                 _f_param2 = 1.7159;
239             break;
240         case GAUSSIAN:
241             max_val = 1.; min_val = 0.05;
242             max_val1 = 1.; min_val1 = 0.02;
243             if (fabs(_f_param1) < FLT_EPSILON)
244                 _f_param1 = 1.;
245             if (fabs(_f_param2) < FLT_EPSILON)
246                 _f_param2 = 1.;
247             break;
248         case RELU:
249             if (fabs(_f_param1) < FLT_EPSILON)
250                 _f_param1 = 1;
251             min_val = max_val = min_val1 = max_val1 = 0.;
252             _f_param2 = 0.;
253             break;
254         case LEAKYRELU:
255             if (fabs(_f_param1) < FLT_EPSILON)
256                 _f_param1 = 0.01;
257             min_val = max_val = min_val1 = max_val1 = 0.;
258             _f_param2 = 0.;
259             break;
260         default:
261             min_val = max_val = min_val1 = max_val1 = 0.;
262             _f_param1 = 1.;
263             _f_param2 = 0.;
264         }
265 
266         f_param1 = _f_param1;
267         f_param2 = _f_param2;
268     }
269 
270 
init_weights()271     void init_weights()
272     {
273         int i, j, k, l_count = layer_count();
274 
275         for( i = 1; i < l_count; i++ )
276         {
277             int n1 = layer_sizes[i-1];
278             int n2 = layer_sizes[i];
279             double val = 0, G = n2 > 2 ? 0.7*pow((double)n1,1./(n2-1)) : 1.;
280             double* w = weights[i].ptr<double>();
281 
282             // initialize weights using Nguyen-Widrow algorithm
283             for( j = 0; j < n2; j++ )
284             {
285                 double s = 0;
286                 for( k = 0; k <= n1; k++ )
287                 {
288                     val = rng.uniform(0., 1.)*2-1.;
289                     w[k*n2 + j] = val;
290                     s += fabs(val);
291                 }
292 
293                 if( i < l_count - 1 )
294                 {
295                     s = 1./(s - fabs(val));
296                     for( k = 0; k <= n1; k++ )
297                         w[k*n2 + j] *= s;
298                     w[n1*n2 + j] *= G*(-1+j*2./n2);
299                 }
300             }
301         }
302     }
303 
getLayerSizes() const304     Mat getLayerSizes() const CV_OVERRIDE
305     {
306         return Mat_<int>(layer_sizes, true);
307     }
308 
setLayerSizes(InputArray _layer_sizes)309     void setLayerSizes( InputArray _layer_sizes ) CV_OVERRIDE
310     {
311         clear();
312 
313         _layer_sizes.copyTo(layer_sizes);
314         int l_count = layer_count();
315 
316         weights.resize(l_count + 2);
317         max_lsize = 0;
318 
319         if( l_count > 0 )
320         {
321             for( int i = 0; i < l_count; i++ )
322             {
323                 int n = layer_sizes[i];
324                 if( n < 1 + (0 < i && i < l_count-1))
325                     CV_Error( CV_StsOutOfRange,
326                              "there should be at least one input and one output "
327                              "and every hidden layer must have more than 1 neuron" );
328                 max_lsize = std::max( max_lsize, n );
329                 if( i > 0 )
330                     weights[i].create(layer_sizes[i-1]+1, n, CV_64F);
331             }
332 
333             int ninputs = layer_sizes.front();
334             int noutputs = layer_sizes.back();
335             weights[0].create(1, ninputs*2, CV_64F);
336             weights[l_count].create(1, noutputs*2, CV_64F);
337             weights[l_count+1].create(1, noutputs*2, CV_64F);
338         }
339     }
340 
predict(InputArray _inputs,OutputArray _outputs,int) const341     float predict( InputArray _inputs, OutputArray _outputs, int ) const CV_OVERRIDE
342     {
343         if( !trained )
344             CV_Error( CV_StsError, "The network has not been trained or loaded" );
345 
346         Mat inputs = _inputs.getMat();
347         int type = inputs.type(), l_count = layer_count();
348         int n = inputs.rows, dn0 = n;
349 
350         CV_Assert( (type == CV_32F || type == CV_64F) && inputs.cols == layer_sizes[0] );
351         int noutputs = layer_sizes[l_count-1];
352         Mat outputs;
353 
354         int min_buf_sz = 2*max_lsize;
355         int buf_sz = n*min_buf_sz;
356 
357         if( buf_sz > max_buf_sz )
358         {
359             dn0 = max_buf_sz/min_buf_sz;
360             dn0 = std::max( dn0, 1 );
361             buf_sz = dn0*min_buf_sz;
362         }
363 
364         cv::AutoBuffer<double> _buf(buf_sz+noutputs);
365         double* buf = _buf.data();
366 
367         if( !_outputs.needed() )
368         {
369             CV_Assert( n == 1 );
370             outputs = Mat(n, noutputs, type, buf + buf_sz);
371         }
372         else
373         {
374             _outputs.create(n, noutputs, type);
375             outputs = _outputs.getMat();
376         }
377 
378         int dn = 0;
379         for( int i = 0; i < n; i += dn )
380         {
381             dn = std::min( dn0, n - i );
382 
383             Mat layer_in = inputs.rowRange(i, i + dn);
384             Mat layer_out( dn, layer_in.cols, CV_64F, buf);
385 
386             scale_input( layer_in, layer_out );
387             layer_in = layer_out;
388 
389             for( int j = 1; j < l_count; j++ )
390             {
391                 double* data = buf + ((j&1) ? max_lsize*dn0 : 0);
392                 int cols = layer_sizes[j];
393 
394                 layer_out = Mat(dn, cols, CV_64F, data);
395                 Mat w = weights[j].rowRange(0, layer_in.cols);
396                 gemm(layer_in, w, 1, noArray(), 0, layer_out);
397                 calc_activ_func( layer_out, weights[j] );
398 
399                 layer_in = layer_out;
400             }
401 
402             layer_out = outputs.rowRange(i, i + dn);
403             scale_output( layer_in, layer_out );
404         }
405 
406         if( n == 1 )
407         {
408             int maxIdx[] = {0, 0};
409             minMaxIdx(outputs, 0, 0, 0, maxIdx);
410             return (float)(maxIdx[0] + maxIdx[1]);
411         }
412 
413         return 0.f;
414     }
415 
scale_input(const Mat & _src,Mat & _dst) const416     void scale_input( const Mat& _src, Mat& _dst ) const
417     {
418         int cols = _src.cols;
419         const double* w = weights[0].ptr<double>();
420 
421         if( _src.type() == CV_32F )
422         {
423             for( int i = 0; i < _src.rows; i++ )
424             {
425                 const float* src = _src.ptr<float>(i);
426                 double* dst = _dst.ptr<double>(i);
427                 for( int j = 0; j < cols; j++ )
428                     dst[j] = src[j]*w[j*2] + w[j*2+1];
429             }
430         }
431         else
432         {
433             for( int i = 0; i < _src.rows; i++ )
434             {
435                 const double* src = _src.ptr<double>(i);
436                 double* dst = _dst.ptr<double>(i);
437                 for( int j = 0; j < cols; j++ )
438                     dst[j] = src[j]*w[j*2] + w[j*2+1];
439             }
440         }
441     }
442 
scale_output(const Mat & _src,Mat & _dst) const443     void scale_output( const Mat& _src, Mat& _dst ) const
444     {
445         int cols = _src.cols;
446         const double* w = weights[layer_count()].ptr<double>();
447 
448         if( _dst.type() == CV_32F )
449         {
450             for( int i = 0; i < _src.rows; i++ )
451             {
452                 const double* src = _src.ptr<double>(i);
453                 float* dst = _dst.ptr<float>(i);
454                 for( int j = 0; j < cols; j++ )
455                     dst[j] = (float)(src[j]*w[j*2] + w[j*2+1]);
456             }
457         }
458         else
459         {
460             for( int i = 0; i < _src.rows; i++ )
461             {
462                 const double* src = _src.ptr<double>(i);
463                 double* dst = _dst.ptr<double>(i);
464                 for( int j = 0; j < cols; j++ )
465                     dst[j] = src[j]*w[j*2] + w[j*2+1];
466             }
467         }
468     }
469 
calc_activ_func(Mat & sums,const Mat & w) const470     void calc_activ_func(Mat& sums, const Mat& w) const
471     {
472         const double* bias = w.ptr<double>(w.rows - 1);
473         int i, j, n = sums.rows, cols = sums.cols;
474         double scale = 0, scale2 = f_param2;
475 
476         switch (activ_func)
477         {
478         case IDENTITY:
479             scale = 1.;
480             break;
481         case SIGMOID_SYM:
482             scale = -f_param1;
483             break;
484         case GAUSSIAN:
485             scale = -f_param1*f_param1;
486             break;
487         case RELU:
488             scale = 1;
489             break;
490         case LEAKYRELU:
491             scale = 1;
492             break;
493         default:
494             ;
495         }
496 
497         CV_Assert(sums.isContinuous());
498 
499         if (activ_func != GAUSSIAN)
500         {
501             for (i = 0; i < n; i++)
502             {
503                 double* data = sums.ptr<double>(i);
504                 for (j = 0; j < cols; j++)
505                 {
506                     data[j] = (data[j] + bias[j])*scale;
507                     if (activ_func == RELU)
508                         if (data[j] < 0)
509                             data[j] = 0;
510                     if (activ_func == LEAKYRELU)
511                         if (data[j] < 0)
512                             data[j] *= f_param1;
513                 }
514             }
515 
516             if (activ_func == IDENTITY || activ_func == RELU || activ_func == LEAKYRELU)
517                 return;
518         }
519         else
520         {
521             for (i = 0; i < n; i++)
522             {
523                 double* data = sums.ptr<double>(i);
524                 for (j = 0; j < cols; j++)
525                 {
526                     double t = data[j] + bias[j];
527                     data[j] = t*t*scale;
528                 }
529             }
530         }
531 
532         exp(sums, sums);
533 
534         if (sums.isContinuous())
535         {
536             cols *= n;
537             n = 1;
538         }
539 
540         switch (activ_func)
541         {
542         case SIGMOID_SYM:
543             for (i = 0; i < n; i++)
544             {
545                 double* data = sums.ptr<double>(i);
546                 for (j = 0; j < cols; j++)
547                 {
548                     if (!cvIsInf(data[j]))
549                     {
550                         double t = scale2*(1. - data[j]) / (1. + data[j]);
551                         data[j] = t;
552                     }
553                     else
554                     {
555                         data[j] = -scale2;
556                     }
557                 }
558             }
559             break;
560 
561         case GAUSSIAN:
562             for (i = 0; i < n; i++)
563             {
564                 double* data = sums.ptr<double>(i);
565                 for (j = 0; j < cols; j++)
566                     data[j] = scale2*data[j];
567             }
568             break;
569 
570         default:
571             ;
572         }
573     }
574 
calc_activ_func_deriv(Mat & _xf,Mat & _df,const Mat & w) const575     void calc_activ_func_deriv(Mat& _xf, Mat& _df, const Mat& w) const
576     {
577         const double* bias = w.ptr<double>(w.rows - 1);
578         int i, j, n = _xf.rows, cols = _xf.cols;
579 
580         if (activ_func == IDENTITY)
581         {
582             for (i = 0; i < n; i++)
583             {
584                 double* xf = _xf.ptr<double>(i);
585                 double* df = _df.ptr<double>(i);
586 
587                 for (j = 0; j < cols; j++)
588                 {
589                     xf[j] += bias[j];
590                     df[j] = 1;
591                 }
592             }
593         }
594         else if (activ_func == RELU)
595         {
596             for (i = 0; i < n; i++)
597             {
598                 double* xf = _xf.ptr<double>(i);
599                 double* df = _df.ptr<double>(i);
600 
601                 for (j = 0; j < cols; j++)
602                 {
603                     xf[j] += bias[j];
604                     if (xf[j] < 0)
605                     {
606                         xf[j] = 0;
607                         df[j] = 0;
608                     }
609                     else
610                         df[j] = 1;
611                 }
612             }
613         }
614         else if (activ_func == LEAKYRELU)
615         {
616             for (i = 0; i < n; i++)
617             {
618                 double* xf = _xf.ptr<double>(i);
619                 double* df = _df.ptr<double>(i);
620 
621                 for (j = 0; j < cols; j++)
622                 {
623                     xf[j] += bias[j];
624                     if (xf[j] < 0)
625                     {
626                         xf[j] = f_param1*xf[j];
627                         df[j] = f_param1;
628                     }
629                     else
630                         df[j] = 1;
631                 }
632             }
633         }
634         else if (activ_func == GAUSSIAN)
635         {
636             double scale = -f_param1*f_param1;
637             double scale2 = scale*f_param2;
638             for (i = 0; i < n; i++)
639             {
640                 double* xf = _xf.ptr<double>(i);
641                 double* df = _df.ptr<double>(i);
642 
643                 for (j = 0; j < cols; j++)
644                 {
645                     double t = xf[j] + bias[j];
646                     df[j] = t * 2 * scale2;
647                     xf[j] = t*t*scale;
648                 }
649             }
650             exp(_xf, _xf);
651 
652             for (i = 0; i < n; i++)
653             {
654                 double* xf = _xf.ptr<double>(i);
655                 double* df = _df.ptr<double>(i);
656 
657                 for (j = 0; j < cols; j++)
658                     df[j] *= xf[j];
659             }
660         }
661         else
662         {
663             double scale = f_param1;
664             double scale2 = f_param2;
665 
666             for (i = 0; i < n; i++)
667             {
668                 double* xf = _xf.ptr<double>(i);
669                 double* df = _df.ptr<double>(i);
670 
671                 for (j = 0; j < cols; j++)
672                 {
673                     xf[j] = (xf[j] + bias[j])*scale;
674                     df[j] = -fabs(xf[j]);
675                 }
676             }
677 
678             exp(_df, _df);
679 
680             // ((1+exp(-ax))^-1)'=a*((1+exp(-ax))^-2)*exp(-ax);
681             // ((1-exp(-ax))/(1+exp(-ax)))'=(a*exp(-ax)*(1+exp(-ax)) + a*exp(-ax)*(1-exp(-ax)))/(1+exp(-ax))^2=
682             // 2*a*exp(-ax)/(1+exp(-ax))^2
683             scale *= 2 * f_param2;
684             for (i = 0; i < n; i++)
685             {
686                 double* xf = _xf.ptr<double>(i);
687                 double* df = _df.ptr<double>(i);
688 
689                 for (j = 0; j < cols; j++)
690                 {
691                     int s0 = xf[j] > 0 ? 1 : -1;
692                     double t0 = 1. / (1. + df[j]);
693                     double t1 = scale*df[j] * t0*t0;
694                     t0 *= scale2*(1. - df[j])*s0;
695                     df[j] = t1;
696                     xf[j] = t0;
697                 }
698             }
699         }
700     }
701 
calc_input_scale(const Mat & inputs,int flags)702     void calc_input_scale( const Mat& inputs, int flags )
703     {
704         bool reset_weights = (flags & UPDATE_WEIGHTS) == 0;
705         bool no_scale = (flags & NO_INPUT_SCALE) != 0;
706         double* scale = weights[0].ptr<double>();
707         int count = inputs.rows;
708 
709         if( reset_weights )
710         {
711             int i, j, vcount = layer_sizes[0];
712             int type = inputs.type();
713             double a = no_scale ? 1. : 0.;
714 
715             for( j = 0; j < vcount; j++ )
716                 scale[2*j] = a, scale[j*2+1] = 0.;
717 
718             if( no_scale )
719                 return;
720 
721             for( i = 0; i < count; i++ )
722             {
723                 const uchar* p = inputs.ptr(i);
724                 const float* f = (const float*)p;
725                 const double* d = (const double*)p;
726                 for( j = 0; j < vcount; j++ )
727                 {
728                     double t = type == CV_32F ? (double)f[j] : d[j];
729                     scale[j*2] += t;
730                     scale[j*2+1] += t*t;
731                 }
732             }
733 
734             for( j = 0; j < vcount; j++ )
735             {
736                 double s = scale[j*2], s2 = scale[j*2+1];
737                 double m = s/count, sigma2 = s2/count - m*m;
738                 scale[j*2] = sigma2 < DBL_EPSILON ? 1 : 1./sqrt(sigma2);
739                 scale[j*2+1] = -m*scale[j*2];
740             }
741         }
742     }
743 
calc_output_scale(const Mat & outputs,int flags)744     void calc_output_scale( const Mat& outputs, int flags )
745     {
746         int i, j, vcount = layer_sizes.back();
747         int type = outputs.type();
748         double m = min_val, M = max_val, m1 = min_val1, M1 = max_val1;
749         bool reset_weights = (flags & UPDATE_WEIGHTS) == 0;
750         bool no_scale = (flags & NO_OUTPUT_SCALE) != 0;
751         int l_count = layer_count();
752         double* scale = weights[l_count].ptr<double>();
753         double* inv_scale = weights[l_count+1].ptr<double>();
754         int count = outputs.rows;
755 
756         if( reset_weights )
757         {
758             double a0 = no_scale ? 1 : DBL_MAX, b0 = no_scale ? 0 : -DBL_MAX;
759 
760             for( j = 0; j < vcount; j++ )
761             {
762                 scale[2*j] = inv_scale[2*j] = a0;
763                 scale[j*2+1] = inv_scale[2*j+1] = b0;
764             }
765 
766             if( no_scale )
767                 return;
768         }
769 
770         for( i = 0; i < count; i++ )
771         {
772             const uchar* p = outputs.ptr(i);
773             const float* f = (const float*)p;
774             const double* d = (const double*)p;
775 
776             for( j = 0; j < vcount; j++ )
777             {
778                 double t = type == CV_32F ? (double)f[j] : d[j];
779 
780                 if( reset_weights )
781                 {
782                     double mj = scale[j*2], Mj = scale[j*2+1];
783                     if( mj > t ) mj = t;
784                     if( Mj < t ) Mj = t;
785 
786                     scale[j*2] = mj;
787                     scale[j*2+1] = Mj;
788                 }
789                 else if( !no_scale )
790                 {
791                     t = t*inv_scale[j*2] + inv_scale[2*j+1];
792                     if( t < m1 || t > M1 )
793                         CV_Error( CV_StsOutOfRange,
794                                  "Some of new output training vector components run exceed the original range too much" );
795                 }
796             }
797         }
798 
799         if( reset_weights )
800             for( j = 0; j < vcount; j++ )
801             {
802                 // map mj..Mj to m..M
803                 double mj = scale[j*2], Mj = scale[j*2+1];
804                 double a, b;
805                 double delta = Mj - mj;
806                 if( delta < DBL_EPSILON )
807                     a = 1, b = (M + m - Mj - mj)*0.5;
808                 else
809                     a = (M - m)/delta, b = m - mj*a;
810                 inv_scale[j*2] = a; inv_scale[j*2+1] = b;
811                 a = 1./a; b = -b*a;
812                 scale[j*2] = a; scale[j*2+1] = b;
813             }
814     }
815 
prepare_to_train(const Mat & inputs,const Mat & outputs,Mat & sample_weights,int flags)816     void prepare_to_train( const Mat& inputs, const Mat& outputs,
817                            Mat& sample_weights, int flags )
818     {
819         if( layer_sizes.empty() )
820             CV_Error( CV_StsError,
821                      "The network has not been created. Use method create or the appropriate constructor" );
822 
823         if( (inputs.type() != CV_32F && inputs.type() != CV_64F) ||
824             inputs.cols != layer_sizes[0] )
825             CV_Error( CV_StsBadArg,
826                      "input training data should be a floating-point matrix with "
827                      "the number of rows equal to the number of training samples and "
828                      "the number of columns equal to the size of 0-th (input) layer" );
829 
830         if( (outputs.type() != CV_32F && outputs.type() != CV_64F) ||
831             outputs.cols != layer_sizes.back() )
832             CV_Error( CV_StsBadArg,
833                      "output training data should be a floating-point matrix with "
834                      "the number of rows equal to the number of training samples and "
835                      "the number of columns equal to the size of last (output) layer" );
836 
837         if( inputs.rows != outputs.rows )
838             CV_Error( CV_StsUnmatchedSizes, "The numbers of input and output samples do not match" );
839 
840         Mat temp;
841         double s = sum(sample_weights)[0];
842         sample_weights.convertTo(temp, CV_64F, 1./s);
843         sample_weights = temp;
844 
845         calc_input_scale( inputs, flags );
846         calc_output_scale( outputs, flags );
847     }
848 
train(const Ptr<TrainData> & trainData,int flags)849     bool train( const Ptr<TrainData>& trainData, int flags ) CV_OVERRIDE
850     {
851         CV_Assert(!trainData.empty());
852         const int MAX_ITER = 1000;
853         const double DEFAULT_EPSILON = FLT_EPSILON;
854 
855         // initialize training data
856         Mat inputs = trainData->getTrainSamples();
857         Mat outputs = trainData->getTrainResponses();
858         Mat sw = trainData->getTrainSampleWeights();
859         prepare_to_train( inputs, outputs, sw, flags );
860 
861         // ... and link weights
862         if( !(flags & UPDATE_WEIGHTS) )
863             init_weights();
864 
865         TermCriteria termcrit;
866         termcrit.type = TermCriteria::COUNT + TermCriteria::EPS;
867         termcrit.maxCount = std::max((params.termCrit.type & CV_TERMCRIT_ITER ? params.termCrit.maxCount : MAX_ITER), 1);
868         termcrit.epsilon = std::max((params.termCrit.type & CV_TERMCRIT_EPS ? params.termCrit.epsilon : DEFAULT_EPSILON), DBL_EPSILON);
869 
870         int iter = 0;
871         switch(params.trainMethod){
872         case ANN_MLP::BACKPROP:
873             iter = train_backprop(inputs, outputs, sw, termcrit);
874             break;
875         case ANN_MLP::RPROP:
876             iter = train_rprop(inputs, outputs, sw, termcrit);
877             break;
878         case ANN_MLP::ANNEAL:
879             iter = train_anneal(trainData);
880             break;
881         }
882         trained = iter > 0;
883         return trained;
884     }
train_anneal(const Ptr<TrainData> & trainData)885     int train_anneal(const Ptr<TrainData>& trainData)
886     {
887         CV_Assert(!trainData.empty());
888         SimulatedAnnealingANN_MLP s(*this, trainData);
889         trained = true; // Enable call to CalcError
890         int iter = simulatedAnnealingSolver(s, params.initialT, params.finalT, params.coolingRatio, params.itePerStep, NULL, params.rEnergy);
891         trained =false;
892         return iter + 1; // ensure that 'train()' call is always successful
893     }
894 
train_backprop(const Mat & inputs,const Mat & outputs,const Mat & _sw,TermCriteria termCrit)895     int train_backprop( const Mat& inputs, const Mat& outputs, const Mat& _sw, TermCriteria termCrit )
896     {
897         int i, j, k;
898         double prev_E = DBL_MAX*0.5, E = 0;
899         int itype = inputs.type(), otype = outputs.type();
900 
901         int count = inputs.rows;
902 
903         int iter = -1, max_iter = termCrit.maxCount*count;
904         double epsilon = (termCrit.type & CV_TERMCRIT_EPS) ? termCrit.epsilon*count : 0;
905 
906         int l_count = layer_count();
907         int ivcount = layer_sizes[0];
908         int ovcount = layer_sizes.back();
909 
910         // allocate buffers
911         vector<vector<double> > x(l_count);
912         vector<vector<double> > df(l_count);
913         vector<Mat> dw(l_count);
914 
915         for( i = 0; i < l_count; i++ )
916         {
917             int n = layer_sizes[i];
918             x[i].resize(n+1);
919             df[i].resize(n);
920             dw[i] = Mat::zeros(weights[i].size(), CV_64F);
921         }
922 
923         Mat _idx_m(1, count, CV_32S);
924         int* _idx = _idx_m.ptr<int>();
925         for( i = 0; i < count; i++ )
926             _idx[i] = i;
927 
928         AutoBuffer<double> _buf(max_lsize*2);
929         double* buf[] = { _buf.data(), _buf.data() + max_lsize };
930 
931         const double* sw = _sw.empty() ? 0 : _sw.ptr<double>();
932 
933         // run back-propagation loop
934         /*
935          y_i = w_i*x_{i-1}
936          x_i = f(y_i)
937          E = 1/2*||u - x_N||^2
938          grad_N = (x_N - u)*f'(y_i)
939          dw_i(t) = momentum*dw_i(t-1) + dw_scale*x_{i-1}*grad_i
940          w_i(t+1) = w_i(t) + dw_i(t)
941          grad_{i-1} = w_i^t*grad_i
942         */
943         for( iter = 0; iter < max_iter; iter++ )
944         {
945             int idx = iter % count;
946             double sweight = sw ? count*sw[idx] : 1.;
947 
948             if( idx == 0 )
949             {
950                 //printf("%d. E = %g\n", iter/count, E);
951                 if( fabs(prev_E - E) < epsilon )
952                     break;
953                 prev_E = E;
954                 E = 0;
955 
956                 // shuffle indices
957                 for( i = 0; i <count; i++ )
958                 {
959                     j = rng.uniform(0, count);
960                     k = rng.uniform(0, count);
961                     std::swap(_idx[j], _idx[k]);
962                 }
963             }
964 
965             idx = _idx[idx];
966 
967             const uchar* x0data_p = inputs.ptr(idx);
968             const float* x0data_f = (const float*)x0data_p;
969             const double* x0data_d = (const double*)x0data_p;
970 
971             double* w = weights[0].ptr<double>();
972             for( j = 0; j < ivcount; j++ )
973                 x[0][j] = (itype == CV_32F ? (double)x0data_f[j] : x0data_d[j])*w[j*2] + w[j*2 + 1];
974 
975             Mat x1( 1, ivcount, CV_64F, &x[0][0] );
976 
977             // forward pass, compute y[i]=w*x[i-1], x[i]=f(y[i]), df[i]=f'(y[i])
978             for( i = 1; i < l_count; i++ )
979             {
980                 int n = layer_sizes[i];
981                 Mat x2(1, n, CV_64F, &x[i][0] );
982                 Mat _w = weights[i].rowRange(0, x1.cols);
983                 gemm(x1, _w, 1, noArray(), 0, x2);
984                 Mat _df(1, n, CV_64F, &df[i][0] );
985                 calc_activ_func_deriv( x2, _df, weights[i] );
986                 x1 = x2;
987             }
988 
989             Mat grad1( 1, ovcount, CV_64F, buf[l_count&1] );
990             w = weights[l_count+1].ptr<double>();
991 
992             // calculate error
993             const uchar* udata_p = outputs.ptr(idx);
994             const float* udata_f = (const float*)udata_p;
995             const double* udata_d = (const double*)udata_p;
996 
997             double* gdata = grad1.ptr<double>();
998             for( k = 0; k < ovcount; k++ )
999             {
1000                 double t = (otype == CV_32F ? (double)udata_f[k] : udata_d[k])*w[k*2] + w[k*2+1] - x[l_count-1][k];
1001                 gdata[k] = t*sweight;
1002                 E += t*t;
1003             }
1004             E *= sweight;
1005 
1006             // backward pass, update weights
1007             for( i = l_count-1; i > 0; i-- )
1008             {
1009                 int n1 = layer_sizes[i-1], n2 = layer_sizes[i];
1010                 Mat _df(1, n2, CV_64F, &df[i][0]);
1011                 multiply( grad1, _df, grad1 );
1012                 Mat _x(n1+1, 1, CV_64F, &x[i-1][0]);
1013                 x[i-1][n1] = 1.;
1014                 gemm( _x, grad1, params.bpDWScale, dw[i], params.bpMomentScale, dw[i] );
1015                 add( weights[i], dw[i], weights[i] );
1016                 if( i > 1 )
1017                 {
1018                     Mat grad2(1, n1, CV_64F, buf[i&1]);
1019                     Mat _w = weights[i].rowRange(0, n1);
1020                     gemm( grad1, _w, 1, noArray(), 0, grad2, GEMM_2_T );
1021                     grad1 = grad2;
1022                 }
1023             }
1024         }
1025 
1026         iter /= count;
1027         return iter;
1028     }
1029 
1030     struct RPropLoop : public ParallelLoopBody
1031     {
RPropLoopcv::ml::CV_FINAL::RPropLoop1032         RPropLoop(ANN_MLPImpl* _ann,
1033                   const Mat& _inputs, const Mat& _outputs, const Mat& _sw,
1034                   int _dcount0, vector<Mat>& _dEdw, double* _E)
1035         {
1036             ann = _ann;
1037             inputs = _inputs;
1038             outputs = _outputs;
1039             sw = _sw.ptr<double>();
1040             dcount0 = _dcount0;
1041             dEdw = &_dEdw;
1042             pE = _E;
1043         }
1044 
1045         ANN_MLPImpl* ann;
1046         vector<Mat>* dEdw;
1047         Mat inputs, outputs;
1048         const double* sw;
1049         int dcount0;
1050         double* pE;
1051 
operator ()cv::ml::CV_FINAL::RPropLoop1052         void operator()(const Range& range) const CV_OVERRIDE
1053         {
1054             double inv_count = 1./inputs.rows;
1055             int ivcount = ann->layer_sizes.front();
1056             int ovcount = ann->layer_sizes.back();
1057             int itype = inputs.type(), otype = outputs.type();
1058             int count = inputs.rows;
1059             int i, j, k, l_count = ann->layer_count();
1060             vector<vector<double> > x(l_count);
1061             vector<vector<double> > df(l_count);
1062             vector<double> _buf(ann->max_lsize*dcount0*2);
1063             double* buf[] = { &_buf[0], &_buf[ann->max_lsize*dcount0] };
1064             double E = 0;
1065 
1066             for( i = 0; i < l_count; i++ )
1067             {
1068                 x[i].resize(ann->layer_sizes[i]*dcount0);
1069                 df[i].resize(ann->layer_sizes[i]*dcount0);
1070             }
1071 
1072             for( int si = range.start; si < range.end; si++ )
1073             {
1074                 int i0 = si*dcount0, i1 = std::min((si + 1)*dcount0, count);
1075                 int dcount = i1 - i0;
1076                 const double* w = ann->weights[0].ptr<double>();
1077 
1078                 // grab and preprocess input data
1079                 for( i = 0; i < dcount; i++ )
1080                 {
1081                     const uchar* x0data_p = inputs.ptr(i0 + i);
1082                     const float* x0data_f = (const float*)x0data_p;
1083                     const double* x0data_d = (const double*)x0data_p;
1084 
1085                     double* xdata = &x[0][i*ivcount];
1086                     for( j = 0; j < ivcount; j++ )
1087                         xdata[j] = (itype == CV_32F ? (double)x0data_f[j] : x0data_d[j])*w[j*2] + w[j*2+1];
1088                 }
1089                 Mat x1(dcount, ivcount, CV_64F, &x[0][0]);
1090 
1091                 // forward pass, compute y[i]=w*x[i-1], x[i]=f(y[i]), df[i]=f'(y[i])
1092                 for( i = 1; i < l_count; i++ )
1093                 {
1094                     Mat x2( dcount, ann->layer_sizes[i], CV_64F, &x[i][0] );
1095                     Mat _w = ann->weights[i].rowRange(0, x1.cols);
1096                     gemm( x1, _w, 1, noArray(), 0, x2 );
1097                     Mat _df( x2.size(), CV_64F, &df[i][0] );
1098                     ann->calc_activ_func_deriv( x2, _df, ann->weights[i] );
1099                     x1 = x2;
1100                 }
1101 
1102                 Mat grad1(dcount, ovcount, CV_64F, buf[l_count & 1]);
1103 
1104                 w = ann->weights[l_count+1].ptr<double>();
1105 
1106                 // calculate error
1107                 for( i = 0; i < dcount; i++ )
1108                 {
1109                     const uchar* udata_p = outputs.ptr(i0+i);
1110                     const float* udata_f = (const float*)udata_p;
1111                     const double* udata_d = (const double*)udata_p;
1112 
1113                     const double* xdata = &x[l_count-1][i*ovcount];
1114                     double* gdata = grad1.ptr<double>(i);
1115                     double sweight = sw ? sw[si+i] : inv_count, E1 = 0;
1116 
1117                     for( j = 0; j < ovcount; j++ )
1118                     {
1119                         double t = (otype == CV_32F ? (double)udata_f[j] : udata_d[j])*w[j*2] + w[j*2+1] - xdata[j];
1120                         gdata[j] = t*sweight;
1121                         E1 += t*t;
1122                     }
1123                     E += sweight*E1;
1124                 }
1125 
1126                 for( i = l_count-1; i > 0; i-- )
1127                 {
1128                     int n1 = ann->layer_sizes[i-1], n2 = ann->layer_sizes[i];
1129                     Mat _df(dcount, n2, CV_64F, &df[i][0]);
1130                     multiply(grad1, _df, grad1);
1131 
1132                     {
1133                         AutoLock lock(ann->mtx);
1134                         Mat _dEdw = dEdw->at(i).rowRange(0, n1);
1135                         x1 = Mat(dcount, n1, CV_64F, &x[i-1][0]);
1136                         gemm(x1, grad1, 1, _dEdw, 1, _dEdw, GEMM_1_T);
1137 
1138                         // update bias part of dEdw
1139                         double* dst = dEdw->at(i).ptr<double>(n1);
1140                         for( k = 0; k < dcount; k++ )
1141                         {
1142                             const double* src = grad1.ptr<double>(k);
1143                             for( j = 0; j < n2; j++ )
1144                                 dst[j] += src[j];
1145                         }
1146                     }
1147 
1148                     Mat grad2( dcount, n1, CV_64F, buf[i&1] );
1149                     if( i > 1 )
1150                     {
1151                         Mat _w = ann->weights[i].rowRange(0, n1);
1152                         gemm(grad1, _w, 1, noArray(), 0, grad2, GEMM_2_T);
1153                     }
1154                     grad1 = grad2;
1155                 }
1156             }
1157             {
1158                 AutoLock lock(ann->mtx);
1159                 *pE += E;
1160             }
1161         }
1162     };
1163 
train_rprop(const Mat & inputs,const Mat & outputs,const Mat & _sw,TermCriteria termCrit)1164     int train_rprop( const Mat& inputs, const Mat& outputs, const Mat& _sw, TermCriteria termCrit )
1165     {
1166         const int max_buf_size = 1 << 16;
1167         int i, iter = -1, count = inputs.rows;
1168 
1169         double prev_E = DBL_MAX*0.5;
1170 
1171         int max_iter = termCrit.maxCount;
1172         double epsilon = termCrit.epsilon;
1173         double dw_plus = params.rpDWPlus;
1174         double dw_minus = params.rpDWMinus;
1175         double dw_min = params.rpDWMin;
1176         double dw_max = params.rpDWMax;
1177 
1178         int l_count = layer_count();
1179 
1180         // allocate buffers
1181         vector<Mat> dw(l_count), dEdw(l_count), prev_dEdw_sign(l_count);
1182 
1183         int total = 0;
1184         for( i = 0; i < l_count; i++ )
1185         {
1186             total += layer_sizes[i];
1187             dw[i].create(weights[i].size(), CV_64F);
1188             dw[i].setTo(Scalar::all(params.rpDW0));
1189             prev_dEdw_sign[i] = Mat::zeros(weights[i].size(), CV_8S);
1190             dEdw[i] = Mat::zeros(weights[i].size(), CV_64F);
1191         }
1192         CV_Assert(total > 0);
1193         int dcount0 = max_buf_size/(2*total);
1194         dcount0 = std::max( dcount0, 1 );
1195         dcount0 = std::min( dcount0, count );
1196         int chunk_count = (count + dcount0 - 1)/dcount0;
1197 
1198         // run rprop loop
1199         /*
1200          y_i(t) = w_i(t)*x_{i-1}(t)
1201          x_i(t) = f(y_i(t))
1202          E = sum_over_all_samples(1/2*||u - x_N||^2)
1203          grad_N = (x_N - u)*f'(y_i)
1204 
1205          std::min(dw_i{jk}(t)*dw_plus, dw_max), if dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) > 0
1206          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
1207          dw_i{jk}(t-1) else
1208 
1209          if (dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) < 0)
1210          dE/dw_i{jk}(t)<-0
1211          else
1212          w_i{jk}(t+1) = w_i{jk}(t) + dw_i{jk}(t)
1213          grad_{i-1}(t) = w_i^t(t)*grad_i(t)
1214          */
1215         for( iter = 0; iter < max_iter; iter++ )
1216         {
1217             double E = 0;
1218 
1219             for( i = 0; i < l_count; i++ )
1220                 dEdw[i].setTo(Scalar::all(0));
1221 
1222             // first, iterate through all the samples and compute dEdw
1223             RPropLoop invoker(this, inputs, outputs, _sw, dcount0, dEdw, &E);
1224             parallel_for_(Range(0, chunk_count), invoker);
1225             //invoker(Range(0, chunk_count));
1226 
1227             // now update weights
1228             for( i = 1; i < l_count; i++ )
1229             {
1230                 int n1 = layer_sizes[i-1], n2 = layer_sizes[i];
1231                 for( int k = 0; k <= n1; k++ )
1232                 {
1233                     CV_Assert(weights[i].size() == Size(n2, n1+1));
1234                     double* wk = weights[i].ptr<double>(k);
1235                     double* dwk = dw[i].ptr<double>(k);
1236                     double* dEdwk = dEdw[i].ptr<double>(k);
1237                     schar* prevEk = prev_dEdw_sign[i].ptr<schar>(k);
1238 
1239                     for( int j = 0; j < n2; j++ )
1240                     {
1241                         double Eval = dEdwk[j];
1242                         double dval = dwk[j];
1243                         double wval = wk[j];
1244                         int s = CV_SIGN(Eval);
1245                         int ss = prevEk[j]*s;
1246                         if( ss > 0 )
1247                         {
1248                             dval *= dw_plus;
1249                             dval = std::min( dval, dw_max );
1250                             dwk[j] = dval;
1251                             wk[j] = wval + dval*s;
1252                         }
1253                         else if( ss < 0 )
1254                         {
1255                             dval *= dw_minus;
1256                             dval = std::max( dval, dw_min );
1257                             prevEk[j] = 0;
1258                             dwk[j] = dval;
1259                             wk[j] = wval + dval*s;
1260                         }
1261                         else
1262                         {
1263                             prevEk[j] = (schar)s;
1264                             wk[j] = wval + dval*s;
1265                         }
1266                         dEdwk[j] = 0.;
1267                     }
1268                 }
1269             }
1270 
1271             //printf("%d. E = %g\n", iter, E);
1272             if( fabs(prev_E - E) < epsilon )
1273                 break;
1274             prev_E = E;
1275         }
1276 
1277         return iter;
1278     }
1279 
write_params(FileStorage & fs) const1280     void write_params( FileStorage& fs ) const
1281     {
1282         const char* activ_func_name = activ_func == IDENTITY ? "IDENTITY" :
1283                                       activ_func == SIGMOID_SYM ? "SIGMOID_SYM" :
1284                                       activ_func == GAUSSIAN ? "GAUSSIAN" :
1285                                       activ_func == RELU ? "RELU" :
1286                                       activ_func == LEAKYRELU ? "LEAKYRELU" : 0;
1287 
1288         if( activ_func_name )
1289             fs << "activation_function" << activ_func_name;
1290         else
1291             fs << "activation_function_id" << activ_func;
1292 
1293         if( activ_func != IDENTITY )
1294         {
1295             fs << "f_param1" << f_param1;
1296             fs << "f_param2" << f_param2;
1297         }
1298 
1299         fs << "min_val" << min_val << "max_val" << max_val << "min_val1" << min_val1 << "max_val1" << max_val1;
1300 
1301         fs << "training_params" << "{";
1302         if( params.trainMethod == ANN_MLP::BACKPROP )
1303         {
1304             fs << "train_method" << "BACKPROP";
1305             fs << "dw_scale" << params.bpDWScale;
1306             fs << "moment_scale" << params.bpMomentScale;
1307         }
1308         else if (params.trainMethod == ANN_MLP::RPROP)
1309         {
1310             fs << "train_method" << "RPROP";
1311             fs << "dw0" << params.rpDW0;
1312             fs << "dw_plus" << params.rpDWPlus;
1313             fs << "dw_minus" << params.rpDWMinus;
1314             fs << "dw_min" << params.rpDWMin;
1315             fs << "dw_max" << params.rpDWMax;
1316         }
1317         else if (params.trainMethod == ANN_MLP::ANNEAL)
1318         {
1319             fs << "train_method" << "ANNEAL";
1320             fs << "initialT" << params.initialT;
1321             fs << "finalT" << params.finalT;
1322             fs << "coolingRatio" << params.coolingRatio;
1323             fs << "itePerStep" << params.itePerStep;
1324         }
1325         else
1326             CV_Error(CV_StsError, "Unknown training method");
1327 
1328         fs << "term_criteria" << "{";
1329         if( params.termCrit.type & TermCriteria::EPS )
1330             fs << "epsilon" << params.termCrit.epsilon;
1331         if( params.termCrit.type & TermCriteria::COUNT )
1332             fs << "iterations" << params.termCrit.maxCount;
1333         fs << "}" << "}";
1334     }
1335 
write(FileStorage & fs) const1336     void write( FileStorage& fs ) const CV_OVERRIDE
1337     {
1338         if( layer_sizes.empty() )
1339             return;
1340         int i, l_count = layer_count();
1341 
1342         writeFormat(fs);
1343         fs << "layer_sizes" << layer_sizes;
1344 
1345         write_params( fs );
1346 
1347         size_t esz = weights[0].elemSize();
1348 
1349         fs << "input_scale" << "[";
1350         fs.writeRaw("d", weights[0].ptr(), weights[0].total()*esz);
1351 
1352         fs << "]" << "output_scale" << "[";
1353         fs.writeRaw("d", weights[l_count].ptr(), weights[l_count].total()*esz);
1354 
1355         fs << "]" << "inv_output_scale" << "[";
1356         fs.writeRaw("d", weights[l_count+1].ptr(), weights[l_count+1].total()*esz);
1357 
1358         fs << "]" << "weights" << "[";
1359         for( i = 1; i < l_count; i++ )
1360         {
1361             fs << "[";
1362             fs.writeRaw("d", weights[i].ptr(), weights[i].total()*esz);
1363             fs << "]";
1364         }
1365         fs << "]";
1366     }
1367 
read_params(const FileNode & fn)1368     void read_params( const FileNode& fn )
1369     {
1370         String activ_func_name = (String)fn["activation_function"];
1371         if( !activ_func_name.empty() )
1372         {
1373             activ_func = activ_func_name == "SIGMOID_SYM" ? SIGMOID_SYM :
1374                          activ_func_name == "IDENTITY" ? IDENTITY :
1375                          activ_func_name == "RELU" ? RELU :
1376                          activ_func_name == "LEAKYRELU" ? LEAKYRELU :
1377                          activ_func_name == "GAUSSIAN" ? GAUSSIAN : -1;
1378             CV_Assert( activ_func >= 0 );
1379         }
1380         else
1381             activ_func = (int)fn["activation_function_id"];
1382 
1383         f_param1 = (double)fn["f_param1"];
1384         f_param2 = (double)fn["f_param2"];
1385 
1386         setActivationFunction( activ_func, f_param1, f_param2);
1387 
1388         min_val = (double)fn["min_val"];
1389         max_val = (double)fn["max_val"];
1390         min_val1 = (double)fn["min_val1"];
1391         max_val1 = (double)fn["max_val1"];
1392 
1393         FileNode tpn = fn["training_params"];
1394         params = AnnParams();
1395 
1396         if( !tpn.empty() )
1397         {
1398             String tmethod_name = (String)tpn["train_method"];
1399 
1400             if( tmethod_name == "BACKPROP" )
1401             {
1402                 params.trainMethod = ANN_MLP::BACKPROP;
1403                 params.bpDWScale = (double)tpn["dw_scale"];
1404                 params.bpMomentScale = (double)tpn["moment_scale"];
1405             }
1406             else if (tmethod_name == "RPROP")
1407             {
1408                 params.trainMethod = ANN_MLP::RPROP;
1409                 params.rpDW0 = (double)tpn["dw0"];
1410                 params.rpDWPlus = (double)tpn["dw_plus"];
1411                 params.rpDWMinus = (double)tpn["dw_minus"];
1412                 params.rpDWMin = (double)tpn["dw_min"];
1413                 params.rpDWMax = (double)tpn["dw_max"];
1414             }
1415             else if (tmethod_name == "ANNEAL")
1416             {
1417                 params.trainMethod = ANN_MLP::ANNEAL;
1418                 params.initialT = (double)tpn["initialT"];
1419                 params.finalT = (double)tpn["finalT"];
1420                 params.coolingRatio = (double)tpn["coolingRatio"];
1421                 params.itePerStep = tpn["itePerStep"];
1422             }
1423             else
1424                 CV_Error(CV_StsParseError, "Unknown training method (should be BACKPROP or RPROP)");
1425 
1426             FileNode tcn = tpn["term_criteria"];
1427             if( !tcn.empty() )
1428             {
1429                 FileNode tcn_e = tcn["epsilon"];
1430                 FileNode tcn_i = tcn["iterations"];
1431                 params.termCrit.type = 0;
1432                 if( !tcn_e.empty() )
1433                 {
1434                     params.termCrit.type |= TermCriteria::EPS;
1435                     params.termCrit.epsilon = (double)tcn_e;
1436                 }
1437                 if( !tcn_i.empty() )
1438                 {
1439                     params.termCrit.type |= TermCriteria::COUNT;
1440                     params.termCrit.maxCount = (int)tcn_i;
1441                 }
1442             }
1443         }
1444     }
1445 
read(const FileNode & fn)1446     void read( const FileNode& fn ) CV_OVERRIDE
1447     {
1448         clear();
1449 
1450         vector<int> _layer_sizes;
1451         readVectorOrMat(fn["layer_sizes"], _layer_sizes);
1452         setLayerSizes( _layer_sizes );
1453 
1454         int i, l_count = layer_count();
1455         read_params(fn);
1456 
1457         size_t esz = weights[0].elemSize();
1458 
1459         FileNode w = fn["input_scale"];
1460         w.readRaw("d", weights[0].ptr(), weights[0].total()*esz);
1461 
1462         w = fn["output_scale"];
1463         w.readRaw("d", weights[l_count].ptr(), weights[l_count].total()*esz);
1464 
1465         w = fn["inv_output_scale"];
1466         w.readRaw("d", weights[l_count+1].ptr(), weights[l_count+1].total()*esz);
1467 
1468         FileNodeIterator w_it = fn["weights"].begin();
1469 
1470         for( i = 1; i < l_count; i++, ++w_it )
1471             (*w_it).readRaw("d", weights[i].ptr(), weights[i].total()*esz);
1472         trained = true;
1473     }
1474 
getWeights(int layerIdx) const1475     Mat getWeights(int layerIdx) const CV_OVERRIDE
1476     {
1477         CV_Assert( 0 <= layerIdx && layerIdx < (int)weights.size() );
1478         return weights[layerIdx];
1479     }
1480 
isTrained() const1481     bool isTrained() const CV_OVERRIDE
1482     {
1483         return trained;
1484     }
1485 
isClassifier() const1486     bool isClassifier() const CV_OVERRIDE
1487     {
1488         return false;
1489     }
1490 
getVarCount() const1491     int getVarCount() const CV_OVERRIDE
1492     {
1493         return layer_sizes.empty() ? 0 : layer_sizes[0];
1494     }
1495 
getDefaultName() const1496     String getDefaultName() const CV_OVERRIDE
1497     {
1498         return "opencv_ml_ann_mlp";
1499     }
1500 
1501     vector<int> layer_sizes;
1502     vector<Mat> weights;
1503     double f_param1, f_param2;
1504     double min_val, max_val, min_val1, max_val1;
1505     int activ_func;
1506     int max_lsize, max_buf_sz;
1507     AnnParams params;
1508     RNG rng;
1509     Mutex mtx;
1510     bool trained;
1511 };
1512 
1513 
1514 
1515 
create()1516 Ptr<ANN_MLP> ANN_MLP::create()
1517 {
1518     return makePtr<ANN_MLPImpl>();
1519 }
1520 
load(const String & filepath)1521 Ptr<ANN_MLP> ANN_MLP::load(const String& filepath)
1522 {
1523     FileStorage fs;
1524     fs.open(filepath, FileStorage::READ);
1525     CV_Assert(fs.isOpened());
1526     Ptr<ANN_MLP> ann = makePtr<ANN_MLPImpl>();
1527     ((ANN_MLPImpl*)ann.get())->read(fs.getFirstTopLevelNode());
1528     return ann;
1529 }
1530 
1531 }}
1532 
1533 /* End of file. */
1534