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 CvANN_MLP_TrainParams::CvANN_MLP_TrainParams()
  44 {
  45     term_crit = cvTermCriteria( CV_TERMCRIT_ITER + CV_TERMCRIT_EPS, 1000, 0.01 );
  46     train_method = RPROP;
  47     bp_dw_scale = bp_moment_scale = 0.1;
  48     rp_dw0 = 0.1; rp_dw_plus = 1.2; rp_dw_minus = 0.5;
  49     rp_dw_min = FLT_EPSILON; rp_dw_max = 50.;
  50 }
  51
  52
  53 CvANN_MLP_TrainParams::CvANN_MLP_TrainParams( CvTermCriteria _term_crit,
  54                                               int _train_method,
  55                                               double _param1, double _param2 )
  56 {
  57     term_crit = _term_crit;
  58     train_method = _train_method;
  59     bp_dw_scale = bp_moment_scale = 0.1;
  60     rp_dw0 = 1.; rp_dw_plus = 1.2; rp_dw_minus = 0.5;
  61     rp_dw_min = FLT_EPSILON; rp_dw_max = 50.;
  62
  63     if( train_method == RPROP )
  64     {
  65         rp_dw0 = _param1;
  66         if( rp_dw0 < FLT_EPSILON )
  67             rp_dw0 = 1.;
  68         rp_dw_min = _param2;
  69         rp_dw_min = MAX( rp_dw_min, 0 );
  70     }
  71     else if( train_method == BACKPROP )
  72     {
  73         bp_dw_scale = _param1;
  74         if( bp_dw_scale <= 0 )
  75             bp_dw_scale = 0.1;
  76         bp_dw_scale = MAX( bp_dw_scale, 1e-3 );
  77         bp_dw_scale = MIN( bp_dw_scale, 1 );
  78         bp_moment_scale = _param2;
  79         if( bp_moment_scale < 0 )
  80             bp_moment_scale = 0.1;
  81         bp_moment_scale = MIN( bp_moment_scale, 1 );
  82     }
  83     else
  84         train_method = RPROP;
  85 }
  86
  87
  88 CvANN_MLP_TrainParams::~CvANN_MLP_TrainParams()
  89 {
  90 }
  91
  92
  93 CvANN_MLP::CvANN_MLP()
  94 {
  95     layer_sizes = wbuf = 0;
  96     min_val = max_val = min_val1 = max_val1 = 0.;
  97     weights = 0;
  98     rng = &cv::theRNG();
  99     default_model_name = "my_nn";
 100     clear();
 101 }
 102
 103
 104 CvANN_MLP::CvANN_MLP( const CvMat* _layer_sizes,
 105                       int _activ_func,
 106                       double _f_param1, double _f_param2 )
 107 {
 108     layer_sizes = wbuf = 0;
 109     min_val = max_val = min_val1 = max_val1 = 0.;
 110     weights = 0;
 111     rng = &cv::theRNG();
 112     default_model_name = "my_nn";
 113     create( _layer_sizes, _activ_func, _f_param1, _f_param2 );
 114 }
 115
 116
 117 CvANN_MLP::~CvANN_MLP()
 118 {
 119     clear();
 120 }
 121
 122
 123 void CvANN_MLP::clear()
 124 {
 125     cvReleaseMat( &layer_sizes );
 126     cvReleaseMat( &wbuf );
 127     cvFree( &weights );
 128     activ_func = SIGMOID_SYM;
 129     f_param1 = f_param2 = 1;
 130     max_buf_sz = 1 << 12;
 131 }
 132
 133
 134 void CvANN_MLP::set_activ_func( int _activ_func, double _f_param1, double _f_param2 )
 135 {
 136     CV_FUNCNAME( "CvANN_MLP::set_activ_func" );
 137
 138     __BEGIN__;
 139
 140     if( _activ_func < 0 || _activ_func > GAUSSIAN )
 141         CV_ERROR( CV_StsOutOfRange, "Unknown activation function" );
 142
 143     activ_func = _activ_func;
 144
 145     switch( activ_func )
 146     {
 147     case SIGMOID_SYM:
 148         max_val = 0.95; min_val = -max_val;
 149         max_val1 = 0.98; min_val1 = -max_val1;
 150         if( fabs(_f_param1) < FLT_EPSILON )
 151             _f_param1 = 2./3;
 152         if( fabs(_f_param2) < FLT_EPSILON )
 153             _f_param2 = 1.7159;
 154         break;
 155     case GAUSSIAN:
 156         max_val = 1.; min_val = 0.05;
 157         max_val1 = 1.; min_val1 = 0.02;
 158         if( fabs(_f_param1) < FLT_EPSILON )
 159             _f_param1 = 1.;
 160         if( fabs(_f_param2) < FLT_EPSILON )
 161             _f_param2 = 1.;
 162         break;
 163     default:
 164         min_val = max_val = min_val1 = max_val1 = 0.;
 165         _f_param1 = 1.;
 166         _f_param2 = 0.;
 167     }
 168
 169     f_param1 = _f_param1;
 170     f_param2 = _f_param2;
 171
 172     __END__;
 173 }
 174
 175
 176 void CvANN_MLP::init_weights()
 177 {
 178     int i, j, k;
 179
 180     for( i = 1; i < layer_sizes->cols; i++ )
 181     {
 182         int n1 = layer_sizes->data.i[i-1];
 183         int n2 = layer_sizes->data.i[i];
 184         double val = 0, G = n2 > 2 ? 0.7*pow((double)n1,1./(n2-1)) : 1.;
 185         double* w = weights[i];
 186
 187         // initialize weights using Nguyen-Widrow algorithm
 188         for( j = 0; j < n2; j++ )
 189         {
 190             double s = 0;
 191             for( k = 0; k <= n1; k++ )
 192             {
 193                 val = rng->uniform(0., 1.)*2-1.;
 194                 w[k*n2 + j] = val;
 195                 s += fabs(val);
 196             }
 197
 198             if( i < layer_sizes->cols - 1 )
 199             {
 200                 s = 1./(s - fabs(val));
 201                 for( k = 0; k <= n1; k++ )
 202                     w[k*n2 + j] *= s;
 203                 w[n1*n2 + j] *= G*(-1+j*2./n2);
 204             }
 205         }
 206     }
 207 }
 208
 209
 210 void CvANN_MLP::create( const CvMat* _layer_sizes, int _activ_func,
 211                         double _f_param1, double _f_param2 )
 212 {
 213     CV_FUNCNAME( "CvANN_MLP::create" );
 214
 215     __BEGIN__;
 216
 217     int i, l_step, l_count, buf_sz = 0;
 218     int *l_src, *l_dst;
 219
 220     clear();
 221
 222     if( !CV_IS_MAT(_layer_sizes) ||
 223         (_layer_sizes->cols != 1 && _layer_sizes->rows != 1) ||
 224         CV_MAT_TYPE(_layer_sizes->type) != CV_32SC1 )
 225         CV_ERROR( CV_StsBadArg,
 226         "The array of layer neuron counters must be an integer vector" );
 227
 228     CV_CALL( set_activ_func( _activ_func, _f_param1, _f_param2 ));
 229
 230     l_count = _layer_sizes->rows + _layer_sizes->cols - 1;
 231     l_src = _layer_sizes->data.i;
 232     l_step = CV_IS_MAT_CONT(_layer_sizes->type) ? 1 :
 233                 _layer_sizes->step / sizeof(l_src[0]);
 234     CV_CALL( layer_sizes = cvCreateMat( 1, l_count, CV_32SC1 ));
 235     l_dst = layer_sizes->data.i;
 236
 237     max_count = 0;
 238     for( i = 0; i < l_count; i++ )
 239     {
 240         int n = l_src[i*l_step];
 241         if( n < 1 + (0 < i && i < l_count-1))
 242             CV_ERROR( CV_StsOutOfRange,
 243             "there should be at least one input and one output "
 244             "and every hidden layer must have more than 1 neuron" );
 245         l_dst[i] = n;
 246         max_count = MAX( max_count, n );
 247         if( i > 0 )
 248             buf_sz += (l_dst[i-1]+1)*n;
 249     }
 250
 251     buf_sz += (l_dst[0] + l_dst[l_count-1]*2)*2;
 252
 253     CV_CALL( wbuf = cvCreateMat( 1, buf_sz, CV_64F ));
 254     CV_CALL( weights = (double**)cvAlloc( (l_count+2)*sizeof(weights[0]) ));
 255
 256     weights[0] = wbuf->data.db;
 257     weights[1] = weights[0] + l_dst[0]*2;
 258     for( i = 1; i < l_count; i++ )
 259         weights[i+1] = weights[i] + (l_dst[i-1] + 1)*l_dst[i];
 260     weights[l_count+1] = weights[l_count] + l_dst[l_count-1]*2;
 261
 262     __END__;
 263 }
 264
 265
 266 float CvANN_MLP::predict( const CvMat* _inputs, CvMat* _outputs ) const
 267 {
 268     int i, j, n, dn = 0, l_count, dn0, buf_sz, min_buf_sz;
 269
 270     if( !layer_sizes )
 271         CV_Error( CV_StsError, "The network has not been initialized" );
 272
 273     if( !CV_IS_MAT(_inputs) || !CV_IS_MAT(_outputs) ||
 274         !CV_ARE_TYPES_EQ(_inputs,_outputs) ||
 275         (CV_MAT_TYPE(_inputs->type) != CV_32FC1 &&
 276         CV_MAT_TYPE(_inputs->type) != CV_64FC1) ||
 277         _inputs->rows != _outputs->rows )
 278         CV_Error( CV_StsBadArg, "Both input and output must be floating-point matrices "
 279                                 "of the same type and have the same number of rows" );
 280
 281     if( _inputs->cols != layer_sizes->data.i[0] )
 282         CV_Error( CV_StsBadSize, "input matrix must have the same number of columns as "
 283                                  "the number of neurons in the input layer" );
 284
 285     if( _outputs->cols != layer_sizes->data.i[layer_sizes->cols - 1] )
 286         CV_Error( CV_StsBadSize, "output matrix must have the same number of columns as "
 287                                  "the number of neurons in the output layer" );
 288     n = dn0 = _inputs->rows;
 289     min_buf_sz = 2*max_count;
 290     buf_sz = n*min_buf_sz;
 291
 292     if( buf_sz > max_buf_sz )
 293     {
 294         dn0 = max_buf_sz/min_buf_sz;
 295         dn0 = MAX( dn0, 1 );
 296         buf_sz = dn0*min_buf_sz;
 297     }
 298
 299     cv::AutoBuffer<double> buf(buf_sz);
 300     l_count = layer_sizes->cols;
 301
 302     for( i = 0; i < n; i += dn )
 303     {
 304         CvMat hdr[2], _w, *layer_in = &hdr[0], *layer_out = &hdr[1], *temp;
 305         dn = MIN( dn0, n - i );
 306
 307         cvGetRows( _inputs, layer_in, i, i + dn );
 308         cvInitMatHeader( layer_out, dn, layer_in->cols, CV_64F, &buf[0] );
 309
 310         scale_input( layer_in, layer_out );
 311         CV_SWAP( layer_in, layer_out, temp );
 312
 313         for( j = 1; j < l_count; j++ )
 314         {
 315             double* data = buf + (j&1 ? max_count*dn0 : 0);
 316             int cols = layer_sizes->data.i[j];
 317
 318             cvInitMatHeader( layer_out, dn, cols, CV_64F, data );
 319             cvInitMatHeader( &_w, layer_in->cols, layer_out->cols, CV_64F, weights[j] );
 320             cvGEMM( layer_in, &_w, 1, 0, 0, layer_out );
 321             calc_activ_func( layer_out, _w.data.db + _w.rows*_w.cols );
 322
 323             CV_SWAP( layer_in, layer_out, temp );
 324         }
 325
 326         cvGetRows( _outputs, layer_out, i, i + dn );
 327         scale_output( layer_in, layer_out );
 328     }
 329
 330     return 0.f;
 331 }
 332
 333
 334 void CvANN_MLP::scale_input( const CvMat* _src, CvMat* _dst ) const
 335 {
 336     int i, j, cols = _src->cols;
 337     double* dst = _dst->data.db;
 338     const double* w = weights[0];
 339     int step = _src->step;
 340
 341     if( CV_MAT_TYPE( _src->type ) == CV_32F )
 342     {
 343         const float* src = _src->data.fl;
 344         step /= sizeof(src[0]);
 345
 346         for( i = 0; i < _src->rows; i++, src += step, dst += cols )
 347             for( j = 0; j < cols; j++ )
 348                 dst[j] = src[j]*w[j*2] + w[j*2+1];
 349     }
 350     else
 351     {
 352         const double* src = _src->data.db;
 353         step /= sizeof(src[0]);
 354
 355         for( i = 0; i < _src->rows; i++, src += step, dst += cols )
 356             for( j = 0; j < cols; j++ )
 357                 dst[j] = src[j]*w[j*2] + w[j*2+1];
 358     }
 359 }
 360
 361
 362 void CvANN_MLP::scale_output( const CvMat* _src, CvMat* _dst ) const
 363 {
 364     int i, j, cols = _src->cols;
 365     const double* src = _src->data.db;
 366     const double* w = weights[layer_sizes->cols];
 367     int step = _dst->step;
 368
 369     if( CV_MAT_TYPE( _dst->type ) == CV_32F )
 370     {
 371         float* dst = _dst->data.fl;
 372         step /= sizeof(dst[0]);
 373
 374         for( i = 0; i < _src->rows; i++, src += cols, dst += step )
 375             for( j = 0; j < cols; j++ )
 376                 dst[j] = (float)(src[j]*w[j*2] + w[j*2+1]);
 377     }
 378     else
 379     {
 380         double* dst = _dst->data.db;
 381         step /= sizeof(dst[0]);
 382
 383         for( i = 0; i < _src->rows; i++, src += cols, dst += step )
 384             for( j = 0; j < cols; j++ )
 385                 dst[j] = src[j]*w[j*2] + w[j*2+1];
 386     }
 387 }
 388
 389
 390 void CvANN_MLP::calc_activ_func( CvMat* sums, const double* bias ) const
 391 {
 392     int i, j, n = sums->rows, cols = sums->cols;
 393     double* data = sums->data.db;
 394     double scale = 0, scale2 = f_param2;
 395
 396     switch( activ_func )
 397     {
 398     case IDENTITY:
 399         scale = 1.;
 400         break;
 401     case SIGMOID_SYM:
 402         scale = -f_param1;
 403         break;
 404     case GAUSSIAN:
 405         scale = -f_param1*f_param1;
 406         break;
 407     default:
 408         ;
 409     }
 410
 411     assert( CV_IS_MAT_CONT(sums->type) );
 412
 413     if( activ_func != GAUSSIAN )
 414     {
 415         for( i = 0; i < n; i++, data += cols )
 416             for( j = 0; j < cols; j++ )
 417                 data[j] = (data[j] + bias[j])*scale;
 418
 419         if( activ_func == IDENTITY )
 420             return;
 421     }
 422     else
 423     {
 424         for( i = 0; i < n; i++, data += cols )
 425             for( j = 0; j < cols; j++ )
 426             {
 427                 double t = data[j] + bias[j];
 428                 data[j] = t*t*scale;
 429             }
 430     }
 431
 432     cvExp( sums, sums );
 433
 434     n *= cols;
 435     data -= n;
 436
 437     switch( activ_func )
 438     {
 439     case SIGMOID_SYM:
 440         for( i = 0; i <= n - 4; i += 4 )
 441         {
 442             double x0 = 1.+data[i], x1 = 1.+data[i+1], x2 = 1.+data[i+2], x3 = 1.+data[i+3];
 443             double a = x0*x1, b = x2*x3, d = scale2/(a*b), t0, t1;
 444             a *= d; b *= d;
 445             t0 = (2 - x0)*b*x1; t1 = (2 - x1)*b*x0;
 446             data[i] = t0; data[i+1] = t1;
 447             t0 = (2 - x2)*a*x3; t1 = (2 - x3)*a*x2;
 448             data[i+2] = t0; data[i+3] = t1;
 449         }
 450
 451         for( ; i < n; i++ )
 452         {
 453             double t = scale2*(1. - data[i])/(1. + data[i]);
 454             data[i] = t;
 455         }
 456         break;
 457
 458     case GAUSSIAN:
 459         for( i = 0; i < n; i++ )
 460             data[i] = scale2*data[i];
 461         break;
 462
 463     default:
 464         ;
 465     }
 466 }
 467
 468
 469 void CvANN_MLP::calc_activ_func_deriv( CvMat* _xf, CvMat* _df,
 470                                        const double* bias ) const
 471 {
 472     int i, j, n = _xf->rows, cols = _xf->cols;
 473     double* xf = _xf->data.db;
 474     double* df = _df->data.db;
 475     double scale, scale2 = f_param2;
 476     assert( CV_IS_MAT_CONT( _xf->type & _df->type ) );
 477
 478     if( activ_func == IDENTITY )
 479     {
 480         for( i = 0; i < n; i++, xf += cols, df += cols )
 481             for( j = 0; j < cols; j++ )
 482             {
 483                 xf[j] += bias[j];
 484                 df[j] = 1;
 485             }
 486         return;
 487     }
 488     else if( activ_func == GAUSSIAN )
 489     {
 490         scale = -f_param1*f_param1;
 491         scale2 *= scale;
 492         for( i = 0; i < n; i++, xf += cols, df += cols )
 493             for( j = 0; j < cols; j++ )
 494             {
 495                 double t = xf[j] + bias[j];
 496                 df[j] = t*2*scale2;
 497                 xf[j] = t*t*scale;
 498             }
 499         cvExp( _xf, _xf );
 500
 501         n *= cols;
 502         xf -= n; df -= n;
 503
 504         for( i = 0; i < n; i++ )
 505             df[i] *= xf[i];
 506     }
 507     else
 508     {
 509         scale = f_param1;
 510         for( i = 0; i < n; i++, xf += cols, df += cols )
 511             for( j = 0; j < cols; j++ )
 512             {
 513                 xf[j] = (xf[j] + bias[j])*scale;
 514                 df[j] = -fabs(xf[j]);
 515             }
 516
 517         cvExp( _df, _df );
 518
 519         n *= cols;
 520         xf -= n; df -= n;
 521
 522         // ((1+exp(-ax))^-1)'=a*((1+exp(-ax))^-2)*exp(-ax);
 523         // ((1-exp(-ax))/(1+exp(-ax)))'=(a*exp(-ax)*(1+exp(-ax)) + a*exp(-ax)*(1-exp(-ax)))/(1+exp(-ax))^2=
 524         // 2*a*exp(-ax)/(1+exp(-ax))^2
 525         scale *= 2*f_param2;
 526         for( i = 0; i < n; i++ )
 527         {
 528             int s0 = xf[i] > 0 ? 1 : -1;
 529             double t0 = 1./(1. + df[i]);
 530             double t1 = scale*df[i]*t0*t0;
 531             t0 *= scale2*(1. - df[i])*s0;
 532             df[i] = t1;
 533             xf[i] = t0;
 534         }
 535     }
 536 }
 537
 538
 539 void CvANN_MLP::calc_input_scale( const CvVectors* vecs, int flags )
 540 {
 541     bool reset_weights = (flags & UPDATE_WEIGHTS) == 0;
 542     bool no_scale = (flags & NO_INPUT_SCALE) != 0;
 543     double* scale = weights[0];
 544     int count = vecs->count;
 545
 546     if( reset_weights )
 547     {
 548         int i, j, vcount = layer_sizes->data.i[0];
 549         int type = vecs->type;
 550         double a = no_scale ? 1. : 0.;
 551
 552         for( j = 0; j < vcount; j++ )
 553             scale[2*j] = a, scale[j*2+1] = 0.;
 554
 555         if( no_scale )
 556             return;
 557
 558         for( i = 0; i < count; i++ )
 559         {
 560             const float* f = vecs->data.fl[i];
 561             const double* d = vecs->data.db[i];
 562             for( j = 0; j < vcount; j++ )
 563             {
 564                 double t = type == CV_32F ? (double)f[j] : d[j];
 565                 scale[j*2] += t;
 566                 scale[j*2+1] += t*t;
 567             }
 568         }
 569
 570         for( j = 0; j < vcount; j++ )
 571         {
 572             double s = scale[j*2], s2 = scale[j*2+1];
 573             double m = s/count, sigma2 = s2/count - m*m;
 574             scale[j*2] = sigma2 < DBL_EPSILON ? 1 : 1./sqrt(sigma2);
 575             scale[j*2+1] = -m*scale[j*2];
 576         }
 577     }
 578 }
 579
 580
 581 void CvANN_MLP::calc_output_scale( const CvVectors* vecs, int flags )
 582 {
 583     int i, j, vcount = layer_sizes->data.i[layer_sizes->cols-1];
 584     int type = vecs->type;
 585     double m = min_val, M = max_val, m1 = min_val1, M1 = max_val1;
 586     bool reset_weights = (flags & UPDATE_WEIGHTS) == 0;
 587     bool no_scale = (flags & NO_OUTPUT_SCALE) != 0;
 588     int l_count = layer_sizes->cols;
 589     double* scale = weights[l_count];
 590     double* inv_scale = weights[l_count+1];
 591     int count = vecs->count;
 592
 593     CV_FUNCNAME( "CvANN_MLP::calc_output_scale" );
 594
 595     __BEGIN__;
 596
 597     if( reset_weights )
 598     {
 599         double a0 = no_scale ? 1 : DBL_MAX, b0 = no_scale ? 0 : -DBL_MAX;
 600
 601         for( j = 0; j < vcount; j++ )
 602         {
 603             scale[2*j] = inv_scale[2*j] = a0;
 604             scale[j*2+1] = inv_scale[2*j+1] = b0;
 605         }
 606
 607         if( no_scale )
 608             EXIT;
 609     }
 610
 611     for( i = 0; i < count; i++ )
 612     {
 613         const float* f = vecs->data.fl[i];
 614         const double* d = vecs->data.db[i];
 615
 616         for( j = 0; j < vcount; j++ )
 617         {
 618             double t = type == CV_32F ? (double)f[j] : d[j];
 619
 620             if( reset_weights )
 621             {
 622                 double mj = scale[j*2], Mj = scale[j*2+1];
 623                 if( mj > t ) mj = t;
 624                 if( Mj < t ) Mj = t;
 625
 626                 scale[j*2] = mj;
 627                 scale[j*2+1] = Mj;
 628             }
 629             else
 630             {
 631                 t = t*inv_scale[j*2] + inv_scale[2*j+1];
 632                 if( t < m1 || t > M1 )
 633                     CV_ERROR( CV_StsOutOfRange,
 634                     "Some of new output training vector components run exceed the original range too much" );
 635             }
 636         }
 637     }
 638
 639     if( reset_weights )
 640         for( j = 0; j < vcount; j++ )
 641         {
 642             // map mj..Mj to m..M
 643             double mj = scale[j*2], Mj = scale[j*2+1];
 644             double a, b;
 645             double delta = Mj - mj;
 646             if( delta < DBL_EPSILON )
 647                 a = 1, b = (M + m - Mj - mj)*0.5;
 648             else
 649                 a = (M - m)/delta, b = m - mj*a;
 650             inv_scale[j*2] = a; inv_scale[j*2+1] = b;
 651             a = 1./a; b = -b*a;
 652             scale[j*2] = a; scale[j*2+1] = b;
 653         }
 654
 655     __END__;
 656 }
 657
 658
 659 bool CvANN_MLP::prepare_to_train( const CvMat* _inputs, const CvMat* _outputs,
 660             const CvMat* _sample_weights, const CvMat* _sample_idx,
 661             CvVectors* _ivecs, CvVectors* _ovecs, double** _sw, int _flags )
 662 {
 663     bool ok = false;
 664     CvMat* sample_idx = 0;
 665     CvVectors ivecs, ovecs;
 666     double* sw = 0;
 667     int count = 0;
 668
 669     CV_FUNCNAME( "CvANN_MLP::prepare_to_train" );
 670
 671     ivecs.data.ptr = ovecs.data.ptr = 0;
 672     assert( _ivecs && _ovecs );
 673
 674     __BEGIN__;
 675
 676     const int* sidx = 0;
 677     int i, sw_type = 0, sw_count = 0;
 678     int sw_step = 0;
 679     double sw_sum = 0;
 680
 681     if( !layer_sizes )
 682         CV_ERROR( CV_StsError,
 683         "The network has not been created. Use method create or the appropriate constructor" );
 684
 685     if( !CV_IS_MAT(_inputs) || (CV_MAT_TYPE(_inputs->type) != CV_32FC1 &&
 686         CV_MAT_TYPE(_inputs->type) != CV_64FC1) || _inputs->cols != layer_sizes->data.i[0] )
 687         CV_ERROR( CV_StsBadArg,
 688         "input training data should be a floating-point matrix with"
 689         "the number of rows equal to the number of training samples and "
 690         "the number of columns equal to the size of 0-th (input) layer" );
 691
 692     if( !CV_IS_MAT(_outputs) || (CV_MAT_TYPE(_outputs->type) != CV_32FC1 &&
 693         CV_MAT_TYPE(_outputs->type) != CV_64FC1) ||
 694         _outputs->cols != layer_sizes->data.i[layer_sizes->cols - 1] )
 695         CV_ERROR( CV_StsBadArg,
 696         "output training data should be a floating-point matrix with"
 697         "the number of rows equal to the number of training samples and "
 698         "the number of columns equal to the size of last (output) layer" );
 699
 700     if( _inputs->rows != _outputs->rows )
 701         CV_ERROR( CV_StsUnmatchedSizes, "The numbers of input and output samples do not match" );
 702
 703     if( _sample_idx )
 704     {
 705         CV_CALL( sample_idx = cvPreprocessIndexArray( _sample_idx, _inputs->rows ));
 706         sidx = sample_idx->data.i;
 707         count = sample_idx->cols + sample_idx->rows - 1;
 708     }
 709     else
 710         count = _inputs->rows;
 711
 712     if( _sample_weights )
 713     {
 714         if( !CV_IS_MAT(_sample_weights) )
 715             CV_ERROR( CV_StsBadArg, "sample_weights (if passed) must be a valid matrix" );
 716
 717         sw_type = CV_MAT_TYPE(_sample_weights->type);
 718         sw_count = _sample_weights->cols + _sample_weights->rows - 1;
 719
 720         if( (sw_type != CV_32FC1 && sw_type != CV_64FC1) ||
 721             (_sample_weights->cols != 1 && _sample_weights->rows != 1) ||
 722             (sw_count != count && sw_count != _inputs->rows) )
 723             CV_ERROR( CV_StsBadArg,
 724             "sample_weights must be 1d floating-point vector containing weights "
 725             "of all or selected training samples" );
 726
 727         sw_step = CV_IS_MAT_CONT(_sample_weights->type) ? 1 :
 728             _sample_weights->step/CV_ELEM_SIZE(sw_type);
 729
 730         CV_CALL( sw = (double*)cvAlloc( count*sizeof(sw[0]) ));
 731     }
 732
 733     CV_CALL( ivecs.data.ptr = (uchar**)cvAlloc( count*sizeof(ivecs.data.ptr[0]) ));
 734     CV_CALL( ovecs.data.ptr = (uchar**)cvAlloc( count*sizeof(ovecs.data.ptr[0]) ));
 735
 736     ivecs.type = CV_MAT_TYPE(_inputs->type);
 737     ovecs.type = CV_MAT_TYPE(_outputs->type);
 738     ivecs.count = ovecs.count = count;
 739
 740     for( i = 0; i < count; i++ )
 741     {
 742         int idx = sidx ? sidx[i] : i;
 743         ivecs.data.ptr[i] = _inputs->data.ptr + idx*_inputs->step;
 744         ovecs.data.ptr[i] = _outputs->data.ptr + idx*_outputs->step;
 745         if( sw )
 746         {
 747             int si = sw_count == count ? i : idx;
 748             double w = sw_type == CV_32FC1 ?
 749                 (double)_sample_weights->data.fl[si*sw_step] :
 750                 _sample_weights->data.db[si*sw_step];
 751             sw[i] = w;
 752             if( w < 0 )
 753                 CV_ERROR( CV_StsOutOfRange, "some of sample weights are negative" );
 754             sw_sum += w;
 755         }
 756     }
 757
 758     // normalize weights
 759     if( sw )
 760     {
 761         sw_sum = sw_sum > DBL_EPSILON ? 1./sw_sum : 0;
 762         for( i = 0; i < count; i++ )
 763             sw[i] *= sw_sum;
 764     }
 765
 766     calc_input_scale( &ivecs, _flags );
 767     CV_CALL( calc_output_scale( &ovecs, _flags ));
 768
 769     ok = true;
 770
 771     __END__;
 772
 773     if( !ok )
 774     {
 775         cvFree( &ivecs.data.ptr );
 776         cvFree( &ovecs.data.ptr );
 777         cvFree( &sw );
 778     }
 779
 780     cvReleaseMat( &sample_idx );
 781     *_ivecs = ivecs;
 782     *_ovecs = ovecs;
 783     *_sw = sw;
 784
 785     return ok;
 786 }
 787
 788
 789 int CvANN_MLP::train( const CvMat* _inputs, const CvMat* _outputs,
 790                       const CvMat* _sample_weights, const CvMat* _sample_idx,
 791                       CvANN_MLP_TrainParams _params, int flags )
 792 {
 793     const int MAX_ITER = 1000;
 794     const double DEFAULT_EPSILON = FLT_EPSILON;
 795
 796     double* sw = 0;
 797     CvVectors x0, u;
 798     int iter = -1;
 799
 800     x0.data.ptr = u.data.ptr = 0;
 801
 802     CV_FUNCNAME( "CvANN_MLP::train" );
 803
 804     __BEGIN__;
 805
 806     int max_iter;
 807     double epsilon;
 808
 809     params = _params;
 810
 811     // initialize training data
 812     CV_CALL( prepare_to_train( _inputs, _outputs, _sample_weights,
 813                                _sample_idx, &x0, &u, &sw, flags ));
 814
 815     // ... and link weights
 816     if( !(flags & UPDATE_WEIGHTS) )
 817         init_weights();
 818
 819     max_iter = params.term_crit.type & CV_TERMCRIT_ITER ? params.term_crit.max_iter : MAX_ITER;
 820     max_iter = MAX( max_iter, 1 );
 821
 822     epsilon = params.term_crit.type & CV_TERMCRIT_EPS ? params.term_crit.epsilon : DEFAULT_EPSILON;
 823     epsilon = MAX(epsilon, DBL_EPSILON);
 824
 825     params.term_crit.type = CV_TERMCRIT_ITER + CV_TERMCRIT_EPS;
 826     params.term_crit.max_iter = max_iter;
 827     params.term_crit.epsilon = epsilon;
 828
 829     if( params.train_method == CvANN_MLP_TrainParams::BACKPROP )
 830     {
 831         CV_CALL( iter = train_backprop( x0, u, sw ));
 832     }
 833     else
 834     {
 835         CV_CALL( iter = train_rprop( x0, u, sw ));
 836     }
 837
 838     __END__;
 839
 840     cvFree( &x0.data.ptr );
 841     cvFree( &u.data.ptr );
 842     cvFree( &sw );
 843
 844     return iter;
 845 }
 846
 847
 848 int CvANN_MLP::train_backprop( CvVectors x0, CvVectors u, const double* sw )
 849 {
 850     CvMat* dw = 0;
 851     CvMat* buf = 0;
 852     double **x = 0, **df = 0;
 853     CvMat* _idx = 0;
 854     int iter = -1, count = x0.count;
 855
 856     CV_FUNCNAME( "CvANN_MLP::train_backprop" );
 857
 858     __BEGIN__;
 859
 860     int i, j, k, ivcount, ovcount, l_count, total = 0, max_iter;
 861     double *buf_ptr;
 862     double prev_E = DBL_MAX*0.5, E = 0, epsilon;
 863
 864     max_iter = params.term_crit.max_iter*count;
 865     epsilon = params.term_crit.epsilon*count;
 866
 867     l_count = layer_sizes->cols;
 868     ivcount = layer_sizes->data.i[0];
 869     ovcount = layer_sizes->data.i[l_count-1];
 870
 871     // allocate buffers
 872     for( i = 0; i < l_count; i++ )
 873         total += layer_sizes->data.i[i] + 1;
 874
 875     CV_CALL( dw = cvCreateMat( wbuf->rows, wbuf->cols, wbuf->type ));
 876     cvZero( dw );
 877     CV_CALL( buf = cvCreateMat( 1, (total + max_count)*2, CV_64F ));
 878     CV_CALL( _idx = cvCreateMat( 1, count, CV_32SC1 ));
 879     for( i = 0; i < count; i++ )
 880         _idx->data.i[i] = i;
 881
 882     CV_CALL( x = (double**)cvAlloc( total*2*sizeof(x[0]) ));
 883     df = x + total;
 884     buf_ptr = buf->data.db;
 885
 886     for( j = 0; j < l_count; j++ )
 887     {
 888         x[j] = buf_ptr;
 889         df[j] = x[j] + layer_sizes->data.i[j];
 890         buf_ptr += (df[j] - x[j])*2;
 891     }
 892
 893     // run back-propagation loop
 894     /*
 895         y_i = w_i*x_{i-1}
 896         x_i = f(y_i)
 897         E = 1/2*||u - x_N||^2
 898         grad_N = (x_N - u)*f'(y_i)
 899         dw_i(t) = momentum*dw_i(t-1) + dw_scale*x_{i-1}*grad_i
 900         w_i(t+1) = w_i(t) + dw_i(t)
 901         grad_{i-1} = w_i^t*grad_i
 902     */
 903     for( iter = 0; iter < max_iter; iter++ )
 904     {
 905         int idx = iter % count;
 906         double* w = weights[0];
 907         double sweight = sw ? count*sw[idx] : 1.;
 908         CvMat _w, _dw, hdr1, hdr2, ghdr1, ghdr2, _df;
 909         CvMat *x1 = &hdr1, *x2 = &hdr2, *grad1 = &ghdr1, *grad2 = &ghdr2, *temp;
 910
 911         if( idx == 0 )
 912         {
 913             //printf("%d. E = %g\n", iter/count, E);
 914             if( fabs(prev_E - E) < epsilon )
 915                 break;
 916             prev_E = E;
 917             E = 0;
 918
 919             // shuffle indices
 920             for( i = 0; i < count; i++ )
 921             {
 922                 int tt;
 923                 j = (*rng)(count);
 924                 k = (*rng)(count);
 925                 CV_SWAP( _idx->data.i[j], _idx->data.i[k], tt );
 926             }
 927         }
 928
 929         idx = _idx->data.i[idx];
 930
 931         if( x0.type == CV_32F )
 932         {
 933             const float* x0data = x0.data.fl[idx];
 934             for( j = 0; j < ivcount; j++ )
 935                 x[0][j] = x0data[j]*w[j*2] + w[j*2 + 1];
 936         }
 937         else
 938         {
 939             const double* x0data = x0.data.db[idx];
 940             for( j = 0; j < ivcount; j++ )
 941                 x[0][j] = x0data[j]*w[j*2] + w[j*2 + 1];
 942         }
 943
 944         cvInitMatHeader( x1, 1, ivcount, CV_64F, x[0] );
 945
 946         // forward pass, compute y[i]=w*x[i-1], x[i]=f(y[i]), df[i]=f'(y[i])
 947         for( i = 1; i < l_count; i++ )
 948         {
 949             cvInitMatHeader( x2, 1, layer_sizes->data.i[i], CV_64F, x[i] );
 950             cvInitMatHeader( &_w, x1->cols, x2->cols, CV_64F, weights[i] );
 951             cvGEMM( x1, &_w, 1, 0, 0, x2 );
 952             _df = *x2;
 953             _df.data.db = df[i];
 954             calc_activ_func_deriv( x2, &_df, _w.data.db + _w.rows*_w.cols );
 955             CV_SWAP( x1, x2, temp );
 956         }
 957
 958         cvInitMatHeader( grad1, 1, ovcount, CV_64F, buf_ptr );
 959         *grad2 = *grad1;
 960         grad2->data.db = buf_ptr + max_count;
 961
 962         w = weights[l_count+1];
 963
 964         // calculate error
 965         if( u.type == CV_32F )
 966         {
 967             const float* udata = u.data.fl[idx];
 968             for( k = 0; k < ovcount; k++ )
 969             {
 970                 double t = udata[k]*w[k*2] + w[k*2+1] - x[l_count-1][k];
 971                 grad1->data.db[k] = t*sweight;
 972                 E += t*t;
 973             }
 974         }
 975         else
 976         {
 977             const double* udata = u.data.db[idx];
 978             for( k = 0; k < ovcount; k++ )
 979             {
 980                 double t = udata[k]*w[k*2] + w[k*2+1] - x[l_count-1][k];
 981                 grad1->data.db[k] = t*sweight;
 982                 E += t*t;
 983             }
 984         }
 985         E *= sweight;
 986
 987         // backward pass, update weights
 988         for( i = l_count-1; i > 0; i-- )
 989         {
 990             int n1 = layer_sizes->data.i[i-1], n2 = layer_sizes->data.i[i];
 991             cvInitMatHeader( &_df, 1, n2, CV_64F, df[i] );
 992             cvMul( grad1, &_df, grad1 );
 993             cvInitMatHeader( &_w, n1+1, n2, CV_64F, weights[i] );
 994             cvInitMatHeader( &_dw, n1+1, n2, CV_64F, dw->data.db + (weights[i] - weights[0]) );
 995             cvInitMatHeader( x1, n1+1, 1, CV_64F, x[i-1] );
 996             x[i-1][n1] = 1.;
 997             cvGEMM( x1, grad1, params.bp_dw_scale, &_dw, params.bp_moment_scale, &_dw );
 998             cvAdd( &_w, &_dw, &_w );
 999             if( i > 1 )
1000             {
1001                 grad2->cols = n1;
1002                 _w.rows = n1;
1003                 cvGEMM( grad1, &_w, 1, 0, 0, grad2, CV_GEMM_B_T );
1004             }
1005             CV_SWAP( grad1, grad2, temp );
1006         }
1007     }
1008
1009     iter /= count;
1010
1011     __END__;
1012
1013     cvReleaseMat( &dw );
1014     cvReleaseMat( &buf );
1015     cvReleaseMat( &_idx );
1016     cvFree( &x );
1017
1018     return iter;
1019 }
1020
1021 struct rprop_loop : cv::ParallelLoopBody {
1022   rprop_loop(const CvANN_MLP* _point, double**& _weights, int& _count, int& _ivcount, CvVectors* _x0,
1023      int& _l_count, CvMat*& _layer_sizes, int& _ovcount, int& _max_count,
1024      CvVectors* _u, const double*& _sw, double& _inv_count, CvMat*& _dEdw, int& _dcount0, double* _E, int _buf_sz)
1025   {
1026     point = _point;
1027     weights = _weights;
1028     count = _count;
1029     ivcount = _ivcount;
1030     x0 = _x0;
1031     l_count = _l_count;
1032     layer_sizes = _layer_sizes;
1033     ovcount = _ovcount;
1034     max_count = _max_count;
1035     u = _u;
1036     sw = _sw;
1037     inv_count = _inv_count;
1038     dEdw = _dEdw;
1039     dcount0 = _dcount0;
1040     E = _E;
1041     buf_sz = _buf_sz;
1042   }
1043
1044   const CvANN_MLP* point;
1045   double** weights;
1046   int count;
1047   int ivcount;
1048   CvVectors* x0;
1049   int l_count;
1050   CvMat* layer_sizes;
1051   int ovcount;
1052   int max_count;
1053   CvVectors* u;
1054   const double* sw;
1055   double inv_count;
1056   CvMat* dEdw;
1057   int dcount0;
1058   double* E;
1059   int buf_sz;
1060
1061
1062   void operator()( const cv::Range& range ) const
1063   {
1064     double* buf_ptr;
1065     double** x = 0;
1066     double **df = 0;
1067     int total = 0;
1068
1069     for(int i = 0; i < l_count; i++ )
1070         total += layer_sizes->data.i[i];
1071     CvMat* buf;
1072     buf = cvCreateMat( 1, buf_sz, CV_64F );
1073     x = (double**)cvAlloc( total*2*sizeof(x[0]) );
1074     df = x + total;
1075     buf_ptr = buf->data.db;
1076     for(int i = 0; i < l_count; i++ )
1077     {
1078         x[i] = buf_ptr;
1079         df[i] = x[i] + layer_sizes->data.i[i]*dcount0;
1080         buf_ptr += (df[i] - x[i])*2;
1081     }
1082
1083     for(int si = range.start; si < range.end; si++ )
1084     {
1085         if (si % dcount0 != 0) continue;
1086         int n1, n2, k;
1087         double* w;
1088         CvMat _w, _dEdw, hdr1, hdr2, ghdr1, ghdr2, _df;
1089         CvMat *x1, *x2, *grad1, *grad2, *temp;
1090         int dcount = 0;
1091
1092         dcount = MIN(count - si , dcount0 );
1093         w = weights[0];
1094         grad1 = &ghdr1; grad2 = &ghdr2;
1095         x1 = &hdr1; x2 = &hdr2;
1096
1097         // grab and preprocess input data
1098         if( x0->type == CV_32F )
1099     {
1100             for(int i = 0; i < dcount; i++ )
1101             {
1102                 const float* x0data = x0->data.fl[si+i];
1103                 double* xdata = x[0]+i*ivcount;
1104                 for(int j = 0; j < ivcount; j++ )
1105                     xdata[j] = x0data[j]*w[j*2] + w[j*2+1];
1106             }
1107     }
1108         else
1109             for(int i = 0; i < dcount; i++ )
1110             {
1111                 const double* x0data = x0->data.db[si+i];
1112                 double* xdata = x[0]+i*ivcount;
1113                 for(int j = 0; j < ivcount; j++ )
1114                     xdata[j] = x0data[j]*w[j*2] + w[j*2+1];
1115             }
1116         cvInitMatHeader( x1, dcount, ivcount, CV_64F, x[0] );
1117
1118         // forward pass, compute y[i]=w*x[i-1], x[i]=f(y[i]), df[i]=f'(y[i])
1119         for(int i = 1; i < l_count; i++ )
1120         {
1121             cvInitMatHeader( x2, dcount, layer_sizes->data.i[i], CV_64F, x[i] );
1122             cvInitMatHeader( &_w, x1->cols, x2->cols, CV_64F, weights[i] );
1123             cvGEMM( x1, &_w, 1, 0, 0, x2 );
1124             _df = *x2;
1125             _df.data.db = df[i];
1126             point->calc_activ_func_deriv( x2, &_df, _w.data.db + _w.rows*_w.cols );
1127             CV_SWAP( x1, x2, temp );
1128         }
1129         cvInitMatHeader( grad1, dcount, ovcount, CV_64F, buf_ptr );
1130
1131         w = weights[l_count+1];
1132         grad2->data.db = buf_ptr + max_count*dcount;
1133
1134         // calculate error
1135         if( u->type == CV_32F )
1136             for(int i = 0; i < dcount; i++ )
1137             {
1138                 const float* udata = u->data.fl[si+i];
1139                 const double* xdata = x[l_count-1] + i*ovcount;
1140                 double* gdata = grad1->data.db + i*ovcount;
1141                 double sweight = sw ? sw[si+i] : inv_count, E1 = 0;
1142
1143                 for(int j = 0; j < ovcount; j++ )
1144                 {
1145                     double t = udata[j]*w[j*2] + w[j*2+1] - xdata[j];
1146                     gdata[j] = t*sweight;
1147                     E1 += t*t;
1148                 }
1149                 *E += sweight*E1;
1150             }
1151         else
1152             for(int i = 0; i < dcount; i++ )
1153             {
1154                 const double* udata = u->data.db[si+i];
1155                 const double* xdata = x[l_count-1] + i*ovcount;
1156                 double* gdata = grad1->data.db + i*ovcount;
1157                 double sweight = sw ? sw[si+i] : inv_count, E1 = 0;
1158
1159                 for(int j = 0; j < ovcount; j++ )
1160                 {
1161                     double t = udata[j]*w[j*2] + w[j*2+1] - xdata[j];
1162                     gdata[j] = t*sweight;
1163                     E1 += t*t;
1164                 }
1165                 *E += sweight*E1;
1166             }
1167
1168         // backward pass, update dEdw
1169         static cv::Mutex mutex;
1170
1171         for(int i = l_count-1; i > 0; i-- )
1172         {
1173             n1 = layer_sizes->data.i[i-1]; n2 = layer_sizes->data.i[i];
1174             cvInitMatHeader( &_df, dcount, n2, CV_64F, df[i] );
1175             cvMul( grad1, &_df, grad1 );
1176
1177             {
1178                 cv::AutoLock lock(mutex);
1179                 cvInitMatHeader( &_dEdw, n1, n2, CV_64F, dEdw->data.db+(weights[i]-weights[0]) );
1180                 cvInitMatHeader( x1, dcount, n1, CV_64F, x[i-1] );
1181                 cvGEMM( x1, grad1, 1, &_dEdw, 1, &_dEdw, CV_GEMM_A_T );
1182
1183                 // update bias part of dEdw
1184                 for( k = 0; k < dcount; k++ )
1185                 {
1186                     double* dst = _dEdw.data.db + n1*n2;
1187                     const double* src = grad1->data.db + k*n2;
1188                     for(int j = 0; j < n2; j++ )
1189                         dst[j] += src[j];
1190                 }
1191
1192                 if (i > 1)
1193                     cvInitMatHeader( &_w, n1, n2, CV_64F, weights[i] );
1194            }
1195
1196            cvInitMatHeader( grad2, dcount, n1, CV_64F, grad2->data.db );
1197            if( i > 1 )
1198                cvGEMM( grad1, &_w, 1, 0, 0, grad2, CV_GEMM_B_T );
1199            CV_SWAP( grad1, grad2, temp );
1200         }
1201     }
1202     cvFree(&x);
1203     cvReleaseMat( &buf );
1204 }
1205
1206 };
1207
1208
1209 int CvANN_MLP::train_rprop( CvVectors x0, CvVectors u, const double* sw )
1210 {
1211     const int max_buf_size = 1 << 16;
1212     CvMat* dw = 0;
1213     CvMat* dEdw = 0;
1214     CvMat* prev_dEdw_sign = 0;
1215     CvMat* buf = 0;
1216     double **x = 0, **df = 0;
1217     int iter = -1, count = x0.count;
1218
1219     CV_FUNCNAME( "CvANN_MLP::train" );
1220
1221     __BEGIN__;
1222
1223     int i, ivcount, ovcount, l_count, total = 0, max_iter, buf_sz, dcount0;
1224     double *buf_ptr;
1225     double prev_E = DBL_MAX*0.5, epsilon;
1226     double dw_plus, dw_minus, dw_min, dw_max;
1227     double inv_count;
1228
1229     max_iter = params.term_crit.max_iter;
1230     epsilon = params.term_crit.epsilon;
1231     dw_plus = params.rp_dw_plus;
1232     dw_minus = params.rp_dw_minus;
1233     dw_min = params.rp_dw_min;
1234     dw_max = params.rp_dw_max;
1235
1236     l_count = layer_sizes->cols;
1237     ivcount = layer_sizes->data.i[0];
1238     ovcount = layer_sizes->data.i[l_count-1];
1239
1240     // allocate buffers
1241     for( i = 0; i < l_count; i++ )
1242         total += layer_sizes->data.i[i];
1243
1244     CV_CALL( dw = cvCreateMat( wbuf->rows, wbuf->cols, wbuf->type ));
1245     cvSet( dw, cvScalarAll(params.rp_dw0) );
1246     CV_CALL( dEdw = cvCreateMat( wbuf->rows, wbuf->cols, wbuf->type ));
1247     cvZero( dEdw );
1248     CV_CALL( prev_dEdw_sign = cvCreateMat( wbuf->rows, wbuf->cols, CV_8SC1 ));
1249     cvZero( prev_dEdw_sign );
1250
1251     inv_count = 1./count;
1252     dcount0 = max_buf_size/(2*total);
1253     dcount0 = MAX( dcount0, 1 );
1254     dcount0 = MIN( dcount0, count );
1255     buf_sz = dcount0*(total + max_count)*2;
1256
1257     CV_CALL( buf = cvCreateMat( 1, buf_sz, CV_64F ));
1258
1259     CV_CALL( x = (double**)cvAlloc( total*2*sizeof(x[0]) ));
1260     df = x + total;
1261     buf_ptr = buf->data.db;
1262
1263     for( i = 0; i < l_count; i++ )
1264     {
1265         x[i] = buf_ptr;
1266         df[i] = x[i] + layer_sizes->data.i[i]*dcount0;
1267         buf_ptr += (df[i] - x[i])*2;
1268     }
1269
1270     // run rprop loop
1271     /*
1272         y_i(t) = w_i(t)*x_{i-1}(t)
1273         x_i(t) = f(y_i(t))
1274         E = sum_over_all_samples(1/2*||u - x_N||^2)
1275         grad_N = (x_N - u)*f'(y_i)
1276
1277                       MIN(dw_i{jk}(t)*dw_plus, dw_max), if dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) > 0
1278         dw_i{jk}(t) = MAX(dw_i{jk}(t)*dw_minus, dw_min), if dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) < 0
1279                       dw_i{jk}(t-1) else
1280
1281         if (dE/dw_i{jk}(t)*dE/dw_i{jk}(t-1) < 0)
1282            dE/dw_i{jk}(t)<-0
1283         else
1284            w_i{jk}(t+1) = w_i{jk}(t) + dw_i{jk}(t)
1285         grad_{i-1}(t) = w_i^t(t)*grad_i(t)
1286     */
1287     for( iter = 0; iter < max_iter; iter++ )
1288     {
1289         int n1, n2, j, k;
1290         double E = 0;
1291
1292         // first, iterate through all the samples and compute dEdw
1293         cv::parallel_for_(cv::Range(0, count),
1294             rprop_loop(this, weights, count, ivcount, &x0, l_count, layer_sizes,
1295                        ovcount, max_count, &u, sw, inv_count, dEdw, dcount0, &E, buf_sz)
1296         );
1297
1298         // now update weights
1299         for( i = 1; i < l_count; i++ )
1300         {
1301             n1 = layer_sizes->data.i[i-1]; n2 = layer_sizes->data.i[i];
1302             for( k = 0; k <= n1; k++ )
1303             {
1304                 double* wk = weights[i]+k*n2;
1305                 size_t delta = wk - weights[0];
1306                 double* dwk = dw->data.db + delta;
1307                 double* dEdwk = dEdw->data.db + delta;
1308                 char* prevEk = (char*)(prev_dEdw_sign->data.ptr + delta);
1309
1310                 for( j = 0; j < n2; j++ )
1311                 {
1312                     double Eval = dEdwk[j];
1313                     double dval = dwk[j];
1314                     double wval = wk[j];
1315                     int s = CV_SIGN(Eval);
1316                     int ss = prevEk[j]*s;
1317                     if( ss > 0 )
1318                     {
1319                         dval *= dw_plus;
1320                         dval = MIN( dval, dw_max );
1321                         dwk[j] = dval;
1322                         wk[j] = wval + dval*s;
1323                     }
1324                     else if( ss < 0 )
1325                     {
1326                         dval *= dw_minus;
1327                         dval = MAX( dval, dw_min );
1328                         prevEk[j] = 0;
1329                         dwk[j] = dval;
1330                         wk[j] = wval + dval*s;
1331                     }
1332                     else
1333                     {
1334                         prevEk[j] = (char)s;
1335                         wk[j] = wval + dval*s;
1336                     }
1337                     dEdwk[j] = 0.;
1338                 }
1339             }
1340         }
1341
1342         //printf("%d. E = %g\n", iter, E);
1343         if( fabs(prev_E - E) < epsilon )
1344             break;
1345         prev_E = E;
1346         E = 0;
1347     }
1348
1349     __END__;
1350
1351     cvReleaseMat( &dw );
1352     cvReleaseMat( &dEdw );
1353     cvReleaseMat( &prev_dEdw_sign );
1354     cvReleaseMat( &buf );
1355     cvFree( &x );
1356
1357     return iter;
1358 }
1359
1360
1361 void CvANN_MLP::write_params( CvFileStorage* fs ) const
1362 {
1363     //CV_FUNCNAME( "CvANN_MLP::write_params" );
1364
1365     __BEGIN__;
1366
1367     const char* activ_func_name = activ_func == IDENTITY ? "IDENTITY" :
1368                             activ_func == SIGMOID_SYM ? "SIGMOID_SYM" :
1369                             activ_func == GAUSSIAN ? "GAUSSIAN" : 0;
1370
1371     if( activ_func_name )
1372         cvWriteString( fs, "activation_function", activ_func_name );
1373     else
1374         cvWriteInt( fs, "activation_function", activ_func );
1375
1376     if( activ_func != IDENTITY )
1377     {
1378         cvWriteReal( fs, "f_param1", f_param1 );
1379         cvWriteReal( fs, "f_param2", f_param2 );
1380     }
1381
1382     cvWriteReal( fs, "min_val", min_val );
1383     cvWriteReal( fs, "max_val", max_val );
1384     cvWriteReal( fs, "min_val1", min_val1 );
1385     cvWriteReal( fs, "max_val1", max_val1 );
1386
1387     cvStartWriteStruct( fs, "training_params", CV_NODE_MAP );
1388     if( params.train_method == CvANN_MLP_TrainParams::BACKPROP )
1389     {
1390         cvWriteString( fs, "train_method", "BACKPROP" );
1391         cvWriteReal( fs, "dw_scale", params.bp_dw_scale );
1392         cvWriteReal( fs, "moment_scale", params.bp_moment_scale );
1393     }
1394     else if( params.train_method == CvANN_MLP_TrainParams::RPROP )
1395     {
1396         cvWriteString( fs, "train_method", "RPROP" );
1397         cvWriteReal( fs, "dw0", params.rp_dw0 );
1398         cvWriteReal( fs, "dw_plus", params.rp_dw_plus );
1399         cvWriteReal( fs, "dw_minus", params.rp_dw_minus );
1400         cvWriteReal( fs, "dw_min", params.rp_dw_min );
1401         cvWriteReal( fs, "dw_max", params.rp_dw_max );
1402     }
1403
1404     cvStartWriteStruct( fs, "term_criteria", CV_NODE_MAP + CV_NODE_FLOW );
1405     if( params.term_crit.type & CV_TERMCRIT_EPS )
1406         cvWriteReal( fs, "epsilon", params.term_crit.epsilon );
1407     if( params.term_crit.type & CV_TERMCRIT_ITER )
1408         cvWriteInt( fs, "iterations", params.term_crit.max_iter );
1409     cvEndWriteStruct( fs );
1410
1411     cvEndWriteStruct( fs );
1412
1413     __END__;
1414 }
1415
1416
1417 void CvANN_MLP::write( CvFileStorage* fs, const char* name ) const
1418 {
1419     CV_FUNCNAME( "CvANN_MLP::write" );
1420
1421     __BEGIN__;
1422
1423     int i, l_count = layer_sizes->cols;
1424
1425     if( !layer_sizes )
1426         CV_ERROR( CV_StsError, "The network has not been initialized" );
1427
1428     cvStartWriteStruct( fs, name, CV_NODE_MAP, CV_TYPE_NAME_ML_ANN_MLP );
1429
1430     cvWrite( fs, "layer_sizes", layer_sizes );
1431
1432     write_params( fs );
1433
1434     cvStartWriteStruct( fs, "input_scale", CV_NODE_SEQ + CV_NODE_FLOW );
1435     cvWriteRawData( fs, weights[0], layer_sizes->data.i[0]*2, "d" );
1436     cvEndWriteStruct( fs );
1437
1438     cvStartWriteStruct( fs, "output_scale", CV_NODE_SEQ + CV_NODE_FLOW );
1439     cvWriteRawData( fs, weights[l_count], layer_sizes->data.i[l_count-1]*2, "d" );
1440     cvEndWriteStruct( fs );
1441
1442     cvStartWriteStruct( fs, "inv_output_scale", CV_NODE_SEQ + CV_NODE_FLOW );
1443     cvWriteRawData( fs, weights[l_count+1], layer_sizes->data.i[l_count-1]*2, "d" );
1444     cvEndWriteStruct( fs );
1445
1446     cvStartWriteStruct( fs, "weights", CV_NODE_SEQ );
1447     for( i = 1; i < l_count; i++ )
1448     {
1449         cvStartWriteStruct( fs, 0, CV_NODE_SEQ + CV_NODE_FLOW );
1450         cvWriteRawData( fs, weights[i], (layer_sizes->data.i[i-1]+1)*layer_sizes->data.i[i], "d" );
1451         cvEndWriteStruct( fs );
1452     }
1453
1454     cvEndWriteStruct( fs );
1455
1456     __END__;
1457 }
1458
1459
1460 void CvANN_MLP::read_params( CvFileStorage* fs, CvFileNode* node )
1461 {
1462     //CV_FUNCNAME( "CvANN_MLP::read_params" );
1463
1464     __BEGIN__;
1465
1466     const char* activ_func_name = cvReadStringByName( fs, node, "activation_function", 0 );
1467     CvFileNode* tparams_node;
1468
1469     if( activ_func_name )
1470         activ_func = strcmp( activ_func_name, "SIGMOID_SYM" ) == 0 ? SIGMOID_SYM :
1471                      strcmp( activ_func_name, "IDENTITY" ) == 0 ? IDENTITY :
1472                      strcmp( activ_func_name, "GAUSSIAN" ) == 0 ? GAUSSIAN : 0;
1473     else
1474         activ_func = cvReadIntByName( fs, node, "activation_function" );
1475
1476     f_param1 = cvReadRealByName( fs, node, "f_param1", 0 );
1477     f_param2 = cvReadRealByName( fs, node, "f_param2", 0 );
1478
1479     set_activ_func( activ_func, f_param1, f_param2 );
1480
1481     min_val = cvReadRealByName( fs, node, "min_val", 0. );
1482     max_val = cvReadRealByName( fs, node, "max_val", 1. );
1483     min_val1 = cvReadRealByName( fs, node, "min_val1", 0. );
1484     max_val1 = cvReadRealByName( fs, node, "max_val1", 1. );
1485
1486     tparams_node = cvGetFileNodeByName( fs, node, "training_params" );
1487     params = CvANN_MLP_TrainParams();
1488
1489     if( tparams_node )
1490     {
1491         const char* tmethod_name = cvReadStringByName( fs, tparams_node, "train_method", "" );
1492         CvFileNode* tcrit_node;
1493
1494         if( strcmp( tmethod_name, "BACKPROP" ) == 0 )
1495         {
1496             params.train_method = CvANN_MLP_TrainParams::BACKPROP;
1497             params.bp_dw_scale = cvReadRealByName( fs, tparams_node, "dw_scale", 0 );
1498             params.bp_moment_scale = cvReadRealByName( fs, tparams_node, "moment_scale", 0 );
1499         }
1500         else if( strcmp( tmethod_name, "RPROP" ) == 0 )
1501         {
1502             params.train_method = CvANN_MLP_TrainParams::RPROP;
1503             params.rp_dw0 = cvReadRealByName( fs, tparams_node, "dw0", 0 );
1504             params.rp_dw_plus = cvReadRealByName( fs, tparams_node, "dw_plus", 0 );
1505             params.rp_dw_minus = cvReadRealByName( fs, tparams_node, "dw_minus", 0 );
1506             params.rp_dw_min = cvReadRealByName( fs, tparams_node, "dw_min", 0 );
1507             params.rp_dw_max = cvReadRealByName( fs, tparams_node, "dw_max", 0 );
1508         }
1509
1510         tcrit_node = cvGetFileNodeByName( fs, tparams_node, "term_criteria" );
1511         if( tcrit_node )
1512         {
1513             params.term_crit.epsilon = cvReadRealByName( fs, tcrit_node, "epsilon", -1 );
1514             params.term_crit.max_iter = cvReadIntByName( fs, tcrit_node, "iterations", -1 );
1515             params.term_crit.type = (params.term_crit.epsilon >= 0 ? CV_TERMCRIT_EPS : 0) +
1516                                    (params.term_crit.max_iter >= 0 ? CV_TERMCRIT_ITER : 0);
1517         }
1518     }
1519
1520     __END__;
1521 }
1522
1523
1524 void CvANN_MLP::read( CvFileStorage* fs, CvFileNode* node )
1525 {
1526     CvMat* _layer_sizes = 0;
1527
1528     CV_FUNCNAME( "CvANN_MLP::read" );
1529
1530     __BEGIN__;
1531
1532     CvFileNode* w;
1533     CvSeqReader reader;
1534     int i, l_count;
1535
1536     _layer_sizes = (CvMat*)cvReadByName( fs, node, "layer_sizes" );
1537     CV_CALL( create( _layer_sizes, SIGMOID_SYM, 0, 0 ));
1538
1539     cvReleaseMat( &_layer_sizes );
1540     _layer_sizes = NULL;
1541
1542     l_count = layer_sizes->cols;
1543
1544     CV_CALL( read_params( fs, node ));
1545
1546     w = cvGetFileNodeByName( fs, node, "input_scale" );
1547     if( !w || CV_NODE_TYPE(w->tag) != CV_NODE_SEQ ||
1548         w->data.seq->total != layer_sizes->data.i[0]*2 )
1549         CV_ERROR( CV_StsParseError, "input_scale tag is not found or is invalid" );
1550
1551     CV_CALL( cvReadRawData( fs, w, weights[0], "d" ));
1552
1553     w = cvGetFileNodeByName( fs, node, "output_scale" );
1554     if( !w || CV_NODE_TYPE(w->tag) != CV_NODE_SEQ ||
1555         w->data.seq->total != layer_sizes->data.i[l_count-1]*2 )
1556         CV_ERROR( CV_StsParseError, "output_scale tag is not found or is invalid" );
1557
1558     CV_CALL( cvReadRawData( fs, w, weights[l_count], "d" ));
1559
1560     w = cvGetFileNodeByName( fs, node, "inv_output_scale" );
1561     if( !w || CV_NODE_TYPE(w->tag) != CV_NODE_SEQ ||
1562         w->data.seq->total != layer_sizes->data.i[l_count-1]*2 )
1563         CV_ERROR( CV_StsParseError, "inv_output_scale tag is not found or is invalid" );
1564
1565     CV_CALL( cvReadRawData( fs, w, weights[l_count+1], "d" ));
1566
1567     w = cvGetFileNodeByName( fs, node, "weights" );
1568     if( !w || CV_NODE_TYPE(w->tag) != CV_NODE_SEQ ||
1569         w->data.seq->total != l_count - 1 )
1570         CV_ERROR( CV_StsParseError, "weights tag is not found or is invalid" );
1571
1572     cvStartReadSeq( w->data.seq, &reader );
1573
1574     for( i = 1; i < l_count; i++ )
1575     {
1576         w = (CvFileNode*)reader.ptr;
1577         CV_CALL( cvReadRawData( fs, w, weights[i], "d" ));
1578         CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );
1579     }
1580
1581     __END__;
1582 }
1583
1584 using namespace cv;
1585
1586 CvANN_MLP::CvANN_MLP( const Mat& _layer_sizes, int _activ_func,
1587                       double _f_param1, double _f_param2 )
1588 {
1589     layer_sizes = wbuf = 0;
1590     min_val = max_val = min_val1 = max_val1 = 0.;
1591     weights = 0;
1592     rng = &cv::theRNG();
1593     default_model_name = "my_nn";
1594     create( _layer_sizes, _activ_func, _f_param1, _f_param2 );
1595 }
1596
1597 void CvANN_MLP::create( const Mat& _layer_sizes, int _activ_func,
1598                        double _f_param1, double _f_param2 )
1599 {
1600     CvMat cvlayer_sizes = _layer_sizes;
1601     create( &cvlayer_sizes, _activ_func, _f_param1, _f_param2 );
1602 }
1603
1604 int CvANN_MLP::train( const Mat& _inputs, const Mat& _outputs,
1605                      const Mat& _sample_weights, const Mat& _sample_idx,
1606                      CvANN_MLP_TrainParams _params, int flags )
1607 {
1608     CvMat inputs = _inputs, outputs = _outputs, sweights = _sample_weights, sidx = _sample_idx;
1609     return train(&inputs, &outputs, sweights.data.ptr ? &sweights : 0,
1610                  sidx.data.ptr ? &sidx : 0, _params, flags);
1611 }
1612
1613 float CvANN_MLP::predict( const Mat& _inputs, Mat& _outputs ) const
1614 {
1615     CV_Assert(layer_sizes != 0);
1616     _outputs.create(_inputs.rows, layer_sizes->data.i[layer_sizes->cols-1], _inputs.type());
1617     CvMat inputs = _inputs, outputs = _outputs;
1618
1619     return predict(&inputs, &outputs);
1620 }
1621
1622 /* End of file. */