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