text: /*M///////////////////////////////////////////////////////////////////////////////////////
  
    IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
  
    By downloading, copying, installing or using the software you agree to this license.
    If you do not agree to this license, do not download, install,
    copy or use the software.
  
  
                          Intel License Agreement
  
   Copyright (C) 2000, Intel Corporation, all rights reserved.
   Third party copyrights are property of their respective owners.
  
   Redistribution and use in source and binary forms, with or without modification,
   are permitted provided that the following conditions are met:
  
     * Redistribution's of source code must retain the above copyright notice,
       this list of conditions and the following disclaimer.
  
     * Redistribution's in binary form must reproduce the above copyright notice,
       this list of conditions and the following disclaimer in the documentation
       and/or other materials provided with the distribution.
  
     * The name of Intel Corporation may not be used to endorse or promote products
       derived from this software without specific prior written permission.
  
   This software is provided by the copyright holders and contributors "as is" and
   any express or implied warranties, including, but not limited to, the implied
   warranties of merchantability and fitness for a particular purpose are disclaimed.
   In no event shall the Intel Corporation or contributors be liable for any direct,
   indirect, incidental, special, exemplary, or consequential damages
   (including, but not limited to, procurement of substitute goods or services;
   loss of use, data, or profits; or business interruption) however caused
   and on any theory of liability, whether in contract, strict liability,
   or tort (including negligence or otherwise) arising in any way out of
   the use of this software, even if advised of the possibility of such damage.
  
  M*/
  
  #include "precomp.hpp"
  #include <kode.h>
  
  #define sprintf	kdSprintfKHR
  
  static const float ord_nan = FLT_MAX*0.5f;
  static const int min_block_size = 1 << 16;
  static const int block_size_delta = 1 << 10;
  
  #define CV_CMP_NUM_PTR(a,b) (*(a) < *(b))
  static CV_IMPLEMENT_QSORT_EX( icvSortIntPtr, int*, CV_CMP_NUM_PTR, int )
  
  #define CV_CMP_PAIRS(a,b) (*((a).i) < *((b).i))
  static CV_IMPLEMENT_QSORT_EX( icvSortPairs, CvPair16u32s, CV_CMP_PAIRS, int )
  
  ///
  
  void CvERTreeTrainData::set_data( const CvMat* _train_data, int _tflag,
      const CvMat* _responses, const CvMat* _var_idx, const CvMat* _sample_idx,
      const CvMat* _var_type, const CvMat* _missing_mask, const CvDTreeParams& _params,
      bool _shared, bool _add_labels, bool _update_data )
  {
      CvMat* sample_indices = 0;
      CvMat* var_type0 = 0;
      CvMat* tmp_map = 0;
      int** int_ptr = 0;
      CvPair16u32s* pair16u32s_ptr = 0;
      CvDTreeTrainData* data = 0;
      float *_fdst = 0;
      int *_idst = 0;
      unsigned short* udst = 0;
      int* idst = 0;
  
      CV_FUNCNAME( "CvERTreeTrainData::set_data" );
  
      __BEGIN__;
      
      int sample_all = 0, r_type, cv_n;
      int total_c_count = 0;
      int tree_block_size, temp_block_size, max_split_size, nv_size, cv_size = 0;
      int ds_step, dv_step, ms_step = 0, mv_step = 0; // {data|mask}{sample|var}_step
      int vi, i, size;
      char err[100];
      const int *sidx = 0, *vidx = 0;
      
      if ( _params.use_surrogates )
          CV_ERROR(CV_StsBadArg, "CvERTrees do not support surrogate splits");
          
      if( _update_data && data_root )
      {
          CV_ERROR(CV_StsBadArg, "CvERTrees do not support data update");
      }
  
      clear();
  
      var_all = 0;
      rng = &cv::theRNG();
  
      CV_CALL( set_params( _params ));
  
      // check parameter types and sizes
      CV_CALL( cvCheckTrainData( _train_data, _tflag, _missing_mask, &var_all, &sample_all ));
  
      train_data = _train_data;
      responses = _responses;
      missing_mask = _missing_mask;
  
      if( _tflag == CV_ROW_SAMPLE )
      {
          ds_step = _train_data->step/CV_ELEM_SIZE(_train_data->type);
          dv_step = 1;
          if( _missing_mask )
              ms_step = _missing_mask->step, mv_step = 1;
      }
      else
      {
          dv_step = _train_data->step/CV_ELEM_SIZE(_train_data->type);
          ds_step = 1;
          if( _missing_mask )
              mv_step = _missing_mask->step, ms_step = 1;
      }
      tflag = _tflag;
  
      sample_count = sample_all;
      var_count = var_all;
  
      if( _sample_idx )
      {
          CV_CALL( sample_indices = cvPreprocessIndexArray( _sample_idx, sample_all ));
          sidx = sample_indices->data.i;
          sample_count = sample_indices->rows + sample_indices->cols - 1;
      }
  
      if( _var_idx )
      {
          CV_CALL( var_idx = cvPreprocessIndexArray( _var_idx, var_all ));
          vidx = var_idx->data.i;
          var_count = var_idx->rows + var_idx->cols - 1;
      }
  
      if( !CV_IS_MAT(_responses) ||
          (CV_MAT_TYPE(_responses->type) != CV_32SC1 &&
           CV_MAT_TYPE(_responses->type) != CV_32FC1) ||
          (_responses->rows != 1 && _responses->cols != 1) ||
          _responses->rows + _responses->cols - 1 != sample_all )
          CV_ERROR( CV_StsBadArg, "The array of _responses must be an integer or "
                    "floating-point vector containing as many elements as "
                    "the total number of samples in the training data matrix" );
     
      is_buf_16u = false;
      if ( sample_count < 65536 )  
          is_buf_16u = true;                                
      
      r_type = CV_VAR_CATEGORICAL;
      if( _var_type )
          CV_CALL( var_type0 = cvPreprocessVarType( _var_type, var_idx, var_count, &r_type ));
  
      CV_CALL( var_type = cvCreateMat( 1, var_count+2, CV_32SC1 ));
      
      cat_var_count = 0;
      ord_var_count = -1;
  
      is_classifier = r_type == CV_VAR_CATEGORICAL;
  
      // step 0. calc the number of categorical vars
      for( vi = 0; vi < var_count; vi++ )
      {
          char vt = var_type0 ? var_type0->data.ptr[vi] : CV_VAR_ORDERED;
          var_type->data.i[vi] = vt == CV_VAR_CATEGORICAL ? cat_var_count++ : ord_var_count--;
      }
  
      ord_var_count = ~ord_var_count;
      cv_n = params.cv_folds;
      // set the two last elements of var_type array to be able
      // to locate responses and cross-validation labels using
      // the corresponding get_* functions.
      var_type->data.i[var_count] = cat_var_count;
      var_type->data.i[var_count+1] = cat_var_count+1;
  
      // in case of single ordered predictor we need dummy cv_labels
      // for safe split_node_data() operation
      have_labels = cv_n > 0 || (ord_var_count == 1 && cat_var_count == 0) || _add_labels;
  
      work_var_count = cat_var_count + (is_classifier ? 1 : 0) + (have_labels ? 1 : 0);
      buf_size = (work_var_count + 1)*sample_count;
      shared = _shared;
      buf_count = shared ? 2 : 1;
      
      if ( is_buf_16u )
      {
          CV_CALL( buf = cvCreateMat( buf_count, buf_size, CV_16UC1 ));
          CV_CALL( pair16u32s_ptr = (CvPair16u32s*)cvAlloc( sample_count*sizeof(pair16u32s_ptr[0]) ));
      }
      else
      {
          CV_CALL( buf = cvCreateMat( buf_count, buf_size, CV_32SC1 ));
          CV_CALL( int_ptr = (int**)cvAlloc( sample_count*sizeof(int_ptr[0]) ));
      }    
  
      size = is_classifier ? cat_var_count+1 : cat_var_count;
      size = !size ? 1 : size;
      CV_CALL( cat_count = cvCreateMat( 1, size, CV_32SC1 ));
      CV_CALL( cat_ofs = cvCreateMat( 1, size, CV_32SC1 ));
      
      size = is_classifier ? (cat_var_count + 1)*params.max_categories : cat_var_count*params.max_categories;
      size = !size ? 1 : size;
      CV_CALL( cat_map = cvCreateMat( 1, size, CV_32SC1 ));
  
      // now calculate the maximum size of split,
      // create memory storage that will keep nodes and splits of the decision tree
      // allocate root node and the buffer for the whole training data
      max_split_size = cvAlign(sizeof(CvDTreeSplit) +
          (MAX(0,sample_count - 33)/32)*sizeof(int),sizeof(void*));
      tree_block_size = MAX((int)sizeof(CvDTreeNode)*8, max_split_size);
      tree_block_size = MAX(tree_block_size + block_size_delta, min_block_size);
      CV_CALL( tree_storage = cvCreateMemStorage( tree_block_size ));
      CV_CALL( node_heap = cvCreateSet( 0, sizeof(*node_heap), sizeof(CvDTreeNode), tree_storage ));
  
      nv_size = var_count*sizeof(int);
      nv_size = cvAlign(MAX( nv_size, (int)sizeof(CvSetElem) ), sizeof(void*));
  
      temp_block_size = nv_size;
  
      if( cv_n )
      {
          if( sample_count < cv_n*MAX(params.min_sample_count,10) )
              CV_ERROR( CV_StsOutOfRange,
                  "The many folds in cross-validation for such a small dataset" );
  
          cv_size = cvAlign( cv_n*(sizeof(int) + sizeof(double)*2), sizeof(double) );
          temp_block_size = MAX(temp_block_size, cv_size);
      }
  
      temp_block_size = MAX( temp_block_size + block_size_delta, min_block_size );
      CV_CALL( temp_storage = cvCreateMemStorage( temp_block_size ));
      CV_CALL( nv_heap = cvCreateSet( 0, sizeof(*nv_heap), nv_size, temp_storage ));
      if( cv_size )
          CV_CALL( cv_heap = cvCreateSet( 0, sizeof(*cv_heap), cv_size, temp_storage ));
  
      CV_CALL( data_root = new_node( 0, sample_count, 0, 0 ));
  
      max_c_count = 1;
  
      _fdst = 0;
      _idst = 0;
      if (ord_var_count)
          _fdst = (float*)cvAlloc(sample_count*sizeof(_fdst[0]));
      if (is_buf_16u && (cat_var_count || is_classifier))
          _idst = (int*)cvAlloc(sample_count*sizeof(_idst[0]));
  
      // transform the training data to convenient representation
      for( vi = 0; vi <= var_count; vi++ )
      {
          int ci;
          const uchar* mask = 0;
          int m_step = 0, step;
          const int* idata = 0;
          const float* fdata = 0;
          int num_valid = 0;
  
          if( vi < var_count ) // analyze i-th input variable
          {
              int vi0 = vidx ? vidx[vi] : vi;
              ci = get_var_type(vi);
              step = ds_step; m_step = ms_step;
              if( CV_MAT_TYPE(_train_data->type) == CV_32SC1 )
                  idata = _train_data->data.i + vi0*dv_step;
              else
                  fdata = _train_data->data.fl + vi0*dv_step;
              if( _missing_mask )
                  mask = _missing_mask->data.ptr + vi0*mv_step;
          }
          else // analyze _responses
          {
              ci = cat_var_count;
              step = CV_IS_MAT_CONT(_responses->type) ?
                  1 : _responses->step / CV_ELEM_SIZE(_responses->type);
              if( CV_MAT_TYPE(_responses->type) == CV_32SC1 )
                  idata = _responses->data.i;
              else
                  fdata = _responses->data.fl;
          }
  
          if( (vi < var_count && ci>=0) ||
              (vi == var_count && is_classifier) ) // process categorical variable or response
          {
              int c_count, prev_label;
              int* c_map;
              
              if (is_buf_16u)
                  udst = (unsigned short*)(buf->data.s + ci*sample_count);
              else
                  idst = buf->data.i + ci*sample_count;
              
              // copy data
              for( i = 0; i < sample_count; i++ )
              {
                  int val = INT_MAX, si = sidx ? sidx[i] : i;
                  if( !mask || !mask[si*m_step] )
                  {
                      if( idata )
                          val = idata[si*step];
                      else
                      {
                          float t = fdata[si*step];
                          val = cvRound(t);
                          if( val != t )
                          {
                              sprintf( err, "%d-th value of %d-th (categorical) "
                                  "variable is not an integer", i, vi );
                              CV_ERROR( CV_StsBadArg, err );
                          }
                      }
  
                      if( val == INT_MAX )
                      {
                          sprintf( err, "%d-th value of %d-th (categorical) "
                              "variable is too large", i, vi );
                          CV_ERROR( CV_StsBadArg, err );
                      }
                      num_valid++;
                  }
                  if (is_buf_16u)
                  {
                      _idst[i] = val;
                      pair16u32s_ptr[i].u = udst + i;
                      pair16u32s_ptr[i].i = _idst + i;
                  }   
                  else
                  {
                      idst[i] = val;
                      int_ptr[i] = idst + i;
                  }
              }
  
              c_count = num_valid > 0;
  
              if (is_buf_16u)
              {
                  icvSortPairs( pair16u32s_ptr, sample_count, 0 );
                  // count the categories
                  for( i = 1; i < num_valid; i++ )
                      if (*pair16u32s_ptr[i].i != *pair16u32s_ptr[i-1].i)
                          c_count ++ ;
              }
              else
              {
                  icvSortIntPtr( int_ptr, sample_count, 0 );
                  // count the categories
                  for( i = 1; i < num_valid; i++ )
                      c_count += *int_ptr[i] != *int_ptr[i-1];
              }
  
              if( vi > 0 )
                  max_c_count = MAX( max_c_count, c_count );
              cat_count->data.i[ci] = c_count;
              cat_ofs->data.i[ci] = total_c_count;
  
              // resize cat_map, if need
              if( cat_map->cols < total_c_count + c_count )
              {
                  tmp_map = cat_map;
                  CV_CALL( cat_map = cvCreateMat( 1,
                      MAX(cat_map->cols*3/2,total_c_count+c_count), CV_32SC1 ));
                  for( i = 0; i < total_c_count; i++ )
                      cat_map->data.i[i] = tmp_map->data.i[i];
                  cvReleaseMat( &tmp_map );
              }
  
              c_map = cat_map->data.i + total_c_count;
              total_c_count += c_count;
  
              c_count = -1;
              if (is_buf_16u)
              {
                  // compact the class indices and build the map
                  prev_label = ~*pair16u32s_ptr[0].i;
                  for( i = 0; i < num_valid; i++ )
                  {
                      int cur_label = *pair16u32s_ptr[i].i;
                      if( cur_label != prev_label )
                          c_map[++c_count] = prev_label = cur_label;
                      *pair16u32s_ptr[i].u = (unsigned short)c_count;
                  }
                  // replace labels for missing values with 65535
                  for( ; i < sample_count; i++ )
                      *pair16u32s_ptr[i].u = 65535;
              }
              else
              {
                  // compact the class indices and build the map
                  prev_label = ~*int_ptr[0];
                  for( i = 0; i < num_valid; i++ )
                  {
                      int cur_label = *int_ptr[i];
                      if( cur_label != prev_label )
                          c_map[++c_count] = prev_label = cur_label;
                      *int_ptr[i] = c_count;
                  }
                  // replace labels for missing values with -1
                  for( ; i < sample_count; i++ )
                      *int_ptr[i] = -1;
              }           
          }
          else if( ci < 0 ) // process ordered variable
          {
              for( i = 0; i < sample_count; i++ )
              {
                  float val = ord_nan;
                  int si = sidx ? sidx[i] : i;
                  if( !mask || !mask[si*m_step] )
                  {
                      if( idata )
                          val = (float)idata[si*step];
                      else
                          val = fdata[si*step];
  
                      if( fabs(val) >= ord_nan )
                      {
                          sprintf( err, "%d-th value of %d-th (ordered) "
                              "variable (=%g) is too large", i, vi, val );
                          CV_ERROR( CV_StsBadArg, err );
                      }
                      num_valid++;
                  }
              }
          }
          if( vi < var_count )
              data_root->set_num_valid(vi, num_valid);
      }
  
      // set sample labels
      if (is_buf_16u)
          udst = (unsigned short*)(buf->data.s + get_work_var_count()*sample_count);
      else
          idst = buf->data.i + get_work_var_count()*sample_count;
  
      for (i = 0; i < sample_count; i++)
      {
          if (udst)
              udst[i] = sidx ? (unsigned short)sidx[i] : (unsigned short)i;
          else
              idst[i] = sidx ? sidx[i] : i;
      }
  
      if( cv_n )
      {
          unsigned short* udst = 0;
          int* idst = 0;
  
          if (is_buf_16u)
          {
              udst = (unsigned short*)(buf->data.s + (get_work_var_count()-1)*sample_count);
              for( i = vi = 0; i < sample_count; i++ )
              {
                  udst[i] = (unsigned short)vi++;
                  vi &= vi < cv_n ? -1 : 0;
              }
  
              for( i = 0; i < sample_count; i++ )
              {
                  int a = (*rng)(sample_count);
                  int b = (*rng)(sample_count);
                  unsigned short unsh = (unsigned short)vi;
                  CV_SWAP( udst[a], udst[b], unsh );
              }
          }
          else
          {
              idst = buf->data.i + (get_work_var_count()-1)*sample_count;
              for( i = vi = 0; i < sample_count; i++ )
              {
                  idst[i] = vi++;
                  vi &= vi < cv_n ? -1 : 0;
              }
  
              for( i = 0; i < sample_count; i++ )
              {
                  int a = (*rng)(sample_count);
                  int b = (*rng)(sample_count);
                  CV_SWAP( idst[a], idst[b], vi );
              }
          }
      }
  
      if ( cat_map ) 
          cat_map->cols = MAX( total_c_count, 1 );
  
      max_split_size = cvAlign(sizeof(CvDTreeSplit) +
          (MAX(0,max_c_count - 33)/32)*sizeof(int),sizeof(void*));
      CV_CALL( split_heap = cvCreateSet( 0, sizeof(*split_heap), max_split_size, tree_storage ));
  
      have_priors = is_classifier && params.priors;
      if( is_classifier )
      {
          int m = get_num_classes();
          double sum = 0;
          CV_CALL( priors = cvCreateMat( 1, m, CV_64F ));
          for( i = 0; i < m; i++ )
          {
              double val = have_priors ? params.priors[i] : 1.;
              if( val <= 0 )
                  CV_ERROR( CV_StsOutOfRange, "Every class weight should be positive" );
              priors->data.db[i] = val;
              sum += val;
          }
  
          // normalize weights
          if( have_priors )
              cvScale( priors, priors, 1./sum );
  
          CV_CALL( priors_mult = cvCloneMat( priors ));
          CV_CALL( counts = cvCreateMat( 1, m, CV_32SC1 ));
      }
  
      CV_CALL( direction = cvCreateMat( 1, sample_count, CV_8UC1 ));
      CV_CALL( split_buf = cvCreateMat( 1, sample_count, CV_32SC1 ));
  
      __END__;
  
      if( data )
          delete data;
  
      if (_fdst)
          cvFree( &_fdst );
      if (_idst)
          cvFree( &_idst );
      cvFree( &int_ptr );
      cvReleaseMat( &var_type0 );
      cvReleaseMat( &sample_indices );
      cvReleaseMat( &tmp_map );
  }
  
  void CvERTreeTrainData::get_ord_var_data( CvDTreeNode* n, int vi, float* ord_values_buf, int* missing_buf,
                                            const float** ord_values, const int** missing, int* sample_indices_buf )
  {
      int vidx = var_idx ? var_idx->data.i[vi] : vi;
      int node_sample_count = n->sample_count; 
      // may use missing_buf as buffer for sample indices!
      const int* sample_indices = get_sample_indices(n, sample_indices_buf ? sample_indices_buf : missing_buf);
  
      int td_step = train_data->step/CV_ELEM_SIZE(train_data->type);
      int m_step = missing_mask ? missing_mask->step/CV_ELEM_SIZE(missing_mask->type) : 1;
      if( tflag == CV_ROW_SAMPLE )
      {
          for( int i = 0; i < node_sample_count; i++ )
          {
              int idx = sample_indices[i];
              missing_buf[i] = missing_mask ? *(missing_mask->data.ptr + idx * m_step + vi) : 0;
              ord_values_buf[i] = *(train_data->data.fl + idx * td_step + vidx);
          }
      }
      else
          for( int i = 0; i < node_sample_count; i++ )
          {
              int idx = sample_indices[i];
              missing_buf[i] = missing_mask ? *(missing_mask->data.ptr + vi* m_step + idx) : 0;
              ord_values_buf[i] = *(train_data->data.fl + vidx* td_step + idx);
          }
      *ord_values = ord_values_buf;
      *missing = missing_buf;
  }
  
  
  const int* CvERTreeTrainData::get_sample_indices( CvDTreeNode* n, int* indices_buf )
  {
      return get_cat_var_data( n, var_count + (is_classifier ? 1 : 0) + (have_labels ? 1 : 0), indices_buf );
  }
  
  
  const int* CvERTreeTrainData::get_cv_labels( CvDTreeNode* n, int* labels_buf )
  {
      if (have_labels)
          return get_cat_var_data( n, var_count + (is_classifier ? 1 : 0), labels_buf );
      return 0;
  }
  
  
  const int* CvERTreeTrainData::get_cat_var_data( CvDTreeNode* n, int vi, int* cat_values_buf )
  {
      int ci = get_var_type( vi);
      const int* cat_values = 0;
      if( !is_buf_16u )
          cat_values = buf->data.i + n->buf_idx*buf->cols + ci*sample_count + n->offset;
      else {
          const unsigned short* short_values = (const unsigned short*)(buf->data.s + n->buf_idx*buf->cols + 
              ci*sample_count + n->offset);
          for( int i = 0; i < n->sample_count; i++ )
              cat_values_buf[i] = short_values[i];
          cat_values = cat_values_buf;
      }
      return cat_values;
  }
  
  void CvERTreeTrainData::get_vectors( const CvMat* _subsample_idx,
                                      float* values, uchar* missing,
                                      float* responses, bool get_class_idx )
  {
      CvMat* subsample_idx = 0;
      CvMat* subsample_co = 0;
  
      cv::AutoBuffer<uchar> inn_buf(sample_count*(sizeof(float) + sizeof(int)));
  
      CV_FUNCNAME( "CvERTreeTrainData::get_vectors" );
  
      __BEGIN__;
  
      int i, vi, total = sample_count, count = total, cur_ofs = 0;
      int* sidx = 0;
      int* co = 0;
  
      if( _subsample_idx )
      {
          CV_CALL( subsample_idx = cvPreprocessIndexArray( _subsample_idx, sample_count ));
          sidx = subsample_idx->data.i;
          CV_CALL( subsample_co = cvCreateMat( 1, sample_count*2, CV_32SC1 ));
          co = subsample_co->data.i;
          cvZero( subsample_co );
          count = subsample_idx->cols + subsample_idx->rows - 1;
          for( i = 0; i < count; i++ )
              co[sidx[i]*2]++;
          for( i = 0; i < total; i++ )
          {
              int count_i = co[i*2];
              if( count_i )
              {
                  co[i*2+1] = cur_ofs*var_count;
                  cur_ofs += count_i;
              }
          }
      }
  
      if( missing )
          memset( missing, 1, count*var_count );
  
      for( vi = 0; vi < var_count; vi++ )
      {
          int ci = get_var_type(vi);
          if( ci >= 0 ) // categorical
          {
              float* dst = values + vi;
              uchar* m = missing ? missing + vi : 0;
              int* lbls_buf = (int*)(uchar*)inn_buf;
              const int* src = get_cat_var_data(data_root, vi, lbls_buf);
  
              for( i = 0; i < count; i++, dst += var_count )
              {
                  int idx = sidx ? sidx[i] : i;
                  int val = src[idx];
                  *dst = (float)val;
                  if( m )
                  {
                      *m = (!is_buf_16u && val < 0) || (is_buf_16u && (val == 65535));
                      m += var_count;
                  }
              }
          }
          else // ordered
          {
              int* mis_buf = (int*)(uchar*)inn_buf;
              const float *dst = 0;
              const int* mis = 0;
              get_ord_var_data(data_root, vi, values + vi, mis_buf, &dst, &mis, 0);
              for (int si = 0; si < total; si++)
                  *(missing + vi + si) = mis[si] == 0 ? 0 : 1;
          }
      }
  
      // copy responses
      if( responses )
      {
          if( is_classifier )
          {
              int* lbls_buf = (int*)(uchar*)inn_buf;
              const int* src = get_class_labels(data_root, lbls_buf);
              for( i = 0; i < count; i++ )
              {
                  int idx = sidx ? sidx[i] : i;
                  int val = get_class_idx ? src[idx] :
                      cat_map->data.i[cat_ofs->data.i[cat_var_count]+src[idx]];
                  responses[i] = (float)val;
              }
          }
          else           
          {
              float* _values_buf = (float*)(uchar*)inn_buf;
              int* sample_idx_buf = (int*)(_values_buf + sample_count);
              const float* _values = get_ord_responses(data_root, _values_buf, sample_idx_buf);
              for( i = 0; i < count; i++ )
              {
                  int idx = sidx ? sidx[i] : i;
                  responses[i] = _values[idx];
              }
          }
      }
  
      __END__;
  
      cvReleaseMat( &subsample_idx );
      cvReleaseMat( &subsample_co );
  }
  
  CvDTreeNode* CvERTreeTrainData::subsample_data( const CvMat* _subsample_idx )
  {
      CvDTreeNode* root = 0;
      
      CV_FUNCNAME( "CvERTreeTrainData::subsample_data" );
  
      __BEGIN__;
  
      if( !data_root )
          CV_ERROR( CV_StsError, "No training data has been set" );
  
      if( !_subsample_idx )
      {
          // make a copy of the root node
          CvDTreeNode temp;
          int i;
          root = new_node( 0, 1, 0, 0 );
          temp = *root;
          *root = *data_root;
          root->num_valid = temp.num_valid;
          if( root->num_valid )
          {
              for( i = 0; i < var_count; i++ )
                  root->num_valid[i] = data_root->num_valid[i];
          }
          root->cv_Tn = temp.cv_Tn;
          root->cv_node_risk = temp.cv_node_risk;
          root->cv_node_error = temp.cv_node_error;
      }
      else
          CV_ERROR( CV_StsError, "_subsample_idx must be null for extra-trees" );
      __END__;
  
      return root;
  }
  
  double CvForestERTree::calc_node_dir( CvDTreeNode* node )
  {
      char* dir = (char*)data->direction->data.ptr;
      int i, n = node->sample_count, vi = node->split->var_idx;
      double L, R;
  
      assert( !node->split->inversed );
  
      if( data->get_var_type(vi) >= 0 ) // split on categorical var
      {
          cv::AutoBuffer<uchar> inn_buf(n*sizeof(int)*(!data->have_priors ? 1 : 2));
          int* labels_buf = (int*)(uchar*)inn_buf;
          const int* labels = data->get_cat_var_data( node, vi, labels_buf );
          const int* subset = node->split->subset;
          if( !data->have_priors )
          {
              int sum = 0, sum_abs = 0;
  
              for( i = 0; i < n; i++ )
              {
                  int idx = labels[i];
                  int d = ( ((idx >= 0)&&(!data->is_buf_16u)) || ((idx != 65535)&&(data->is_buf_16u)) ) ?
                      CV_DTREE_CAT_DIR(idx,subset) : 0;
                  sum += d; sum_abs += d & 1;
                  dir[i] = (char)d;
              }
  
              R = (sum_abs + sum) >> 1;
              L = (sum_abs - sum) >> 1;
          }
          else
          {
              const double* priors = data->priors_mult->data.db;
              double sum = 0, sum_abs = 0;
              int *responses_buf = labels_buf + n;
              const int* responses = data->get_class_labels(node, responses_buf);
  
              for( i = 0; i < n; i++ )
              {
                  int idx = labels[i];
                  double w = priors[responses[i]];
                  int d = idx >= 0 ? CV_DTREE_CAT_DIR(idx,subset) : 0;
                  sum += d*w; sum_abs += (d & 1)*w;
                  dir[i] = (char)d;
              }
  
              R = (sum_abs + sum) * 0.5;
              L = (sum_abs - sum) * 0.5;
          }
      }
      else // split on ordered var
      {
          float split_val = node->split->ord.c;
          cv::AutoBuffer<uchar> inn_buf(n*(sizeof(int)*(!data->have_priors ? 1 : 2) + sizeof(float)));
          float* val_buf = (float*)(uchar*)inn_buf;
          int* missing_buf = (int*)(val_buf + n);
          const float* val = 0;
          const int* missing = 0;
          data->get_ord_var_data( node, vi, val_buf, missing_buf, &val, &missing, 0 );
  
          if( !data->have_priors )
          {
              L = R = 0;
              for( i = 0; i < n; i++ )
              {
                  if ( missing[i] )
                      dir[i] = (char)0;
                  else
                  {
                      if ( val[i] < split_val)
                      {
                          dir[i] = (char)-1;
                          L++;
                      }
                      else
                      {
                          dir[i] = (char)1;
                          R++;
                      }
                  }
              }
          }
          else
          {
              const double* priors = data->priors_mult->data.db;
              int* responses_buf = missing_buf + n;
              const int* responses = data->get_class_labels(node, responses_buf);
              L = R = 0;
              for( i = 0; i < n; i++ )
              {
                  if ( missing[i] )
                      dir[i] = (char)0;
                  else
                  {
                      double w = priors[responses[i]];
                      if ( val[i] < split_val)
                      {
                          dir[i] = (char)-1;
                           L += w;
                      }
                      else
                      {
                          dir[i] = (char)1;
                          R += w;
                      }
                  }
              }
          }
      }
  
      node->maxlr = MAX( L, R );
      return node->split->quality/(L + R);
  }
  
  CvDTreeSplit* CvForestERTree::find_split_ord_class( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split,
                                                      uchar* _ext_buf )
  {
      const float epsilon = FLT_EPSILON*2;
      const float split_delta = (1 + FLT_EPSILON) * FLT_EPSILON;
  
      int n = node->sample_count, i;
      int m = data->get_num_classes();
  
      cv::AutoBuffer<uchar> inn_buf;
      if( !_ext_buf )
          inn_buf.allocate(n*(2*sizeof(int) + sizeof(float)));
      uchar* ext_buf = _ext_buf ? _ext_buf : (uchar*)inn_buf;
      float* values_buf = (float*)ext_buf;
      int* missing_buf = (int*)(values_buf + n);
      const float* values = 0;
      const int* missing = 0;
      data->get_ord_var_data( node, vi, values_buf, missing_buf, &values, &missing, 0 );
      int* responses_buf = missing_buf + n;
      const int* responses = data->get_class_labels( node, responses_buf );
  
      double lbest_val = 0, rbest_val = 0, best_val = init_quality, split_val = 0;
      const double* priors = data->have_priors ? data->priors_mult->data.db : 0;
      bool is_find_split = false;
      float pmin, pmax;
      int smpi = 0;
      while ( missing[smpi] && (smpi < n) )
          smpi++;
      assert(smpi < n);
  
      pmin = values[smpi];
      pmax = pmin;
      for (; smpi < n; smpi++)
      {
          float ptemp = values[smpi];
          int m = missing[smpi];
          if (m) continue;
          if ( ptemp < pmin)
              pmin = ptemp;
          if ( ptemp > pmax)
              pmax = ptemp;
      }
      float fdiff = pmax-pmin;
      if (fdiff > epsilon)
      {
          is_find_split = true;
          cv::RNG* rng = data->rng;
          split_val = pmin + rng->uniform(0.f, 1.f) * fdiff ;
          if (split_val - pmin <= FLT_EPSILON)
              split_val = pmin + split_delta;
          if (pmax - split_val <= FLT_EPSILON)
              split_val = pmax - split_delta;       
  
          // calculate Gini index
          if ( !priors )
          {
              cv::AutoBuffer<int> lrc(m*2);
              int *lc = lrc, *rc = lc + m;
              int L = 0, R = 0;
      
              // init arrays of class instance counters on both sides of the split
              for( i = 0; i < m; i++ )
              {
                  lc[i] = 0;
                  rc[i] = 0;
              }
              for( int si = 0; si < n; si++ )
              {
                  int r = responses[si];
                  float val = values[si];
                  int m = missing[si];
                  if (m) continue;
                  if ( val < split_val )
                  {
                      lc[r]++;
                      L++;
                  }
                  else
                  {
                      rc[r]++;
                      R++;
                  }
              }
              for (int i = 0; i < m; i++)
              {
                  lbest_val += lc[i]*lc[i];
                  rbest_val += rc[i]*rc[i];
              }
              best_val = (lbest_val*R + rbest_val*L) / ((double)(L*R));
          }
          else
          {
              cv::AutoBuffer<double> lrc(m*2);
              double *lc = lrc, *rc = lc + m;
              double L = 0, R = 0;
      
              // init arrays of class instance counters on both sides of the split
              for( i = 0; i < m; i++ )
              {
                  lc[i] = 0;
                  rc[i] = 0;
              }
              for( int si = 0; si < n; si++ )
              {
                  int r = responses[si];
                  float val = values[si];
                  int m = missing[si];
                  double p = priors[r];
                  if (m) continue;
                  if ( val < split_val )
                  {
                      lc[r] += p;
                      L += p;
                  }
                  else
                  {
                      rc[r] += p;
                      R += p;
                  }
              }
              for (int i = 0; i < m; i++)
              {
                  lbest_val += lc[i]*lc[i];
                  rbest_val += rc[i]*rc[i];
              }
              best_val = (lbest_val*R + rbest_val*L) / (L*R);
          }
          
      }
  
      CvDTreeSplit* split = 0;
      if( is_find_split )
      {
          split = _split ? _split : data->new_split_ord( 0, 0.0f, 0, 0, 0.0f );
          split->var_idx = vi;
          split->ord.c = (float)split_val;
          split->ord.split_point = -1;
          split->inversed = 0;
          split->quality = (float)best_val;
      }
      return split;
  }
  
  CvDTreeSplit* CvForestERTree::find_split_cat_class( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split,
                                                      uchar* _ext_buf )
  {
      int ci = data->get_var_type(vi);
      int n = node->sample_count;
      int cm = data->get_num_classes(); 
      int vm = data->cat_count->data.i[ci];
      double best_val = init_quality;
      CvDTreeSplit *split = 0;
  
      if ( vm > 1 )
      {
          cv::AutoBuffer<int> inn_buf;
          if( !_ext_buf )
              inn_buf.allocate(2*n);
          int* ext_buf = _ext_buf ? (int*)_ext_buf : (int*)inn_buf;
  
          const int* labels = data->get_cat_var_data( node, vi, ext_buf );
          const int* responses = data->get_class_labels( node, ext_buf + n );
      
          const double* priors = data->have_priors ? data->priors_mult->data.db : 0;       
  
          // create random class mask
          cv::AutoBuffer<int> valid_cidx(vm);
          for (int i = 0; i < vm; i++)
          {
              valid_cidx[i] = -1;
          }
          for (int si = 0; si < n; si++)
          {
              int c = labels[si];
              if ( ((c == 65535) && data->is_buf_16u) || ((c<0) && (!data->is_buf_16u)) )
                  continue;
              valid_cidx[c]++;
          }
  
          int valid_ccount = 0;
          for (int i = 0; i < vm; i++)
              if (valid_cidx[i] >= 0)
              {
                  valid_cidx[i] = valid_ccount;
                  valid_ccount++;
              }
          if (valid_ccount > 1)
          {
              CvRNG* rng = forest->get_rng();
              int l_cval_count = 1 + cvRandInt(rng) % (valid_ccount-1);
  
              CvMat* var_class_mask = cvCreateMat( 1, valid_ccount, CV_8UC1 );
              CvMat submask;
              memset(var_class_mask->data.ptr, 0, valid_ccount*CV_ELEM_SIZE(var_class_mask->type));
              cvGetCols( var_class_mask, &submask, 0, l_cval_count );
              cvSet( &submask, cvScalar(1) );
              for (int i = 0; i < valid_ccount; i++)
              {
                  uchar temp;
                  int i1 =  cvRandInt( rng ) % valid_ccount;
                  int i2 = cvRandInt( rng ) % valid_ccount;
                  CV_SWAP( var_class_mask->data.ptr[i1], var_class_mask->data.ptr[i2], temp );
              }
  
              split = _split ? _split : data->new_split_cat( 0, -1.0f );
              split->var_idx = vi;
              memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));
  
              // calculate Gini index
              double lbest_val = 0, rbest_val = 0;
              if( !priors )
              {
                  cv::AutoBuffer<int> lrc(cm*2);
                  int *lc = lrc, *rc = lc + cm;
                  int L = 0, R = 0;
                  // init arrays of class instance counters on both sides of the split
                  for(int i = 0; i < cm; i++ )
                  {
                      lc[i] = 0;
                      rc[i] = 0;
                  }
                  for( int si = 0; si < n; si++ )
                  {
                      int r = responses[si];
                      int var_class_idx = labels[si];
                      if ( ((var_class_idx == 65535) && data->is_buf_16u) || ((var_class_idx<0) && (!data->is_buf_16u)) )
                          continue;
                      int mask_class_idx = valid_cidx[var_class_idx];
                      if (var_class_mask->data.ptr[mask_class_idx])
                      {
                          lc[r]++;
                          L++;                 
                          split->subset[var_class_idx >> 5] |= 1 << (var_class_idx & 31);
                      }
                      else
                      {
                          rc[r]++;
                          R++;
                      }
                  }
                  for (int i = 0; i < cm; i++)
                  {
                      lbest_val += lc[i]*lc[i];
                      rbest_val += rc[i]*rc[i];
                  }                
                  best_val = (lbest_val*R + rbest_val*L) / ((double)(L*R));
              }
              else
              {
                  cv::AutoBuffer<int> lrc(cm*2);
                  int *lc = lrc, *rc = lc + cm;
                  double L = 0, R = 0;
                  // init arrays of class instance counters on both sides of the split
                  for(int i = 0; i < cm; i++ )
                  {
                      lc[i] = 0;
                      rc[i] = 0;
                  }
                  for( int si = 0; si < n; si++ )
                  {
                      int r = responses[si];
                      int var_class_idx = labels[si];
                      if ( ((var_class_idx == 65535) && data->is_buf_16u) || ((var_class_idx<0) && (!data->is_buf_16u)) )
                          continue;
                      double p = priors[si];
                      int mask_class_idx = valid_cidx[var_class_idx];
                      
                      if (var_class_mask->data.ptr[mask_class_idx])
                      {
                          lc[r]+=(int)p;
                          L+=p;                 
                          split->subset[var_class_idx >> 5] |= 1 << (var_class_idx & 31);
                      }
                      else
                      {
                          rc[r]+=(int)p;
                          R+=p;
                      }
                  }
                  for (int i = 0; i < cm; i++)
                  {
                      lbest_val += lc[i]*lc[i];
                      rbest_val += rc[i]*rc[i];
                  }
                  best_val = (lbest_val*R + rbest_val*L) / (L*R);
              }
              split->quality = (float)best_val;
  
              cvReleaseMat(&var_class_mask);
          }   
      }  
  
      return split;
  }
  
  CvDTreeSplit* CvForestERTree::find_split_ord_reg( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split,
                                                    uchar* _ext_buf )
  {
      const float epsilon = FLT_EPSILON*2;
      const float split_delta = (1 + FLT_EPSILON) * FLT_EPSILON;
      int n = node->sample_count;
      cv::AutoBuffer<uchar> inn_buf;
      if( !_ext_buf )
          inn_buf.allocate(n*(2*sizeof(int) + 2*sizeof(float)));
      uchar* ext_buf = _ext_buf ? _ext_buf : (uchar*)inn_buf;
      float* values_buf = (float*)ext_buf;
      int* missing_buf = (int*)(values_buf + n);
      const float* values = 0;
      const int* missing = 0;
      data->get_ord_var_data( node, vi, values_buf, missing_buf, &values, &missing, 0 );
      float* responses_buf =  (float*)(missing_buf + n);
      int* sample_indices_buf =  (int*)(responses_buf + n);
      const float* responses = data->get_ord_responses( node, responses_buf, sample_indices_buf );
  
      double best_val = init_quality, split_val = 0, lsum = 0, rsum = 0;
      int L = 0, R = 0;
  
      bool is_find_split = false;
      float pmin, pmax;
      int smpi = 0;
      while ( missing[smpi] && (smpi < n) )
          smpi++;
  
      assert(smpi < n);
  
      pmin = values[smpi];
      pmax = pmin;
      for (; smpi < n; smpi++)
      {
          float ptemp = values[smpi];
          int m = missing[smpi];
          if (m) continue;
          if ( ptemp < pmin)
              pmin = ptemp;
          if ( ptemp > pmax)
              pmax = ptemp;
      }
      float fdiff = pmax-pmin;
      if (fdiff > epsilon)
      {
          is_find_split = true;
          cv::RNG* rng = data->rng;
          split_val = pmin + rng->uniform(0.f, 1.f) * fdiff ;
          if (split_val - pmin <= FLT_EPSILON)
              split_val = pmin + split_delta;
          if (pmax - split_val <= FLT_EPSILON)
              split_val = pmax - split_delta;    
  
          for (int si = 0; si < n; si++)
          {
              float r = responses[si];
              float val = values[si];
              int m = missing[si];
              if (m) continue;
              if (val < split_val)
              {
                  lsum += r;
                  L++;
              }
              else
              {
                  rsum += r;
                  R++;            
              }
          }
          best_val = (lsum*lsum*R + rsum*rsum*L)/((double)L*R);
      }
  
      CvDTreeSplit* split = 0;
      if( is_find_split )
      {
          split = _split ? _split : data->new_split_ord( 0, 0.0f, 0, 0, 0.0f );
          split->var_idx = vi;
          split->ord.c = (float)split_val;
          split->ord.split_point = -1;
          split->inversed = 0;
          split->quality = (float)best_val;
      }
      return split;
  }
  
  CvDTreeSplit* CvForestERTree::find_split_cat_reg( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split,
                                                    uchar* _ext_buf )
  {
      int ci = data->get_var_type(vi);
      int n = node->sample_count;
      int vm = data->cat_count->data.i[ci];
      double best_val = init_quality;
      CvDTreeSplit *split = 0;
      float lsum = 0, rsum = 0;
  
      if ( vm > 1 )
      {
          int base_size =  vm*sizeof(int);
          cv::AutoBuffer<uchar> inn_buf(base_size);
          if( !_ext_buf )
              inn_buf.allocate(base_size + n*(2*sizeof(int) + sizeof(float)));
          uchar* base_buf = (uchar*)inn_buf;
          uchar* ext_buf = _ext_buf ? _ext_buf : base_buf + base_size;
          int* labels_buf = (int*)ext_buf;
          const int* labels = data->get_cat_var_data( node, vi, labels_buf );
          float* responses_buf =  (float*)(labels_buf + n);
          int* sample_indices_buf = (int*)(responses_buf + n);
          const float* responses = data->get_ord_responses( node, responses_buf, sample_indices_buf );
  
          // create random class mask
          int *valid_cidx = (int*)base_buf;
          for (int i = 0; i < vm; i++)
          {
              valid_cidx[i] = -1;
          }
          for (int si = 0; si < n; si++)
          {
              int c = labels[si];
              if ( ((c == 65535) && data->is_buf_16u) || ((c<0) && (!data->is_buf_16u)) )
                          continue;
              valid_cidx[c]++;
          }
  
          int valid_ccount = 0;
          for (int i = 0; i < vm; i++)
              if (valid_cidx[i] >= 0)
              {
                  valid_cidx[i] = valid_ccount;
                  valid_ccount++;
              }
          if (valid_ccount > 1)
          {
              CvRNG* rng = forest->get_rng();
              int l_cval_count = 1 + cvRandInt(rng) % (valid_ccount-1);
  
              CvMat* var_class_mask = cvCreateMat( 1, valid_ccount, CV_8UC1 );
              CvMat submask;
              memset(var_class_mask->data.ptr, 0, valid_ccount*CV_ELEM_SIZE(var_class_mask->type));
              cvGetCols( var_class_mask, &submask, 0, l_cval_count );
              cvSet( &submask, cvScalar(1) );
              for (int i = 0; i < valid_ccount; i++)
              {
                  uchar temp;
                  int i1 = cvRandInt( rng ) % valid_ccount;
                  int i2 = cvRandInt( rng ) % valid_ccount;
                  CV_SWAP( var_class_mask->data.ptr[i1], var_class_mask->data.ptr[i2], temp );
              }
  
              split = _split ? _split : data->new_split_cat( 0, -1.0f);
              split->var_idx = vi;
              memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));
  
              int L = 0, R = 0;
              for( int si = 0; si < n; si++ )
              {
                  float r = responses[si];
                  int var_class_idx = labels[si];
                  if ( ((var_class_idx == 65535) && data->is_buf_16u) || ((var_class_idx<0) && (!data->is_buf_16u)) )
                          continue;
                  int mask_class_idx = valid_cidx[var_class_idx];
                  if (var_class_mask->data.ptr[mask_class_idx])
                  {
                      lsum += r;
                      L++;                 
                      split->subset[var_class_idx >> 5] |= 1 << (var_class_idx & 31);
                  }
                  else
                  {
                      rsum += r;
                      R++;
                  }
              }
              best_val = (lsum*lsum*R + rsum*rsum*L)/((double)L*R);
  
              split->quality = (float)best_val;
  
              cvReleaseMat(&var_class_mask);
          }   
      }  
  
      return split;
  }
  
  void CvForestERTree::split_node_data( CvDTreeNode* node )
  {
      int vi, i, n = node->sample_count, nl, nr, scount = data->sample_count;
      char* dir = (char*)data->direction->data.ptr;
      CvDTreeNode *left = 0, *right = 0;
      int new_buf_idx = data->get_child_buf_idx( node );
      CvMat* buf = data->buf;
      cv::AutoBuffer<int> temp_buf(n);
  
      complete_node_dir(node);
  
      for( i = nl = nr = 0; i < n; i++ )
      {
          int d = dir[i];
          nr += d;
          nl += d^1;
      }
  
      bool split_input_data;
      node->left = left = data->new_node( node, nl, new_buf_idx, node->offset );
      node->right = right = data->new_node( node, nr, new_buf_idx, node->offset + nl );
  
      split_input_data = node->depth + 1 < data->params.max_depth &&
          (node->left->sample_count > data->params.min_sample_count ||
          node->right->sample_count > data->params.min_sample_count);
  
      cv::AutoBuffer<uchar> inn_buf(n*(sizeof(int)+sizeof(float)));
      // split ordered vars
      for( vi = 0; vi < data->var_count; vi++ )
      {
          int ci = data->get_var_type(vi);
          if (ci >= 0) continue;
          
          int n1 = node->get_num_valid(vi), nr1 = 0;
          float* values_buf = (float*)(uchar*)inn_buf;
          int* missing_buf = (int*)(values_buf + n);
          const float* values = 0;
          const int* missing = 0;
          data->get_ord_var_data( node, vi, values_buf, missing_buf, &values, &missing, 0 );
  
          for( i = 0; i < n; i++ )
              nr1 += (!missing[i] & dir[i]);
          left->set_num_valid(vi, n1 - nr1);
          right->set_num_valid(vi, nr1);                
      }
      // split categorical vars, responses and cv_labels using new_idx relocation table
      for( vi = 0; vi < data->get_work_var_count() + data->ord_var_count; vi++ )
      {
          int ci = data->get_var_type(vi);
          if (ci < 0) continue;
  
          int n1 = node->get_num_valid(vi), nr1 = 0;
          const int* src_lbls = data->get_cat_var_data(node, vi, (int*)(uchar*)inn_buf);
  
          for(i = 0; i < n; i++)
              temp_buf[i] = src_lbls[i];
  
          if (data->is_buf_16u)
          {
              unsigned short *ldst = (unsigned short *)(buf->data.s + left->buf_idx*buf->cols + 
                  ci*scount + left->offset);
              unsigned short *rdst = (unsigned short *)(buf->data.s + right->buf_idx*buf->cols + 
                  ci*scount + right->offset);
              
              for( i = 0; i < n; i++ )
              {
                  int d = dir[i];
                  int idx = temp_buf[i];
                  if (d)
                  {
                      *rdst = (unsigned short)idx;
                      rdst++;
                      nr1 += (idx != 65535);
                  }
                  else
                  {
                      *ldst = (unsigned short)idx;
                      ldst++;
                  }
              }
  
              if( vi < data->var_count )
              {
                  left->set_num_valid(vi, n1 - nr1);
                  right->set_num_valid(vi, nr1);
              }
          }
          else
          {
              int *ldst = buf->data.i + left->buf_idx*buf->cols + 
                  ci*scount + left->offset;
              int *rdst = buf->data.i + right->buf_idx*buf->cols + 
                  ci*scount + right->offset;
              
              for( i = 0; i < n; i++ )
              {
                  int d = dir[i];
                  int idx = temp_buf[i];
                  if (d)
                  {
                      *rdst = idx;
                      rdst++;
                      nr1 += (idx >= 0);
                  }
                  else
                  {
                      *ldst = idx;
                      ldst++;
                  }
                  
              }
  
              if( vi < data->var_count )
              {
                  left->set_num_valid(vi, n1 - nr1);
                  right->set_num_valid(vi, nr1);
              }
          }        
      }
  
      // split sample indices
      int *sample_idx_src_buf = (int*)(uchar*)inn_buf;
      const int* sample_idx_src = 0;
      if (split_input_data)
      {
          sample_idx_src = data->get_sample_indices(node, sample_idx_src_buf);
  
          for(i = 0; i < n; i++)
              temp_buf[i] = sample_idx_src[i];
  
          int pos = data->get_work_var_count();
         
          if (data->is_buf_16u)
          {
              unsigned short* ldst = (unsigned short*)(buf->data.s + left->buf_idx*buf->cols + 
                  pos*scount + left->offset);
              unsigned short* rdst = (unsigned short*)(buf->data.s + right->buf_idx*buf->cols + 
                  pos*scount + right->offset);
              
              for (i = 0; i < n; i++)
              {
                  int d = dir[i];
                  unsigned short idx = (unsigned short)temp_buf[i];
                  if (d)
                  {
                      *rdst = idx;
                      rdst++;
                  }
                  else
                  {
                      *ldst = idx;
                      ldst++;
                  }
              }
          }
          else
          {
              int* ldst = buf->data.i + left->buf_idx*buf->cols + 
                  pos*scount + left->offset;
              int* rdst = buf->data.i + right->buf_idx*buf->cols + 
                  pos*scount + right->offset;
              for (i = 0; i < n; i++)
              {
                  int d = dir[i];
                  int idx = temp_buf[i];
                  if (d)
                  {
                      *rdst = idx;
                      rdst++;
                  }
                  else
                  {
                      *ldst = idx;
                      ldst++;
                  }
              }
          }
      }
      
      // deallocate the parent node data that is not needed anymore
      data->free_node_data(node);    
  }
  
  CvERTrees::CvERTrees()
  {
  }
  
  CvERTrees::~CvERTrees()
  {
  }
  
  bool CvERTrees::train( const CvMat* _train_data, int _tflag,
                          const CvMat* _responses, const CvMat* _var_idx,
                          const CvMat* _sample_idx, const CvMat* _var_type,
                          const CvMat* _missing_mask, CvRTParams params )
  {
      bool result = false;
  
      CV_FUNCNAME("CvERTrees::train");
      __BEGIN__
      int var_count = 0;
  
      clear();
  
      CvDTreeParams tree_params( params.max_depth, params.min_sample_count,
          params.regression_accuracy, params.use_surrogates, params.max_categories,
          params.cv_folds, params.use_1se_rule, false, params.priors );
  
      data = new CvERTreeTrainData();
      CV_CALL(data->set_data( _train_data, _tflag, _responses, _var_idx,
          _sample_idx, _var_type, _missing_mask, tree_params, true));
  
      var_count = data->var_count;
      if( params.nactive_vars > var_count )
          params.nactive_vars = var_count;
      else if( params.nactive_vars == 0 )
          params.nactive_vars = (int)sqrt((double)var_count);
      else if( params.nactive_vars < 0 )
          CV_ERROR( CV_StsBadArg, "<nactive_vars> must be non-negative" );
  
      // Create mask of active variables at the tree nodes
      CV_CALL(active_var_mask = cvCreateMat( 1, var_count, CV_8UC1 ));
      if( params.calc_var_importance )
      {
          CV_CALL(var_importance  = cvCreateMat( 1, var_count, CV_32FC1 ));
          cvZero(var_importance);
      }
      { // initialize active variables mask
          CvMat submask1, submask2;
          CV_Assert( (active_var_mask->cols >= 1) && (params.nactive_vars > 0) && (params.nactive_vars <= active_var_mask->cols) );
          cvGetCols( active_var_mask, &submask1, 0, params.nactive_vars );
          cvSet( &submask1, cvScalar(1) );
          if( params.nactive_vars < active_var_mask->cols )
          {
              cvGetCols( active_var_mask, &submask2, params.nactive_vars, var_count );
              cvZero( &submask2 );
          }
      }
  
      CV_CALL(result = grow_forest( params.term_crit ));
  
      result = true;
  
      __END__
      return result;
      
  }
  
  bool CvERTrees::train( CvMLData* data, CvRTParams params)
  {
     bool result = false;
  
      CV_FUNCNAME( "CvERTrees::train" );
  
      __BEGIN__;
  
      CV_CALL( result = CvRTrees::train( data, params) );
  
      __END__;
  
      return result;
  }
  
  bool CvERTrees::grow_forest( const CvTermCriteria term_crit )
  {
      bool result = false;
  
      CvMat* sample_idx_for_tree      = 0;
  
      CV_FUNCNAME("CvERTrees::grow_forest");
      __BEGIN__;
  
      const int max_ntrees = term_crit.max_iter;
      const double max_oob_err = term_crit.epsilon;
  
      const int dims = data->var_count;
      float maximal_response = 0;
  
      CvMat* oob_sample_votes	   = 0;
      CvMat* oob_responses       = 0;
  
      float* oob_samples_perm_ptr= 0;
  
      float* samples_ptr     = 0;
      uchar* missing_ptr     = 0;
      float* true_resp_ptr   = 0;
      bool is_oob_or_vimportance = ((max_oob_err > 0) && (term_crit.type != CV_TERMCRIT_ITER)) || var_importance;
  
      // oob_predictions_sum[i] = sum of predicted values for the i-th sample
      // oob_num_of_predictions[i] = number of summands
      //                            (number of predictions for the i-th sample)
      // initialize these variable to avoid warning C4701
      CvMat oob_predictions_sum = cvMat( 1, 1, CV_32FC1 );
      CvMat oob_num_of_predictions = cvMat( 1, 1, CV_32FC1 );
       
      nsamples = data->sample_count;
      nclasses = data->get_num_classes();
  
      if ( is_oob_or_vimportance )
      {
          if( data->is_classifier )
          {
              CV_CALL(oob_sample_votes = cvCreateMat( nsamples, nclasses, CV_32SC1 ));
              cvZero(oob_sample_votes);
          }
          else
          {
              // oob_responses[0,i] = oob_predictions_sum[i]
              //    = sum of predicted values for the i-th sample
              // oob_responses[1,i] = oob_num_of_predictions[i]
              //    = number of summands (number of predictions for the i-th sample)
              CV_CALL(oob_responses = cvCreateMat( 2, nsamples, CV_32FC1 ));
              cvZero(oob_responses);
              cvGetRow( oob_responses, &oob_predictions_sum, 0 );
              cvGetRow( oob_responses, &oob_num_of_predictions, 1 );
          }
          
          CV_CALL(oob_samples_perm_ptr     = (float*)cvAlloc( sizeof(float)*nsamples*dims ));
          CV_CALL(samples_ptr              = (float*)cvAlloc( sizeof(float)*nsamples*dims ));
          CV_CALL(missing_ptr              = (uchar*)cvAlloc( sizeof(uchar)*nsamples*dims ));
          CV_CALL(true_resp_ptr            = (float*)cvAlloc( sizeof(float)*nsamples ));            
  
          CV_CALL(data->get_vectors( 0, samples_ptr, missing_ptr, true_resp_ptr ));
          {
              double minval, maxval;
              CvMat responses = cvMat(1, nsamples, CV_32FC1, true_resp_ptr);
              cvMinMaxLoc( &responses, &minval, &maxval );
              maximal_response = (float)MAX( MAX( fabs(minval), fabs(maxval) ), 0 );
          }
      }
     
      trees = (CvForestTree**)cvAlloc( sizeof(trees[0])*max_ntrees );
      memset( trees, 0, sizeof(trees[0])*max_ntrees );
  
      CV_CALL(sample_idx_for_tree = cvCreateMat( 1, nsamples, CV_32SC1 ));
  
      for (int i = 0; i < nsamples; i++)
          sample_idx_for_tree->data.i[i] = i;
      ntrees = 0;
      while( ntrees < max_ntrees )
      {
          int i, oob_samples_count = 0;
          double ncorrect_responses = 0; // used for estimation of variable importance
          CvForestTree* tree = 0;
  
          trees[ntrees] = new CvForestERTree();
          tree = (CvForestERTree*)trees[ntrees];
          CV_CALL(tree->train( data, 0, this ));
  
          if ( is_oob_or_vimportance )
          {
              CvMat sample, missing;
              // form array of OOB samples indices and get these samples
              sample   = cvMat( 1, dims, CV_32FC1, samples_ptr );
              missing  = cvMat( 1, dims, CV_8UC1,  missing_ptr );
  
              oob_error = 0;
              for( i = 0; i < nsamples; i++,
                  sample.data.fl += dims, missing.data.ptr += dims )
              {
                  CvDTreeNode* predicted_node = 0;
                  
                  // predict oob samples
                  if( !predicted_node )
                      CV_CALL(predicted_node = tree->predict(&sample, &missing, true));
  
                  if( !data->is_classifier ) //regression
                  {
                      double avg_resp, resp = predicted_node->value;
                      oob_predictions_sum.data.fl[i] += (float)resp;
                      oob_num_of_predictions.data.fl[i] += 1;
  
                      // compute oob error
                      avg_resp = oob_predictions_sum.data.fl[i]/oob_num_of_predictions.data.fl[i];
                      avg_resp -= true_resp_ptr[i];
                      oob_error += avg_resp*avg_resp;
                      resp = (resp - true_resp_ptr[i])/maximal_response;
                      ncorrect_responses += exp( -resp*resp );
                  }
                  else //classification
                  {
                      double prdct_resp;
                      CvPoint max_loc;
                      CvMat votes;
  
                      cvGetRow(oob_sample_votes, &votes, i);
                      votes.data.i[predicted_node->class_idx]++;
  
                      // compute oob error
                      cvMinMaxLoc( &votes, 0, 0, 0, &max_loc );
  
                      prdct_resp = data->cat_map->data.i[max_loc.x];
                      oob_error += (fabs(prdct_resp - true_resp_ptr[i]) < FLT_EPSILON) ? 0 : 1;
  
                      ncorrect_responses += cvRound(predicted_node->value - true_resp_ptr[i]) == 0;
                  }
                  oob_samples_count++;
              }
              if( oob_samples_count > 0 )
                  oob_error /= (double)oob_samples_count;
  
              // estimate variable importance
              if( var_importance && oob_samples_count > 0 )
              {
                  int m;
  
                  memcpy( oob_samples_perm_ptr, samples_ptr, dims*nsamples*sizeof(float));
                  for( m = 0; m < dims; m++ )
                  {
                      double ncorrect_responses_permuted = 0;
                      // randomly permute values of the m-th variable in the oob samples
                      float* mth_var_ptr = oob_samples_perm_ptr + m;
  
                      for( i = 0; i < nsamples; i++ )
                      {
                          int i1, i2;
                          float temp;
  
                          i1 = (*rng)(nsamples);
                          i2 = (*rng)(nsamples);
                          CV_SWAP( mth_var_ptr[i1*dims], mth_var_ptr[i2*dims], temp );
  
                          // turn values of (m-1)-th variable, that were permuted
                          // at the previous iteration, untouched
                          if( m > 1 )
                              oob_samples_perm_ptr[i*dims+m-1] = samples_ptr[i*dims+m-1];
                      }
  
                      // predict "permuted" cases and calculate the number of votes for the
                      // correct class in the variable-m-permuted oob data
                      sample  = cvMat( 1, dims, CV_32FC1, oob_samples_perm_ptr );
                      missing = cvMat( 1, dims, CV_8UC1, missing_ptr );
                      for( i = 0; i < nsamples; i++,
                          sample.data.fl += dims, missing.data.ptr += dims )
                      {
                          double predct_resp, true_resp;
  
                          predct_resp = tree->predict(&sample, &missing, true)->value;
                          true_resp   = true_resp_ptr[i];
                          if( data->is_classifier )
                              ncorrect_responses_permuted += cvRound(true_resp - predct_resp) == 0;
                          else
                          {
                              true_resp = (true_resp - predct_resp)/maximal_response;
                              ncorrect_responses_permuted += exp( -true_resp*true_resp );
                          }
                      }
                      var_importance->data.fl[m] += (float)(ncorrect_responses
                          - ncorrect_responses_permuted);
                  }
              }
          }
          ntrees++;
          if( term_crit.type != CV_TERMCRIT_ITER && oob_error < max_oob_err )
              break;
      }
      if( var_importance )
      {
          for ( int vi = 0; vi < var_importance->cols; vi++ )
                  var_importance->data.fl[vi] = ( var_importance->data.fl[vi] > 0 ) ?
                      var_importance->data.fl[vi] : 0;
          cvNormalize( var_importance, var_importance, 1., 0, CV_L1 );
      }
  
      result = true;
      
      cvFree( &oob_samples_perm_ptr );
      cvFree( &samples_ptr );
      cvFree( &missing_ptr );
      cvFree( &true_resp_ptr );
      
      cvReleaseMat( &sample_idx_for_tree );
  
      cvReleaseMat( &oob_sample_votes );
      cvReleaseMat( &oob_responses );
  
      __END__;
  
      return result;
  }
  
  using namespace cv;
  
  bool CvERTrees::train( const Mat& _train_data, int _tflag,
                        const Mat& _responses, const Mat& _var_idx,
                        const Mat& _sample_idx, const Mat& _var_type,
                        const Mat& _missing_mask, CvRTParams params )
  {
      CvMat tdata = _train_data, responses = _responses, vidx = _var_idx,
      sidx = _sample_idx, vtype = _var_type, mmask = _missing_mask;
      return train(&tdata, _tflag, &responses, vidx.data.ptr ? &vidx : 0,
                   sidx.data.ptr ? &sidx : 0, vtype.data.ptr ? &vtype : 0,
                   mmask.data.ptr ? &mmask : 0, params);
  }
  
  // End of file.
  
  