-/*M///////////////////////////////////////////////////////////////////////////////////////\r
-//\r
-// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
-//\r
-// By downloading, copying, installing or using the software you agree to this license.\r
-// If you do not agree to this license, do not download, install,\r
-// copy or use the software.\r
-//\r
-//\r
-// Intel License Agreement\r
-//\r
-// Copyright (C) 2000, Intel Corporation, all rights reserved.\r
-// Third party copyrights are property of their respective owners.\r
-//\r
-// Redistribution and use in source and binary forms, with or without modification,\r
-// are permitted provided that the following conditions are met:\r
-//\r
-// * Redistribution's of source code must retain the above copyright notice,\r
-// this list of conditions and the following disclaimer.\r
-//\r
-// * Redistribution's in binary form must reproduce the above copyright notice,\r
-// this list of conditions and the following disclaimer in the documentation\r
-// and/or other materials provided with the distribution.\r
-//\r
-// * The name of Intel Corporation may not be used to endorse or promote products\r
-// derived from this software without specific prior written permission.\r
-//\r
-// This software is provided by the copyright holders and contributors "as is" and\r
-// any express or implied warranties, including, but not limited to, the implied\r
-// warranties of merchantability and fitness for a particular purpose are disclaimed.\r
-// In no event shall the Intel Corporation or contributors be liable for any direct,\r
-// indirect, incidental, special, exemplary, or consequential damages\r
-// (including, but not limited to, procurement of substitute goods or services;\r
-// loss of use, data, or profits; or business interruption) however caused\r
-// and on any theory of liability, whether in contract, strict liability,\r
-// or tort (including negligence or otherwise) arising in any way out of\r
-// the use of this software, even if advised of the possibility of such damage.\r
-//\r
-//M*/\r
-\r
-#include "_ml.h"\r
-#include <ctype.h>\r
-\r
-static const float ord_nan = FLT_MAX*0.5f;\r
-static const int min_block_size = 1 << 16;\r
-static const int block_size_delta = 1 << 10;\r
-\r
-CvDTreeTrainData::CvDTreeTrainData()\r
-{\r
- var_idx = var_type = cat_count = cat_ofs = cat_map =\r
- priors = priors_mult = counts = buf = direction = split_buf = responses_copy = 0;\r
- tree_storage = temp_storage = 0;\r
-\r
- clear();\r
-}\r
-\r
-\r
-CvDTreeTrainData::CvDTreeTrainData( const CvMat* _train_data, int _tflag,\r
- const CvMat* _responses, const CvMat* _var_idx,\r
- const CvMat* _sample_idx, const CvMat* _var_type,\r
- const CvMat* _missing_mask, const CvDTreeParams& _params,\r
- bool _shared, bool _add_labels )\r
-{\r
- var_idx = var_type = cat_count = cat_ofs = cat_map =\r
- priors = priors_mult = counts = buf = direction = split_buf = responses_copy = 0;\r
-\r
- tree_storage = temp_storage = 0;\r
-\r
- set_data( _train_data, _tflag, _responses, _var_idx, _sample_idx,\r
- _var_type, _missing_mask, _params, _shared, _add_labels );\r
-}\r
-\r
-\r
-CvDTreeTrainData::~CvDTreeTrainData()\r
-{\r
- clear();\r
-}\r
-\r
-\r
-bool CvDTreeTrainData::set_params( const CvDTreeParams& _params )\r
-{\r
- bool ok = false;\r
-\r
- CV_FUNCNAME( "CvDTreeTrainData::set_params" );\r
-\r
- __BEGIN__;\r
-\r
- // set parameters\r
- params = _params;\r
-\r
- if( params.max_categories < 2 )\r
- CV_ERROR( CV_StsOutOfRange, "params.max_categories should be >= 2" );\r
- params.max_categories = MIN( params.max_categories, 15 );\r
-\r
- if( params.max_depth < 0 )\r
- CV_ERROR( CV_StsOutOfRange, "params.max_depth should be >= 0" );\r
- params.max_depth = MIN( params.max_depth, 25 );\r
-\r
- params.min_sample_count = MAX(params.min_sample_count,1);\r
-\r
- if( params.cv_folds < 0 )\r
- CV_ERROR( CV_StsOutOfRange,\r
- "params.cv_folds should be =0 (the tree is not pruned) "\r
- "or n>0 (tree is pruned using n-fold cross-validation)" );\r
-\r
- if( params.cv_folds == 1 )\r
- params.cv_folds = 0;\r
-\r
- if( params.regression_accuracy < 0 )\r
- CV_ERROR( CV_StsOutOfRange, "params.regression_accuracy should be >= 0" );\r
-\r
- ok = true;\r
-\r
- __END__;\r
-\r
- return ok;\r
-}\r
-\r
-#define CV_CMP_NUM_PTR(a,b) (*(a) < *(b))\r
-static CV_IMPLEMENT_QSORT_EX( icvSortIntPtr, int*, CV_CMP_NUM_PTR, int )\r
-static CV_IMPLEMENT_QSORT_EX( icvSortDblPtr, double*, CV_CMP_NUM_PTR, int )\r
-\r
-#define CV_CMP_NUM_IDX(i,j) (aux[i] < aux[j])\r
-static CV_IMPLEMENT_QSORT_EX( icvSortIntAux, int, CV_CMP_NUM_IDX, const float* )\r
-static CV_IMPLEMENT_QSORT_EX( icvSortUShAux, unsigned short, CV_CMP_NUM_IDX, const float* )\r
-\r
-#define CV_CMP_PAIRS(a,b) (*((a).i) < *((b).i))\r
-static CV_IMPLEMENT_QSORT_EX( icvSortPairs, CvPair16u32s, CV_CMP_PAIRS, int )\r
-\r
-void CvDTreeTrainData::set_data( const CvMat* _train_data, int _tflag,\r
- const CvMat* _responses, const CvMat* _var_idx, const CvMat* _sample_idx,\r
- const CvMat* _var_type, const CvMat* _missing_mask, const CvDTreeParams& _params,\r
- bool _shared, bool _add_labels, bool _update_data )\r
-{\r
- CvMat* sample_indices = 0;\r
- CvMat* var_type0 = 0;\r
- CvMat* tmp_map = 0;\r
- int** int_ptr = 0;\r
- CvPair16u32s* pair16u32s_ptr = 0;\r
- CvDTreeTrainData* data = 0;\r
- float *_fdst = 0;\r
- int *_idst = 0;\r
- unsigned short* udst = 0;\r
- int* idst = 0;\r
-\r
- CV_FUNCNAME( "CvDTreeTrainData::set_data" );\r
-\r
- __BEGIN__;\r
-\r
- int sample_all = 0, r_type = 0, cv_n;\r
- int total_c_count = 0;\r
- int tree_block_size, temp_block_size, max_split_size, nv_size, cv_size = 0;\r
- int ds_step, dv_step, ms_step = 0, mv_step = 0; // {data|mask}{sample|var}_step\r
- int vi, i, size;\r
- char err[100];\r
- const int *sidx = 0, *vidx = 0;\r
- \r
- if( _update_data && data_root )\r
- {\r
- data = new CvDTreeTrainData( _train_data, _tflag, _responses, _var_idx,\r
- _sample_idx, _var_type, _missing_mask, _params, _shared, _add_labels );\r
-\r
- // compare new and old train data\r
- if( !(data->var_count == var_count &&\r
- cvNorm( data->var_type, var_type, CV_C ) < FLT_EPSILON &&\r
- cvNorm( data->cat_count, cat_count, CV_C ) < FLT_EPSILON &&\r
- cvNorm( data->cat_map, cat_map, CV_C ) < FLT_EPSILON) )\r
- CV_ERROR( CV_StsBadArg,\r
- "The new training data must have the same types and the input and output variables "\r
- "and the same categories for categorical variables" );\r
-\r
- cvReleaseMat( &priors );\r
- cvReleaseMat( &priors_mult );\r
- cvReleaseMat( &buf );\r
- cvReleaseMat( &direction );\r
- cvReleaseMat( &split_buf );\r
- cvReleaseMemStorage( &temp_storage );\r
-\r
- priors = data->priors; data->priors = 0;\r
- priors_mult = data->priors_mult; data->priors_mult = 0;\r
- buf = data->buf; data->buf = 0;\r
- buf_count = data->buf_count; buf_size = data->buf_size;\r
- sample_count = data->sample_count;\r
-\r
- direction = data->direction; data->direction = 0;\r
- split_buf = data->split_buf; data->split_buf = 0;\r
- temp_storage = data->temp_storage; data->temp_storage = 0;\r
- nv_heap = data->nv_heap; cv_heap = data->cv_heap;\r
-\r
- data_root = new_node( 0, sample_count, 0, 0 );\r
- EXIT;\r
- }\r
-\r
- clear();\r
-\r
- var_all = 0;\r
- rng = cvRNG(-1);\r
-\r
- CV_CALL( set_params( _params ));\r
-\r
- // check parameter types and sizes\r
- CV_CALL( cvCheckTrainData( _train_data, _tflag, _missing_mask, &var_all, &sample_all ));\r
-\r
- train_data = _train_data;\r
- responses = _responses;\r
-\r
- if( _tflag == CV_ROW_SAMPLE )\r
- {\r
- ds_step = _train_data->step/CV_ELEM_SIZE(_train_data->type);\r
- dv_step = 1;\r
- if( _missing_mask )\r
- ms_step = _missing_mask->step, mv_step = 1;\r
- }\r
- else\r
- {\r
- dv_step = _train_data->step/CV_ELEM_SIZE(_train_data->type);\r
- ds_step = 1;\r
- if( _missing_mask )\r
- mv_step = _missing_mask->step, ms_step = 1;\r
- }\r
- tflag = _tflag;\r
-\r
- sample_count = sample_all;\r
- var_count = var_all;\r
- \r
- if( _sample_idx )\r
- {\r
- CV_CALL( sample_indices = cvPreprocessIndexArray( _sample_idx, sample_all ));\r
- sidx = sample_indices->data.i;\r
- sample_count = sample_indices->rows + sample_indices->cols - 1;\r
- }\r
-\r
- if( _var_idx )\r
- {\r
- CV_CALL( var_idx = cvPreprocessIndexArray( _var_idx, var_all ));\r
- vidx = var_idx->data.i;\r
- var_count = var_idx->rows + var_idx->cols - 1;\r
- }\r
-\r
- is_buf_16u = false; \r
- if ( sample_count < 65536 ) \r
- is_buf_16u = true; \r
- \r
- if( !CV_IS_MAT(_responses) ||\r
- (CV_MAT_TYPE(_responses->type) != CV_32SC1 &&\r
- CV_MAT_TYPE(_responses->type) != CV_32FC1) ||\r
- (_responses->rows != 1 && _responses->cols != 1) ||\r
- _responses->rows + _responses->cols - 1 != sample_all )\r
- CV_ERROR( CV_StsBadArg, "The array of _responses must be an integer or "\r
- "floating-point vector containing as many elements as "\r
- "the total number of samples in the training data matrix" );\r
- \r
- \r
- CV_CALL( var_type0 = cvPreprocessVarType( _var_type, var_idx, var_count, &r_type ));\r
-\r
- CV_CALL( var_type = cvCreateMat( 1, var_count+2, CV_32SC1 ));\r
- \r
- \r
- cat_var_count = 0;\r
- ord_var_count = -1;\r
-\r
- is_classifier = r_type == CV_VAR_CATEGORICAL;\r
-\r
- // step 0. calc the number of categorical vars\r
- for( vi = 0; vi < var_count; vi++ )\r
- {\r
- var_type->data.i[vi] = var_type0->data.ptr[vi] == CV_VAR_CATEGORICAL ?\r
- cat_var_count++ : ord_var_count--;\r
- }\r
-\r
- ord_var_count = ~ord_var_count;\r
- cv_n = params.cv_folds;\r
- // set the two last elements of var_type array to be able\r
- // to locate responses and cross-validation labels using\r
- // the corresponding get_* functions.\r
- var_type->data.i[var_count] = cat_var_count;\r
- var_type->data.i[var_count+1] = cat_var_count+1;\r
-\r
- // in case of single ordered predictor we need dummy cv_labels\r
- // for safe split_node_data() operation\r
- have_labels = cv_n > 0 || (ord_var_count == 1 && cat_var_count == 0) || _add_labels;\r
-\r
- work_var_count = var_count + (is_classifier ? 1 : 0) + (have_labels ? 1 : 0);\r
- buf_size = (work_var_count + 1)*sample_count;\r
- shared = _shared;\r
- buf_count = shared ? 2 : 1;\r
- \r
- if ( is_buf_16u )\r
- {\r
- CV_CALL( buf = cvCreateMat( buf_count, buf_size, CV_16UC1 ));\r
- CV_CALL( pair16u32s_ptr = (CvPair16u32s*)cvAlloc( sample_count*sizeof(pair16u32s_ptr[0]) ));\r
- }\r
- else\r
- {\r
- CV_CALL( buf = cvCreateMat( buf_count, buf_size, CV_32SC1 ));\r
- CV_CALL( int_ptr = (int**)cvAlloc( sample_count*sizeof(int_ptr[0]) ));\r
- } \r
-\r
- size = is_classifier ? (cat_var_count+1) : cat_var_count;\r
- size = !size ? 1 : size;\r
- CV_CALL( cat_count = cvCreateMat( 1, size, CV_32SC1 ));\r
- CV_CALL( cat_ofs = cvCreateMat( 1, size, CV_32SC1 ));\r
- \r
- size = is_classifier ? (cat_var_count + 1)*params.max_categories : cat_var_count*params.max_categories;\r
- size = !size ? 1 : size;\r
- CV_CALL( cat_map = cvCreateMat( 1, size, CV_32SC1 ));\r
-\r
- // now calculate the maximum size of split,\r
- // create memory storage that will keep nodes and splits of the decision tree\r
- // allocate root node and the buffer for the whole training data\r
- max_split_size = cvAlign(sizeof(CvDTreeSplit) +\r
- (MAX(0,sample_count - 33)/32)*sizeof(int),sizeof(void*));\r
- tree_block_size = MAX((int)sizeof(CvDTreeNode)*8, max_split_size);\r
- tree_block_size = MAX(tree_block_size + block_size_delta, min_block_size);\r
- CV_CALL( tree_storage = cvCreateMemStorage( tree_block_size ));\r
- CV_CALL( node_heap = cvCreateSet( 0, sizeof(*node_heap), sizeof(CvDTreeNode), tree_storage ));\r
-\r
- nv_size = var_count*sizeof(int);\r
- nv_size = cvAlign(MAX( nv_size, (int)sizeof(CvSetElem) ), sizeof(void*));\r
-\r
- temp_block_size = nv_size;\r
-\r
- if( cv_n )\r
- {\r
- if( sample_count < cv_n*MAX(params.min_sample_count,10) )\r
- CV_ERROR( CV_StsOutOfRange,\r
- "The many folds in cross-validation for such a small dataset" );\r
-\r
- cv_size = cvAlign( cv_n*(sizeof(int) + sizeof(double)*2), sizeof(double) );\r
- temp_block_size = MAX(temp_block_size, cv_size);\r
- }\r
-\r
- temp_block_size = MAX( temp_block_size + block_size_delta, min_block_size );\r
- CV_CALL( temp_storage = cvCreateMemStorage( temp_block_size ));\r
- CV_CALL( nv_heap = cvCreateSet( 0, sizeof(*nv_heap), nv_size, temp_storage ));\r
- if( cv_size )\r
- CV_CALL( cv_heap = cvCreateSet( 0, sizeof(*cv_heap), cv_size, temp_storage ));\r
-\r
- CV_CALL( data_root = new_node( 0, sample_count, 0, 0 ));\r
-\r
- max_c_count = 1;\r
-\r
- _fdst = 0;\r
- _idst = 0;\r
- if (ord_var_count)\r
- _fdst = (float*)cvAlloc(sample_count*sizeof(_fdst[0]));\r
- if (is_buf_16u && (cat_var_count || is_classifier))\r
- _idst = (int*)cvAlloc(sample_count*sizeof(_idst[0]));\r
-\r
- // transform the training data to convenient representation\r
- for( vi = 0; vi <= var_count; vi++ )\r
- {\r
- int ci;\r
- const uchar* mask = 0;\r
- int m_step = 0, step;\r
- const int* idata = 0;\r
- const float* fdata = 0;\r
- int num_valid = 0;\r
-\r
- if( vi < var_count ) // analyze i-th input variable\r
- {\r
- int vi0 = vidx ? vidx[vi] : vi;\r
- ci = get_var_type(vi);\r
- step = ds_step; m_step = ms_step;\r
- if( CV_MAT_TYPE(_train_data->type) == CV_32SC1 )\r
- idata = _train_data->data.i + vi0*dv_step;\r
- else\r
- fdata = _train_data->data.fl + vi0*dv_step;\r
- if( _missing_mask )\r
- mask = _missing_mask->data.ptr + vi0*mv_step;\r
- }\r
- else // analyze _responses\r
- {\r
- ci = cat_var_count;\r
- step = CV_IS_MAT_CONT(_responses->type) ?\r
- 1 : _responses->step / CV_ELEM_SIZE(_responses->type);\r
- if( CV_MAT_TYPE(_responses->type) == CV_32SC1 )\r
- idata = _responses->data.i;\r
- else\r
- fdata = _responses->data.fl;\r
- }\r
-\r
- if( (vi < var_count && ci>=0) ||\r
- (vi == var_count && is_classifier) ) // process categorical variable or response\r
- {\r
- int c_count, prev_label;\r
- int* c_map;\r
- \r
- if (is_buf_16u)\r
- udst = (unsigned short*)(buf->data.s + vi*sample_count);\r
- else\r
- idst = buf->data.i + vi*sample_count;\r
- \r
- // copy data\r
- for( i = 0; i < sample_count; i++ )\r
- {\r
- int val = INT_MAX, si = sidx ? sidx[i] : i;\r
- if( !mask || !mask[si*m_step] )\r
- {\r
- if( idata )\r
- val = idata[si*step];\r
- else\r
- {\r
- float t = fdata[si*step];\r
- val = cvRound(t);\r
- if( val != t )\r
- {\r
- sprintf( err, "%d-th value of %d-th (categorical) "\r
- "variable is not an integer", i, vi );\r
- CV_ERROR( CV_StsBadArg, err );\r
- }\r
- }\r
-\r
- if( val == INT_MAX )\r
- {\r
- sprintf( err, "%d-th value of %d-th (categorical) "\r
- "variable is too large", i, vi );\r
- CV_ERROR( CV_StsBadArg, err );\r
- }\r
- num_valid++;\r
- }\r
- if (is_buf_16u)\r
- {\r
- _idst[i] = val;\r
- pair16u32s_ptr[i].u = udst + i;\r
- pair16u32s_ptr[i].i = _idst + i;\r
- } \r
- else\r
- {\r
- idst[i] = val;\r
- int_ptr[i] = idst + i;\r
- }\r
- }\r
-\r
- c_count = num_valid > 0;\r
- if (is_buf_16u)\r
- {\r
- icvSortPairs( pair16u32s_ptr, sample_count, 0 );\r
- // count the categories\r
- for( i = 1; i < num_valid; i++ )\r
- if (*pair16u32s_ptr[i].i != *pair16u32s_ptr[i-1].i)\r
- c_count ++ ;\r
- }\r
- else\r
- {\r
- icvSortIntPtr( int_ptr, sample_count, 0 );\r
- // count the categories\r
- for( i = 1; i < num_valid; i++ )\r
- c_count += *int_ptr[i] != *int_ptr[i-1];\r
- }\r
-\r
- if( vi > 0 )\r
- max_c_count = MAX( max_c_count, c_count );\r
- cat_count->data.i[ci] = c_count;\r
- cat_ofs->data.i[ci] = total_c_count;\r
-\r
- // resize cat_map, if need\r
- if( cat_map->cols < total_c_count + c_count )\r
- {\r
- tmp_map = cat_map;\r
- CV_CALL( cat_map = cvCreateMat( 1,\r
- MAX(cat_map->cols*3/2,total_c_count+c_count), CV_32SC1 ));\r
- for( i = 0; i < total_c_count; i++ )\r
- cat_map->data.i[i] = tmp_map->data.i[i];\r
- cvReleaseMat( &tmp_map );\r
- }\r
-\r
- c_map = cat_map->data.i + total_c_count;\r
- total_c_count += c_count;\r
-\r
- c_count = -1;\r
- if (is_buf_16u)\r
- {\r
- // compact the class indices and build the map\r
- prev_label = ~*pair16u32s_ptr[0].i;\r
- for( i = 0; i < num_valid; i++ )\r
- {\r
- int cur_label = *pair16u32s_ptr[i].i;\r
- if( cur_label != prev_label )\r
- c_map[++c_count] = prev_label = cur_label;\r
- *pair16u32s_ptr[i].u = (unsigned short)c_count;\r
- }\r
- // replace labels for missing values with -1\r
- for( ; i < sample_count; i++ )\r
- *pair16u32s_ptr[i].u = 65535;\r
- }\r
- else\r
- {\r
- // compact the class indices and build the map\r
- prev_label = ~*int_ptr[0];\r
- for( i = 0; i < num_valid; i++ )\r
- {\r
- int cur_label = *int_ptr[i];\r
- if( cur_label != prev_label )\r
- c_map[++c_count] = prev_label = cur_label;\r
- *int_ptr[i] = c_count;\r
- }\r
- // replace labels for missing values with -1\r
- for( ; i < sample_count; i++ )\r
- *int_ptr[i] = -1;\r
- } \r
- }\r
- else if( ci < 0 ) // process ordered variable\r
- {\r
- if (is_buf_16u)\r
- udst = (unsigned short*)(buf->data.s + vi*sample_count);\r
- else\r
- idst = buf->data.i + vi*sample_count;\r
-\r
- for( i = 0; i < sample_count; i++ )\r
- {\r
- float val = ord_nan;\r
- int si = sidx ? sidx[i] : i;\r
- if( !mask || !mask[si*m_step] )\r
- {\r
- if( idata )\r
- val = (float)idata[si*step];\r
- else\r
- val = fdata[si*step];\r
-\r
- if( fabs(val) >= ord_nan )\r
- {\r
- sprintf( err, "%d-th value of %d-th (ordered) "\r
- "variable (=%g) is too large", i, vi, val );\r
- CV_ERROR( CV_StsBadArg, err );\r
- }\r
- }\r
- num_valid++;\r
- if (is_buf_16u)\r
- udst[i] = (unsigned short)i;\r
- else\r
- idst[i] = i; // ïåðåíåñòè âûøå â if( idata )\r
- _fdst[i] = val;\r
- \r
- }\r
- if (is_buf_16u)\r
- icvSortUShAux( udst, num_valid, _fdst);\r
- else\r
- icvSortIntAux( idst, /*or num_valid?\*/ sample_count, _fdst );\r
- }\r
- \r
- if( vi < var_count )\r
- data_root->set_num_valid(vi, num_valid);\r
- }\r
-\r
- // set sample labels\r
- if (is_buf_16u)\r
- udst = (unsigned short*)(buf->data.s + work_var_count*sample_count);\r
- else\r
- idst = buf->data.i + work_var_count*sample_count;\r
-\r
- for (i = 0; i < sample_count; i++)\r
- {\r
- if (udst)\r
- udst[i] = sidx ? (unsigned short)sidx[i] : (unsigned short)i;\r
- else\r
- idst[i] = sidx ? sidx[i] : i;\r
- }\r
-\r
- if( cv_n )\r
- {\r
- unsigned short* udst = 0;\r
- int* idst = 0;\r
- CvRNG* r = &rng;\r
-\r
- if (is_buf_16u)\r
- {\r
- udst = (unsigned short*)(buf->data.s + (get_work_var_count()-1)*sample_count);\r
- for( i = vi = 0; i < sample_count; i++ )\r
- {\r
- udst[i] = (unsigned short)vi++;\r
- vi &= vi < cv_n ? -1 : 0;\r
- }\r
-\r
- for( i = 0; i < sample_count; i++ )\r
- {\r
- int a = cvRandInt(r) % sample_count;\r
- int b = cvRandInt(r) % sample_count;\r
- unsigned short unsh = (unsigned short)vi;\r
- CV_SWAP( udst[a], udst[b], unsh );\r
- }\r
- }\r
- else\r
- {\r
- idst = buf->data.i + (get_work_var_count()-1)*sample_count;\r
- for( i = vi = 0; i < sample_count; i++ )\r
- {\r
- idst[i] = vi++;\r
- vi &= vi < cv_n ? -1 : 0;\r
- }\r
-\r
- for( i = 0; i < sample_count; i++ )\r
- {\r
- int a = cvRandInt(r) % sample_count;\r
- int b = cvRandInt(r) % sample_count;\r
- CV_SWAP( idst[a], idst[b], vi );\r
- }\r
- }\r
- }\r
-\r
- if ( cat_map ) \r
- cat_map->cols = MAX( total_c_count, 1 );\r
-\r
- max_split_size = cvAlign(sizeof(CvDTreeSplit) +\r
- (MAX(0,max_c_count - 33)/32)*sizeof(int),sizeof(void*));\r
- CV_CALL( split_heap = cvCreateSet( 0, sizeof(*split_heap), max_split_size, tree_storage ));\r
-\r
- have_priors = is_classifier && params.priors;\r
- if( is_classifier )\r
- {\r
- int m = get_num_classes();\r
- double sum = 0;\r
- CV_CALL( priors = cvCreateMat( 1, m, CV_64F ));\r
- for( i = 0; i < m; i++ )\r
- {\r
- double val = have_priors ? params.priors[i] : 1.;\r
- if( val <= 0 )\r
- CV_ERROR( CV_StsOutOfRange, "Every class weight should be positive" );\r
- priors->data.db[i] = val;\r
- sum += val;\r
- }\r
-\r
- // normalize weights\r
- if( have_priors )\r
- cvScale( priors, priors, 1./sum );\r
-\r
- CV_CALL( priors_mult = cvCloneMat( priors ));\r
- CV_CALL( counts = cvCreateMat( 1, m, CV_32SC1 ));\r
- }\r
-\r
-\r
- CV_CALL( direction = cvCreateMat( 1, sample_count, CV_8UC1 ));\r
- CV_CALL( split_buf = cvCreateMat( 1, sample_count, CV_32SC1 ));\r
-\r
- {\r
- int maxNumThreads = 1;\r
-#ifdef _OPENMP\r
- maxNumThreads = cv::getNumThreads();\r
-#endif\r
- pred_float_buf.resize(maxNumThreads);\r
- pred_int_buf.resize(maxNumThreads);\r
- resp_float_buf.resize(maxNumThreads);\r
- resp_int_buf.resize(maxNumThreads);\r
- cv_lables_buf.resize(maxNumThreads);\r
- sample_idx_buf.resize(maxNumThreads);\r
- for( int ti = 0; ti < maxNumThreads; ti++ )\r
- {\r
- pred_float_buf[ti].resize(sample_count);\r
- pred_int_buf[ti].resize(sample_count);\r
- resp_float_buf[ti].resize(sample_count);\r
- resp_int_buf[ti].resize(sample_count);\r
- cv_lables_buf[ti].resize(sample_count);\r
- sample_idx_buf[ti].resize(sample_count);\r
- }\r
- }\r
-\r
- __END__;\r
-\r
- if( data )\r
- delete data;\r
-\r
- if (_fdst)\r
- cvFree( &_fdst );\r
- if (_idst)\r
- cvFree( &_idst );\r
- cvFree( &int_ptr );\r
- cvReleaseMat( &var_type0 );\r
- cvReleaseMat( &sample_indices );\r
- cvReleaseMat( &tmp_map );\r
-}\r
-\r
-void CvDTreeTrainData::do_responses_copy()\r
-{\r
- responses_copy = cvCreateMat( responses->rows, responses->cols, responses->type );\r
- cvCopy( responses, responses_copy);\r
- responses = responses_copy;\r
-}\r
-\r
-CvDTreeNode* CvDTreeTrainData::subsample_data( const CvMat* _subsample_idx )\r
-{\r
- CvDTreeNode* root = 0;\r
- CvMat* isubsample_idx = 0;\r
- CvMat* subsample_co = 0;\r
-\r
- CV_FUNCNAME( "CvDTreeTrainData::subsample_data" );\r
-\r
- __BEGIN__;\r
-\r
- if( !data_root )\r
- CV_ERROR( CV_StsError, "No training data has been set" );\r
-\r
- if( _subsample_idx )\r
- CV_CALL( isubsample_idx = cvPreprocessIndexArray( _subsample_idx, sample_count ));\r
-\r
- if( !isubsample_idx )\r
- {\r
- // make a copy of the root node\r
- CvDTreeNode temp;\r
- int i;\r
- root = new_node( 0, 1, 0, 0 );\r
- temp = *root;\r
- *root = *data_root;\r
- root->num_valid = temp.num_valid;\r
- if( root->num_valid )\r
- {\r
- for( i = 0; i < var_count; i++ )\r
- root->num_valid[i] = data_root->num_valid[i];\r
- }\r
- root->cv_Tn = temp.cv_Tn;\r
- root->cv_node_risk = temp.cv_node_risk;\r
- root->cv_node_error = temp.cv_node_error;\r
- }\r
- else\r
- {\r
- int* sidx = isubsample_idx->data.i;\r
- // co - array of count/offset pairs (to handle duplicated values in _subsample_idx)\r
- int* co, cur_ofs = 0;\r
- int vi, i;\r
- int work_var_count = get_work_var_count();\r
- int count = isubsample_idx->rows + isubsample_idx->cols - 1;\r
-\r
- root = new_node( 0, count, 1, 0 );\r
-\r
- CV_CALL( subsample_co = cvCreateMat( 1, sample_count*2, CV_32SC1 ));\r
- cvZero( subsample_co );\r
- co = subsample_co->data.i;\r
- for( i = 0; i < count; i++ )\r
- co[sidx[i]*2]++;\r
- for( i = 0; i < sample_count; i++ )\r
- {\r
- if( co[i*2] )\r
- {\r
- co[i*2+1] = cur_ofs;\r
- cur_ofs += co[i*2];\r
- }\r
- else\r
- co[i*2+1] = -1;\r
- }\r
-\r
- for( vi = 0; vi < work_var_count; vi++ )\r
- {\r
- int ci = get_var_type(vi);\r
-\r
- if( ci >= 0 || vi >= var_count )\r
- {\r
- int* src_buf = get_pred_int_buf();\r
- const int* src = 0;\r
- int num_valid = 0;\r
- \r
- get_cat_var_data( data_root, vi, src_buf, &src );\r
-\r
- if (is_buf_16u)\r
- {\r
- unsigned short* udst = (unsigned short*)(buf->data.s + root->buf_idx*buf->cols + \r
- vi*sample_count + root->offset);\r
- for( i = 0; i < count; i++ )\r
- {\r
- int val = src[sidx[i]];\r
- udst[i] = (unsigned short)val;\r
- num_valid += val >= 0;\r
- }\r
- }\r
- else\r
- {\r
- int* idst = buf->data.i + root->buf_idx*buf->cols + \r
- vi*sample_count + root->offset;\r
- for( i = 0; i < count; i++ )\r
- {\r
- int val = src[sidx[i]];\r
- idst[i] = val;\r
- num_valid += val >= 0;\r
- }\r
- }\r
-\r
- if( vi < var_count )\r
- root->set_num_valid(vi, num_valid);\r
- }\r
- else\r
- {\r
- int *src_idx_buf = get_pred_int_buf();\r
- const int* src_idx = 0;\r
- float *src_val_buf = get_pred_float_buf();\r
- const float* src_val = 0;\r
- int j = 0, idx, count_i;\r
- int num_valid = data_root->get_num_valid(vi);\r
-\r
- get_ord_var_data( data_root, vi, src_val_buf, src_idx_buf, &src_val, &src_idx );\r
- if (is_buf_16u)\r
- {\r
- unsigned short* udst_idx = (unsigned short*)(buf->data.s + root->buf_idx*buf->cols + \r
- vi*sample_count + data_root->offset);\r
- for( i = 0; i < num_valid; i++ )\r
- {\r
- idx = src_idx[i];\r
- count_i = co[idx*2];\r
- if( count_i )\r
- for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )\r
- udst_idx[j] = (unsigned short)cur_ofs;\r
- }\r
-\r
- root->set_num_valid(vi, j);\r
-\r
- for( ; i < sample_count; i++ )\r
- {\r
- idx = src_idx[i];\r
- count_i = co[idx*2];\r
- if( count_i )\r
- for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )\r
- udst_idx[j] = (unsigned short)cur_ofs;\r
- }\r
- }\r
- else\r
- {\r
- int* idst_idx = buf->data.i + root->buf_idx*buf->cols + \r
- vi*sample_count + root->offset;\r
- for( i = 0; i < num_valid; i++ )\r
- {\r
- idx = src_idx[i];\r
- count_i = co[idx*2];\r
- if( count_i )\r
- for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )\r
- idst_idx[j] = cur_ofs;\r
- }\r
-\r
- root->set_num_valid(vi, j);\r
-\r
- for( ; i < sample_count; i++ )\r
- {\r
- idx = src_idx[i];\r
- count_i = co[idx*2];\r
- if( count_i )\r
- for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )\r
- idst_idx[j] = cur_ofs;\r
- }\r
- }\r
- }\r
- }\r
- // sample indices subsampling\r
- int* sample_idx_src_buf = get_sample_idx_buf();\r
- const int* sample_idx_src = 0;\r
- get_sample_indices(data_root, sample_idx_src_buf, &sample_idx_src);\r
- if (is_buf_16u)\r
- {\r
- unsigned short* sample_idx_dst = (unsigned short*)(buf->data.s + root->buf_idx*buf->cols + \r
- get_work_var_count()*sample_count + root->offset); \r
- for (i = 0; i < count; i++)\r
- sample_idx_dst[i] = (unsigned short)sample_idx_src[sidx[i]];\r
- }\r
- else\r
- {\r
- int* sample_idx_dst = buf->data.i + root->buf_idx*buf->cols + \r
- get_work_var_count()*sample_count + root->offset; \r
- for (i = 0; i < count; i++)\r
- sample_idx_dst[i] = sample_idx_src[sidx[i]];\r
- }\r
- }\r
-\r
- __END__;\r
-\r
- cvReleaseMat( &isubsample_idx );\r
- cvReleaseMat( &subsample_co );\r
-\r
- return root;\r
-}\r
-\r
-\r
-void CvDTreeTrainData::get_vectors( const CvMat* _subsample_idx,\r
- float* values, uchar* missing,\r
- float* responses, bool get_class_idx )\r
-{\r
- CvMat* subsample_idx = 0;\r
- CvMat* subsample_co = 0;\r
-\r
- CV_FUNCNAME( "CvDTreeTrainData::get_vectors" );\r
-\r
- __BEGIN__;\r
-\r
- int i, vi, total = sample_count, count = total, cur_ofs = 0;\r
- int* sidx = 0;\r
- int* co = 0;\r
-\r
- if( _subsample_idx )\r
- {\r
- CV_CALL( subsample_idx = cvPreprocessIndexArray( _subsample_idx, sample_count ));\r
- sidx = subsample_idx->data.i;\r
- CV_CALL( subsample_co = cvCreateMat( 1, sample_count*2, CV_32SC1 ));\r
- co = subsample_co->data.i;\r
- cvZero( subsample_co );\r
- count = subsample_idx->cols + subsample_idx->rows - 1;\r
- for( i = 0; i < count; i++ )\r
- co[sidx[i]*2]++;\r
- for( i = 0; i < total; i++ )\r
- {\r
- int count_i = co[i*2];\r
- if( count_i )\r
- {\r
- co[i*2+1] = cur_ofs*var_count;\r
- cur_ofs += count_i;\r
- }\r
- }\r
- }\r
-\r
- if( missing )\r
- memset( missing, 1, count*var_count );\r
-\r
- for( vi = 0; vi < var_count; vi++ )\r
- {\r
- int ci = get_var_type(vi);\r
- if( ci >= 0 ) // categorical\r
- {\r
- float* dst = values + vi;\r
- uchar* m = missing ? missing + vi : 0;\r
- int* src_buf = get_pred_int_buf();\r
- const int* src = 0; \r
- get_cat_var_data(data_root, vi, src_buf, &src);\r
-\r
- for( i = 0; i < count; i++, dst += var_count )\r
- {\r
- int idx = sidx ? sidx[i] : i;\r
- int val = src[idx];\r
- *dst = (float)val;\r
- if( m )\r
- {\r
- *m = (!is_buf_16u && val < 0) || (is_buf_16u && (val == 65535));\r
- m += var_count;\r
- }\r
- }\r
- }\r
- else // ordered\r
- {\r
- float* dst = values + vi;\r
- uchar* m = missing ? missing + vi : 0;\r
- int count1 = data_root->get_num_valid(vi);\r
- float *src_val_buf = get_pred_float_buf();\r
- const float *src_val = 0;\r
- int* src_idx_buf = get_pred_int_buf();\r
- const int* src_idx = 0;\r
- get_ord_var_data(data_root, vi, src_val_buf, src_idx_buf, &src_val, &src_idx);\r
-\r
- for( i = 0; i < count1; i++ )\r
- {\r
- int idx = src_idx[i];\r
- int count_i = 1;\r
- if( co )\r
- {\r
- count_i = co[idx*2];\r
- cur_ofs = co[idx*2+1];\r
- }\r
- else\r
- cur_ofs = idx*var_count;\r
- if( count_i )\r
- {\r
- float val = src_val[i];\r
- for( ; count_i > 0; count_i--, cur_ofs += var_count )\r
- {\r
- dst[cur_ofs] = val;\r
- if( m )\r
- m[cur_ofs] = 0;\r
- }\r
- }\r
- }\r
- }\r
- }\r
-\r
- // copy responses\r
- if( responses )\r
- {\r
- if( is_classifier )\r
- {\r
- int* src_buf = get_resp_int_buf();\r
- const int* src = 0;\r
- get_class_labels(data_root, src_buf, &src);\r
- for( i = 0; i < count; i++ )\r
- {\r
- int idx = sidx ? sidx[i] : i;\r
- int val = get_class_idx ? src[idx] :\r
- cat_map->data.i[cat_ofs->data.i[cat_var_count]+src[idx]];\r
- responses[i] = (float)val;\r
- }\r
- }\r
- else\r
- {\r
- float *_values_buf = get_resp_float_buf();\r
- const float* _values = 0;\r
- get_ord_responses(data_root, _values_buf, &_values);\r
- for( i = 0; i < count; i++ )\r
- {\r
- int idx = sidx ? sidx[i] : i;\r
- responses[i] = _values[idx];\r
- }\r
- }\r
- }\r
-\r
- __END__;\r
-\r
- cvReleaseMat( &subsample_idx );\r
- cvReleaseMat( &subsample_co );\r
-}\r
-\r
-\r
-CvDTreeNode* CvDTreeTrainData::new_node( CvDTreeNode* parent, int count,\r
- int storage_idx, int offset )\r
-{\r
- CvDTreeNode* node = (CvDTreeNode*)cvSetNew( node_heap );\r
-\r
- node->sample_count = count;\r
- node->depth = parent ? parent->depth + 1 : 0;\r
- node->parent = parent;\r
- node->left = node->right = 0;\r
- node->split = 0;\r
- node->value = 0;\r
- node->class_idx = 0;\r
- node->maxlr = 0.;\r
-\r
- node->buf_idx = storage_idx;\r
- node->offset = offset;\r
- if( nv_heap )\r
- node->num_valid = (int*)cvSetNew( nv_heap );\r
- else\r
- node->num_valid = 0;\r
- node->alpha = node->node_risk = node->tree_risk = node->tree_error = 0.;\r
- node->complexity = 0;\r
-\r
- if( params.cv_folds > 0 && cv_heap )\r
- {\r
- int cv_n = params.cv_folds;\r
- node->Tn = INT_MAX;\r
- node->cv_Tn = (int*)cvSetNew( cv_heap );\r
- node->cv_node_risk = (double*)cvAlignPtr(node->cv_Tn + cv_n, sizeof(double));\r
- node->cv_node_error = node->cv_node_risk + cv_n;\r
- }\r
- else\r
- {\r
- node->Tn = 0;\r
- node->cv_Tn = 0;\r
- node->cv_node_risk = 0;\r
- node->cv_node_error = 0;\r
- }\r
-\r
- return node;\r
-}\r
-\r
-\r
-CvDTreeSplit* CvDTreeTrainData::new_split_ord( int vi, float cmp_val,\r
- int split_point, int inversed, float quality )\r
-{\r
- CvDTreeSplit* split = (CvDTreeSplit*)cvSetNew( split_heap );\r
- split->var_idx = vi;\r
- split->condensed_idx = INT_MIN;\r
- split->ord.c = cmp_val;\r
- split->ord.split_point = split_point;\r
- split->inversed = inversed;\r
- split->quality = quality;\r
- split->next = 0;\r
-\r
- return split;\r
-}\r
-\r
-\r
-CvDTreeSplit* CvDTreeTrainData::new_split_cat( int vi, float quality )\r
-{\r
- CvDTreeSplit* split = (CvDTreeSplit*)cvSetNew( split_heap );\r
- int i, n = (max_c_count + 31)/32;\r
-\r
- split->var_idx = vi;\r
- split->condensed_idx = INT_MIN;\r
- split->inversed = 0;\r
- split->quality = quality;\r
- for( i = 0; i < n; i++ )\r
- split->subset[i] = 0;\r
- split->next = 0;\r
-\r
- return split;\r
-}\r
-\r
-\r
-void CvDTreeTrainData::free_node( CvDTreeNode* node )\r
-{\r
- CvDTreeSplit* split = node->split;\r
- free_node_data( node );\r
- while( split )\r
- {\r
- CvDTreeSplit* next = split->next;\r
- cvSetRemoveByPtr( split_heap, split );\r
- split = next;\r
- }\r
- node->split = 0;\r
- cvSetRemoveByPtr( node_heap, node );\r
-}\r
-\r
-\r
-void CvDTreeTrainData::free_node_data( CvDTreeNode* node )\r
-{\r
- if( node->num_valid )\r
- {\r
- cvSetRemoveByPtr( nv_heap, node->num_valid );\r
- node->num_valid = 0;\r
- }\r
- // do not free cv_* fields, as all the cross-validation related data is released at once.\r
-}\r
-\r
-\r
-void CvDTreeTrainData::free_train_data()\r
-{\r
- cvReleaseMat( &counts );\r
- cvReleaseMat( &buf );\r
- cvReleaseMat( &direction );\r
- cvReleaseMat( &split_buf );\r
- cvReleaseMemStorage( &temp_storage );\r
- cvReleaseMat( &responses_copy );\r
- pred_float_buf.clear();\r
- pred_int_buf.clear();\r
- resp_float_buf.clear();\r
- resp_int_buf.clear();\r
- cv_lables_buf.clear();\r
- sample_idx_buf.clear();\r
-\r
- cv_heap = nv_heap = 0;\r
-}\r
-\r
-\r
-void CvDTreeTrainData::clear()\r
-{\r
- free_train_data();\r
-\r
- cvReleaseMemStorage( &tree_storage );\r
-\r
- cvReleaseMat( &var_idx );\r
- cvReleaseMat( &var_type );\r
- cvReleaseMat( &cat_count );\r
- cvReleaseMat( &cat_ofs );\r
- cvReleaseMat( &cat_map );\r
- cvReleaseMat( &priors );\r
- cvReleaseMat( &priors_mult );\r
- \r
- node_heap = split_heap = 0;\r
-\r
- sample_count = var_all = var_count = max_c_count = ord_var_count = cat_var_count = 0;\r
- have_labels = have_priors = is_classifier = false;\r
-\r
- buf_count = buf_size = 0;\r
- shared = false;\r
- \r
- data_root = 0;\r
-\r
- rng = cvRNG(-1);\r
-}\r
-\r
-\r
-int CvDTreeTrainData::get_num_classes() const\r
-{\r
- return is_classifier ? cat_count->data.i[cat_var_count] : 0;\r
-}\r
-\r
-\r
-int CvDTreeTrainData::get_var_type(int vi) const\r
-{\r
- return var_type->data.i[vi];\r
-}\r
-\r
-int CvDTreeTrainData::get_ord_var_data( CvDTreeNode* n, int vi, float* ord_values_buf, int* indices_buf, const float** ord_values, const int** indices )\r
-{\r
- int vidx = var_idx ? var_idx->data.i[vi] : vi;\r
- int node_sample_count = n->sample_count; \r
- int* sample_indices_buf = get_sample_idx_buf();\r
- const int* sample_indices = 0;\r
- int td_step = train_data->step/CV_ELEM_SIZE(train_data->type);\r
-\r
- get_sample_indices(n, sample_indices_buf, &sample_indices);\r
-\r
- if( !is_buf_16u )\r
- *indices = buf->data.i + n->buf_idx*buf->cols + \r
- vi*sample_count + n->offset;\r
- else {\r
- const unsigned short* short_indices = (const unsigned short*)(buf->data.s + n->buf_idx*buf->cols + \r
- vi*sample_count + n->offset );\r
- for( int i = 0; i < node_sample_count; i++ )\r
- indices_buf[i] = short_indices[i];\r
- *indices = indices_buf;\r
- }\r
- \r
- if( tflag == CV_ROW_SAMPLE )\r
- {\r
- for( int i = 0; i < node_sample_count && \r
- ((((*indices)[i] >= 0) && !is_buf_16u) || (((*indices)[i] != 65535) && is_buf_16u)); i++ )\r
- {\r
- int idx = (*indices)[i];\r
- idx = sample_indices[idx];\r
- ord_values_buf[i] = *(train_data->data.fl + idx * td_step + vidx);\r
- }\r
- }\r
- else\r
- for( int i = 0; i < node_sample_count && \r
- ((((*indices)[i] >= 0) && !is_buf_16u) || (((*indices)[i] != 65535) && is_buf_16u)); i++ )\r
- {\r
- int idx = (*indices)[i];\r
- idx = sample_indices[idx];\r
- ord_values_buf[i] = *(train_data->data.fl + vidx* td_step + idx);\r
- }\r
- \r
- *ord_values = ord_values_buf;\r
- return 0; //TODO: return the number of non-missing values\r
-}\r
-\r
-\r
-void CvDTreeTrainData::get_class_labels( CvDTreeNode* n, int* labels_buf, const int** labels )\r
-{\r
- if (is_classifier)\r
- get_cat_var_data( n, var_count, labels_buf, labels );\r
-}\r
-\r
-void CvDTreeTrainData::get_sample_indices( CvDTreeNode* n, int* indices_buf, const int** indices )\r
-{\r
- get_cat_var_data( n, get_work_var_count(), indices_buf, indices );\r
-}\r
-\r
-void CvDTreeTrainData::get_ord_responses( CvDTreeNode* n, float* values_buf, const float** values)\r
-{\r
- int sample_count = n->sample_count;\r
- int* indices_buf = get_sample_idx_buf();\r
- const int* indices = 0;\r
-\r
- int r_step = responses->step/CV_ELEM_SIZE(responses->type);\r
-\r
- get_sample_indices(n, indices_buf, &indices);\r
-\r
- \r
- for( int i = 0; i < sample_count && \r
- (((indices[i] >= 0) && !is_buf_16u) || ((indices[i] != 65535) && is_buf_16u)); i++ )\r
- {\r
- int idx = indices[i];\r
- values_buf[i] = *(responses->data.fl + idx * r_step);\r
- }\r
- \r
- *values = values_buf; \r
-}\r
-\r
-\r
-void CvDTreeTrainData::get_cv_labels( CvDTreeNode* n, int* labels_buf, const int** labels )\r
-{\r
- if (have_labels)\r
- get_cat_var_data( n, get_work_var_count()- 1, labels_buf, labels );\r
-}\r
-\r
-\r
-int CvDTreeTrainData::get_cat_var_data( CvDTreeNode* n, int vi, int* cat_values_buf, const int** cat_values )\r
-{\r
- if( !is_buf_16u )\r
- *cat_values = buf->data.i + n->buf_idx*buf->cols + \r
- vi*sample_count + n->offset;\r
- else {\r
- const unsigned short* short_values = (const unsigned short*)(buf->data.s + n->buf_idx*buf->cols + \r
- vi*sample_count + n->offset);\r
- for( int i = 0; i < n->sample_count; i++ )\r
- cat_values_buf[i] = short_values[i];\r
- *cat_values = cat_values_buf;\r
- }\r
-\r
- return 0; //TODO: return the number of non-missing values\r
-}\r
-\r
-\r
-int CvDTreeTrainData::get_child_buf_idx( CvDTreeNode* n )\r
-{\r
- int idx = n->buf_idx + 1;\r
- if( idx >= buf_count )\r
- idx = shared ? 1 : 0;\r
- return idx;\r
-}\r
-\r
-\r
-void CvDTreeTrainData::write_params( CvFileStorage* fs )\r
-{\r
- CV_FUNCNAME( "CvDTreeTrainData::write_params" );\r
-\r
- __BEGIN__;\r
-\r
- int vi, vcount = var_count;\r
-\r
- cvWriteInt( fs, "is_classifier", is_classifier ? 1 : 0 );\r
- cvWriteInt( fs, "var_all", var_all );\r
- cvWriteInt( fs, "var_count", var_count );\r
- cvWriteInt( fs, "ord_var_count", ord_var_count );\r
- cvWriteInt( fs, "cat_var_count", cat_var_count );\r
-\r
- cvStartWriteStruct( fs, "training_params", CV_NODE_MAP );\r
- cvWriteInt( fs, "use_surrogates", params.use_surrogates ? 1 : 0 );\r
-\r
- if( is_classifier )\r
- {\r
- cvWriteInt( fs, "max_categories", params.max_categories );\r
- }\r
- else\r
- {\r
- cvWriteReal( fs, "regression_accuracy", params.regression_accuracy );\r
- }\r
-\r
- cvWriteInt( fs, "max_depth", params.max_depth );\r
- cvWriteInt( fs, "min_sample_count", params.min_sample_count );\r
- cvWriteInt( fs, "cross_validation_folds", params.cv_folds );\r
-\r
- if( params.cv_folds > 1 )\r
- {\r
- cvWriteInt( fs, "use_1se_rule", params.use_1se_rule ? 1 : 0 );\r
- cvWriteInt( fs, "truncate_pruned_tree", params.truncate_pruned_tree ? 1 : 0 );\r
- }\r
-\r
- if( priors )\r
- cvWrite( fs, "priors", priors );\r
-\r
- cvEndWriteStruct( fs );\r
-\r
- if( var_idx )\r
- cvWrite( fs, "var_idx", var_idx );\r
-\r
- cvStartWriteStruct( fs, "var_type", CV_NODE_SEQ+CV_NODE_FLOW );\r
-\r
- for( vi = 0; vi < vcount; vi++ )\r
- cvWriteInt( fs, 0, var_type->data.i[vi] >= 0 );\r
-\r
- cvEndWriteStruct( fs );\r
-\r
- if( cat_count && (cat_var_count > 0 || is_classifier) )\r
- {\r
- CV_ASSERT( cat_count != 0 );\r
- cvWrite( fs, "cat_count", cat_count );\r
- cvWrite( fs, "cat_map", cat_map );\r
- }\r
-\r
- __END__;\r
-}\r
-\r
-\r
-void CvDTreeTrainData::read_params( CvFileStorage* fs, CvFileNode* node )\r
-{\r
- CV_FUNCNAME( "CvDTreeTrainData::read_params" );\r
-\r
- __BEGIN__;\r
-\r
- CvFileNode *tparams_node, *vartype_node;\r
- CvSeqReader reader;\r
- int vi, max_split_size, tree_block_size;\r
-\r
- is_classifier = (cvReadIntByName( fs, node, "is_classifier" ) != 0);\r
- var_all = cvReadIntByName( fs, node, "var_all" );\r
- var_count = cvReadIntByName( fs, node, "var_count", var_all );\r
- cat_var_count = cvReadIntByName( fs, node, "cat_var_count" );\r
- ord_var_count = cvReadIntByName( fs, node, "ord_var_count" );\r
-\r
- tparams_node = cvGetFileNodeByName( fs, node, "training_params" );\r
-\r
- if( tparams_node ) // training parameters are not necessary\r
- {\r
- params.use_surrogates = cvReadIntByName( fs, tparams_node, "use_surrogates", 1 ) != 0;\r
-\r
- if( is_classifier )\r
- {\r
- params.max_categories = cvReadIntByName( fs, tparams_node, "max_categories" );\r
- }\r
- else\r
- {\r
- params.regression_accuracy =\r
- (float)cvReadRealByName( fs, tparams_node, "regression_accuracy" );\r
- }\r
-\r
- params.max_depth = cvReadIntByName( fs, tparams_node, "max_depth" );\r
- params.min_sample_count = cvReadIntByName( fs, tparams_node, "min_sample_count" );\r
- params.cv_folds = cvReadIntByName( fs, tparams_node, "cross_validation_folds" );\r
-\r
- if( params.cv_folds > 1 )\r
- {\r
- params.use_1se_rule = cvReadIntByName( fs, tparams_node, "use_1se_rule" ) != 0;\r
- params.truncate_pruned_tree =\r
- cvReadIntByName( fs, tparams_node, "truncate_pruned_tree" ) != 0;\r
- }\r
-\r
- priors = (CvMat*)cvReadByName( fs, tparams_node, "priors" );\r
- if( priors )\r
- {\r
- if( !CV_IS_MAT(priors) )\r
- CV_ERROR( CV_StsParseError, "priors must stored as a matrix" );\r
- priors_mult = cvCloneMat( priors );\r
- }\r
- }\r
-\r
- CV_CALL( var_idx = (CvMat*)cvReadByName( fs, node, "var_idx" ));\r
- if( var_idx )\r
- {\r
- if( !CV_IS_MAT(var_idx) ||\r
- (var_idx->cols != 1 && var_idx->rows != 1) ||\r
- var_idx->cols + var_idx->rows - 1 != var_count ||\r
- CV_MAT_TYPE(var_idx->type) != CV_32SC1 )\r
- CV_ERROR( CV_StsParseError,\r
- "var_idx (if exist) must be valid 1d integer vector containing <var_count> elements" );\r
-\r
- for( vi = 0; vi < var_count; vi++ )\r
- if( (unsigned)var_idx->data.i[vi] >= (unsigned)var_all )\r
- CV_ERROR( CV_StsOutOfRange, "some of var_idx elements are out of range" );\r
- }\r
-\r
- ////// read var type\r
- CV_CALL( var_type = cvCreateMat( 1, var_count + 2, CV_32SC1 ));\r
-\r
- cat_var_count = 0;\r
- ord_var_count = -1;\r
- vartype_node = cvGetFileNodeByName( fs, node, "var_type" );\r
-\r
- if( vartype_node && CV_NODE_TYPE(vartype_node->tag) == CV_NODE_INT && var_count == 1 )\r
- var_type->data.i[0] = vartype_node->data.i ? cat_var_count++ : ord_var_count--;\r
- else\r
- {\r
- if( !vartype_node || CV_NODE_TYPE(vartype_node->tag) != CV_NODE_SEQ ||\r
- vartype_node->data.seq->total != var_count )\r
- CV_ERROR( CV_StsParseError, "var_type must exist and be a sequence of 0's and 1's" );\r
-\r
- cvStartReadSeq( vartype_node->data.seq, &reader );\r
-\r
- for( vi = 0; vi < var_count; vi++ )\r
- {\r
- CvFileNode* n = (CvFileNode*)reader.ptr;\r
- if( CV_NODE_TYPE(n->tag) != CV_NODE_INT || (n->data.i & ~1) )\r
- CV_ERROR( CV_StsParseError, "var_type must exist and be a sequence of 0's and 1's" );\r
- var_type->data.i[vi] = n->data.i ? cat_var_count++ : ord_var_count--;\r
- CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );\r
- }\r
- }\r
- var_type->data.i[var_count] = cat_var_count;\r
-\r
- ord_var_count = ~ord_var_count;\r
- if( cat_var_count != cat_var_count || ord_var_count != ord_var_count )\r
- CV_ERROR( CV_StsParseError, "var_type is inconsistent with cat_var_count and ord_var_count" );\r
- //////\r
-\r
- if( cat_var_count > 0 || is_classifier )\r
- {\r
- int ccount, total_c_count = 0;\r
- CV_CALL( cat_count = (CvMat*)cvReadByName( fs, node, "cat_count" ));\r
- CV_CALL( cat_map = (CvMat*)cvReadByName( fs, node, "cat_map" ));\r
-\r
- if( !CV_IS_MAT(cat_count) || !CV_IS_MAT(cat_map) ||\r
- (cat_count->cols != 1 && cat_count->rows != 1) ||\r
- CV_MAT_TYPE(cat_count->type) != CV_32SC1 ||\r
- cat_count->cols + cat_count->rows - 1 != cat_var_count + is_classifier ||\r
- (cat_map->cols != 1 && cat_map->rows != 1) ||\r
- CV_MAT_TYPE(cat_map->type) != CV_32SC1 )\r
- CV_ERROR( CV_StsParseError,\r
- "Both cat_count and cat_map must exist and be valid 1d integer vectors of an appropriate size" );\r
-\r
- ccount = cat_var_count + is_classifier;\r
-\r
- CV_CALL( cat_ofs = cvCreateMat( 1, ccount + 1, CV_32SC1 ));\r
- cat_ofs->data.i[0] = 0;\r
- max_c_count = 1;\r
-\r
- for( vi = 0; vi < ccount; vi++ )\r
- {\r
- int val = cat_count->data.i[vi];\r
- if( val <= 0 )\r
- CV_ERROR( CV_StsOutOfRange, "some of cat_count elements are out of range" );\r
- max_c_count = MAX( max_c_count, val );\r
- cat_ofs->data.i[vi+1] = total_c_count += val;\r
- }\r
-\r
- if( cat_map->cols + cat_map->rows - 1 != total_c_count )\r
- CV_ERROR( CV_StsBadSize,\r
- "cat_map vector length is not equal to the total number of categories in all categorical vars" );\r
- }\r
-\r
- max_split_size = cvAlign(sizeof(CvDTreeSplit) +\r
- (MAX(0,max_c_count - 33)/32)*sizeof(int),sizeof(void*));\r
-\r
- tree_block_size = MAX((int)sizeof(CvDTreeNode)*8, max_split_size);\r
- tree_block_size = MAX(tree_block_size + block_size_delta, min_block_size);\r
- CV_CALL( tree_storage = cvCreateMemStorage( tree_block_size ));\r
- CV_CALL( node_heap = cvCreateSet( 0, sizeof(node_heap[0]),\r
- sizeof(CvDTreeNode), tree_storage ));\r
- CV_CALL( split_heap = cvCreateSet( 0, sizeof(split_heap[0]),\r
- max_split_size, tree_storage ));\r
-\r
- __END__;\r
-}\r
-\r
-float* CvDTreeTrainData::get_pred_float_buf()\r
-{\r
- return &pred_float_buf[cv::getThreadNum()][0];\r
-}\r
-int* CvDTreeTrainData::get_pred_int_buf()\r
-{\r
- return &pred_int_buf[cv::getThreadNum()][0];\r
-}\r
-float* CvDTreeTrainData::get_resp_float_buf()\r
-{\r
- return &resp_float_buf[cv::getThreadNum()][0];\r
-}\r
-int* CvDTreeTrainData::get_resp_int_buf()\r
-{\r
- return &resp_int_buf[cv::getThreadNum()][0];\r
-}\r
-int* CvDTreeTrainData::get_cv_lables_buf()\r
-{\r
- return &cv_lables_buf[cv::getThreadNum()][0];\r
-}\r
-int* CvDTreeTrainData::get_sample_idx_buf()\r
-{\r
- return &sample_idx_buf[cv::getThreadNum()][0];\r
-}\r
-\r
-/////////////////////// Decision Tree /////////////////////////\r
-\r
-CvDTree::CvDTree()\r
-{\r
- data = 0;\r
- var_importance = 0;\r
- default_model_name = "my_tree";\r
-\r
- clear();\r
-}\r
-\r
-\r
-void CvDTree::clear()\r
-{\r
- cvReleaseMat( &var_importance );\r
- if( data )\r
- {\r
- if( !data->shared )\r
- delete data;\r
- else\r
- free_tree();\r
- data = 0;\r
- }\r
- root = 0;\r
- pruned_tree_idx = -1;\r
-}\r
-\r
-\r
-CvDTree::~CvDTree()\r
-{\r
- clear();\r
-}\r
-\r
-\r
-const CvDTreeNode* CvDTree::get_root() const\r
-{\r
- return root;\r
-}\r
-\r
-\r
-int CvDTree::get_pruned_tree_idx() const\r
-{\r
- return pruned_tree_idx;\r
-}\r
-\r
-\r
-CvDTreeTrainData* CvDTree::get_data()\r
-{\r
- return data;\r
-}\r
-\r
-\r
-bool CvDTree::train( const CvMat* _train_data, int _tflag,\r
- const CvMat* _responses, const CvMat* _var_idx,\r
- const CvMat* _sample_idx, const CvMat* _var_type,\r
- const CvMat* _missing_mask, CvDTreeParams _params )\r
-{\r
- bool result = false;\r
-\r
- CV_FUNCNAME( "CvDTree::train" );\r
-\r
- __BEGIN__;\r
-\r
- clear();\r
- data = new CvDTreeTrainData( _train_data, _tflag, _responses,\r
- _var_idx, _sample_idx, _var_type,\r
- _missing_mask, _params, false );\r
- CV_CALL( result = do_train(0) );\r
-\r
- __END__;\r
-\r
- return result;\r
-}\r
-\r
-bool CvDTree::train( CvMLData* _data, CvDTreeParams _params )\r
-{\r
- bool result = false;\r
-\r
- CV_FUNCNAME( "CvDTree::train" );\r
-\r
- __BEGIN__;\r
-\r
- const CvMat* values = _data->get_values();\r
- const CvMat* response = _data->get_response();\r
- const CvMat* missing = _data->get_missing();\r
- const CvMat* var_types = _data->get_var_types();\r
- const CvMat* train_sidx = _data->get_train_sample_idx();\r
- const CvMat* var_idx = _data->get_var_idx();\r
-\r
- CV_CALL( result = train( values, CV_ROW_SAMPLE, response, var_idx,\r
- train_sidx, var_types, missing, _params ) );\r
-\r
- __END__;\r
-\r
- return result;\r
-}\r
-\r
-bool CvDTree::train( CvDTreeTrainData* _data, const CvMat* _subsample_idx )\r
-{\r
- bool result = false;\r
-\r
- CV_FUNCNAME( "CvDTree::train" );\r
-\r
- __BEGIN__;\r
-\r
- clear();\r
- data = _data;\r
- data->shared = true;\r
- CV_CALL( result = do_train(_subsample_idx));\r
-\r
- __END__;\r
-\r
- return result;\r
-}\r
-\r
-\r
-bool CvDTree::do_train( const CvMat* _subsample_idx )\r
-{\r
- bool result = false;\r
-\r
- CV_FUNCNAME( "CvDTree::do_train" );\r
-\r
- __BEGIN__;\r
-\r
- root = data->subsample_data( _subsample_idx );\r
-\r
- CV_CALL( try_split_node(root));\r
-\r
- if( data->params.cv_folds > 0 )\r
- CV_CALL( prune_cv());\r
-\r
- if( !data->shared )\r
- data->free_train_data();\r
-\r
- result = true;\r
-\r
- __END__;\r
-\r
- return result;\r
-}\r
-\r
-\r
-void CvDTree::try_split_node( CvDTreeNode* node )\r
-{\r
- CvDTreeSplit* best_split = 0;\r
- int i, n = node->sample_count, vi;\r
- bool can_split = true;\r
- double quality_scale;\r
-\r
- calc_node_value( node );\r
-\r
- if( node->sample_count <= data->params.min_sample_count ||\r
- node->depth >= data->params.max_depth )\r
- can_split = false;\r
-\r
- if( can_split && data->is_classifier )\r
- {\r
- // check if we have a "pure" node,\r
- // we assume that cls_count is filled by calc_node_value()\r
- int* cls_count = data->counts->data.i;\r
- int nz = 0, m = data->get_num_classes();\r
- for( i = 0; i < m; i++ )\r
- nz += cls_count[i] != 0;\r
- if( nz == 1 ) // there is only one class\r
- can_split = false;\r
- }\r
- else if( can_split )\r
- {\r
- if( sqrt(node->node_risk)/n < data->params.regression_accuracy )\r
- can_split = false;\r
- }\r
-\r
- if( can_split )\r
- {\r
- best_split = find_best_split(node);\r
- // TODO: check the split quality ...\r
- node->split = best_split;\r
- }\r
- if( !can_split || !best_split )\r
- {\r
- data->free_node_data(node);\r
- return;\r
- }\r
-\r
- quality_scale = calc_node_dir( node );\r
- if( data->params.use_surrogates )\r
- {\r
- // find all the surrogate splits\r
- // and sort them by their similarity to the primary one\r
- for( vi = 0; vi < data->var_count; vi++ )\r
- {\r
- CvDTreeSplit* split;\r
- int ci = data->get_var_type(vi);\r
-\r
- if( vi == best_split->var_idx )\r
- continue;\r
-\r
- if( ci >= 0 )\r
- split = find_surrogate_split_cat( node, vi );\r
- else\r
- split = find_surrogate_split_ord( node, vi );\r
-\r
- if( split )\r
- {\r
- // insert the split\r
- CvDTreeSplit* prev_split = node->split;\r
- split->quality = (float)(split->quality*quality_scale);\r
-\r
- while( prev_split->next &&\r
- prev_split->next->quality > split->quality )\r
- prev_split = prev_split->next;\r
- split->next = prev_split->next;\r
- prev_split->next = split;\r
- }\r
- }\r
- }\r
- split_node_data( node );\r
- try_split_node( node->left );\r
- try_split_node( node->right );\r
-}\r
-\r
-\r
-// calculate direction (left(-1),right(1),missing(0))\r
-// for each sample using the best split\r
-// the function returns scale coefficients for surrogate split quality factors.\r
-// the scale is applied to normalize surrogate split quality relatively to the\r
-// best (primary) split quality. That is, if a surrogate split is absolutely\r
-// identical to the primary split, its quality will be set to the maximum value =\r
-// quality of the primary split; otherwise, it will be lower.\r
-// besides, the function compute node->maxlr,\r
-// minimum possible quality (w/o considering the above mentioned scale)\r
-// for a surrogate split. Surrogate splits with quality less than node->maxlr\r
-// are not discarded.\r
-double CvDTree::calc_node_dir( CvDTreeNode* node )\r
-{\r
- char* dir = (char*)data->direction->data.ptr;\r
- int i, n = node->sample_count, vi = node->split->var_idx;\r
- double L, R;\r
-\r
- assert( !node->split->inversed );\r
-\r
- if( data->get_var_type(vi) >= 0 ) // split on categorical var\r
- {\r
- int* labels_buf = data->get_pred_int_buf();\r
- const int* labels = 0;\r
- const int* subset = node->split->subset;\r
- data->get_cat_var_data( node, vi, labels_buf, &labels );\r
- if( !data->have_priors )\r
- {\r
- int sum = 0, sum_abs = 0;\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- int idx = labels[i];\r
- int d = ( ((idx >= 0)&&(!data->is_buf_16u)) || ((idx != 65535)&&(data->is_buf_16u)) ) ?\r
- CV_DTREE_CAT_DIR(idx,subset) : 0;\r
- sum += d; sum_abs += d & 1;\r
- dir[i] = (char)d;\r
- }\r
-\r
- R = (sum_abs + sum) >> 1;\r
- L = (sum_abs - sum) >> 1;\r
- }\r
- else\r
- {\r
- const double* priors = data->priors_mult->data.db;\r
- double sum = 0, sum_abs = 0;\r
- int *responses_buf = data->get_resp_int_buf();\r
- const int* responses;\r
- data->get_class_labels(node, responses_buf, &responses);\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- int idx = labels[i];\r
- double w = priors[responses[i]];\r
- int d = idx >= 0 ? CV_DTREE_CAT_DIR(idx,subset) : 0;\r
- sum += d*w; sum_abs += (d & 1)*w;\r
- dir[i] = (char)d;\r
- }\r
-\r
- R = (sum_abs + sum) * 0.5;\r
- L = (sum_abs - sum) * 0.5;\r
- }\r
- }\r
- else // split on ordered var\r
- {\r
- int split_point = node->split->ord.split_point;\r
- int n1 = node->get_num_valid(vi);\r
- \r
- float* val_buf = data->get_pred_float_buf();\r
- const float* val = 0;\r
- int* sorted_buf = data->get_pred_int_buf();\r
- const int* sorted = 0;\r
- data->get_ord_var_data( node, vi, val_buf, sorted_buf, &val, &sorted);\r
- \r
- assert( 0 <= split_point && split_point < n1-1 );\r
-\r
- if( !data->have_priors )\r
- {\r
- for( i = 0; i <= split_point; i++ )\r
- dir[sorted[i]] = (char)-1;\r
- for( ; i < n1; i++ )\r
- dir[sorted[i]] = (char)1;\r
- for( ; i < n; i++ )\r
- dir[sorted[i]] = (char)0;\r
-\r
- L = split_point-1;\r
- R = n1 - split_point + 1;\r
- }\r
- else\r
- {\r
- const double* priors = data->priors_mult->data.db;\r
- int* responses_buf = data->get_resp_int_buf();\r
- const int* responses = 0;\r
- data->get_class_labels(node, responses_buf, &responses);\r
- L = R = 0;\r
-\r
- for( i = 0; i <= split_point; i++ )\r
- {\r
- int idx = sorted[i];\r
- double w = priors[responses[idx]];\r
- dir[idx] = (char)-1;\r
- L += w;\r
- }\r
-\r
- for( ; i < n1; i++ )\r
- {\r
- int idx = sorted[i];\r
- double w = priors[responses[idx]];\r
- dir[idx] = (char)1;\r
- R += w;\r
- }\r
-\r
- for( ; i < n; i++ )\r
- dir[sorted[i]] = (char)0;\r
- }\r
- }\r
- node->maxlr = MAX( L, R );\r
- return node->split->quality/(L + R);\r
-}\r
-\r
-CvDTreeSplit* CvDTree::find_best_split( CvDTreeNode* node )\r
-{\r
- int vi;\r
- CvDTreeSplit *bestSplit = 0;\r
- int maxNumThreads = 1;\r
-#ifdef _OPENMP\r
- maxNumThreads = cv::getNumThreads();\r
-#endif\r
- vector<CvDTreeSplit*> splits(maxNumThreads);\r
- vector<CvDTreeSplit*> bestSplits(maxNumThreads);\r
- for (int i = 0; i < maxNumThreads; i++)\r
- {\r
- splits[i] = data->new_split_cat( 0, -1.0f );\r
- bestSplits[i] = data->new_split_cat( 0, -1.0f );\r
- }\r
-\r
- bool can_split = false;\r
-#ifdef _OPENMP\r
-#pragma omp parallel for num_threads(maxNumThreads) schedule(dynamic)\r
-#endif\r
- for( vi = 0; vi < data->var_count; vi++ )\r
- {\r
- CvDTreeSplit *res, *t;\r
- int threadIdx = cv::getThreadNum();\r
- int ci = data->get_var_type(vi);\r
- if( node->get_num_valid(vi) <= 1 )\r
- continue;\r
-\r
- if( data->is_classifier )\r
- {\r
- if( ci >= 0 )\r
- res = find_split_cat_class( node, vi, bestSplits[threadIdx]->quality, splits[threadIdx] );\r
- else\r
- res = find_split_ord_class( node, vi, bestSplits[threadIdx]->quality, splits[threadIdx] );\r
- }\r
- else\r
- {\r
- if( ci >= 0 )\r
- res = find_split_cat_reg( node, vi, bestSplits[threadIdx]->quality, splits[threadIdx] );\r
- else\r
- res = find_split_ord_reg( node, vi, bestSplits[threadIdx]->quality, splits[threadIdx] );\r
- }\r
-\r
- if( res )\r
- {\r
- can_split = true;\r
- if( bestSplits[threadIdx]->quality < splits[threadIdx]->quality )\r
- CV_SWAP( bestSplits[threadIdx], splits[threadIdx], t );\r
- }\r
- }\r
- if ( can_split )\r
- {\r
- bestSplit = bestSplits[0];\r
- for(int i = 1; i < maxNumThreads; i++)\r
- {\r
- if( bestSplit->quality < bestSplits[i]->quality )\r
- bestSplit = bestSplits[i];\r
- }\r
- }\r
- for(int i = 0; i < maxNumThreads; i++)\r
- {\r
- cvSetRemoveByPtr( data->split_heap, splits[i] );\r
- if( bestSplits[i] != bestSplit )\r
- cvSetRemoveByPtr( data->split_heap, bestSplits[i] );\r
- }\r
- return bestSplit;\r
-}\r
-\r
-CvDTreeSplit* CvDTree::find_split_ord_class( CvDTreeNode* node, int vi,\r
- float init_quality, CvDTreeSplit* _split )\r
-{\r
- const float epsilon = FLT_EPSILON*2;\r
- int n = node->sample_count;\r
- int n1 = node->get_num_valid(vi);\r
- int m = data->get_num_classes();\r
-\r
- float* values_buf = data->get_pred_float_buf();\r
- const float* values = 0;\r
- int* indices_buf = data->get_pred_int_buf();\r
- const int* indices = 0;\r
- data->get_ord_var_data( node, vi, values_buf, indices_buf, &values, &indices );\r
- int* responses_buf = data->get_resp_int_buf();\r
- const int* responses = 0;\r
- data->get_class_labels( node, responses_buf, &responses );\r
-\r
- const int* rc0 = data->counts->data.i;\r
- int* lc = (int*)cvStackAlloc(m*sizeof(lc[0]));\r
- int* rc = (int*)cvStackAlloc(m*sizeof(rc[0]));\r
- int i, best_i = -1;\r
- double lsum2 = 0, rsum2 = 0, best_val = init_quality;\r
- const double* priors = data->have_priors ? data->priors_mult->data.db : 0;\r
-\r
- // init arrays of class instance counters on both sides of the split\r
- for( i = 0; i < m; i++ )\r
- {\r
- lc[i] = 0;\r
- rc[i] = rc0[i];\r
- }\r
-\r
- // compensate for missing values\r
- for( i = n1; i < n; i++ )\r
- {\r
- rc[responses[indices[i]]]--;\r
- }\r
-\r
- if( !priors )\r
- {\r
- int L = 0, R = n1;\r
-\r
- for( i = 0; i < m; i++ )\r
- rsum2 += (double)rc[i]*rc[i];\r
-\r
- for( i = 0; i < n1 - 1; i++ )\r
- {\r
- int idx = responses[indices[i]];\r
- int lv, rv;\r
- L++; R--;\r
- lv = lc[idx]; rv = rc[idx];\r
- lsum2 += lv*2 + 1;\r
- rsum2 -= rv*2 - 1;\r
- lc[idx] = lv + 1; rc[idx] = rv - 1;\r
-\r
- if( values[i] + epsilon < values[i+1] )\r
- {\r
- double val = (lsum2*R + rsum2*L)/((double)L*R);\r
- if( best_val < val )\r
- {\r
- best_val = val;\r
- best_i = i;\r
- }\r
- }\r
- }\r
- }\r
- else\r
- {\r
- double L = 0, R = 0;\r
- for( i = 0; i < m; i++ )\r
- {\r
- double wv = rc[i]*priors[i];\r
- R += wv;\r
- rsum2 += wv*wv;\r
- }\r
-\r
- for( i = 0; i < n1 - 1; i++ )\r
- {\r
- int idx = responses[indices[i]];\r
- int lv, rv;\r
- double p = priors[idx], p2 = p*p;\r
- L += p; R -= p;\r
- lv = lc[idx]; rv = rc[idx];\r
- lsum2 += p2*(lv*2 + 1);\r
- rsum2 -= p2*(rv*2 - 1);\r
- lc[idx] = lv + 1; rc[idx] = rv - 1;\r
-\r
- if( values[i] + epsilon < values[i+1] )\r
- {\r
- double val = (lsum2*R + rsum2*L)/((double)L*R);\r
- if( best_val < val )\r
- {\r
- best_val = val;\r
- best_i = i;\r
- }\r
- }\r
- }\r
- }\r
-\r
- CvDTreeSplit* split = 0;\r
- if( best_i >= 0 )\r
- {\r
- split = _split ? _split : data->new_split_ord( 0, 0.0f, 0, 0, 0.0f );\r
- split->var_idx = vi;\r
- split->ord.c = (values[best_i] + values[best_i+1])*0.5f;\r
- split->ord.split_point = best_i;\r
- split->inversed = 0;\r
- split->quality = (float)best_val;\r
- }\r
- return split;\r
-}\r
-\r
-\r
-void CvDTree::cluster_categories( const int* vectors, int n, int m,\r
- int* csums, int k, int* labels )\r
-{\r
- // TODO: consider adding priors (class weights) and sample weights to the clustering algorithm\r
- int iters = 0, max_iters = 100;\r
- int i, j, idx;\r
- double* buf = (double*)cvStackAlloc( (n + k)*sizeof(buf[0]) );\r
- double *v_weights = buf, *c_weights = buf + k;\r
- bool modified = true;\r
- CvRNG* r = &data->rng;\r
-\r
- // assign labels randomly\r
- for( i = idx = 0; i < n; i++ )\r
- {\r
- int sum = 0;\r
- const int* v = vectors + i*m;\r
- labels[i] = idx++;\r
- idx &= idx < k ? -1 : 0;\r
-\r
- // compute weight of each vector\r
- for( j = 0; j < m; j++ )\r
- sum += v[j];\r
- v_weights[i] = sum ? 1./sum : 0.;\r
- }\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- int i1 = cvRandInt(r) % n;\r
- int i2 = cvRandInt(r) % n;\r
- CV_SWAP( labels[i1], labels[i2], j );\r
- }\r
-\r
- for( iters = 0; iters <= max_iters; iters++ )\r
- {\r
- // calculate csums\r
- for( i = 0; i < k; i++ )\r
- {\r
- for( j = 0; j < m; j++ )\r
- csums[i*m + j] = 0;\r
- }\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- const int* v = vectors + i*m;\r
- int* s = csums + labels[i]*m;\r
- for( j = 0; j < m; j++ )\r
- s[j] += v[j];\r
- }\r
-\r
- // exit the loop here, when we have up-to-date csums\r
- if( iters == max_iters || !modified )\r
- break;\r
-\r
- modified = false;\r
-\r
- // calculate weight of each cluster\r
- for( i = 0; i < k; i++ )\r
- {\r
- const int* s = csums + i*m;\r
- int sum = 0;\r
- for( j = 0; j < m; j++ )\r
- sum += s[j];\r
- c_weights[i] = sum ? 1./sum : 0;\r
- }\r
-\r
- // now for each vector determine the closest cluster\r
- for( i = 0; i < n; i++ )\r
- {\r
- const int* v = vectors + i*m;\r
- double alpha = v_weights[i];\r
- double min_dist2 = DBL_MAX;\r
- int min_idx = -1;\r
-\r
- for( idx = 0; idx < k; idx++ )\r
- {\r
- const int* s = csums + idx*m;\r
- double dist2 = 0., beta = c_weights[idx];\r
- for( j = 0; j < m; j++ )\r
- {\r
- double t = v[j]*alpha - s[j]*beta;\r
- dist2 += t*t;\r
- }\r
- if( min_dist2 > dist2 )\r
- {\r
- min_dist2 = dist2;\r
- min_idx = idx;\r
- }\r
- }\r
-\r
- if( min_idx != labels[i] )\r
- modified = true;\r
- labels[i] = min_idx;\r
- }\r
- }\r
-}\r
-\r
-\r
-CvDTreeSplit* CvDTree::find_split_cat_class( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split )\r
-{\r
- int ci = data->get_var_type(vi);\r
- int n = node->sample_count;\r
- int m = data->get_num_classes();\r
- int _mi = data->cat_count->data.i[ci], mi = _mi;\r
-\r
- int* labels_buf = data->get_pred_int_buf();\r
- const int* labels = 0;\r
- data->get_cat_var_data(node, vi, labels_buf, &labels);\r
- int *responses_buf = data->get_resp_int_buf();\r
- const int* responses = 0;\r
- data->get_class_labels(node, responses_buf, &responses);\r
-\r
- int* lc = (int*)cvStackAlloc(m*sizeof(lc[0]));\r
- int* rc = (int*)cvStackAlloc(m*sizeof(rc[0]));\r
- int* _cjk = (int*)cvStackAlloc(m*(mi+1)*sizeof(_cjk[0]))+m, *cjk = _cjk;\r
- double* c_weights = (double*)cvStackAlloc( mi*sizeof(c_weights[0]) );\r
- int* cluster_labels = 0;\r
- int** int_ptr = 0;\r
- int i, j, k, idx;\r
- double L = 0, R = 0;\r
- double best_val = init_quality;\r
- int prevcode = 0, best_subset = -1, subset_i, subset_n, subtract = 0;\r
- const double* priors = data->priors_mult->data.db;\r
-\r
- // init array of counters:\r
- // c_{jk} - number of samples that have vi-th input variable = j and response = k.\r
- for( j = -1; j < mi; j++ )\r
- for( k = 0; k < m; k++ )\r
- cjk[j*m + k] = 0;\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- j = ( labels[i] == 65535 && data->is_buf_16u) ? -1 : labels[i];\r
- k = responses[i];\r
- cjk[j*m + k]++;\r
- }\r
-\r
- if( m > 2 )\r
- {\r
- if( mi > data->params.max_categories )\r
- {\r
- mi = MIN(data->params.max_categories, n);\r
- cjk += _mi*m;\r
- cluster_labels = (int*)cvStackAlloc(mi*sizeof(cluster_labels[0]));\r
- cluster_categories( _cjk, _mi, m, cjk, mi, cluster_labels );\r
- }\r
- subset_i = 1;\r
- subset_n = 1 << mi;\r
- }\r
- else\r
- {\r
- assert( m == 2 );\r
- int_ptr = (int**)cvStackAlloc( mi*sizeof(int_ptr[0]) );\r
- for( j = 0; j < mi; j++ )\r
- int_ptr[j] = cjk + j*2 + 1;\r
- icvSortIntPtr( int_ptr, mi, 0 );\r
- subset_i = 0;\r
- subset_n = mi;\r
- }\r
-\r
- for( k = 0; k < m; k++ )\r
- {\r
- int sum = 0;\r
- for( j = 0; j < mi; j++ )\r
- sum += cjk[j*m + k];\r
- rc[k] = sum;\r
- lc[k] = 0;\r
- }\r
-\r
- for( j = 0; j < mi; j++ )\r
- {\r
- double sum = 0;\r
- for( k = 0; k < m; k++ )\r
- sum += cjk[j*m + k]*priors[k];\r
- c_weights[j] = sum;\r
- R += c_weights[j];\r
- }\r
-\r
- for( ; subset_i < subset_n; subset_i++ )\r
- {\r
- double weight;\r
- int* crow;\r
- double lsum2 = 0, rsum2 = 0;\r
-\r
- if( m == 2 )\r
- idx = (int)(int_ptr[subset_i] - cjk)/2;\r
- else\r
- {\r
- int graycode = (subset_i>>1)^subset_i;\r
- int diff = graycode ^ prevcode;\r
-\r
- // determine index of the changed bit.\r
- Cv32suf u;\r
- idx = diff >= (1 << 16) ? 16 : 0;\r
- u.f = (float)(((diff >> 16) | diff) & 65535);\r
- idx += (u.i >> 23) - 127;\r
- subtract = graycode < prevcode;\r
- prevcode = graycode;\r
- }\r
-\r
- crow = cjk + idx*m;\r
- weight = c_weights[idx];\r
- if( weight < FLT_EPSILON )\r
- continue;\r
-\r
- if( !subtract )\r
- {\r
- for( k = 0; k < m; k++ )\r
- {\r
- int t = crow[k];\r
- int lval = lc[k] + t;\r
- int rval = rc[k] - t;\r
- double p = priors[k], p2 = p*p;\r
- lsum2 += p2*lval*lval;\r
- rsum2 += p2*rval*rval;\r
- lc[k] = lval; rc[k] = rval;\r
- }\r
- L += weight;\r
- R -= weight;\r
- }\r
- else\r
- {\r
- for( k = 0; k < m; k++ )\r
- {\r
- int t = crow[k];\r
- int lval = lc[k] - t;\r
- int rval = rc[k] + t;\r
- double p = priors[k], p2 = p*p;\r
- lsum2 += p2*lval*lval;\r
- rsum2 += p2*rval*rval;\r
- lc[k] = lval; rc[k] = rval;\r
- }\r
- L -= weight;\r
- R += weight;\r
- }\r
-\r
- if( L > FLT_EPSILON && R > FLT_EPSILON )\r
- {\r
- double val = (lsum2*R + rsum2*L)/((double)L*R);\r
- if( best_val < val )\r
- {\r
- best_val = val;\r
- best_subset = subset_i;\r
- }\r
- }\r
- }\r
-\r
- CvDTreeSplit* split = 0;\r
- if( best_subset >= 0 ) \r
- {\r
- split = _split ? _split : data->new_split_cat( 0, -1.0f );\r
- split->var_idx = vi;\r
- split->quality = (float)best_val;\r
- memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));\r
- if( m == 2 )\r
- {\r
- for( i = 0; i <= best_subset; i++ )\r
- {\r
- idx = (int)(int_ptr[i] - cjk) >> 1;\r
- split->subset[idx >> 5] |= 1 << (idx & 31);\r
- }\r
- }\r
- else\r
- {\r
- for( i = 0; i < _mi; i++ )\r
- {\r
- idx = cluster_labels ? cluster_labels[i] : i;\r
- if( best_subset & (1 << idx) )\r
- split->subset[i >> 5] |= 1 << (i & 31);\r
- }\r
- }\r
- }\r
- return split;\r
-}\r
-\r
-\r
-CvDTreeSplit* CvDTree::find_split_ord_reg( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split )\r
-{\r
- const float epsilon = FLT_EPSILON*2;\r
- int n = node->sample_count;\r
- int n1 = node->get_num_valid(vi);\r
-\r
- float* values_buf = data->get_pred_float_buf();\r
- const float* values = 0;\r
- int* indices_buf = data->get_pred_int_buf();\r
- const int* indices = 0;\r
- data->get_ord_var_data( node, vi, values_buf, indices_buf, &values, &indices );\r
- float* responses_buf = data->get_resp_float_buf();\r
- const float* responses = 0;\r
- data->get_ord_responses( node, responses_buf, &responses );\r
-\r
- int i, best_i = -1;\r
- double best_val = init_quality, lsum = 0, rsum = node->value*n;\r
- int L = 0, R = n1;\r
-\r
- // compensate for missing values\r
- for( i = n1; i < n; i++ )\r
- rsum -= responses[indices[i]];\r
-\r
- // find the optimal split\r
- for( i = 0; i < n1 - 1; i++ )\r
- {\r
- float t = responses[indices[i]];\r
- L++; R--;\r
- lsum += t;\r
- rsum -= t;\r
-\r
- if( values[i] + epsilon < values[i+1] )\r
- {\r
- double val = (lsum*lsum*R + rsum*rsum*L)/((double)L*R);\r
- if( best_val < val )\r
- {\r
- best_val = val;\r
- best_i = i;\r
- }\r
- }\r
- }\r
-\r
- CvDTreeSplit* split = 0;\r
- if( best_i >= 0 )\r
- {\r
- split = _split ? _split : data->new_split_ord( 0, 0.0f, 0, 0, 0.0f );\r
- split->var_idx = vi;\r
- split->ord.c = (values[best_i] + values[best_i+1])*0.5f;\r
- split->ord.split_point = best_i;\r
- split->inversed = 0;\r
- split->quality = (float)best_val;\r
- }\r
- return split;\r
-}\r
-\r
-CvDTreeSplit* CvDTree::find_split_cat_reg( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split )\r
-{\r
- int ci = data->get_var_type(vi);\r
- int n = node->sample_count;\r
- int mi = data->cat_count->data.i[ci];\r
- int* labels_buf = data->get_pred_int_buf();\r
- const int* labels = 0;\r
- float* responses_buf = data->get_resp_float_buf();\r
- const float* responses = 0;\r
- data->get_cat_var_data(node, vi, labels_buf, &labels);\r
- data->get_ord_responses(node, responses_buf, &responses);\r
-\r
- double* sum = (double*)cvStackAlloc( (mi+1)*sizeof(sum[0]) ) + 1;\r
- int* counts = (int*)cvStackAlloc( (mi+1)*sizeof(counts[0]) ) + 1;\r
- double** sum_ptr = (double**)cvStackAlloc( (mi+1)*sizeof(sum_ptr[0]) );\r
- int i, L = 0, R = 0;\r
- double best_val = init_quality, lsum = 0, rsum = 0;\r
- int best_subset = -1, subset_i;\r
-\r
- for( i = -1; i < mi; i++ )\r
- sum[i] = counts[i] = 0;\r
-\r
- // calculate sum response and weight of each category of the input var\r
- for( i = 0; i < n; i++ )\r
- {\r
- int idx = ( (labels[i] == 65535) && data->is_buf_16u ) ? -1 : labels[i];\r
- double s = sum[idx] + responses[i];\r
- int nc = counts[idx] + 1;\r
- sum[idx] = s;\r
- counts[idx] = nc;\r
- }\r
-\r
- // calculate average response in each category\r
- for( i = 0; i < mi; i++ )\r
- {\r
- R += counts[i];\r
- rsum += sum[i];\r
- sum[i] /= MAX(counts[i],1);\r
- sum_ptr[i] = sum + i;\r
- }\r
-\r
- icvSortDblPtr( sum_ptr, mi, 0 );\r
-\r
- // revert back to unnormalized sums\r
- // (there should be a very little loss of accuracy)\r
- for( i = 0; i < mi; i++ )\r
- sum[i] *= counts[i];\r
-\r
- for( subset_i = 0; subset_i < mi-1; subset_i++ )\r
- {\r
- int idx = (int)(sum_ptr[subset_i] - sum);\r
- int ni = counts[idx];\r
-\r
- if( ni )\r
- {\r
- double s = sum[idx];\r
- lsum += s; L += ni;\r
- rsum -= s; R -= ni;\r
-\r
- if( L && R )\r
- {\r
- double val = (lsum*lsum*R + rsum*rsum*L)/((double)L*R);\r
- if( best_val < val )\r
- {\r
- best_val = val;\r
- best_subset = subset_i;\r
- }\r
- }\r
- }\r
- }\r
-\r
- CvDTreeSplit* split = 0;\r
- if( best_subset >= 0 )\r
- {\r
- split = _split ? _split : data->new_split_cat( 0, -1.0f);\r
- split->var_idx = vi;\r
- split->quality = (float)best_val;\r
- memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));\r
- for( i = 0; i <= best_subset; i++ )\r
- {\r
- int idx = (int)(sum_ptr[i] - sum);\r
- split->subset[idx >> 5] |= 1 << (idx & 31);\r
- }\r
- }\r
- return split;\r
-}\r
-\r
-CvDTreeSplit* CvDTree::find_surrogate_split_ord( CvDTreeNode* node, int vi )\r
-{\r
- const float epsilon = FLT_EPSILON*2;\r
- const char* dir = (char*)data->direction->data.ptr;\r
- int n1 = node->get_num_valid(vi);\r
- float* values_buf = data->get_pred_float_buf();\r
- const float* values = 0;\r
- int* indices_buf = data->get_pred_int_buf();\r
- const int* indices = 0;\r
- data->get_ord_var_data( node, vi, values_buf, indices_buf, &values, &indices );\r
- // LL - number of samples that both the primary and the surrogate splits send to the left\r
- // LR - ... primary split sends to the left and the surrogate split sends to the right\r
- // RL - ... primary split sends to the right and the surrogate split sends to the left\r
- // RR - ... both send to the right\r
- int i, best_i = -1, best_inversed = 0;\r
- double best_val;\r
-\r
- if( !data->have_priors )\r
- {\r
- int LL = 0, RL = 0, LR, RR;\r
- int worst_val = cvFloor(node->maxlr), _best_val = worst_val;\r
- int sum = 0, sum_abs = 0;\r
-\r
- for( i = 0; i < n1; i++ )\r
- {\r
- int d = dir[indices[i]];\r
- sum += d; sum_abs += d & 1;\r
- }\r
-\r
- // sum_abs = R + L; sum = R - L\r
- RR = (sum_abs + sum) >> 1;\r
- LR = (sum_abs - sum) >> 1;\r
-\r
- // initially all the samples are sent to the right by the surrogate split,\r
- // LR of them are sent to the left by primary split, and RR - to the right.\r
- // now iteratively compute LL, LR, RL and RR for every possible surrogate split value.\r
- for( i = 0; i < n1 - 1; i++ )\r
- {\r
- int d = dir[indices[i]];\r
-\r
- if( d < 0 )\r
- {\r
- LL++; LR--;\r
- if( LL + RR > _best_val && values[i] + epsilon < values[i+1] )\r
- {\r
- best_val = LL + RR;\r
- best_i = i; best_inversed = 0;\r
- }\r
- }\r
- else if( d > 0 )\r
- {\r
- RL++; RR--;\r
- if( RL + LR > _best_val && values[i] + epsilon < values[i+1] )\r
- {\r
- best_val = RL + LR;\r
- best_i = i; best_inversed = 1;\r
- }\r
- }\r
- }\r
- best_val = _best_val;\r
- }\r
- else\r
- {\r
- double LL = 0, RL = 0, LR, RR;\r
- double worst_val = node->maxlr;\r
- double sum = 0, sum_abs = 0;\r
- const double* priors = data->priors_mult->data.db;\r
- int* responses_buf = data->get_resp_int_buf();\r
- const int* responses = 0;\r
- data->get_class_labels(node, responses_buf, &responses);\r
- best_val = worst_val;\r
-\r
- for( i = 0; i < n1; i++ )\r
- {\r
- int idx = indices[i];\r
- double w = priors[responses[idx]];\r
- int d = dir[idx];\r
- sum += d*w; sum_abs += (d & 1)*w;\r
- }\r
-\r
- // sum_abs = R + L; sum = R - L\r
- RR = (sum_abs + sum)*0.5;\r
- LR = (sum_abs - sum)*0.5;\r
-\r
- // initially all the samples are sent to the right by the surrogate split,\r
- // LR of them are sent to the left by primary split, and RR - to the right.\r
- // now iteratively compute LL, LR, RL and RR for every possible surrogate split value.\r
- for( i = 0; i < n1 - 1; i++ )\r
- {\r
- int idx = indices[i];\r
- double w = priors[responses[idx]];\r
- int d = dir[idx];\r
-\r
- if( d < 0 )\r
- {\r
- LL += w; LR -= w;\r
- if( LL + RR > best_val && values[i] + epsilon < values[i+1] )\r
- {\r
- best_val = LL + RR;\r
- best_i = i; best_inversed = 0;\r
- }\r
- }\r
- else if( d > 0 )\r
- {\r
- RL += w; RR -= w;\r
- if( RL + LR > best_val && values[i] + epsilon < values[i+1] )\r
- {\r
- best_val = RL + LR;\r
- best_i = i; best_inversed = 1;\r
- }\r
- }\r
- }\r
- }\r
- return best_i >= 0 && best_val > node->maxlr ? data->new_split_ord( vi,\r
- (values[best_i] + values[best_i+1])*0.5f, best_i, best_inversed, (float)best_val ) : 0;\r
-}\r
-\r
-\r
-CvDTreeSplit* CvDTree::find_surrogate_split_cat( CvDTreeNode* node, int vi )\r
-{\r
- const char* dir = (char*)data->direction->data.ptr;\r
- int n = node->sample_count;\r
- int* labels_buf = data->get_pred_int_buf();\r
- const int* labels = 0;\r
- data->get_cat_var_data(node, vi, labels_buf, &labels);\r
- // LL - number of samples that both the primary and the surrogate splits send to the left\r
- // LR - ... primary split sends to the left and the surrogate split sends to the right\r
- // RL - ... primary split sends to the right and the surrogate split sends to the left\r
- // RR - ... both send to the right\r
- CvDTreeSplit* split = data->new_split_cat( vi, 0 );\r
- int i, mi = data->cat_count->data.i[data->get_var_type(vi)], l_win = 0;\r
- double best_val = 0;\r
- double* lc = (double*)cvStackAlloc( (mi+1)*2*sizeof(lc[0]) ) + 1;\r
- double* rc = lc + mi + 1;\r
-\r
- for( i = -1; i < mi; i++ )\r
- lc[i] = rc[i] = 0;\r
-\r
- // for each category calculate the weight of samples\r
- // sent to the left (lc) and to the right (rc) by the primary split\r
- if( !data->have_priors )\r
- {\r
- int* _lc = (int*)cvStackAlloc((mi+2)*2*sizeof(_lc[0])) + 1;\r
- int* _rc = _lc + mi + 1;\r
-\r
- for( i = -1; i < mi; i++ )\r
- _lc[i] = _rc[i] = 0;\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- int idx = ( (labels[i] == 65535) && (data->is_buf_16u) ) ? -1 : labels[i];\r
- int d = dir[i];\r
- int sum = _lc[idx] + d;\r
- int sum_abs = _rc[idx] + (d & 1);\r
- _lc[idx] = sum; _rc[idx] = sum_abs;\r
- }\r
-\r
- for( i = 0; i < mi; i++ )\r
- {\r
- int sum = _lc[i];\r
- int sum_abs = _rc[i];\r
- lc[i] = (sum_abs - sum) >> 1;\r
- rc[i] = (sum_abs + sum) >> 1;\r
- }\r
- }\r
- else\r
- {\r
- const double* priors = data->priors_mult->data.db;\r
- int* responses_buf = data->get_resp_int_buf();\r
- const int* responses = 0;\r
- data->get_class_labels(node, responses_buf, &responses);\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- int idx = ( (labels[i] == 65535) && (data->is_buf_16u) ) ? -1 : labels[i];\r
- double w = priors[responses[i]];\r
- int d = dir[i];\r
- double sum = lc[idx] + d*w;\r
- double sum_abs = rc[idx] + (d & 1)*w;\r
- lc[idx] = sum; rc[idx] = sum_abs;\r
- }\r
-\r
- for( i = 0; i < mi; i++ )\r
- {\r
- double sum = lc[i];\r
- double sum_abs = rc[i];\r
- lc[i] = (sum_abs - sum) * 0.5;\r
- rc[i] = (sum_abs + sum) * 0.5;\r
- }\r
- }\r
-\r
- // 2. now form the split.\r
- // in each category send all the samples to the same direction as majority\r
- for( i = 0; i < mi; i++ )\r
- {\r
- double lval = lc[i], rval = rc[i];\r
- if( lval > rval )\r
- {\r
- split->subset[i >> 5] |= 1 << (i & 31);\r
- best_val += lval;\r
- l_win++;\r
- }\r
- else\r
- best_val += rval;\r
- }\r
-\r
- split->quality = (float)best_val;\r
- if( split->quality <= node->maxlr || l_win == 0 || l_win == mi )\r
- cvSetRemoveByPtr( data->split_heap, split ), split = 0;\r
-\r
- return split;\r
-}\r
-\r
-\r
-void CvDTree::calc_node_value( CvDTreeNode* node )\r
-{\r
- int i, j, k, n = node->sample_count, cv_n = data->params.cv_folds;\r
- int* cv_labels_buf = data->get_cv_lables_buf();\r
- const int* cv_labels = 0;\r
- data->get_cv_labels(node, cv_labels_buf, &cv_labels);\r
-\r
- if( data->is_classifier )\r
- {\r
- // in case of classification tree:\r
- // * node value is the label of the class that has the largest weight in the node.\r
- // * node risk is the weighted number of misclassified samples,\r
- // * j-th cross-validation fold value and risk are calculated as above,\r
- // but using the samples with cv_labels(*)!=j.\r
- // * j-th cross-validation fold error is calculated as the weighted number of\r
- // misclassified samples with cv_labels(*)==j.\r
-\r
- // compute the number of instances of each class\r
- int* cls_count = data->counts->data.i;\r
- int* responses_buf = data->get_resp_int_buf();\r
- const int* responses = 0;\r
- data->get_class_labels(node, responses_buf, &responses);\r
- int m = data->get_num_classes();\r
- int* cv_cls_count = (int*)cvStackAlloc(m*cv_n*sizeof(cv_cls_count[0]));\r
- double max_val = -1, total_weight = 0;\r
- int max_k = -1;\r
- double* priors = data->priors_mult->data.db;\r
-\r
- for( k = 0; k < m; k++ )\r
- cls_count[k] = 0;\r
-\r
- if( cv_n == 0 )\r
- {\r
- for( i = 0; i < n; i++ )\r
- cls_count[responses[i]]++;\r
- }\r
- else\r
- {\r
- for( j = 0; j < cv_n; j++ )\r
- for( k = 0; k < m; k++ )\r
- cv_cls_count[j*m + k] = 0;\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- j = cv_labels[i]; k = responses[i];\r
- cv_cls_count[j*m + k]++;\r
- }\r
-\r
- for( j = 0; j < cv_n; j++ )\r
- for( k = 0; k < m; k++ )\r
- cls_count[k] += cv_cls_count[j*m + k];\r
- }\r
-\r
- if( data->have_priors && node->parent == 0 )\r
- {\r
- // compute priors_mult from priors, take the sample ratio into account.\r
- double sum = 0;\r
- for( k = 0; k < m; k++ )\r
- {\r
- int n_k = cls_count[k];\r
- priors[k] = data->priors->data.db[k]*(n_k ? 1./n_k : 0.);\r
- sum += priors[k];\r
- }\r
- sum = 1./sum;\r
- for( k = 0; k < m; k++ )\r
- priors[k] *= sum;\r
- }\r
-\r
- for( k = 0; k < m; k++ )\r
- {\r
- double val = cls_count[k]*priors[k];\r
- total_weight += val;\r
- if( max_val < val )\r
- {\r
- max_val = val;\r
- max_k = k;\r
- }\r
- }\r
-\r
- node->class_idx = max_k;\r
- node->value = data->cat_map->data.i[\r
- data->cat_ofs->data.i[data->cat_var_count] + max_k];\r
- node->node_risk = total_weight - max_val;\r
-\r
- for( j = 0; j < cv_n; j++ )\r
- {\r
- double sum_k = 0, sum = 0, max_val_k = 0;\r
- max_val = -1; max_k = -1;\r
-\r
- for( k = 0; k < m; k++ )\r
- {\r
- double w = priors[k];\r
- double val_k = cv_cls_count[j*m + k]*w;\r
- double val = cls_count[k]*w - val_k;\r
- sum_k += val_k;\r
- sum += val;\r
- if( max_val < val )\r
- {\r
- max_val = val;\r
- max_val_k = val_k;\r
- max_k = k;\r
- }\r
- }\r
-\r
- node->cv_Tn[j] = INT_MAX;\r
- node->cv_node_risk[j] = sum - max_val;\r
- node->cv_node_error[j] = sum_k - max_val_k;\r
- }\r
- }\r
- else\r
- {\r
- // in case of regression tree:\r
- // * node value is 1/n*sum_i(Y_i), where Y_i is i-th response,\r
- // n is the number of samples in the node.\r
- // * node risk is the sum of squared errors: sum_i((Y_i - <node_value>)^2)\r
- // * j-th cross-validation fold value and risk are calculated as above,\r
- // but using the samples with cv_labels(*)!=j.\r
- // * j-th cross-validation fold error is calculated\r
- // using samples with cv_labels(*)==j as the test subset:\r
- // error_j = sum_(i,cv_labels(i)==j)((Y_i - <node_value_j>)^2),\r
- // where node_value_j is the node value calculated\r
- // as described in the previous bullet, and summation is done\r
- // over the samples with cv_labels(*)==j.\r
-\r
- double sum = 0, sum2 = 0;\r
- float* values_buf = data->get_resp_float_buf();\r
- const float* values = 0;\r
- data->get_ord_responses(node, values_buf, &values);\r
- double *cv_sum = 0, *cv_sum2 = 0;\r
- int* cv_count = 0;\r
-\r
- if( cv_n == 0 )\r
- {\r
- for( i = 0; i < n; i++ )\r
- {\r
- double t = values[i];\r
- sum += t;\r
- sum2 += t*t;\r
- }\r
- }\r
- else\r
- {\r
- cv_sum = (double*)cvStackAlloc( cv_n*sizeof(cv_sum[0]) );\r
- cv_sum2 = (double*)cvStackAlloc( cv_n*sizeof(cv_sum2[0]) );\r
- cv_count = (int*)cvStackAlloc( cv_n*sizeof(cv_count[0]) );\r
-\r
- for( j = 0; j < cv_n; j++ )\r
- {\r
- cv_sum[j] = cv_sum2[j] = 0.;\r
- cv_count[j] = 0;\r
- }\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- j = cv_labels[i];\r
- double t = values[i];\r
- double s = cv_sum[j] + t;\r
- double s2 = cv_sum2[j] + t*t;\r
- int nc = cv_count[j] + 1;\r
- cv_sum[j] = s;\r
- cv_sum2[j] = s2;\r
- cv_count[j] = nc;\r
- }\r
-\r
- for( j = 0; j < cv_n; j++ )\r
- {\r
- sum += cv_sum[j];\r
- sum2 += cv_sum2[j];\r
- }\r
- }\r
-\r
- node->node_risk = sum2 - (sum/n)*sum;\r
- node->value = sum/n;\r
-\r
- for( j = 0; j < cv_n; j++ )\r
- {\r
- double s = cv_sum[j], si = sum - s;\r
- double s2 = cv_sum2[j], s2i = sum2 - s2;\r
- int c = cv_count[j], ci = n - c;\r
- double r = si/MAX(ci,1);\r
- node->cv_node_risk[j] = s2i - r*r*ci;\r
- node->cv_node_error[j] = s2 - 2*r*s + c*r*r;\r
- node->cv_Tn[j] = INT_MAX;\r
- }\r
- }\r
-}\r
-\r
-\r
-void CvDTree::complete_node_dir( CvDTreeNode* node )\r
-{\r
- int vi, i, n = node->sample_count, nl, nr, d0 = 0, d1 = -1;\r
- int nz = n - node->get_num_valid(node->split->var_idx);\r
- char* dir = (char*)data->direction->data.ptr;\r
-\r
- // try to complete direction using surrogate splits\r
- if( nz && data->params.use_surrogates )\r
- {\r
- CvDTreeSplit* split = node->split->next;\r
- for( ; split != 0 && nz; split = split->next )\r
- {\r
- int inversed_mask = split->inversed ? -1 : 0;\r
- vi = split->var_idx;\r
-\r
- if( data->get_var_type(vi) >= 0 ) // split on categorical var\r
- {\r
- int* labels_buf = data->get_pred_int_buf();\r
- const int* labels = 0;\r
- data->get_cat_var_data(node, vi, labels_buf, &labels);\r
- const int* subset = split->subset;\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- int idx = labels[i];\r
- if( !dir[i] && ( ((idx >= 0)&&(!data->is_buf_16u)) || ((idx != 65535)&&(data->is_buf_16u)) ))\r
- \r
- {\r
- int d = CV_DTREE_CAT_DIR(idx,subset);\r
- dir[i] = (char)((d ^ inversed_mask) - inversed_mask);\r
- if( --nz )\r
- break;\r
- }\r
- }\r
- }\r
- else // split on ordered var\r
- {\r
- float* values_buf = data->get_pred_float_buf();\r
- const float* values = 0;\r
- int* indices_buf = data->get_pred_int_buf();\r
- const int* indices = 0;\r
- data->get_ord_var_data( node, vi, values_buf, indices_buf, &values, &indices );\r
- int split_point = split->ord.split_point;\r
- int n1 = node->get_num_valid(vi);\r
-\r
- assert( 0 <= split_point && split_point < n-1 );\r
-\r
- for( i = 0; i < n1; i++ )\r
- {\r
- int idx = indices[i];\r
- if( !dir[idx] )\r
- {\r
- int d = i <= split_point ? -1 : 1;\r
- dir[idx] = (char)((d ^ inversed_mask) - inversed_mask);\r
- if( --nz )\r
- break;\r
- }\r
- }\r
- }\r
- }\r
- }\r
-\r
- // find the default direction for the rest\r
- if( nz )\r
- {\r
- for( i = nr = 0; i < n; i++ )\r
- nr += dir[i] > 0;\r
- nl = n - nr - nz;\r
- d0 = nl > nr ? -1 : nr > nl;\r
- }\r
-\r
- // make sure that every sample is directed either to the left or to the right\r
- for( i = 0; i < n; i++ )\r
- {\r
- int d = dir[i];\r
- if( !d )\r
- {\r
- d = d0;\r
- if( !d )\r
- d = d1, d1 = -d1;\r
- }\r
- d = d > 0;\r
- dir[i] = (char)d; // remap (-1,1) to (0,1)\r
- }\r
-}\r
-\r
-\r
-void CvDTree::split_node_data( CvDTreeNode* node )\r
-{\r
- int vi, i, n = node->sample_count, nl, nr, scount = data->sample_count;\r
- char* dir = (char*)data->direction->data.ptr;\r
- CvDTreeNode *left = 0, *right = 0;\r
- int* new_idx = data->split_buf->data.i;\r
- int new_buf_idx = data->get_child_buf_idx( node );\r
- int work_var_count = data->get_work_var_count();\r
- CvMat* buf = data->buf;\r
- int* temp_buf = (int*)cvStackAlloc(n*sizeof(temp_buf[0]));\r
-\r
- complete_node_dir(node);\r
-\r
- for( i = nl = nr = 0; i < n; i++ )\r
- {\r
- int d = dir[i];\r
- // initialize new indices for splitting ordered variables\r
- new_idx[i] = (nl & (d-1)) | (nr & -d); // d ? ri : li\r
- nr += d;\r
- nl += d^1;\r
- }\r
-\r
-\r
- bool split_input_data;\r
- node->left = left = data->new_node( node, nl, new_buf_idx, node->offset );\r
- node->right = right = data->new_node( node, nr, new_buf_idx, node->offset + nl );\r
-\r
- split_input_data = node->depth + 1 < data->params.max_depth &&\r
- (node->left->sample_count > data->params.min_sample_count ||\r
- node->right->sample_count > data->params.min_sample_count);\r
-\r
- // split ordered variables, keep both halves sorted.\r
- for( vi = 0; vi < data->var_count; vi++ )\r
- {\r
- int ci = data->get_var_type(vi);\r
- int n1 = node->get_num_valid(vi);\r
- int *src_idx_buf = data->get_pred_int_buf();\r
- const int* src_idx = 0;\r
- float *src_val_buf = data->get_pred_float_buf();\r
- const float* src_val = 0;\r
- \r
- if( ci >= 0 || !split_input_data )\r
- continue;\r
-\r
- data->get_ord_var_data(node, vi, src_val_buf, src_idx_buf, &src_val, &src_idx);\r
-\r
- for(i = 0; i < n; i++)\r
- temp_buf[i] = src_idx[i];\r
-\r
- if (data->is_buf_16u)\r
- {\r
- unsigned short *ldst, *rdst, *ldst0, *rdst0;\r
- //unsigned short tl, tr;\r
- ldst0 = ldst = (unsigned short*)(buf->data.s + left->buf_idx*buf->cols + \r
- vi*scount + left->offset);\r
- rdst0 = rdst = (unsigned short*)(ldst + nl);\r
-\r
- // split sorted\r
- for( i = 0; i < n1; i++ )\r
- {\r
- int idx = temp_buf[i];\r
- int d = dir[idx];\r
- idx = new_idx[idx];\r
- if (d)\r
- {\r
- *rdst = (unsigned short)idx;\r
- rdst++;\r
- }\r
- else\r
- {\r
- *ldst = (unsigned short)idx;\r
- ldst++;\r
- }\r
- }\r
-\r
- left->set_num_valid(vi, (int)(ldst - ldst0));\r
- right->set_num_valid(vi, (int)(rdst - rdst0));\r
-\r
- // split missing\r
- for( ; i < n; i++ )\r
- {\r
- int idx = temp_buf[i];\r
- int d = dir[idx];\r
- idx = new_idx[idx];\r
- if (d)\r
- {\r
- *rdst = (unsigned short)idx;\r
- rdst++;\r
- }\r
- else\r
- {\r
- *ldst = (unsigned short)idx;\r
- ldst++;\r
- }\r
- }\r
- }\r
- else\r
- {\r
- int *ldst0, *ldst, *rdst0, *rdst;\r
- ldst0 = ldst = buf->data.i + left->buf_idx*buf->cols + \r
- vi*scount + left->offset;\r
- rdst0 = rdst = buf->data.i + right->buf_idx*buf->cols + \r
- vi*scount + right->offset;\r
-\r
- // split sorted\r
- for( i = 0; i < n1; i++ )\r
- {\r
- int idx = temp_buf[i];\r
- int d = dir[idx];\r
- idx = new_idx[idx];\r
- if (d)\r
- {\r
- *rdst = idx;\r
- rdst++;\r
- }\r
- else\r
- {\r
- *ldst = idx;\r
- ldst++;\r
- }\r
- }\r
-\r
- left->set_num_valid(vi, (int)(ldst - ldst0));\r
- right->set_num_valid(vi, (int)(rdst - rdst0));\r
-\r
- // split missing\r
- for( ; i < n; i++ )\r
- {\r
- int idx = temp_buf[i];\r
- int d = dir[idx];\r
- idx = new_idx[idx];\r
- if (d)\r
- {\r
- *rdst = idx;\r
- rdst++;\r
- }\r
- else\r
- {\r
- *ldst = idx;\r
- ldst++;\r
- }\r
- }\r
- }\r
- }\r
-\r
- // split categorical vars, responses and cv_labels using new_idx relocation table\r
- for( vi = 0; vi < work_var_count; vi++ )\r
- {\r
- int ci = data->get_var_type(vi);\r
- int n1 = node->get_num_valid(vi), nr1 = 0;\r
- \r
- if( ci < 0 || (vi < data->var_count && !split_input_data) )\r
- continue;\r
-\r
- int *src_lbls_buf = data->get_pred_int_buf();\r
- const int* src_lbls = 0;\r
- data->get_cat_var_data(node, vi, src_lbls_buf, &src_lbls);\r
-\r
- for(i = 0; i < n; i++)\r
- temp_buf[i] = src_lbls[i];\r
-\r
- if (data->is_buf_16u)\r
- {\r
- unsigned short *ldst = (unsigned short *)(buf->data.s + left->buf_idx*buf->cols + \r
- vi*scount + left->offset);\r
- unsigned short *rdst = (unsigned short *)(buf->data.s + right->buf_idx*buf->cols + \r
- vi*scount + right->offset);\r
- \r
- for( i = 0; i < n; i++ )\r
- {\r
- int d = dir[i];\r
- int idx = temp_buf[i];\r
- if (d)\r
- {\r
- *rdst = (unsigned short)idx;\r
- rdst++;\r
- nr1 += (idx != 65535 )&d;\r
- }\r
- else\r
- {\r
- *ldst = (unsigned short)idx;\r
- ldst++;\r
- }\r
- }\r
-\r
- if( vi < data->var_count )\r
- {\r
- left->set_num_valid(vi, n1 - nr1);\r
- right->set_num_valid(vi, nr1);\r
- }\r
- }\r
- else\r
- {\r
- int *ldst = buf->data.i + left->buf_idx*buf->cols + \r
- vi*scount + left->offset;\r
- int *rdst = buf->data.i + right->buf_idx*buf->cols + \r
- vi*scount + right->offset;\r
- \r
- for( i = 0; i < n; i++ )\r
- {\r
- int d = dir[i];\r
- int idx = temp_buf[i];\r
- if (d)\r
- {\r
- *rdst = idx;\r
- rdst++;\r
- nr1 += (idx >= 0)&d;\r
- }\r
- else\r
- {\r
- *ldst = idx;\r
- ldst++;\r
- }\r
- \r
- }\r
-\r
- if( vi < data->var_count )\r
- {\r
- left->set_num_valid(vi, n1 - nr1);\r
- right->set_num_valid(vi, nr1);\r
- }\r
- } \r
- }\r
-\r
-\r
- // split sample indices\r
- int *sample_idx_src_buf = data->get_sample_idx_buf();\r
- const int* sample_idx_src = 0;\r
- data->get_sample_indices(node, sample_idx_src_buf, &sample_idx_src);\r
-\r
- for(i = 0; i < n; i++)\r
- temp_buf[i] = sample_idx_src[i];\r
-\r
- int pos = data->get_work_var_count();\r
- if (data->is_buf_16u)\r
- {\r
- unsigned short* ldst = (unsigned short*)(buf->data.s + left->buf_idx*buf->cols + \r
- pos*scount + left->offset);\r
- unsigned short* rdst = (unsigned short*)(buf->data.s + right->buf_idx*buf->cols + \r
- pos*scount + right->offset);\r
- for (i = 0; i < n; i++)\r
- {\r
- int d = dir[i];\r
- unsigned short idx = (unsigned short)temp_buf[i];\r
- if (d)\r
- {\r
- *rdst = idx;\r
- rdst++;\r
- }\r
- else\r
- {\r
- *ldst = idx;\r
- ldst++;\r
- }\r
- }\r
- }\r
- else\r
- {\r
- int* ldst = buf->data.i + left->buf_idx*buf->cols + \r
- pos*scount + left->offset;\r
- int* rdst = buf->data.i + right->buf_idx*buf->cols + \r
- pos*scount + right->offset;\r
- for (i = 0; i < n; i++)\r
- {\r
- int d = dir[i];\r
- int idx = temp_buf[i];\r
- if (d)\r
- {\r
- *rdst = idx;\r
- rdst++;\r
- }\r
- else\r
- {\r
- *ldst = idx;\r
- ldst++;\r
- }\r
- }\r
- }\r
- \r
- // deallocate the parent node data that is not needed anymore\r
- data->free_node_data(node); \r
-}\r
-\r
-float CvDTree::calc_error( CvMLData* _data, int type )\r
-{\r
- float err = 0;\r
- const CvMat* values = _data->get_values();\r
- const CvMat* response = _data->get_response();\r
- const CvMat* missing = _data->get_missing();\r
- const CvMat* sample_idx = (type == CV_TEST_ERROR) ? _data->get_test_sample_idx() : _data->get_train_sample_idx();\r
- const CvMat* var_types = _data->get_var_types();\r
- int* sidx = sample_idx ? sample_idx->data.i : 0;\r
- int r_step = CV_IS_MAT_CONT(response->type) ?\r
- 1 : response->step / CV_ELEM_SIZE(response->type);\r
- bool is_classifier = var_types->data.ptr[var_types->cols-1] == CV_VAR_CATEGORICAL;\r
- int sample_count = sample_idx ? sample_idx->cols : 0;\r
- sample_count = (type == CV_TRAIN_ERROR && sample_count == 0) ? values->rows : sample_count;\r
- if ( is_classifier )\r
- {\r
- for( int i = 0; i < sample_count; i++ )\r
- {\r
- CvMat sample, miss;\r
- int si = sidx ? sidx[i] : i;\r
- cvGetRow( values, &sample, si ); \r
- if( missing ) \r
- cvGetRow( missing, &miss, si ); \r
- float r = (float)predict( &sample, missing ? &miss : 0 )->value;\r
- int d = fabs((double)r - response->data.fl[si*r_step]) <= FLT_EPSILON ? 0 : 1;\r
- err += d;\r
- }\r
- err = sample_count ? err / (float)sample_count * 100 : -FLT_MAX;\r
- }\r
- else\r
- {\r
- for( int i = 0; i < sample_count; i++ )\r
- {\r
- CvMat sample, miss;\r
- int si = sidx ? sidx[i] : i;\r
- cvGetRow( values, &sample, si ); \r
- if( missing ) \r
- cvGetRow( missing, &miss, si ); \r
- float r = (float)predict( &sample, missing ? &miss : 0 )->value;\r
- float d = r - response->data.fl[si*r_step];\r
- err += d*d;\r
- }\r
- err = sample_count ? err / (float)sample_count : -FLT_MAX; \r
- }\r
- return err;\r
-}\r
-\r
-void CvDTree::prune_cv()\r
-{\r
- CvMat* ab = 0;\r
- CvMat* temp = 0;\r
- CvMat* err_jk = 0;\r
-\r
- // 1. build tree sequence for each cv fold, calculate error_{Tj,beta_k}.\r
- // 2. choose the best tree index (if need, apply 1SE rule).\r
- // 3. store the best index and cut the branches.\r
-\r
- CV_FUNCNAME( "CvDTree::prune_cv" );\r
-\r
- __BEGIN__;\r
-\r
- int ti, j, tree_count = 0, cv_n = data->params.cv_folds, n = root->sample_count;\r
- // currently, 1SE for regression is not implemented\r
- bool use_1se = data->params.use_1se_rule != 0 && data->is_classifier;\r
- double* err;\r
- double min_err = 0, min_err_se = 0;\r
- int min_idx = -1;\r
-\r
- CV_CALL( ab = cvCreateMat( 1, 256, CV_64F ));\r
-\r
- // build the main tree sequence, calculate alpha's\r
- for(;;tree_count++)\r
- {\r
- double min_alpha = update_tree_rnc(tree_count, -1);\r
- if( cut_tree(tree_count, -1, min_alpha) )\r
- break;\r
-\r
- if( ab->cols <= tree_count )\r
- {\r
- CV_CALL( temp = cvCreateMat( 1, ab->cols*3/2, CV_64F ));\r
- for( ti = 0; ti < ab->cols; ti++ )\r
- temp->data.db[ti] = ab->data.db[ti];\r
- cvReleaseMat( &ab );\r
- ab = temp;\r
- temp = 0;\r
- }\r
-\r
- ab->data.db[tree_count] = min_alpha;\r
- }\r
-\r
- ab->data.db[0] = 0.;\r
-\r
- if( tree_count > 0 )\r
- {\r
- for( ti = 1; ti < tree_count-1; ti++ )\r
- ab->data.db[ti] = sqrt(ab->data.db[ti]*ab->data.db[ti+1]);\r
- ab->data.db[tree_count-1] = DBL_MAX*0.5;\r
-\r
- CV_CALL( err_jk = cvCreateMat( cv_n, tree_count, CV_64F ));\r
- err = err_jk->data.db;\r
-\r
- for( j = 0; j < cv_n; j++ )\r
- {\r
- int tj = 0, tk = 0;\r
- for( ; tk < tree_count; tj++ )\r
- {\r
- double min_alpha = update_tree_rnc(tj, j);\r
- if( cut_tree(tj, j, min_alpha) )\r
- min_alpha = DBL_MAX;\r
-\r
- for( ; tk < tree_count; tk++ )\r
- {\r
- if( ab->data.db[tk] > min_alpha )\r
- break;\r
- err[j*tree_count + tk] = root->tree_error;\r
- }\r
- }\r
- }\r
-\r
- for( ti = 0; ti < tree_count; ti++ )\r
- {\r
- double sum_err = 0;\r
- for( j = 0; j < cv_n; j++ )\r
- sum_err += err[j*tree_count + ti];\r
- if( ti == 0 || sum_err < min_err )\r
- {\r
- min_err = sum_err;\r
- min_idx = ti;\r
- if( use_1se )\r
- min_err_se = sqrt( sum_err*(n - sum_err) );\r
- }\r
- else if( sum_err < min_err + min_err_se )\r
- min_idx = ti;\r
- }\r
- }\r
-\r
- pruned_tree_idx = min_idx;\r
- free_prune_data(data->params.truncate_pruned_tree != 0);\r
-\r
- __END__;\r
-\r
- cvReleaseMat( &err_jk );\r
- cvReleaseMat( &ab );\r
- cvReleaseMat( &temp );\r
-}\r
-\r
-\r
-double CvDTree::update_tree_rnc( int T, int fold )\r
-{\r
- CvDTreeNode* node = root;\r
- double min_alpha = DBL_MAX;\r
-\r
- for(;;)\r
- {\r
- CvDTreeNode* parent;\r
- for(;;)\r
- {\r
- int t = fold >= 0 ? node->cv_Tn[fold] : node->Tn;\r
- if( t <= T || !node->left )\r
- {\r
- node->complexity = 1;\r
- node->tree_risk = node->node_risk;\r
- node->tree_error = 0.;\r
- if( fold >= 0 )\r
- {\r
- node->tree_risk = node->cv_node_risk[fold];\r
- node->tree_error = node->cv_node_error[fold];\r
- }\r
- break;\r
- }\r
- node = node->left;\r
- }\r
-\r
- for( parent = node->parent; parent && parent->right == node;\r
- node = parent, parent = parent->parent )\r
- {\r
- parent->complexity += node->complexity;\r
- parent->tree_risk += node->tree_risk;\r
- parent->tree_error += node->tree_error;\r
-\r
- parent->alpha = ((fold >= 0 ? parent->cv_node_risk[fold] : parent->node_risk)\r
- - parent->tree_risk)/(parent->complexity - 1);\r
- min_alpha = MIN( min_alpha, parent->alpha );\r
- }\r
-\r
- if( !parent )\r
- break;\r
-\r
- parent->complexity = node->complexity;\r
- parent->tree_risk = node->tree_risk;\r
- parent->tree_error = node->tree_error;\r
- node = parent->right;\r
- }\r
-\r
- return min_alpha;\r
-}\r
-\r
-\r
-int CvDTree::cut_tree( int T, int fold, double min_alpha )\r
-{\r
- CvDTreeNode* node = root;\r
- if( !node->left )\r
- return 1;\r
-\r
- for(;;)\r
- {\r
- CvDTreeNode* parent;\r
- for(;;)\r
- {\r
- int t = fold >= 0 ? node->cv_Tn[fold] : node->Tn;\r
- if( t <= T || !node->left )\r
- break;\r
- if( node->alpha <= min_alpha + FLT_EPSILON )\r
- {\r
- if( fold >= 0 )\r
- node->cv_Tn[fold] = T;\r
- else\r
- node->Tn = T;\r
- if( node == root )\r
- return 1;\r
- break;\r
- }\r
- node = node->left;\r
- }\r
-\r
- for( parent = node->parent; parent && parent->right == node;\r
- node = parent, parent = parent->parent )\r
- ;\r
-\r
- if( !parent )\r
- break;\r
-\r
- node = parent->right;\r
- }\r
-\r
- return 0;\r
-}\r
-\r
-\r
-void CvDTree::free_prune_data(bool cut_tree)\r
-{\r
- CvDTreeNode* node = root;\r
-\r
- for(;;)\r
- {\r
- CvDTreeNode* parent;\r
- for(;;)\r
- {\r
- // do not call cvSetRemoveByPtr( cv_heap, node->cv_Tn )\r
- // as we will clear the whole cross-validation heap at the end\r
- node->cv_Tn = 0;\r
- node->cv_node_error = node->cv_node_risk = 0;\r
- if( !node->left )\r
- break;\r
- node = node->left;\r
- }\r
-\r
- for( parent = node->parent; parent && parent->right == node;\r
- node = parent, parent = parent->parent )\r
- {\r
- if( cut_tree && parent->Tn <= pruned_tree_idx )\r
- {\r
- data->free_node( parent->left );\r
- data->free_node( parent->right );\r
- parent->left = parent->right = 0;\r
- }\r
- }\r
-\r
- if( !parent )\r
- break;\r
-\r
- node = parent->right;\r
- }\r
-\r
- if( data->cv_heap )\r
- cvClearSet( data->cv_heap );\r
-}\r
-\r
-\r
-void CvDTree::free_tree()\r
-{\r
- if( root && data && data->shared )\r
- {\r
- pruned_tree_idx = INT_MIN;\r
- free_prune_data(true);\r
- data->free_node(root);\r
- root = 0;\r
- }\r
-}\r
-\r
-CvDTreeNode* CvDTree::predict( const CvMat* _sample,\r
- const CvMat* _missing, bool preprocessed_input ) const\r
-{\r
- CvDTreeNode* result = 0;\r
- int* catbuf = 0;\r
-\r
- CV_FUNCNAME( "CvDTree::predict" );\r
-\r
- __BEGIN__;\r
-\r
- int i, step, mstep = 0;\r
- const float* sample;\r
- const uchar* m = 0;\r
- CvDTreeNode* node = root;\r
- const int* vtype;\r
- const int* vidx;\r
- const int* cmap;\r
- const int* cofs;\r
-\r
- if( !node )\r
- CV_ERROR( CV_StsError, "The tree has not been trained yet" );\r
-\r
- if( !CV_IS_MAT(_sample) || CV_MAT_TYPE(_sample->type) != CV_32FC1 ||\r
- (_sample->cols != 1 && _sample->rows != 1) ||\r
- (_sample->cols + _sample->rows - 1 != data->var_all && !preprocessed_input) ||\r
- (_sample->cols + _sample->rows - 1 != data->var_count && preprocessed_input) )\r
- CV_ERROR( CV_StsBadArg,\r
- "the input sample must be 1d floating-point vector with the same "\r
- "number of elements as the total number of variables used for training" );\r
-\r
- sample = _sample->data.fl;\r
- step = CV_IS_MAT_CONT(_sample->type) ? 1 : _sample->step/sizeof(sample[0]);\r
-\r
- if( data->cat_count && !preprocessed_input ) // cache for categorical variables\r
- {\r
- int n = data->cat_count->cols;\r
- catbuf = (int*)cvStackAlloc(n*sizeof(catbuf[0]));\r
- for( i = 0; i < n; i++ )\r
- catbuf[i] = -1;\r
- }\r
-\r
- if( _missing )\r
- {\r
- if( !CV_IS_MAT(_missing) || !CV_IS_MASK_ARR(_missing) ||\r
- !CV_ARE_SIZES_EQ(_missing, _sample) )\r
- CV_ERROR( CV_StsBadArg,\r
- "the missing data mask must be 8-bit vector of the same size as input sample" );\r
- m = _missing->data.ptr;\r
- mstep = CV_IS_MAT_CONT(_missing->type) ? 1 : _missing->step/sizeof(m[0]);\r
- }\r
-\r
- vtype = data->var_type->data.i;\r
- vidx = data->var_idx && !preprocessed_input ? data->var_idx->data.i : 0;\r
- cmap = data->cat_map ? data->cat_map->data.i : 0;\r
- cofs = data->cat_ofs ? data->cat_ofs->data.i : 0;\r
-\r
- while( node->Tn > pruned_tree_idx && node->left )\r
- {\r
- CvDTreeSplit* split = node->split;\r
- int dir = 0;\r
- for( ; !dir && split != 0; split = split->next )\r
- {\r
- int vi = split->var_idx;\r
- int ci = vtype[vi];\r
- i = vidx ? vidx[vi] : vi;\r
- float val = sample[i*step];\r
- if( m && m[i*mstep] )\r
- continue;\r
- if( ci < 0 ) // ordered\r
- dir = val <= split->ord.c ? -1 : 1;\r
- else // categorical\r
- {\r
- int c;\r
- if( preprocessed_input )\r
- c = cvRound(val);\r
- else\r
- {\r
- c = catbuf[ci];\r
- if( c < 0 )\r
- {\r
- int a = c = cofs[ci];\r
- int b = (ci+1 >= data->cat_ofs->cols) ? data->cat_map->cols : cofs[ci+1];\r
- \r
- int ival = cvRound(val);\r
- if( ival != val )\r
- CV_ERROR( CV_StsBadArg,\r
- "one of input categorical variable is not an integer" );\r
- \r
- int sh = 0;\r
- while( a < b )\r
- {\r
- sh++;\r
- c = (a + b) >> 1;\r
- if( ival < cmap[c] )\r
- b = c;\r
- else if( ival > cmap[c] )\r
- a = c+1;\r
- else\r
- break;\r
- }\r
-\r
- if( c < 0 || ival != cmap[c] )\r
- continue;\r
-\r
- catbuf[ci] = c -= cofs[ci];\r
- }\r
- }\r
- c = ( (c == 65535) && data->is_buf_16u ) ? -1 : c;\r
- dir = CV_DTREE_CAT_DIR(c, split->subset);\r
- }\r
-\r
- if( split->inversed )\r
- dir = -dir;\r
- }\r
-\r
- if( !dir )\r
- {\r
- double diff = node->right->sample_count - node->left->sample_count;\r
- dir = diff < 0 ? -1 : 1;\r
- }\r
- node = dir < 0 ? node->left : node->right;\r
- }\r
-\r
- result = node;\r
-\r
- __END__;\r
-\r
- return result;\r
-}\r
-\r
-\r
-const CvMat* CvDTree::get_var_importance()\r
-{\r
- if( !var_importance )\r
- {\r
- CvDTreeNode* node = root;\r
- double* importance;\r
- if( !node )\r
- return 0;\r
- var_importance = cvCreateMat( 1, data->var_count, CV_64F );\r
- cvZero( var_importance );\r
- importance = var_importance->data.db;\r
-\r
- for(;;)\r
- {\r
- CvDTreeNode* parent;\r
- for( ;; node = node->left )\r
- {\r
- CvDTreeSplit* split = node->split;\r
-\r
- if( !node->left || node->Tn <= pruned_tree_idx )\r
- break;\r
-\r
- for( ; split != 0; split = split->next )\r
- importance[split->var_idx] += split->quality;\r
- }\r
-\r
- for( parent = node->parent; parent && parent->right == node;\r
- node = parent, parent = parent->parent )\r
- ;\r
-\r
- if( !parent )\r
- break;\r
-\r
- node = parent->right;\r
- }\r
-\r
- cvNormalize( var_importance, var_importance, 1., 0, CV_L1 );\r
- }\r
-\r
- return var_importance;\r
-}\r
-\r
-\r
-void CvDTree::write_split( CvFileStorage* fs, CvDTreeSplit* split )\r
-{\r
- int ci;\r
-\r
- cvStartWriteStruct( fs, 0, CV_NODE_MAP + CV_NODE_FLOW );\r
- cvWriteInt( fs, "var", split->var_idx );\r
- cvWriteReal( fs, "quality", split->quality );\r
-\r
- ci = data->get_var_type(split->var_idx);\r
- if( ci >= 0 ) // split on a categorical var\r
- {\r
- int i, n = data->cat_count->data.i[ci], to_right = 0, default_dir;\r
- for( i = 0; i < n; i++ )\r
- to_right += CV_DTREE_CAT_DIR(i,split->subset) > 0;\r
-\r
- // ad-hoc rule when to use inverse categorical split notation\r
- // to achieve more compact and clear representation\r
- default_dir = to_right <= 1 || to_right <= MIN(3, n/2) || to_right <= n/3 ? -1 : 1;\r
-\r
- cvStartWriteStruct( fs, default_dir*(split->inversed ? -1 : 1) > 0 ?\r
- "in" : "not_in", CV_NODE_SEQ+CV_NODE_FLOW );\r
-\r
- for( i = 0; i < n; i++ )\r
- {\r
- int dir = CV_DTREE_CAT_DIR(i,split->subset);\r
- if( dir*default_dir < 0 )\r
- cvWriteInt( fs, 0, i );\r
- }\r
- cvEndWriteStruct( fs );\r
- }\r
- else\r
- cvWriteReal( fs, !split->inversed ? "le" : "gt", split->ord.c );\r
-\r
- cvEndWriteStruct( fs );\r
-}\r
-\r
-\r
-void CvDTree::write_node( CvFileStorage* fs, CvDTreeNode* node )\r
-{\r
- CvDTreeSplit* split;\r
-\r
- cvStartWriteStruct( fs, 0, CV_NODE_MAP );\r
-\r
- cvWriteInt( fs, "depth", node->depth );\r
- cvWriteInt( fs, "sample_count", node->sample_count );\r
- cvWriteReal( fs, "value", node->value );\r
-\r
- if( data->is_classifier )\r
- cvWriteInt( fs, "norm_class_idx", node->class_idx );\r
-\r
- cvWriteInt( fs, "Tn", node->Tn );\r
- cvWriteInt( fs, "complexity", node->complexity );\r
- cvWriteReal( fs, "alpha", node->alpha );\r
- cvWriteReal( fs, "node_risk", node->node_risk );\r
- cvWriteReal( fs, "tree_risk", node->tree_risk );\r
- cvWriteReal( fs, "tree_error", node->tree_error );\r
-\r
- if( node->left )\r
- {\r
- cvStartWriteStruct( fs, "splits", CV_NODE_SEQ );\r
-\r
- for( split = node->split; split != 0; split = split->next )\r
- write_split( fs, split );\r
-\r
- cvEndWriteStruct( fs );\r
- }\r
-\r
- cvEndWriteStruct( fs );\r
-}\r
-\r
-\r
-void CvDTree::write_tree_nodes( CvFileStorage* fs )\r
-{\r
- //CV_FUNCNAME( "CvDTree::write_tree_nodes" );\r
-\r
- __BEGIN__;\r
-\r
- CvDTreeNode* node = root;\r
-\r
- // traverse the tree and save all the nodes in depth-first order\r
- for(;;)\r
- {\r
- CvDTreeNode* parent;\r
- for(;;)\r
- {\r
- write_node( fs, node );\r
- if( !node->left )\r
- break;\r
- node = node->left;\r
- }\r
-\r
- for( parent = node->parent; parent && parent->right == node;\r
- node = parent, parent = parent->parent )\r
- ;\r
-\r
- if( !parent )\r
- break;\r
-\r
- node = parent->right;\r
- }\r
-\r
- __END__;\r
-}\r
-\r
-\r
-void CvDTree::write( CvFileStorage* fs, const char* name )\r
-{\r
- //CV_FUNCNAME( "CvDTree::write" );\r
-\r
- __BEGIN__;\r
-\r
- cvStartWriteStruct( fs, name, CV_NODE_MAP, CV_TYPE_NAME_ML_TREE );\r
-\r
- get_var_importance();\r
- data->write_params( fs );\r
- if( var_importance )\r
- cvWrite( fs, "var_importance", var_importance );\r
- write( fs );\r
-\r
- cvEndWriteStruct( fs );\r
-\r
- __END__;\r
-}\r
-\r
-\r
-void CvDTree::write( CvFileStorage* fs )\r
-{\r
- //CV_FUNCNAME( "CvDTree::write" );\r
-\r
- __BEGIN__;\r
-\r
- cvWriteInt( fs, "best_tree_idx", pruned_tree_idx );\r
-\r
- cvStartWriteStruct( fs, "nodes", CV_NODE_SEQ );\r
- write_tree_nodes( fs );\r
- cvEndWriteStruct( fs );\r
-\r
- __END__;\r
-}\r
-\r
-\r
-CvDTreeSplit* CvDTree::read_split( CvFileStorage* fs, CvFileNode* fnode )\r
-{\r
- CvDTreeSplit* split = 0;\r
-\r
- CV_FUNCNAME( "CvDTree::read_split" );\r
-\r
- __BEGIN__;\r
-\r
- int vi, ci;\r
-\r
- if( !fnode || CV_NODE_TYPE(fnode->tag) != CV_NODE_MAP )\r
- CV_ERROR( CV_StsParseError, "some of the splits are not stored properly" );\r
-\r
- vi = cvReadIntByName( fs, fnode, "var", -1 );\r
- if( (unsigned)vi >= (unsigned)data->var_count )\r
- CV_ERROR( CV_StsOutOfRange, "Split variable index is out of range" );\r
-\r
- ci = data->get_var_type(vi);\r
- if( ci >= 0 ) // split on categorical var\r
- {\r
- int i, n = data->cat_count->data.i[ci], inversed = 0, val;\r
- CvSeqReader reader;\r
- CvFileNode* inseq;\r
- split = data->new_split_cat( vi, 0 );\r
- inseq = cvGetFileNodeByName( fs, fnode, "in" );\r
- if( !inseq )\r
- {\r
- inseq = cvGetFileNodeByName( fs, fnode, "not_in" );\r
- inversed = 1;\r
- }\r
- if( !inseq ||\r
- (CV_NODE_TYPE(inseq->tag) != CV_NODE_SEQ && CV_NODE_TYPE(inseq->tag) != CV_NODE_INT))\r
- CV_ERROR( CV_StsParseError,\r
- "Either 'in' or 'not_in' tags should be inside a categorical split data" );\r
-\r
- if( CV_NODE_TYPE(inseq->tag) == CV_NODE_INT )\r
- {\r
- val = inseq->data.i;\r
- if( (unsigned)val >= (unsigned)n )\r
- CV_ERROR( CV_StsOutOfRange, "some of in/not_in elements are out of range" );\r
-\r
- split->subset[val >> 5] |= 1 << (val & 31);\r
- }\r
- else\r
- {\r
- cvStartReadSeq( inseq->data.seq, &reader );\r
-\r
- for( i = 0; i < reader.seq->total; i++ )\r
- {\r
- CvFileNode* inode = (CvFileNode*)reader.ptr;\r
- val = inode->data.i;\r
- if( CV_NODE_TYPE(inode->tag) != CV_NODE_INT || (unsigned)val >= (unsigned)n )\r
- CV_ERROR( CV_StsOutOfRange, "some of in/not_in elements are out of range" );\r
-\r
- split->subset[val >> 5] |= 1 << (val & 31);\r
- CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );\r
- }\r
- }\r
-\r
- // for categorical splits we do not use inversed splits,\r
- // instead we inverse the variable set in the split\r
- if( inversed )\r
- for( i = 0; i < (n + 31) >> 5; i++ )\r
- split->subset[i] ^= -1;\r
- }\r
- else\r
- {\r
- CvFileNode* cmp_node;\r
- split = data->new_split_ord( vi, 0, 0, 0, 0 );\r
-\r
- cmp_node = cvGetFileNodeByName( fs, fnode, "le" );\r
- if( !cmp_node )\r
- {\r
- cmp_node = cvGetFileNodeByName( fs, fnode, "gt" );\r
- split->inversed = 1;\r
- }\r
-\r
- split->ord.c = (float)cvReadReal( cmp_node );\r
- }\r
-\r
- split->quality = (float)cvReadRealByName( fs, fnode, "quality" );\r
-\r
- __END__;\r
-\r
- return split;\r
-}\r
-\r
-\r
-CvDTreeNode* CvDTree::read_node( CvFileStorage* fs, CvFileNode* fnode, CvDTreeNode* parent )\r
-{\r
- CvDTreeNode* node = 0;\r
-\r
- CV_FUNCNAME( "CvDTree::read_node" );\r
-\r
- __BEGIN__;\r
-\r
- CvFileNode* splits;\r
- int i, depth;\r
-\r
- if( !fnode || CV_NODE_TYPE(fnode->tag) != CV_NODE_MAP )\r
- CV_ERROR( CV_StsParseError, "some of the tree elements are not stored properly" );\r
-\r
- CV_CALL( node = data->new_node( parent, 0, 0, 0 ));\r
- depth = cvReadIntByName( fs, fnode, "depth", -1 );\r
- if( depth != node->depth )\r
- CV_ERROR( CV_StsParseError, "incorrect node depth" );\r
-\r
- node->sample_count = cvReadIntByName( fs, fnode, "sample_count" );\r
- node->value = cvReadRealByName( fs, fnode, "value" );\r
- if( data->is_classifier )\r
- node->class_idx = cvReadIntByName( fs, fnode, "norm_class_idx" );\r
-\r
- node->Tn = cvReadIntByName( fs, fnode, "Tn" );\r
- node->complexity = cvReadIntByName( fs, fnode, "complexity" );\r
- node->alpha = cvReadRealByName( fs, fnode, "alpha" );\r
- node->node_risk = cvReadRealByName( fs, fnode, "node_risk" );\r
- node->tree_risk = cvReadRealByName( fs, fnode, "tree_risk" );\r
- node->tree_error = cvReadRealByName( fs, fnode, "tree_error" );\r
-\r
- splits = cvGetFileNodeByName( fs, fnode, "splits" );\r
- if( splits )\r
- {\r
- CvSeqReader reader;\r
- CvDTreeSplit* last_split = 0;\r
-\r
- if( CV_NODE_TYPE(splits->tag) != CV_NODE_SEQ )\r
- CV_ERROR( CV_StsParseError, "splits tag must stored as a sequence" );\r
-\r
- cvStartReadSeq( splits->data.seq, &reader );\r
- for( i = 0; i < reader.seq->total; i++ )\r
- {\r
- CvDTreeSplit* split;\r
- CV_CALL( split = read_split( fs, (CvFileNode*)reader.ptr ));\r
- if( !last_split )\r
- node->split = last_split = split;\r
- else\r
- last_split = last_split->next = split;\r
-\r
- CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );\r
- }\r
- }\r
-\r
- __END__;\r
-\r
- return node;\r
-}\r
-\r
-\r
-void CvDTree::read_tree_nodes( CvFileStorage* fs, CvFileNode* fnode )\r
-{\r
- CV_FUNCNAME( "CvDTree::read_tree_nodes" );\r
-\r
- __BEGIN__;\r
-\r
- CvSeqReader reader;\r
- CvDTreeNode _root;\r
- CvDTreeNode* parent = &_root;\r
- int i;\r
- parent->left = parent->right = parent->parent = 0;\r
-\r
- cvStartReadSeq( fnode->data.seq, &reader );\r
-\r
- for( i = 0; i < reader.seq->total; i++ )\r
- {\r
- CvDTreeNode* node;\r
-\r
- CV_CALL( node = read_node( fs, (CvFileNode*)reader.ptr, parent != &_root ? parent : 0 ));\r
- if( !parent->left )\r
- parent->left = node;\r
- else\r
- parent->right = node;\r
- if( node->split )\r
- parent = node;\r
- else\r
- {\r
- while( parent && parent->right )\r
- parent = parent->parent;\r
- }\r
-\r
- CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );\r
- }\r
-\r
- root = _root.left;\r
-\r
- __END__;\r
-}\r
-\r
-\r
-void CvDTree::read( CvFileStorage* fs, CvFileNode* fnode )\r
-{\r
- CvDTreeTrainData* _data = new CvDTreeTrainData();\r
- _data->read_params( fs, fnode );\r
-\r
- read( fs, fnode, _data );\r
- get_var_importance();\r
-}\r
-\r
-\r
-// a special entry point for reading weak decision trees from the tree ensembles\r
-void CvDTree::read( CvFileStorage* fs, CvFileNode* node, CvDTreeTrainData* _data )\r
-{\r
- CV_FUNCNAME( "CvDTree::read" );\r
-\r
- __BEGIN__;\r
-\r
- CvFileNode* tree_nodes;\r
-\r
- clear();\r
- data = _data;\r
-\r
- tree_nodes = cvGetFileNodeByName( fs, node, "nodes" );\r
- if( !tree_nodes || CV_NODE_TYPE(tree_nodes->tag) != CV_NODE_SEQ )\r
- CV_ERROR( CV_StsParseError, "nodes tag is missing" );\r
-\r
- pruned_tree_idx = cvReadIntByName( fs, node, "best_tree_idx", -1 );\r
- read_tree_nodes( fs, tree_nodes );\r
-\r
- __END__;\r
-}\r
-\r
-/* End of file. */\r
+/*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 "_ml.h"
+#include <ctype.h>
+
+using namespace cv;
+
+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;
+
+CvDTreeTrainData::CvDTreeTrainData()
+{
+ var_idx = var_type = cat_count = cat_ofs = cat_map =
+ priors = priors_mult = counts = buf = direction = split_buf = responses_copy = 0;
+ tree_storage = temp_storage = 0;
+
+ clear();
+}
+
+
+CvDTreeTrainData::CvDTreeTrainData( 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 )
+{
+ var_idx = var_type = cat_count = cat_ofs = cat_map =
+ priors = priors_mult = counts = buf = direction = split_buf = responses_copy = 0;
+
+ tree_storage = temp_storage = 0;
+
+ set_data( _train_data, _tflag, _responses, _var_idx, _sample_idx,
+ _var_type, _missing_mask, _params, _shared, _add_labels );
+}
+
+
+CvDTreeTrainData::~CvDTreeTrainData()
+{
+ clear();
+}
+
+
+bool CvDTreeTrainData::set_params( const CvDTreeParams& _params )
+{
+ bool ok = false;
+
+ CV_FUNCNAME( "CvDTreeTrainData::set_params" );
+
+ __BEGIN__;
+
+ // set parameters
+ params = _params;
+
+ if( params.max_categories < 2 )
+ CV_ERROR( CV_StsOutOfRange, "params.max_categories should be >= 2" );
+ params.max_categories = MIN( params.max_categories, 15 );
+
+ if( params.max_depth < 0 )
+ CV_ERROR( CV_StsOutOfRange, "params.max_depth should be >= 0" );
+ params.max_depth = MIN( params.max_depth, 25 );
+
+ params.min_sample_count = MAX(params.min_sample_count,1);
+
+ if( params.cv_folds < 0 )
+ CV_ERROR( CV_StsOutOfRange,
+ "params.cv_folds should be =0 (the tree is not pruned) "
+ "or n>0 (tree is pruned using n-fold cross-validation)" );
+
+ if( params.cv_folds == 1 )
+ params.cv_folds = 0;
+
+ if( params.regression_accuracy < 0 )
+ CV_ERROR( CV_StsOutOfRange, "params.regression_accuracy should be >= 0" );
+
+ ok = true;
+
+ __END__;
+
+ return ok;
+}
+
+#define CV_CMP_NUM_PTR(a,b) (*(a) < *(b))
+static CV_IMPLEMENT_QSORT_EX( icvSortIntPtr, int*, CV_CMP_NUM_PTR, int )
+static CV_IMPLEMENT_QSORT_EX( icvSortDblPtr, double*, CV_CMP_NUM_PTR, int )
+
+#define CV_CMP_NUM_IDX(i,j) (aux[i] < aux[j])
+static CV_IMPLEMENT_QSORT_EX( icvSortIntAux, int, CV_CMP_NUM_IDX, const float* )
+static CV_IMPLEMENT_QSORT_EX( icvSortUShAux, unsigned short, CV_CMP_NUM_IDX, const float* )
+
+#define CV_CMP_PAIRS(a,b) (*((a).i) < *((b).i))
+static CV_IMPLEMENT_QSORT_EX( icvSortPairs, CvPair16u32s, CV_CMP_PAIRS, int )
+
+void CvDTreeTrainData::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( "CvDTreeTrainData::set_data" );
+
+ __BEGIN__;
+
+ int sample_all = 0, r_type = 0, 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( _update_data && data_root )
+ {
+ data = new CvDTreeTrainData( _train_data, _tflag, _responses, _var_idx,
+ _sample_idx, _var_type, _missing_mask, _params, _shared, _add_labels );
+
+ // compare new and old train data
+ if( !(data->var_count == var_count &&
+ cvNorm( data->var_type, var_type, CV_C ) < FLT_EPSILON &&
+ cvNorm( data->cat_count, cat_count, CV_C ) < FLT_EPSILON &&
+ cvNorm( data->cat_map, cat_map, CV_C ) < FLT_EPSILON) )
+ CV_ERROR( CV_StsBadArg,
+ "The new training data must have the same types and the input and output variables "
+ "and the same categories for categorical variables" );
+
+ cvReleaseMat( &priors );
+ cvReleaseMat( &priors_mult );
+ cvReleaseMat( &buf );
+ cvReleaseMat( &direction );
+ cvReleaseMat( &split_buf );
+ cvReleaseMemStorage( &temp_storage );
+
+ priors = data->priors; data->priors = 0;
+ priors_mult = data->priors_mult; data->priors_mult = 0;
+ buf = data->buf; data->buf = 0;
+ buf_count = data->buf_count; buf_size = data->buf_size;
+ sample_count = data->sample_count;
+
+ direction = data->direction; data->direction = 0;
+ split_buf = data->split_buf; data->split_buf = 0;
+ temp_storage = data->temp_storage; data->temp_storage = 0;
+ nv_heap = data->nv_heap; cv_heap = data->cv_heap;
+
+ data_root = new_node( 0, sample_count, 0, 0 );
+ EXIT;
+ }
+
+ clear();
+
+ var_all = 0;
+ rng = cvRNG(-1);
+
+ 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;
+
+ 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;
+ }
+
+ is_buf_16u = false;
+ if ( sample_count < 65536 )
+ is_buf_16u = true;
+
+ 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" );
+
+
+ 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++ )
+ {
+ var_type->data.i[vi] = var_type0->data.ptr[vi] == 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 = var_count + (is_classifier ? 1 : 0) // for responses class_labels
+ + (have_labels ? 1 : 0); // for cv_labels
+
+ buf_size = (work_var_count + 1 /*for sample_indices*/) * 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 + vi*sample_count);
+ else
+ idst = buf->data.i + vi*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( fabs(t - val) > FLT_EPSILON )
+ {
+ 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 -1
+ 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
+ {
+ if (is_buf_16u)
+ udst = (unsigned short*)(buf->data.s + vi*sample_count);
+ else
+ idst = buf->data.i + vi*sample_count;
+
+ 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 (is_buf_16u)
+ udst[i] = (unsigned short)i;
+ else
+ idst[i] = i;
+ _fdst[i] = val;
+
+ }
+ if (is_buf_16u)
+ icvSortUShAux( udst, num_valid, _fdst);
+ else
+ icvSortIntAux( idst, /*or num_valid?\*/ sample_count, _fdst );
+ }
+
+ 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 + work_var_count*sample_count);
+ else
+ idst = buf->data.i + 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;
+ CvRNG* r = &rng;
+
+ 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 = cvRandInt(r) % sample_count;
+ int b = cvRandInt(r) % 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 = cvRandInt(r) % sample_count;
+ int b = cvRandInt(r) % 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 );
+ cvFree( &pair16u32s_ptr);
+ cvReleaseMat( &var_type0 );
+ cvReleaseMat( &sample_indices );
+ cvReleaseMat( &tmp_map );
+}
+
+void CvDTreeTrainData::do_responses_copy()
+{
+ responses_copy = cvCreateMat( responses->rows, responses->cols, responses->type );
+ cvCopy( responses, responses_copy);
+ responses = responses_copy;
+}
+
+CvDTreeNode* CvDTreeTrainData::subsample_data( const CvMat* _subsample_idx )
+{
+ CvDTreeNode* root = 0;
+ CvMat* isubsample_idx = 0;
+ CvMat* subsample_co = 0;
+
+ CV_FUNCNAME( "CvDTreeTrainData::subsample_data" );
+
+ __BEGIN__;
+
+ if( !data_root )
+ CV_ERROR( CV_StsError, "No training data has been set" );
+
+ if( _subsample_idx )
+ CV_CALL( isubsample_idx = cvPreprocessIndexArray( _subsample_idx, sample_count ));
+
+ if( !isubsample_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
+ {
+ int* sidx = isubsample_idx->data.i;
+ // co - array of count/offset pairs (to handle duplicated values in _subsample_idx)
+ int* co, cur_ofs = 0;
+ int vi, i;
+ int work_var_count = get_work_var_count();
+ int count = isubsample_idx->rows + isubsample_idx->cols - 1;
+
+ root = new_node( 0, count, 1, 0 );
+
+ CV_CALL( subsample_co = cvCreateMat( 1, sample_count*2, CV_32SC1 ));
+ cvZero( subsample_co );
+ co = subsample_co->data.i;
+ for( i = 0; i < count; i++ )
+ co[sidx[i]*2]++;
+ for( i = 0; i < sample_count; i++ )
+ {
+ if( co[i*2] )
+ {
+ co[i*2+1] = cur_ofs;
+ cur_ofs += co[i*2];
+ }
+ else
+ co[i*2+1] = -1;
+ }
+
+ cv::AutoBuffer<uchar> inn_buf(sample_count*(2*sizeof(int) + sizeof(float)));
+ for( vi = 0; vi < work_var_count; vi++ )
+ {
+ int ci = get_var_type(vi);
+
+ if( ci >= 0 || vi >= var_count )
+ {
+ int num_valid = 0;
+ const int* src = get_cat_var_data( data_root, vi, (int*)(uchar*)inn_buf );
+
+ if (is_buf_16u)
+ {
+ unsigned short* udst = (unsigned short*)(buf->data.s + root->buf_idx*buf->cols +
+ vi*sample_count + root->offset);
+ for( i = 0; i < count; i++ )
+ {
+ int val = src[sidx[i]];
+ udst[i] = (unsigned short)val;
+ num_valid += val >= 0;
+ }
+ }
+ else
+ {
+ int* idst = buf->data.i + root->buf_idx*buf->cols +
+ vi*sample_count + root->offset;
+ for( i = 0; i < count; i++ )
+ {
+ int val = src[sidx[i]];
+ idst[i] = val;
+ num_valid += val >= 0;
+ }
+ }
+
+ if( vi < var_count )
+ root->set_num_valid(vi, num_valid);
+ }
+ else
+ {
+ int *src_idx_buf = (int*)(uchar*)inn_buf;
+ float *src_val_buf = (float*)(src_idx_buf + sample_count);
+ int* sample_indices_buf = (int*)(src_val_buf + sample_count);
+ const int* src_idx = 0;
+ const float* src_val = 0;
+ get_ord_var_data( data_root, vi, src_val_buf, src_idx_buf, &src_val, &src_idx, sample_indices_buf );
+ int j = 0, idx, count_i;
+ int num_valid = data_root->get_num_valid(vi);
+
+ if (is_buf_16u)
+ {
+ unsigned short* udst_idx = (unsigned short*)(buf->data.s + root->buf_idx*buf->cols +
+ vi*sample_count + data_root->offset);
+ for( i = 0; i < num_valid; i++ )
+ {
+ idx = src_idx[i];
+ count_i = co[idx*2];
+ if( count_i )
+ for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )
+ udst_idx[j] = (unsigned short)cur_ofs;
+ }
+
+ root->set_num_valid(vi, j);
+
+ for( ; i < sample_count; i++ )
+ {
+ idx = src_idx[i];
+ count_i = co[idx*2];
+ if( count_i )
+ for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )
+ udst_idx[j] = (unsigned short)cur_ofs;
+ }
+ }
+ else
+ {
+ int* idst_idx = buf->data.i + root->buf_idx*buf->cols +
+ vi*sample_count + root->offset;
+ for( i = 0; i < num_valid; i++ )
+ {
+ idx = src_idx[i];
+ count_i = co[idx*2];
+ if( count_i )
+ for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )
+ idst_idx[j] = cur_ofs;
+ }
+
+ root->set_num_valid(vi, j);
+
+ for( ; i < sample_count; i++ )
+ {
+ idx = src_idx[i];
+ count_i = co[idx*2];
+ if( count_i )
+ for( cur_ofs = co[idx*2+1]; count_i > 0; count_i--, j++, cur_ofs++ )
+ idst_idx[j] = cur_ofs;
+ }
+ }
+ }
+ }
+ // sample indices subsampling
+ const int* sample_idx_src = get_sample_indices(data_root, (int*)(uchar*)inn_buf);
+ if (is_buf_16u)
+ {
+ unsigned short* sample_idx_dst = (unsigned short*)(buf->data.s + root->buf_idx*buf->cols +
+ get_work_var_count()*sample_count + root->offset);
+ for (i = 0; i < count; i++)
+ sample_idx_dst[i] = (unsigned short)sample_idx_src[sidx[i]];
+ }
+ else
+ {
+ int* sample_idx_dst = buf->data.i + root->buf_idx*buf->cols +
+ get_work_var_count()*sample_count + root->offset;
+ for (i = 0; i < count; i++)
+ sample_idx_dst[i] = sample_idx_src[sidx[i]];
+ }
+ }
+
+ __END__;
+
+ cvReleaseMat( &isubsample_idx );
+ cvReleaseMat( &subsample_co );
+
+ return root;
+}
+
+
+void CvDTreeTrainData::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_FUNCNAME( "CvDTreeTrainData::get_vectors" );
+
+ __BEGIN__;
+
+ int i, vi, total = sample_count, count = total, cur_ofs = 0;
+ int* sidx = 0;
+ int* co = 0;
+
+ cv::AutoBuffer<uchar> inn_buf(sample_count*(2*sizeof(int) + sizeof(float)));
+ 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;
+ const int* src = get_cat_var_data(data_root, vi, (int*)(uchar*)inn_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
+ {
+ float* dst = values + vi;
+ uchar* m = missing ? missing + vi : 0;
+ int count1 = data_root->get_num_valid(vi);
+ float *src_val_buf = (float*)(uchar*)inn_buf;
+ int* src_idx_buf = (int*)(src_val_buf + sample_count);
+ int* sample_indices_buf = src_idx_buf + sample_count;
+ const float *src_val = 0;
+ const int* src_idx = 0;
+ get_ord_var_data(data_root, vi, src_val_buf, src_idx_buf, &src_val, &src_idx, sample_indices_buf);
+
+ for( i = 0; i < count1; i++ )
+ {
+ int idx = src_idx[i];
+ int count_i = 1;
+ if( co )
+ {
+ count_i = co[idx*2];
+ cur_ofs = co[idx*2+1];
+ }
+ else
+ cur_ofs = idx*var_count;
+ if( count_i )
+ {
+ float val = src_val[i];
+ for( ; count_i > 0; count_i--, cur_ofs += var_count )
+ {
+ dst[cur_ofs] = val;
+ if( m )
+ m[cur_ofs] = 0;
+ }
+ }
+ }
+ }
+ }
+
+ // copy responses
+ if( responses )
+ {
+ if( is_classifier )
+ {
+ const int* src = get_class_labels(data_root, (int*)(uchar*)inn_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* val_buf = (float*)(uchar*)inn_buf;
+ int* sample_idx_buf = (int*)(val_buf + sample_count);
+ const float* _values = get_ord_responses(data_root, val_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* CvDTreeTrainData::new_node( CvDTreeNode* parent, int count,
+ int storage_idx, int offset )
+{
+ CvDTreeNode* node = (CvDTreeNode*)cvSetNew( node_heap );
+
+ node->sample_count = count;
+ node->depth = parent ? parent->depth + 1 : 0;
+ node->parent = parent;
+ node->left = node->right = 0;
+ node->split = 0;
+ node->value = 0;
+ node->class_idx = 0;
+ node->maxlr = 0.;
+
+ node->buf_idx = storage_idx;
+ node->offset = offset;
+ if( nv_heap )
+ node->num_valid = (int*)cvSetNew( nv_heap );
+ else
+ node->num_valid = 0;
+ node->alpha = node->node_risk = node->tree_risk = node->tree_error = 0.;
+ node->complexity = 0;
+
+ if( params.cv_folds > 0 && cv_heap )
+ {
+ int cv_n = params.cv_folds;
+ node->Tn = INT_MAX;
+ node->cv_Tn = (int*)cvSetNew( cv_heap );
+ node->cv_node_risk = (double*)cvAlignPtr(node->cv_Tn + cv_n, sizeof(double));
+ node->cv_node_error = node->cv_node_risk + cv_n;
+ }
+ else
+ {
+ node->Tn = 0;
+ node->cv_Tn = 0;
+ node->cv_node_risk = 0;
+ node->cv_node_error = 0;
+ }
+
+ return node;
+}
+
+
+CvDTreeSplit* CvDTreeTrainData::new_split_ord( int vi, float cmp_val,
+ int split_point, int inversed, float quality )
+{
+ CvDTreeSplit* split = (CvDTreeSplit*)cvSetNew( split_heap );
+ split->var_idx = vi;
+ split->condensed_idx = INT_MIN;
+ split->ord.c = cmp_val;
+ split->ord.split_point = split_point;
+ split->inversed = inversed;
+ split->quality = quality;
+ split->next = 0;
+
+ return split;
+}
+
+
+CvDTreeSplit* CvDTreeTrainData::new_split_cat( int vi, float quality )
+{
+ CvDTreeSplit* split = (CvDTreeSplit*)cvSetNew( split_heap );
+ int i, n = (max_c_count + 31)/32;
+
+ split->var_idx = vi;
+ split->condensed_idx = INT_MIN;
+ split->inversed = 0;
+ split->quality = quality;
+ for( i = 0; i < n; i++ )
+ split->subset[i] = 0;
+ split->next = 0;
+
+ return split;
+}
+
+
+void CvDTreeTrainData::free_node( CvDTreeNode* node )
+{
+ CvDTreeSplit* split = node->split;
+ free_node_data( node );
+ while( split )
+ {
+ CvDTreeSplit* next = split->next;
+ cvSetRemoveByPtr( split_heap, split );
+ split = next;
+ }
+ node->split = 0;
+ cvSetRemoveByPtr( node_heap, node );
+}
+
+
+void CvDTreeTrainData::free_node_data( CvDTreeNode* node )
+{
+ if( node->num_valid )
+ {
+ cvSetRemoveByPtr( nv_heap, node->num_valid );
+ node->num_valid = 0;
+ }
+ // do not free cv_* fields, as all the cross-validation related data is released at once.
+}
+
+
+void CvDTreeTrainData::free_train_data()
+{
+ cvReleaseMat( &counts );
+ cvReleaseMat( &buf );
+ cvReleaseMat( &direction );
+ cvReleaseMat( &split_buf );
+ cvReleaseMemStorage( &temp_storage );
+ cvReleaseMat( &responses_copy );
+ cv_heap = nv_heap = 0;
+}
+
+
+void CvDTreeTrainData::clear()
+{
+ free_train_data();
+
+ cvReleaseMemStorage( &tree_storage );
+
+ cvReleaseMat( &var_idx );
+ cvReleaseMat( &var_type );
+ cvReleaseMat( &cat_count );
+ cvReleaseMat( &cat_ofs );
+ cvReleaseMat( &cat_map );
+ cvReleaseMat( &priors );
+ cvReleaseMat( &priors_mult );
+
+ node_heap = split_heap = 0;
+
+ sample_count = var_all = var_count = max_c_count = ord_var_count = cat_var_count = 0;
+ have_labels = have_priors = is_classifier = false;
+
+ buf_count = buf_size = 0;
+ shared = false;
+
+ data_root = 0;
+
+ rng = cvRNG(-1);
+}
+
+
+int CvDTreeTrainData::get_num_classes() const
+{
+ return is_classifier ? cat_count->data.i[cat_var_count] : 0;
+}
+
+
+int CvDTreeTrainData::get_var_type(int vi) const
+{
+ return var_type->data.i[vi];
+}
+
+void CvDTreeTrainData::get_ord_var_data( CvDTreeNode* n, int vi, float* ord_values_buf, int* sorted_indices_buf,
+ const float** ord_values, const int** sorted_indices, int* sample_indices_buf )
+{
+ int vidx = var_idx ? var_idx->data.i[vi] : vi;
+ int node_sample_count = n->sample_count;
+ int td_step = train_data->step/CV_ELEM_SIZE(train_data->type);
+
+ const int* sample_indices = get_sample_indices(n, sample_indices_buf);
+
+ if( !is_buf_16u )
+ *sorted_indices = buf->data.i + n->buf_idx*buf->cols +
+ vi*sample_count + n->offset;
+ else {
+ const unsigned short* short_indices = (const unsigned short*)(buf->data.s + n->buf_idx*buf->cols +
+ vi*sample_count + n->offset );
+ for( int i = 0; i < node_sample_count; i++ )
+ sorted_indices_buf[i] = short_indices[i];
+ *sorted_indices = sorted_indices_buf;
+ }
+
+ if( tflag == CV_ROW_SAMPLE )
+ {
+ for( int i = 0; i < node_sample_count &&
+ ((((*sorted_indices)[i] >= 0) && !is_buf_16u) || (((*sorted_indices)[i] != 65535) && is_buf_16u)); i++ )
+ {
+ int idx = (*sorted_indices)[i];
+ idx = sample_indices[idx];
+ ord_values_buf[i] = *(train_data->data.fl + idx * td_step + vidx);
+ }
+ }
+ else
+ for( int i = 0; i < node_sample_count &&
+ ((((*sorted_indices)[i] >= 0) && !is_buf_16u) || (((*sorted_indices)[i] != 65535) && is_buf_16u)); i++ )
+ {
+ int idx = (*sorted_indices)[i];
+ idx = sample_indices[idx];
+ ord_values_buf[i] = *(train_data->data.fl + vidx* td_step + idx);
+ }
+
+ *ord_values = ord_values_buf;
+}
+
+
+const int* CvDTreeTrainData::get_class_labels( CvDTreeNode* n, int* labels_buf )
+{
+ if (is_classifier)
+ return get_cat_var_data( n, var_count, labels_buf);
+ return 0;
+}
+
+const int* CvDTreeTrainData::get_sample_indices( CvDTreeNode* n, int* indices_buf )
+{
+ return get_cat_var_data( n, get_work_var_count(), indices_buf );
+}
+
+const float* CvDTreeTrainData::get_ord_responses( CvDTreeNode* n, float* values_buf, int*sample_indices_buf )
+{
+ int sample_count = n->sample_count;
+ int r_step = CV_IS_MAT_CONT(responses->type) ? 1 : responses->step/CV_ELEM_SIZE(responses->type);
+ const int* indices = get_sample_indices(n, sample_indices_buf);
+
+ for( int i = 0; i < sample_count &&
+ (((indices[i] >= 0) && !is_buf_16u) || ((indices[i] != 65535) && is_buf_16u)); i++ )
+ {
+ int idx = indices[i];
+ values_buf[i] = *(responses->data.fl + idx * r_step);
+ }
+
+ return values_buf;
+}
+
+
+const int* CvDTreeTrainData::get_cv_labels( CvDTreeNode* n, int* labels_buf )
+{
+ if (have_labels)
+ return get_cat_var_data( n, get_work_var_count()- 1, labels_buf);
+ return 0;
+}
+
+
+const int* CvDTreeTrainData::get_cat_var_data( CvDTreeNode* n, int vi, int* cat_values_buf)
+{
+ const int* cat_values = 0;
+ if( !is_buf_16u )
+ cat_values = buf->data.i + n->buf_idx*buf->cols +
+ vi*sample_count + n->offset;
+ else {
+ const unsigned short* short_values = (const unsigned short*)(buf->data.s + n->buf_idx*buf->cols +
+ vi*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;
+}
+
+
+int CvDTreeTrainData::get_child_buf_idx( CvDTreeNode* n )
+{
+ int idx = n->buf_idx + 1;
+ if( idx >= buf_count )
+ idx = shared ? 1 : 0;
+ return idx;
+}
+
+
+void CvDTreeTrainData::write_params( CvFileStorage* fs ) const
+{
+ CV_FUNCNAME( "CvDTreeTrainData::write_params" );
+
+ __BEGIN__;
+
+ int vi, vcount = var_count;
+
+ cvWriteInt( fs, "is_classifier", is_classifier ? 1 : 0 );
+ cvWriteInt( fs, "var_all", var_all );
+ cvWriteInt( fs, "var_count", var_count );
+ cvWriteInt( fs, "ord_var_count", ord_var_count );
+ cvWriteInt( fs, "cat_var_count", cat_var_count );
+
+ cvStartWriteStruct( fs, "training_params", CV_NODE_MAP );
+ cvWriteInt( fs, "use_surrogates", params.use_surrogates ? 1 : 0 );
+
+ if( is_classifier )
+ {
+ cvWriteInt( fs, "max_categories", params.max_categories );
+ }
+ else
+ {
+ cvWriteReal( fs, "regression_accuracy", params.regression_accuracy );
+ }
+
+ cvWriteInt( fs, "max_depth", params.max_depth );
+ cvWriteInt( fs, "min_sample_count", params.min_sample_count );
+ cvWriteInt( fs, "cross_validation_folds", params.cv_folds );
+
+ if( params.cv_folds > 1 )
+ {
+ cvWriteInt( fs, "use_1se_rule", params.use_1se_rule ? 1 : 0 );
+ cvWriteInt( fs, "truncate_pruned_tree", params.truncate_pruned_tree ? 1 : 0 );
+ }
+
+ if( priors )
+ cvWrite( fs, "priors", priors );
+
+ cvEndWriteStruct( fs );
+
+ if( var_idx )
+ cvWrite( fs, "var_idx", var_idx );
+
+ cvStartWriteStruct( fs, "var_type", CV_NODE_SEQ+CV_NODE_FLOW );
+
+ for( vi = 0; vi < vcount; vi++ )
+ cvWriteInt( fs, 0, var_type->data.i[vi] >= 0 );
+
+ cvEndWriteStruct( fs );
+
+ if( cat_count && (cat_var_count > 0 || is_classifier) )
+ {
+ CV_ASSERT( cat_count != 0 );
+ cvWrite( fs, "cat_count", cat_count );
+ cvWrite( fs, "cat_map", cat_map );
+ }
+
+ __END__;
+}
+
+
+void CvDTreeTrainData::read_params( CvFileStorage* fs, CvFileNode* node )
+{
+ CV_FUNCNAME( "CvDTreeTrainData::read_params" );
+
+ __BEGIN__;
+
+ CvFileNode *tparams_node, *vartype_node;
+ CvSeqReader reader;
+ int vi, max_split_size, tree_block_size;
+
+ is_classifier = (cvReadIntByName( fs, node, "is_classifier" ) != 0);
+ var_all = cvReadIntByName( fs, node, "var_all" );
+ var_count = cvReadIntByName( fs, node, "var_count", var_all );
+ cat_var_count = cvReadIntByName( fs, node, "cat_var_count" );
+ ord_var_count = cvReadIntByName( fs, node, "ord_var_count" );
+
+ tparams_node = cvGetFileNodeByName( fs, node, "training_params" );
+
+ if( tparams_node ) // training parameters are not necessary
+ {
+ params.use_surrogates = cvReadIntByName( fs, tparams_node, "use_surrogates", 1 ) != 0;
+
+ if( is_classifier )
+ {
+ params.max_categories = cvReadIntByName( fs, tparams_node, "max_categories" );
+ }
+ else
+ {
+ params.regression_accuracy =
+ (float)cvReadRealByName( fs, tparams_node, "regression_accuracy" );
+ }
+
+ params.max_depth = cvReadIntByName( fs, tparams_node, "max_depth" );
+ params.min_sample_count = cvReadIntByName( fs, tparams_node, "min_sample_count" );
+ params.cv_folds = cvReadIntByName( fs, tparams_node, "cross_validation_folds" );
+
+ if( params.cv_folds > 1 )
+ {
+ params.use_1se_rule = cvReadIntByName( fs, tparams_node, "use_1se_rule" ) != 0;
+ params.truncate_pruned_tree =
+ cvReadIntByName( fs, tparams_node, "truncate_pruned_tree" ) != 0;
+ }
+
+ priors = (CvMat*)cvReadByName( fs, tparams_node, "priors" );
+ if( priors )
+ {
+ if( !CV_IS_MAT(priors) )
+ CV_ERROR( CV_StsParseError, "priors must stored as a matrix" );
+ priors_mult = cvCloneMat( priors );
+ }
+ }
+
+ CV_CALL( var_idx = (CvMat*)cvReadByName( fs, node, "var_idx" ));
+ if( var_idx )
+ {
+ if( !CV_IS_MAT(var_idx) ||
+ (var_idx->cols != 1 && var_idx->rows != 1) ||
+ var_idx->cols + var_idx->rows - 1 != var_count ||
+ CV_MAT_TYPE(var_idx->type) != CV_32SC1 )
+ CV_ERROR( CV_StsParseError,
+ "var_idx (if exist) must be valid 1d integer vector containing <var_count> elements" );
+
+ for( vi = 0; vi < var_count; vi++ )
+ if( (unsigned)var_idx->data.i[vi] >= (unsigned)var_all )
+ CV_ERROR( CV_StsOutOfRange, "some of var_idx elements are out of range" );
+ }
+
+ ////// read var type
+ CV_CALL( var_type = cvCreateMat( 1, var_count + 2, CV_32SC1 ));
+
+ cat_var_count = 0;
+ ord_var_count = -1;
+ vartype_node = cvGetFileNodeByName( fs, node, "var_type" );
+
+ if( vartype_node && CV_NODE_TYPE(vartype_node->tag) == CV_NODE_INT && var_count == 1 )
+ var_type->data.i[0] = vartype_node->data.i ? cat_var_count++ : ord_var_count--;
+ else
+ {
+ if( !vartype_node || CV_NODE_TYPE(vartype_node->tag) != CV_NODE_SEQ ||
+ vartype_node->data.seq->total != var_count )
+ CV_ERROR( CV_StsParseError, "var_type must exist and be a sequence of 0's and 1's" );
+
+ cvStartReadSeq( vartype_node->data.seq, &reader );
+
+ for( vi = 0; vi < var_count; vi++ )
+ {
+ CvFileNode* n = (CvFileNode*)reader.ptr;
+ if( CV_NODE_TYPE(n->tag) != CV_NODE_INT || (n->data.i & ~1) )
+ CV_ERROR( CV_StsParseError, "var_type must exist and be a sequence of 0's and 1's" );
+ var_type->data.i[vi] = n->data.i ? cat_var_count++ : ord_var_count--;
+ CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );
+ }
+ }
+ var_type->data.i[var_count] = cat_var_count;
+
+ ord_var_count = ~ord_var_count;
+ if( cat_var_count != cat_var_count || ord_var_count != ord_var_count )
+ CV_ERROR( CV_StsParseError, "var_type is inconsistent with cat_var_count and ord_var_count" );
+ //////
+
+ if( cat_var_count > 0 || is_classifier )
+ {
+ int ccount, total_c_count = 0;
+ CV_CALL( cat_count = (CvMat*)cvReadByName( fs, node, "cat_count" ));
+ CV_CALL( cat_map = (CvMat*)cvReadByName( fs, node, "cat_map" ));
+
+ if( !CV_IS_MAT(cat_count) || !CV_IS_MAT(cat_map) ||
+ (cat_count->cols != 1 && cat_count->rows != 1) ||
+ CV_MAT_TYPE(cat_count->type) != CV_32SC1 ||
+ cat_count->cols + cat_count->rows - 1 != cat_var_count + is_classifier ||
+ (cat_map->cols != 1 && cat_map->rows != 1) ||
+ CV_MAT_TYPE(cat_map->type) != CV_32SC1 )
+ CV_ERROR( CV_StsParseError,
+ "Both cat_count and cat_map must exist and be valid 1d integer vectors of an appropriate size" );
+
+ ccount = cat_var_count + is_classifier;
+
+ CV_CALL( cat_ofs = cvCreateMat( 1, ccount + 1, CV_32SC1 ));
+ cat_ofs->data.i[0] = 0;
+ max_c_count = 1;
+
+ for( vi = 0; vi < ccount; vi++ )
+ {
+ int val = cat_count->data.i[vi];
+ if( val <= 0 )
+ CV_ERROR( CV_StsOutOfRange, "some of cat_count elements are out of range" );
+ max_c_count = MAX( max_c_count, val );
+ cat_ofs->data.i[vi+1] = total_c_count += val;
+ }
+
+ if( cat_map->cols + cat_map->rows - 1 != total_c_count )
+ CV_ERROR( CV_StsBadSize,
+ "cat_map vector length is not equal to the total number of categories in all categorical vars" );
+ }
+
+ max_split_size = cvAlign(sizeof(CvDTreeSplit) +
+ (MAX(0,max_c_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[0]),
+ sizeof(CvDTreeNode), tree_storage ));
+ CV_CALL( split_heap = cvCreateSet( 0, sizeof(split_heap[0]),
+ max_split_size, tree_storage ));
+
+ __END__;
+}
+
+/////////////////////// Decision Tree /////////////////////////
+
+CvDTree::CvDTree()
+{
+ data = 0;
+ var_importance = 0;
+ default_model_name = "my_tree";
+
+ clear();
+}
+
+
+void CvDTree::clear()
+{
+ cvReleaseMat( &var_importance );
+ if( data )
+ {
+ if( !data->shared )
+ delete data;
+ else
+ free_tree();
+ data = 0;
+ }
+ root = 0;
+ pruned_tree_idx = -1;
+}
+
+
+CvDTree::~CvDTree()
+{
+ clear();
+}
+
+
+const CvDTreeNode* CvDTree::get_root() const
+{
+ return root;
+}
+
+
+int CvDTree::get_pruned_tree_idx() const
+{
+ return pruned_tree_idx;
+}
+
+
+CvDTreeTrainData* CvDTree::get_data()
+{
+ return data;
+}
+
+
+bool CvDTree::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, CvDTreeParams _params )
+{
+ bool result = false;
+
+ CV_FUNCNAME( "CvDTree::train" );
+
+ __BEGIN__;
+
+ clear();
+ data = new CvDTreeTrainData( _train_data, _tflag, _responses,
+ _var_idx, _sample_idx, _var_type,
+ _missing_mask, _params, false );
+ CV_CALL( result = do_train(0) );
+
+ __END__;
+
+ return result;
+}
+
+bool CvDTree::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, CvDTreeParams _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);
+}
+
+
+bool CvDTree::train( CvMLData* _data, CvDTreeParams _params )
+{
+ bool result = false;
+
+ CV_FUNCNAME( "CvDTree::train" );
+
+ __BEGIN__;
+
+ const CvMat* values = _data->get_values();
+ const CvMat* response = _data->get_responses();
+ const CvMat* missing = _data->get_missing();
+ const CvMat* var_types = _data->get_var_types();
+ const CvMat* train_sidx = _data->get_train_sample_idx();
+ const CvMat* var_idx = _data->get_var_idx();
+
+ CV_CALL( result = train( values, CV_ROW_SAMPLE, response, var_idx,
+ train_sidx, var_types, missing, _params ) );
+
+ __END__;
+
+ return result;
+}
+
+bool CvDTree::train( CvDTreeTrainData* _data, const CvMat* _subsample_idx )
+{
+ bool result = false;
+
+ CV_FUNCNAME( "CvDTree::train" );
+
+ __BEGIN__;
+
+ clear();
+ data = _data;
+ data->shared = true;
+ CV_CALL( result = do_train(_subsample_idx));
+
+ __END__;
+
+ return result;
+}
+
+
+bool CvDTree::do_train( const CvMat* _subsample_idx )
+{
+ bool result = false;
+
+ CV_FUNCNAME( "CvDTree::do_train" );
+
+ __BEGIN__;
+
+ root = data->subsample_data( _subsample_idx );
+
+ CV_CALL( try_split_node(root));
+
+ if( data->params.cv_folds > 0 )
+ CV_CALL( prune_cv());
+
+ if( !data->shared )
+ data->free_train_data();
+
+ result = true;
+
+ __END__;
+
+ return result;
+}
+
+
+void CvDTree::try_split_node( CvDTreeNode* node )
+{
+ CvDTreeSplit* best_split = 0;
+ int i, n = node->sample_count, vi;
+ bool can_split = true;
+ double quality_scale;
+
+ calc_node_value( node );
+
+ if( node->sample_count <= data->params.min_sample_count ||
+ node->depth >= data->params.max_depth )
+ can_split = false;
+
+ if( can_split && data->is_classifier )
+ {
+ // check if we have a "pure" node,
+ // we assume that cls_count is filled by calc_node_value()
+ int* cls_count = data->counts->data.i;
+ int nz = 0, m = data->get_num_classes();
+ for( i = 0; i < m; i++ )
+ nz += cls_count[i] != 0;
+ if( nz == 1 ) // there is only one class
+ can_split = false;
+ }
+ else if( can_split )
+ {
+ if( sqrt(node->node_risk)/n < data->params.regression_accuracy )
+ can_split = false;
+ }
+
+ if( can_split )
+ {
+ best_split = find_best_split(node);
+ // TODO: check the split quality ...
+ node->split = best_split;
+ }
+ if( !can_split || !best_split )
+ {
+ data->free_node_data(node);
+ return;
+ }
+
+ quality_scale = calc_node_dir( node );
+ if( data->params.use_surrogates )
+ {
+ // find all the surrogate splits
+ // and sort them by their similarity to the primary one
+ for( vi = 0; vi < data->var_count; vi++ )
+ {
+ CvDTreeSplit* split;
+ int ci = data->get_var_type(vi);
+
+ if( vi == best_split->var_idx )
+ continue;
+
+ if( ci >= 0 )
+ split = find_surrogate_split_cat( node, vi );
+ else
+ split = find_surrogate_split_ord( node, vi );
+
+ if( split )
+ {
+ // insert the split
+ CvDTreeSplit* prev_split = node->split;
+ split->quality = (float)(split->quality*quality_scale);
+
+ while( prev_split->next &&
+ prev_split->next->quality > split->quality )
+ prev_split = prev_split->next;
+ split->next = prev_split->next;
+ prev_split->next = split;
+ }
+ }
+ }
+ split_node_data( node );
+ try_split_node( node->left );
+ try_split_node( node->right );
+}
+
+
+// calculate direction (left(-1),right(1),missing(0))
+// for each sample using the best split
+// the function returns scale coefficients for surrogate split quality factors.
+// the scale is applied to normalize surrogate split quality relatively to the
+// best (primary) split quality. That is, if a surrogate split is absolutely
+// identical to the primary split, its quality will be set to the maximum value =
+// quality of the primary split; otherwise, it will be lower.
+// besides, the function compute node->maxlr,
+// minimum possible quality (w/o considering the above mentioned scale)
+// for a surrogate split. Surrogate splits with quality less than node->maxlr
+// are not discarded.
+double CvDTree::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<int> inn_buf(n*(!data->have_priors ? 1 : 2));
+ int* labels_buf = (int*)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
+ {
+ int split_point = node->split->ord.split_point;
+ int n1 = node->get_num_valid(vi);
+ cv::AutoBuffer<uchar> inn_buf(n*(sizeof(int)*(data->have_priors ? 3 : 2) + sizeof(float)));
+ float* val_buf = (float*)(uchar*)inn_buf;
+ int* sorted_buf = (int*)(val_buf + n);
+ int* sample_idx_buf = sorted_buf + n;
+ const float* val = 0;
+ const int* sorted = 0;
+ data->get_ord_var_data( node, vi, val_buf, sorted_buf, &val, &sorted, sample_idx_buf);
+
+ assert( 0 <= split_point && split_point < n1-1 );
+
+ if( !data->have_priors )
+ {
+ for( i = 0; i <= split_point; i++ )
+ dir[sorted[i]] = (char)-1;
+ for( ; i < n1; i++ )
+ dir[sorted[i]] = (char)1;
+ for( ; i < n; i++ )
+ dir[sorted[i]] = (char)0;
+
+ L = split_point-1;
+ R = n1 - split_point + 1;
+ }
+ else
+ {
+ const double* priors = data->priors_mult->data.db;
+ int* responses_buf = sample_idx_buf + n;
+ const int* responses = data->get_class_labels(node, responses_buf);
+ L = R = 0;
+
+ for( i = 0; i <= split_point; i++ )
+ {
+ int idx = sorted[i];
+ double w = priors[responses[idx]];
+ dir[idx] = (char)-1;
+ L += w;
+ }
+
+ for( ; i < n1; i++ )
+ {
+ int idx = sorted[i];
+ double w = priors[responses[idx]];
+ dir[idx] = (char)1;
+ R += w;
+ }
+
+ for( ; i < n; i++ )
+ dir[sorted[i]] = (char)0;
+ }
+ }
+ node->maxlr = MAX( L, R );
+ return node->split->quality/(L + R);
+}
+
+
+namespace cv
+{
+
+DTreeBestSplitFinder::DTreeBestSplitFinder( CvDTree* _tree, CvDTreeNode* _node)
+{
+ tree = _tree;
+ node = _node;
+ splitSize = tree->get_data()->split_heap->elem_size;
+
+ bestSplit = (CvDTreeSplit*)(new char[splitSize]);
+ memset((CvDTreeSplit*)bestSplit, 0, splitSize);
+ bestSplit->quality = -1;
+ bestSplit->condensed_idx = INT_MIN;
+ split = (CvDTreeSplit*)(new char[splitSize]);
+ memset((CvDTreeSplit*)split, 0, splitSize);
+ //haveSplit = false;
+}
+
+DTreeBestSplitFinder::DTreeBestSplitFinder( const DTreeBestSplitFinder& finder, Split )
+{
+ tree = finder.tree;
+ node = finder.node;
+ splitSize = tree->get_data()->split_heap->elem_size;
+
+ bestSplit = (CvDTreeSplit*)(new char[splitSize]);
+ memcpy((CvDTreeSplit*)(bestSplit), (const CvDTreeSplit*)finder.bestSplit, splitSize);
+ split = (CvDTreeSplit*)(new char[splitSize]);
+ memset((CvDTreeSplit*)split, 0, splitSize);
+}
+
+void DTreeBestSplitFinder::operator()(const BlockedRange& range)
+{
+ int vi, vi1 = range.begin(), vi2 = range.end();
+ int n = node->sample_count;
+ CvDTreeTrainData* data = tree->get_data();
+ AutoBuffer<uchar> inn_buf(2*n*(sizeof(int) + sizeof(float)));
+
+ for( vi = vi1; vi < vi2; vi++ )
+ {
+ CvDTreeSplit *res;
+ int ci = data->get_var_type(vi);
+ if( node->get_num_valid(vi) <= 1 )
+ continue;
+
+ if( data->is_classifier )
+ {
+ if( ci >= 0 )
+ res = tree->find_split_cat_class( node, vi, bestSplit->quality, split, (uchar*)inn_buf );
+ else
+ res = tree->find_split_ord_class( node, vi, bestSplit->quality, split, (uchar*)inn_buf );
+ }
+ else
+ {
+ if( ci >= 0 )
+ res = tree->find_split_cat_reg( node, vi, bestSplit->quality, split, (uchar*)inn_buf );
+ else
+ res = tree->find_split_ord_reg( node, vi, bestSplit->quality, split, (uchar*)inn_buf );
+ }
+
+ if( res && bestSplit->quality < split->quality )
+ memcpy( (CvDTreeSplit*)bestSplit, (CvDTreeSplit*)split, splitSize );
+ }
+}
+
+void DTreeBestSplitFinder::join( DTreeBestSplitFinder& rhs )
+{
+ if( bestSplit->quality < rhs.bestSplit->quality )
+ memcpy( (CvDTreeSplit*)bestSplit, (CvDTreeSplit*)rhs.bestSplit, splitSize );
+}
+}
+
+
+CvDTreeSplit* CvDTree::find_best_split( CvDTreeNode* node )
+{
+ DTreeBestSplitFinder finder( this, node );
+
+ cv::parallel_reduce(cv::BlockedRange(0, data->var_count), finder);
+
+ CvDTreeSplit *bestSplit = data->new_split_cat( 0, -1.0f );
+ memcpy( bestSplit, finder.bestSplit, finder.splitSize );
+
+ return bestSplit;
+}
+
+CvDTreeSplit* CvDTree::find_split_ord_class( CvDTreeNode* node, int vi,
+ float init_quality, CvDTreeSplit* _split, uchar* _ext_buf )
+{
+ const float epsilon = FLT_EPSILON*2;
+ int n = node->sample_count;
+ int n1 = node->get_num_valid(vi);
+ int m = data->get_num_classes();
+
+ int base_size = 2*m*sizeof(int);
+ cv::AutoBuffer<uchar> inn_buf(base_size);
+ if( !_ext_buf )
+ inn_buf.allocate(base_size + n*(3*sizeof(int)+sizeof(float)));
+ uchar* base_buf = (uchar*)inn_buf;
+ uchar* ext_buf = _ext_buf ? _ext_buf : base_buf + base_size;
+ float* values_buf = (float*)ext_buf;
+ int* sorted_indices_buf = (int*)(values_buf + n);
+ int* sample_indices_buf = sorted_indices_buf + n;
+ const float* values = 0;
+ const int* sorted_indices = 0;
+ data->get_ord_var_data( node, vi, values_buf, sorted_indices_buf, &values,
+ &sorted_indices, sample_indices_buf );
+ int* responses_buf = sample_indices_buf + n;
+ const int* responses = data->get_class_labels( node, responses_buf );
+
+ const int* rc0 = data->counts->data.i;
+ int* lc = (int*)base_buf;
+ int* rc = lc + m;
+ int i, best_i = -1;
+ double lsum2 = 0, rsum2 = 0, best_val = init_quality;
+ const double* priors = data->have_priors ? data->priors_mult->data.db : 0;
+
+ // init arrays of class instance counters on both sides of the split
+ for( i = 0; i < m; i++ )
+ {
+ lc[i] = 0;
+ rc[i] = rc0[i];
+ }
+
+ // compensate for missing values
+ for( i = n1; i < n; i++ )
+ {
+ rc[responses[sorted_indices[i]]]--;
+ }
+
+ if( !priors )
+ {
+ int L = 0, R = n1;
+
+ for( i = 0; i < m; i++ )
+ rsum2 += (double)rc[i]*rc[i];
+
+ for( i = 0; i < n1 - 1; i++ )
+ {
+ int idx = responses[sorted_indices[i]];
+ int lv, rv;
+ L++; R--;
+ lv = lc[idx]; rv = rc[idx];
+ lsum2 += lv*2 + 1;
+ rsum2 -= rv*2 - 1;
+ lc[idx] = lv + 1; rc[idx] = rv - 1;
+
+ if( values[i] + epsilon < values[i+1] )
+ {
+ double val = (lsum2*R + rsum2*L)/((double)L*R);
+ if( best_val < val )
+ {
+ best_val = val;
+ best_i = i;
+ }
+ }
+ }
+ }
+ else
+ {
+ double L = 0, R = 0;
+ for( i = 0; i < m; i++ )
+ {
+ double wv = rc[i]*priors[i];
+ R += wv;
+ rsum2 += wv*wv;
+ }
+
+ for( i = 0; i < n1 - 1; i++ )
+ {
+ int idx = responses[sorted_indices[i]];
+ int lv, rv;
+ double p = priors[idx], p2 = p*p;
+ L += p; R -= p;
+ lv = lc[idx]; rv = rc[idx];
+ lsum2 += p2*(lv*2 + 1);
+ rsum2 -= p2*(rv*2 - 1);
+ lc[idx] = lv + 1; rc[idx] = rv - 1;
+
+ if( values[i] + epsilon < values[i+1] )
+ {
+ double val = (lsum2*R + rsum2*L)/((double)L*R);
+ if( best_val < val )
+ {
+ best_val = val;
+ best_i = i;
+ }
+ }
+ }
+ }
+
+ CvDTreeSplit* split = 0;
+ if( best_i >= 0 )
+ {
+ split = _split ? _split : data->new_split_ord( 0, 0.0f, 0, 0, 0.0f );
+ split->var_idx = vi;
+ split->ord.c = (values[best_i] + values[best_i+1])*0.5f;
+ split->ord.split_point = best_i;
+ split->inversed = 0;
+ split->quality = (float)best_val;
+ }
+ return split;
+}
+
+
+void CvDTree::cluster_categories( const int* vectors, int n, int m,
+ int* csums, int k, int* labels )
+{
+ // TODO: consider adding priors (class weights) and sample weights to the clustering algorithm
+ int iters = 0, max_iters = 100;
+ int i, j, idx;
+ double* buf = (double*)cvStackAlloc( (n + k)*sizeof(buf[0]) );
+ double *v_weights = buf, *c_weights = buf + n;
+ bool modified = true;
+ CvRNG* r = &data->rng;
+
+ // assign labels randomly
+ for( i = 0; i < n; i++ )
+ {
+ int sum = 0;
+ const int* v = vectors + i*m;
+ labels[i] = i < k ? i : (cvRandInt(r) % k);
+
+ // compute weight of each vector
+ for( j = 0; j < m; j++ )
+ sum += v[j];
+ v_weights[i] = sum ? 1./sum : 0.;
+ }
+
+ for( i = 0; i < n; i++ )
+ {
+ int i1 = cvRandInt(r) % n;
+ int i2 = cvRandInt(r) % n;
+ CV_SWAP( labels[i1], labels[i2], j );
+ }
+
+ for( iters = 0; iters <= max_iters; iters++ )
+ {
+ // calculate csums
+ for( i = 0; i < k; i++ )
+ {
+ for( j = 0; j < m; j++ )
+ csums[i*m + j] = 0;
+ }
+
+ for( i = 0; i < n; i++ )
+ {
+ const int* v = vectors + i*m;
+ int* s = csums + labels[i]*m;
+ for( j = 0; j < m; j++ )
+ s[j] += v[j];
+ }
+
+ // exit the loop here, when we have up-to-date csums
+ if( iters == max_iters || !modified )
+ break;
+
+ modified = false;
+
+ // calculate weight of each cluster
+ for( i = 0; i < k; i++ )
+ {
+ const int* s = csums + i*m;
+ int sum = 0;
+ for( j = 0; j < m; j++ )
+ sum += s[j];
+ c_weights[i] = sum ? 1./sum : 0;
+ }
+
+ // now for each vector determine the closest cluster
+ for( i = 0; i < n; i++ )
+ {
+ const int* v = vectors + i*m;
+ double alpha = v_weights[i];
+ double min_dist2 = DBL_MAX;
+ int min_idx = -1;
+
+ for( idx = 0; idx < k; idx++ )
+ {
+ const int* s = csums + idx*m;
+ double dist2 = 0., beta = c_weights[idx];
+ for( j = 0; j < m; j++ )
+ {
+ double t = v[j]*alpha - s[j]*beta;
+ dist2 += t*t;
+ }
+ if( min_dist2 > dist2 )
+ {
+ min_dist2 = dist2;
+ min_idx = idx;
+ }
+ }
+
+ if( min_idx != labels[i] )
+ modified = true;
+ labels[i] = min_idx;
+ }
+ }
+}
+
+
+CvDTreeSplit* CvDTree::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 m = data->get_num_classes();
+ int _mi = data->cat_count->data.i[ci], mi = _mi;
+
+ int base_size = m*(3 + mi)*sizeof(int) + (mi+1)*sizeof(double);
+ if( m > 2 && mi > data->params.max_categories )
+ base_size += (m*min(data->params.max_categories, n) + mi)*sizeof(int);
+ else
+ base_size += mi*sizeof(int*);
+ cv::AutoBuffer<uchar> inn_buf(base_size);
+ if( !_ext_buf )
+ inn_buf.allocate(base_size + 2*n*sizeof(int));
+ uchar* base_buf = (uchar*)inn_buf;
+ uchar* ext_buf = _ext_buf ? _ext_buf : base_buf + base_size;
+
+ int* lc = (int*)base_buf;
+ int* rc = lc + m;
+ int* _cjk = rc + m*2, *cjk = _cjk;
+ double* c_weights = (double*)alignPtr(cjk + m*mi, sizeof(double));
+
+ int* labels_buf = (int*)ext_buf;
+ const int* labels = data->get_cat_var_data(node, vi, labels_buf);
+ int* responses_buf = labels_buf + n;
+ const int* responses = data->get_class_labels(node, responses_buf);
+
+ int* cluster_labels = 0;
+ int** int_ptr = 0;
+ int i, j, k, idx;
+ double L = 0, R = 0;
+ double best_val = init_quality;
+ int prevcode = 0, best_subset = -1, subset_i, subset_n, subtract = 0;
+ const double* priors = data->priors_mult->data.db;
+
+ // init array of counters:
+ // c_{jk} - number of samples that have vi-th input variable = j and response = k.
+ for( j = -1; j < mi; j++ )
+ for( k = 0; k < m; k++ )
+ cjk[j*m + k] = 0;
+
+ for( i = 0; i < n; i++ )
+ {
+ j = ( labels[i] == 65535 && data->is_buf_16u) ? -1 : labels[i];
+ k = responses[i];
+ cjk[j*m + k]++;
+ }
+
+ if( m > 2 )
+ {
+ if( mi > data->params.max_categories )
+ {
+ mi = MIN(data->params.max_categories, n);
+ cjk = (int*)(c_weights + _mi);
+ cluster_labels = cjk + m*mi;
+ cluster_categories( _cjk, _mi, m, cjk, mi, cluster_labels );
+ }
+ subset_i = 1;
+ subset_n = 1 << mi;
+ }
+ else
+ {
+ assert( m == 2 );
+ int_ptr = (int**)(c_weights + _mi);
+ for( j = 0; j < mi; j++ )
+ int_ptr[j] = cjk + j*2 + 1;
+ icvSortIntPtr( int_ptr, mi, 0 );
+ subset_i = 0;
+ subset_n = mi;
+ }
+
+ for( k = 0; k < m; k++ )
+ {
+ int sum = 0;
+ for( j = 0; j < mi; j++ )
+ sum += cjk[j*m + k];
+ rc[k] = sum;
+ lc[k] = 0;
+ }
+
+ for( j = 0; j < mi; j++ )
+ {
+ double sum = 0;
+ for( k = 0; k < m; k++ )
+ sum += cjk[j*m + k]*priors[k];
+ c_weights[j] = sum;
+ R += c_weights[j];
+ }
+
+ for( ; subset_i < subset_n; subset_i++ )
+ {
+ double weight;
+ int* crow;
+ double lsum2 = 0, rsum2 = 0;
+
+ if( m == 2 )
+ idx = (int)(int_ptr[subset_i] - cjk)/2;
+ else
+ {
+ int graycode = (subset_i>>1)^subset_i;
+ int diff = graycode ^ prevcode;
+
+ // determine index of the changed bit.
+ Cv32suf u;
+ idx = diff >= (1 << 16) ? 16 : 0;
+ u.f = (float)(((diff >> 16) | diff) & 65535);
+ idx += (u.i >> 23) - 127;
+ subtract = graycode < prevcode;
+ prevcode = graycode;
+ }
+
+ crow = cjk + idx*m;
+ weight = c_weights[idx];
+ if( weight < FLT_EPSILON )
+ continue;
+
+ if( !subtract )
+ {
+ for( k = 0; k < m; k++ )
+ {
+ int t = crow[k];
+ int lval = lc[k] + t;
+ int rval = rc[k] - t;
+ double p = priors[k], p2 = p*p;
+ lsum2 += p2*lval*lval;
+ rsum2 += p2*rval*rval;
+ lc[k] = lval; rc[k] = rval;
+ }
+ L += weight;
+ R -= weight;
+ }
+ else
+ {
+ for( k = 0; k < m; k++ )
+ {
+ int t = crow[k];
+ int lval = lc[k] - t;
+ int rval = rc[k] + t;
+ double p = priors[k], p2 = p*p;
+ lsum2 += p2*lval*lval;
+ rsum2 += p2*rval*rval;
+ lc[k] = lval; rc[k] = rval;
+ }
+ L -= weight;
+ R += weight;
+ }
+
+ if( L > FLT_EPSILON && R > FLT_EPSILON )
+ {
+ double val = (lsum2*R + rsum2*L)/((double)L*R);
+ if( best_val < val )
+ {
+ best_val = val;
+ best_subset = subset_i;
+ }
+ }
+ }
+
+ CvDTreeSplit* split = 0;
+ if( best_subset >= 0 )
+ {
+ split = _split ? _split : data->new_split_cat( 0, -1.0f );
+ split->var_idx = vi;
+ split->quality = (float)best_val;
+ memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));
+ if( m == 2 )
+ {
+ for( i = 0; i <= best_subset; i++ )
+ {
+ idx = (int)(int_ptr[i] - cjk) >> 1;
+ split->subset[idx >> 5] |= 1 << (idx & 31);
+ }
+ }
+ else
+ {
+ for( i = 0; i < _mi; i++ )
+ {
+ idx = cluster_labels ? cluster_labels[i] : i;
+ if( best_subset & (1 << idx) )
+ split->subset[i >> 5] |= 1 << (i & 31);
+ }
+ }
+ }
+ return split;
+}
+
+
+CvDTreeSplit* CvDTree::find_split_ord_reg( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split, uchar* _ext_buf )
+{
+ const float epsilon = FLT_EPSILON*2;
+ int n = node->sample_count;
+ int n1 = node->get_num_valid(vi);
+
+ cv::AutoBuffer<uchar> inn_buf;
+ if( !_ext_buf )
+ inn_buf.allocate(2*n*(sizeof(int) + sizeof(float)));
+ uchar* ext_buf = _ext_buf ? _ext_buf : (uchar*)inn_buf;
+ float* values_buf = (float*)ext_buf;
+ int* sorted_indices_buf = (int*)(values_buf + n);
+ int* sample_indices_buf = sorted_indices_buf + n;
+ const float* values = 0;
+ const int* sorted_indices = 0;
+ data->get_ord_var_data( node, vi, values_buf, sorted_indices_buf, &values, &sorted_indices, sample_indices_buf );
+ float* responses_buf = (float*)(sample_indices_buf + n);
+ const float* responses = data->get_ord_responses( node, responses_buf, sample_indices_buf );
+
+ int i, best_i = -1;
+ double best_val = init_quality, lsum = 0, rsum = node->value*n;
+ int L = 0, R = n1;
+
+ // compensate for missing values
+ for( i = n1; i < n; i++ )
+ rsum -= responses[sorted_indices[i]];
+
+ // find the optimal split
+ for( i = 0; i < n1 - 1; i++ )
+ {
+ float t = responses[sorted_indices[i]];
+ L++; R--;
+ lsum += t;
+ rsum -= t;
+
+ if( values[i] + epsilon < values[i+1] )
+ {
+ double val = (lsum*lsum*R + rsum*rsum*L)/((double)L*R);
+ if( best_val < val )
+ {
+ best_val = val;
+ best_i = i;
+ }
+ }
+ }
+
+ CvDTreeSplit* split = 0;
+ if( best_i >= 0 )
+ {
+ split = _split ? _split : data->new_split_ord( 0, 0.0f, 0, 0, 0.0f );
+ split->var_idx = vi;
+ split->ord.c = (values[best_i] + values[best_i+1])*0.5f;
+ split->ord.split_point = best_i;
+ split->inversed = 0;
+ split->quality = (float)best_val;
+ }
+ return split;
+}
+
+CvDTreeSplit* CvDTree::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 mi = data->cat_count->data.i[ci];
+
+ int base_size = (mi+2)*sizeof(double) + (mi+1)*(sizeof(int) + sizeof(double*));
+ 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);
+
+ double* sum = (double*)cv::alignPtr(base_buf,sizeof(double)) + 1;
+ int* counts = (int*)(sum + mi) + 1;
+ double** sum_ptr = (double**)(counts + mi);
+ int i, L = 0, R = 0;
+ double best_val = init_quality, lsum = 0, rsum = 0;
+ int best_subset = -1, subset_i;
+
+ for( i = -1; i < mi; i++ )
+ sum[i] = counts[i] = 0;
+
+ // calculate sum response and weight of each category of the input var
+ for( i = 0; i < n; i++ )
+ {
+ int idx = ( (labels[i] == 65535) && data->is_buf_16u ) ? -1 : labels[i];
+ double s = sum[idx] + responses[i];
+ int nc = counts[idx] + 1;
+ sum[idx] = s;
+ counts[idx] = nc;
+ }
+
+ // calculate average response in each category
+ for( i = 0; i < mi; i++ )
+ {
+ R += counts[i];
+ rsum += sum[i];
+ sum[i] /= MAX(counts[i],1);
+ sum_ptr[i] = sum + i;
+ }
+
+ icvSortDblPtr( sum_ptr, mi, 0 );
+
+ // revert back to unnormalized sums
+ // (there should be a very little loss of accuracy)
+ for( i = 0; i < mi; i++ )
+ sum[i] *= counts[i];
+
+ for( subset_i = 0; subset_i < mi-1; subset_i++ )
+ {
+ int idx = (int)(sum_ptr[subset_i] - sum);
+ int ni = counts[idx];
+
+ if( ni )
+ {
+ double s = sum[idx];
+ lsum += s; L += ni;
+ rsum -= s; R -= ni;
+
+ if( L && R )
+ {
+ double val = (lsum*lsum*R + rsum*rsum*L)/((double)L*R);
+ if( best_val < val )
+ {
+ best_val = val;
+ best_subset = subset_i;
+ }
+ }
+ }
+ }
+
+ CvDTreeSplit* split = 0;
+ if( best_subset >= 0 )
+ {
+ split = _split ? _split : data->new_split_cat( 0, -1.0f);
+ split->var_idx = vi;
+ split->quality = (float)best_val;
+ memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));
+ for( i = 0; i <= best_subset; i++ )
+ {
+ int idx = (int)(sum_ptr[i] - sum);
+ split->subset[idx >> 5] |= 1 << (idx & 31);
+ }
+ }
+ return split;
+}
+
+CvDTreeSplit* CvDTree::find_surrogate_split_ord( CvDTreeNode* node, int vi, uchar* _ext_buf )
+{
+ const float epsilon = FLT_EPSILON*2;
+ const char* dir = (char*)data->direction->data.ptr;
+ int n = node->sample_count, n1 = node->get_num_valid(vi);
+ cv::AutoBuffer<uchar> inn_buf;
+ if( !_ext_buf )
+ inn_buf.allocate( n*(sizeof(int)*(data->have_priors ? 3 : 2) + sizeof(float)) );
+ uchar* ext_buf = _ext_buf ? _ext_buf : (uchar*)inn_buf;
+ float* values_buf = (float*)ext_buf;
+ int* sorted_indices_buf = (int*)(values_buf + n);
+ int* sample_indices_buf = sorted_indices_buf + n;
+ const float* values = 0;
+ const int* sorted_indices = 0;
+ data->get_ord_var_data( node, vi, values_buf, sorted_indices_buf, &values, &sorted_indices, sample_indices_buf );
+ // LL - number of samples that both the primary and the surrogate splits send to the left
+ // LR - ... primary split sends to the left and the surrogate split sends to the right
+ // RL - ... primary split sends to the right and the surrogate split sends to the left
+ // RR - ... both send to the right
+ int i, best_i = -1, best_inversed = 0;
+ double best_val;
+
+ if( !data->have_priors )
+ {
+ int LL = 0, RL = 0, LR, RR;
+ int worst_val = cvFloor(node->maxlr), _best_val = worst_val;
+ int sum = 0, sum_abs = 0;
+
+ for( i = 0; i < n1; i++ )
+ {
+ int d = dir[sorted_indices[i]];
+ sum += d; sum_abs += d & 1;
+ }
+
+ // sum_abs = R + L; sum = R - L
+ RR = (sum_abs + sum) >> 1;
+ LR = (sum_abs - sum) >> 1;
+
+ // initially all the samples are sent to the right by the surrogate split,
+ // LR of them are sent to the left by primary split, and RR - to the right.
+ // now iteratively compute LL, LR, RL and RR for every possible surrogate split value.
+ for( i = 0; i < n1 - 1; i++ )
+ {
+ int d = dir[sorted_indices[i]];
+
+ if( d < 0 )
+ {
+ LL++; LR--;
+ if( LL + RR > _best_val && values[i] + epsilon < values[i+1] )
+ {
+ best_val = LL + RR;
+ best_i = i; best_inversed = 0;
+ }
+ }
+ else if( d > 0 )
+ {
+ RL++; RR--;
+ if( RL + LR > _best_val && values[i] + epsilon < values[i+1] )
+ {
+ best_val = RL + LR;
+ best_i = i; best_inversed = 1;
+ }
+ }
+ }
+ best_val = _best_val;
+ }
+ else
+ {
+ double LL = 0, RL = 0, LR, RR;
+ double worst_val = node->maxlr;
+ double sum = 0, sum_abs = 0;
+ const double* priors = data->priors_mult->data.db;
+ int* responses_buf = sample_indices_buf + n;
+ const int* responses = data->get_class_labels(node, responses_buf);
+ best_val = worst_val;
+
+ for( i = 0; i < n1; i++ )
+ {
+ int idx = sorted_indices[i];
+ double w = priors[responses[idx]];
+ int d = dir[idx];
+ sum += d*w; sum_abs += (d & 1)*w;
+ }
+
+ // sum_abs = R + L; sum = R - L
+ RR = (sum_abs + sum)*0.5;
+ LR = (sum_abs - sum)*0.5;
+
+ // initially all the samples are sent to the right by the surrogate split,
+ // LR of them are sent to the left by primary split, and RR - to the right.
+ // now iteratively compute LL, LR, RL and RR for every possible surrogate split value.
+ for( i = 0; i < n1 - 1; i++ )
+ {
+ int idx = sorted_indices[i];
+ double w = priors[responses[idx]];
+ int d = dir[idx];
+
+ if( d < 0 )
+ {
+ LL += w; LR -= w;
+ if( LL + RR > best_val && values[i] + epsilon < values[i+1] )
+ {
+ best_val = LL + RR;
+ best_i = i; best_inversed = 0;
+ }
+ }
+ else if( d > 0 )
+ {
+ RL += w; RR -= w;
+ if( RL + LR > best_val && values[i] + epsilon < values[i+1] )
+ {
+ best_val = RL + LR;
+ best_i = i; best_inversed = 1;
+ }
+ }
+ }
+ }
+ return best_i >= 0 && best_val > node->maxlr ? data->new_split_ord( vi,
+ (values[best_i] + values[best_i+1])*0.5f, best_i, best_inversed, (float)best_val ) : 0;
+}
+
+
+CvDTreeSplit* CvDTree::find_surrogate_split_cat( CvDTreeNode* node, int vi, uchar* _ext_buf )
+{
+ const char* dir = (char*)data->direction->data.ptr;
+ int n = node->sample_count;
+ int i, mi = data->cat_count->data.i[data->get_var_type(vi)], l_win = 0;
+
+ int base_size = (2*(mi+1)+1)*sizeof(double) + (!data->have_priors ? 2*(mi+1)*sizeof(int) : 0);
+ cv::AutoBuffer<uchar> inn_buf(base_size);
+ if( !_ext_buf )
+ inn_buf.allocate(base_size + n*(sizeof(int) + (data->have_priors ? sizeof(int) : 0)));
+ 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);
+ // LL - number of samples that both the primary and the surrogate splits send to the left
+ // LR - ... primary split sends to the left and the surrogate split sends to the right
+ // RL - ... primary split sends to the right and the surrogate split sends to the left
+ // RR - ... both send to the right
+ CvDTreeSplit* split = data->new_split_cat( vi, 0 );
+ double best_val = 0;
+ double* lc = (double*)cv::alignPtr(base_buf,sizeof(double)) + 1;
+ double* rc = lc + mi + 1;
+
+ for( i = -1; i < mi; i++ )
+ lc[i] = rc[i] = 0;
+
+ // for each category calculate the weight of samples
+ // sent to the left (lc) and to the right (rc) by the primary split
+ if( !data->have_priors )
+ {
+ int* _lc = (int*)rc + 1;
+ int* _rc = _lc + mi + 1;
+
+ for( i = -1; i < mi; i++ )
+ _lc[i] = _rc[i] = 0;
+
+ for( i = 0; i < n; i++ )
+ {
+ int idx = ( (labels[i] == 65535) && (data->is_buf_16u) ) ? -1 : labels[i];
+ int d = dir[i];
+ int sum = _lc[idx] + d;
+ int sum_abs = _rc[idx] + (d & 1);
+ _lc[idx] = sum; _rc[idx] = sum_abs;
+ }
+
+ for( i = 0; i < mi; i++ )
+ {
+ int sum = _lc[i];
+ int sum_abs = _rc[i];
+ lc[i] = (sum_abs - sum) >> 1;
+ rc[i] = (sum_abs + sum) >> 1;
+ }
+ }
+ else
+ {
+ const double* priors = data->priors_mult->data.db;
+ 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] == 65535) && (data->is_buf_16u) ) ? -1 : labels[i];
+ double w = priors[responses[i]];
+ int d = dir[i];
+ double sum = lc[idx] + d*w;
+ double sum_abs = rc[idx] + (d & 1)*w;
+ lc[idx] = sum; rc[idx] = sum_abs;
+ }
+
+ for( i = 0; i < mi; i++ )
+ {
+ double sum = lc[i];
+ double sum_abs = rc[i];
+ lc[i] = (sum_abs - sum) * 0.5;
+ rc[i] = (sum_abs + sum) * 0.5;
+ }
+ }
+
+ // 2. now form the split.
+ // in each category send all the samples to the same direction as majority
+ for( i = 0; i < mi; i++ )
+ {
+ double lval = lc[i], rval = rc[i];
+ if( lval > rval )
+ {
+ split->subset[i >> 5] |= 1 << (i & 31);
+ best_val += lval;
+ l_win++;
+ }
+ else
+ best_val += rval;
+ }
+
+ split->quality = (float)best_val;
+ if( split->quality <= node->maxlr || l_win == 0 || l_win == mi )
+ cvSetRemoveByPtr( data->split_heap, split ), split = 0;
+
+ return split;
+}
+
+
+void CvDTree::calc_node_value( CvDTreeNode* node )
+{
+ int i, j, k, n = node->sample_count, cv_n = data->params.cv_folds;
+ int m = data->get_num_classes();
+
+ int base_size = data->is_classifier ? m*cv_n*sizeof(int) : 2*cv_n*sizeof(double)+cv_n*sizeof(int);
+ int ext_size = n*(sizeof(int) + (data->is_classifier ? sizeof(int) : sizeof(int)+sizeof(float)));
+ cv::AutoBuffer<uchar> inn_buf(base_size + ext_size);
+ uchar* base_buf = (uchar*)inn_buf;
+ uchar* ext_buf = base_buf + base_size;
+
+ int* cv_labels_buf = (int*)ext_buf;
+ const int* cv_labels = data->get_cv_labels(node, cv_labels_buf);
+
+ if( data->is_classifier )
+ {
+ // in case of classification tree:
+ // * node value is the label of the class that has the largest weight in the node.
+ // * node risk is the weighted number of misclassified samples,
+ // * j-th cross-validation fold value and risk are calculated as above,
+ // but using the samples with cv_labels(*)!=j.
+ // * j-th cross-validation fold error is calculated as the weighted number of
+ // misclassified samples with cv_labels(*)==j.
+
+ // compute the number of instances of each class
+ int* cls_count = data->counts->data.i;
+ int* responses_buf = cv_labels_buf + n;
+ const int* responses = data->get_class_labels(node, responses_buf);
+ int* cv_cls_count = (int*)base_buf;
+ double max_val = -1, total_weight = 0;
+ int max_k = -1;
+ double* priors = data->priors_mult->data.db;
+
+ for( k = 0; k < m; k++ )
+ cls_count[k] = 0;
+
+ if( cv_n == 0 )
+ {
+ for( i = 0; i < n; i++ )
+ cls_count[responses[i]]++;
+ }
+ else
+ {
+ for( j = 0; j < cv_n; j++ )
+ for( k = 0; k < m; k++ )
+ cv_cls_count[j*m + k] = 0;
+
+ for( i = 0; i < n; i++ )
+ {
+ j = cv_labels[i]; k = responses[i];
+ cv_cls_count[j*m + k]++;
+ }
+
+ for( j = 0; j < cv_n; j++ )
+ for( k = 0; k < m; k++ )
+ cls_count[k] += cv_cls_count[j*m + k];
+ }
+
+ if( data->have_priors && node->parent == 0 )
+ {
+ // compute priors_mult from priors, take the sample ratio into account.
+ double sum = 0;
+ for( k = 0; k < m; k++ )
+ {
+ int n_k = cls_count[k];
+ priors[k] = data->priors->data.db[k]*(n_k ? 1./n_k : 0.);
+ sum += priors[k];
+ }
+ sum = 1./sum;
+ for( k = 0; k < m; k++ )
+ priors[k] *= sum;
+ }
+
+ for( k = 0; k < m; k++ )
+ {
+ double val = cls_count[k]*priors[k];
+ total_weight += val;
+ if( max_val < val )
+ {
+ max_val = val;
+ max_k = k;
+ }
+ }
+
+ node->class_idx = max_k;
+ node->value = data->cat_map->data.i[
+ data->cat_ofs->data.i[data->cat_var_count] + max_k];
+ node->node_risk = total_weight - max_val;
+
+ for( j = 0; j < cv_n; j++ )
+ {
+ double sum_k = 0, sum = 0, max_val_k = 0;
+ max_val = -1; max_k = -1;
+
+ for( k = 0; k < m; k++ )
+ {
+ double w = priors[k];
+ double val_k = cv_cls_count[j*m + k]*w;
+ double val = cls_count[k]*w - val_k;
+ sum_k += val_k;
+ sum += val;
+ if( max_val < val )
+ {
+ max_val = val;
+ max_val_k = val_k;
+ max_k = k;
+ }
+ }
+
+ node->cv_Tn[j] = INT_MAX;
+ node->cv_node_risk[j] = sum - max_val;
+ node->cv_node_error[j] = sum_k - max_val_k;
+ }
+ }
+ else
+ {
+ // in case of regression tree:
+ // * node value is 1/n*sum_i(Y_i), where Y_i is i-th response,
+ // n is the number of samples in the node.
+ // * node risk is the sum of squared errors: sum_i((Y_i - <node_value>)^2)
+ // * j-th cross-validation fold value and risk are calculated as above,
+ // but using the samples with cv_labels(*)!=j.
+ // * j-th cross-validation fold error is calculated
+ // using samples with cv_labels(*)==j as the test subset:
+ // error_j = sum_(i,cv_labels(i)==j)((Y_i - <node_value_j>)^2),
+ // where node_value_j is the node value calculated
+ // as described in the previous bullet, and summation is done
+ // over the samples with cv_labels(*)==j.
+
+ double sum = 0, sum2 = 0;
+ float* values_buf = (float*)(cv_labels_buf + n);
+ int* sample_indices_buf = (int*)(values_buf + n);
+ const float* values = data->get_ord_responses(node, values_buf, sample_indices_buf);
+ double *cv_sum = 0, *cv_sum2 = 0;
+ int* cv_count = 0;
+
+ if( cv_n == 0 )
+ {
+ for( i = 0; i < n; i++ )
+ {
+ double t = values[i];
+ sum += t;
+ sum2 += t*t;
+ }
+ }
+ else
+ {
+ cv_sum = (double*)base_buf;
+ cv_sum2 = cv_sum + cv_n;
+ cv_count = (int*)(cv_sum2 + cv_n);
+
+ for( j = 0; j < cv_n; j++ )
+ {
+ cv_sum[j] = cv_sum2[j] = 0.;
+ cv_count[j] = 0;
+ }
+
+ for( i = 0; i < n; i++ )
+ {
+ j = cv_labels[i];
+ double t = values[i];
+ double s = cv_sum[j] + t;
+ double s2 = cv_sum2[j] + t*t;
+ int nc = cv_count[j] + 1;
+ cv_sum[j] = s;
+ cv_sum2[j] = s2;
+ cv_count[j] = nc;
+ }
+
+ for( j = 0; j < cv_n; j++ )
+ {
+ sum += cv_sum[j];
+ sum2 += cv_sum2[j];
+ }
+ }
+
+ node->node_risk = sum2 - (sum/n)*sum;
+ node->value = sum/n;
+
+ for( j = 0; j < cv_n; j++ )
+ {
+ double s = cv_sum[j], si = sum - s;
+ double s2 = cv_sum2[j], s2i = sum2 - s2;
+ int c = cv_count[j], ci = n - c;
+ double r = si/MAX(ci,1);
+ node->cv_node_risk[j] = s2i - r*r*ci;
+ node->cv_node_error[j] = s2 - 2*r*s + c*r*r;
+ node->cv_Tn[j] = INT_MAX;
+ }
+ }
+}
+
+
+void CvDTree::complete_node_dir( CvDTreeNode* node )
+{
+ int vi, i, n = node->sample_count, nl, nr, d0 = 0, d1 = -1;
+ int nz = n - node->get_num_valid(node->split->var_idx);
+ char* dir = (char*)data->direction->data.ptr;
+
+ // try to complete direction using surrogate splits
+ if( nz && data->params.use_surrogates )
+ {
+ cv::AutoBuffer<uchar> inn_buf(n*(2*sizeof(int)+sizeof(float)));
+ CvDTreeSplit* split = node->split->next;
+ for( ; split != 0 && nz; split = split->next )
+ {
+ int inversed_mask = split->inversed ? -1 : 0;
+ vi = split->var_idx;
+
+ if( data->get_var_type(vi) >= 0 ) // split on categorical var
+ {
+ int* labels_buf = (int*)(uchar*)inn_buf;
+ const int* labels = data->get_cat_var_data(node, vi, labels_buf);
+ const int* subset = split->subset;
+
+ for( i = 0; i < n; i++ )
+ {
+ int idx = labels[i];
+ if( !dir[i] && ( ((idx >= 0)&&(!data->is_buf_16u)) || ((idx != 65535)&&(data->is_buf_16u)) ))
+
+ {
+ int d = CV_DTREE_CAT_DIR(idx,subset);
+ dir[i] = (char)((d ^ inversed_mask) - inversed_mask);
+ if( --nz )
+ break;
+ }
+ }
+ }
+ else // split on ordered var
+ {
+ float* values_buf = (float*)(uchar*)inn_buf;
+ int* sorted_indices_buf = (int*)(values_buf + n);
+ int* sample_indices_buf = sorted_indices_buf + n;
+ const float* values = 0;
+ const int* sorted_indices = 0;
+ data->get_ord_var_data( node, vi, values_buf, sorted_indices_buf, &values, &sorted_indices, sample_indices_buf );
+ int split_point = split->ord.split_point;
+ int n1 = node->get_num_valid(vi);
+
+ assert( 0 <= split_point && split_point < n-1 );
+
+ for( i = 0; i < n1; i++ )
+ {
+ int idx = sorted_indices[i];
+ if( !dir[idx] )
+ {
+ int d = i <= split_point ? -1 : 1;
+ dir[idx] = (char)((d ^ inversed_mask) - inversed_mask);
+ if( --nz )
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ // find the default direction for the rest
+ if( nz )
+ {
+ for( i = nr = 0; i < n; i++ )
+ nr += dir[i] > 0;
+ nl = n - nr - nz;
+ d0 = nl > nr ? -1 : nr > nl;
+ }
+
+ // make sure that every sample is directed either to the left or to the right
+ for( i = 0; i < n; i++ )
+ {
+ int d = dir[i];
+ if( !d )
+ {
+ d = d0;
+ if( !d )
+ d = d1, d1 = -d1;
+ }
+ d = d > 0;
+ dir[i] = (char)d; // remap (-1,1) to (0,1)
+ }
+}
+
+
+void CvDTree::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_idx = data->split_buf->data.i;
+ int new_buf_idx = data->get_child_buf_idx( node );
+ int work_var_count = data->get_work_var_count();
+ CvMat* buf = data->buf;
+ cv::AutoBuffer<uchar> inn_buf(n*(3*sizeof(int) + sizeof(float)));
+ int* temp_buf = (int*)(uchar*)inn_buf;
+
+ complete_node_dir(node);
+
+ for( i = nl = nr = 0; i < n; i++ )
+ {
+ int d = dir[i];
+ // initialize new indices for splitting ordered variables
+ new_idx[i] = (nl & (d-1)) | (nr & -d); // d ? ri : li
+ 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);
+
+ // split ordered variables, keep both halves sorted.
+ for( vi = 0; vi < data->var_count; vi++ )
+ {
+ int ci = data->get_var_type(vi);
+
+ if( ci >= 0 || !split_input_data )
+ continue;
+
+ int n1 = node->get_num_valid(vi);
+ float* src_val_buf = (float*)(uchar*)(temp_buf + n);
+ int* src_sorted_idx_buf = (int*)(src_val_buf + n);
+ int* src_sample_idx_buf = src_sorted_idx_buf + n;
+ const float* src_val = 0;
+ const int* src_sorted_idx = 0;
+ data->get_ord_var_data(node, vi, src_val_buf, src_sorted_idx_buf, &src_val, &src_sorted_idx, src_sample_idx_buf);
+
+ for(i = 0; i < n; i++)
+ temp_buf[i] = src_sorted_idx[i];
+
+ if (data->is_buf_16u)
+ {
+ unsigned short *ldst, *rdst, *ldst0, *rdst0;
+ //unsigned short tl, tr;
+ ldst0 = ldst = (unsigned short*)(buf->data.s + left->buf_idx*buf->cols +
+ vi*scount + left->offset);
+ rdst0 = rdst = (unsigned short*)(ldst + nl);
+
+ // split sorted
+ for( i = 0; i < n1; i++ )
+ {
+ int idx = temp_buf[i];
+ int d = dir[idx];
+ idx = new_idx[idx];
+ if (d)
+ {
+ *rdst = (unsigned short)idx;
+ rdst++;
+ }
+ else
+ {
+ *ldst = (unsigned short)idx;
+ ldst++;
+ }
+ }
+
+ left->set_num_valid(vi, (int)(ldst - ldst0));
+ right->set_num_valid(vi, (int)(rdst - rdst0));
+
+ // split missing
+ for( ; i < n; i++ )
+ {
+ int idx = temp_buf[i];
+ int d = dir[idx];
+ idx = new_idx[idx];
+ if (d)
+ {
+ *rdst = (unsigned short)idx;
+ rdst++;
+ }
+ else
+ {
+ *ldst = (unsigned short)idx;
+ ldst++;
+ }
+ }
+ }
+ else
+ {
+ int *ldst0, *ldst, *rdst0, *rdst;
+ ldst0 = ldst = buf->data.i + left->buf_idx*buf->cols +
+ vi*scount + left->offset;
+ rdst0 = rdst = buf->data.i + right->buf_idx*buf->cols +
+ vi*scount + right->offset;
+
+ // split sorted
+ for( i = 0; i < n1; i++ )
+ {
+ int idx = temp_buf[i];
+ int d = dir[idx];
+ idx = new_idx[idx];
+ if (d)
+ {
+ *rdst = idx;
+ rdst++;
+ }
+ else
+ {
+ *ldst = idx;
+ ldst++;
+ }
+ }
+
+ left->set_num_valid(vi, (int)(ldst - ldst0));
+ right->set_num_valid(vi, (int)(rdst - rdst0));
+
+ // split missing
+ for( ; i < n; i++ )
+ {
+ int idx = temp_buf[i];
+ int d = dir[idx];
+ idx = new_idx[idx];
+ if (d)
+ {
+ *rdst = idx;
+ rdst++;
+ }
+ else
+ {
+ *ldst = idx;
+ ldst++;
+ }
+ }
+ }
+ }
+
+ // split categorical vars, responses and cv_labels using new_idx relocation table
+ for( vi = 0; vi < work_var_count; vi++ )
+ {
+ int ci = data->get_var_type(vi);
+ int n1 = node->get_num_valid(vi), nr1 = 0;
+
+ if( ci < 0 || (vi < data->var_count && !split_input_data) )
+ continue;
+
+ int *src_lbls_buf = temp_buf + n;
+ const int* src_lbls = data->get_cat_var_data(node, vi, src_lbls_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 +
+ vi*scount + left->offset);
+ unsigned short *rdst = (unsigned short *)(buf->data.s + right->buf_idx*buf->cols +
+ vi*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 )&d;
+ }
+ 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 +
+ vi*scount + left->offset;
+ int *rdst = buf->data.i + right->buf_idx*buf->cols +
+ vi*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)&d;
+ }
+ 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 = temp_buf + n;
+ const int* 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);
+}
+
+float CvDTree::calc_error( CvMLData* _data, int type, vector<float> *resp )
+{
+ float err = 0;
+ const CvMat* values = _data->get_values();
+ const CvMat* response = _data->get_responses();
+ const CvMat* missing = _data->get_missing();
+ const CvMat* sample_idx = (type == CV_TEST_ERROR) ? _data->get_test_sample_idx() : _data->get_train_sample_idx();
+ const CvMat* var_types = _data->get_var_types();
+ int* sidx = sample_idx ? sample_idx->data.i : 0;
+ int r_step = CV_IS_MAT_CONT(response->type) ?
+ 1 : response->step / CV_ELEM_SIZE(response->type);
+ bool is_classifier = var_types->data.ptr[var_types->cols-1] == CV_VAR_CATEGORICAL;
+ int sample_count = sample_idx ? sample_idx->cols : 0;
+ sample_count = (type == CV_TRAIN_ERROR && sample_count == 0) ? values->rows : sample_count;
+ float* pred_resp = 0;
+ if( resp && (sample_count > 0) )
+ {
+ resp->resize( sample_count );
+ pred_resp = &((*resp)[0]);
+ }
+
+ if ( is_classifier )
+ {
+ for( int i = 0; i < sample_count; i++ )
+ {
+ CvMat sample, miss;
+ int si = sidx ? sidx[i] : i;
+ cvGetRow( values, &sample, si );
+ if( missing )
+ cvGetRow( missing, &miss, si );
+ float r = (float)predict( &sample, missing ? &miss : 0 )->value;
+ if( pred_resp )
+ pred_resp[i] = r;
+ int d = fabs((double)r - response->data.fl[si*r_step]) <= FLT_EPSILON ? 0 : 1;
+ err += d;
+ }
+ err = sample_count ? err / (float)sample_count * 100 : -FLT_MAX;
+ }
+ else
+ {
+ for( int i = 0; i < sample_count; i++ )
+ {
+ CvMat sample, miss;
+ int si = sidx ? sidx[i] : i;
+ cvGetRow( values, &sample, si );
+ if( missing )
+ cvGetRow( missing, &miss, si );
+ float r = (float)predict( &sample, missing ? &miss : 0 )->value;
+ if( pred_resp )
+ pred_resp[i] = r;
+ float d = r - response->data.fl[si*r_step];
+ err += d*d;
+ }
+ err = sample_count ? err / (float)sample_count : -FLT_MAX;
+ }
+ return err;
+}
+
+void CvDTree::prune_cv()
+{
+ CvMat* ab = 0;
+ CvMat* temp = 0;
+ CvMat* err_jk = 0;
+
+ // 1. build tree sequence for each cv fold, calculate error_{Tj,beta_k}.
+ // 2. choose the best tree index (if need, apply 1SE rule).
+ // 3. store the best index and cut the branches.
+
+ CV_FUNCNAME( "CvDTree::prune_cv" );
+
+ __BEGIN__;
+
+ int ti, j, tree_count = 0, cv_n = data->params.cv_folds, n = root->sample_count;
+ // currently, 1SE for regression is not implemented
+ bool use_1se = data->params.use_1se_rule != 0 && data->is_classifier;
+ double* err;
+ double min_err = 0, min_err_se = 0;
+ int min_idx = -1;
+
+ CV_CALL( ab = cvCreateMat( 1, 256, CV_64F ));
+
+ // build the main tree sequence, calculate alpha's
+ for(;;tree_count++)
+ {
+ double min_alpha = update_tree_rnc(tree_count, -1);
+ if( cut_tree(tree_count, -1, min_alpha) )
+ break;
+
+ if( ab->cols <= tree_count )
+ {
+ CV_CALL( temp = cvCreateMat( 1, ab->cols*3/2, CV_64F ));
+ for( ti = 0; ti < ab->cols; ti++ )
+ temp->data.db[ti] = ab->data.db[ti];
+ cvReleaseMat( &ab );
+ ab = temp;
+ temp = 0;
+ }
+
+ ab->data.db[tree_count] = min_alpha;
+ }
+
+ ab->data.db[0] = 0.;
+
+ if( tree_count > 0 )
+ {
+ for( ti = 1; ti < tree_count-1; ti++ )
+ ab->data.db[ti] = sqrt(ab->data.db[ti]*ab->data.db[ti+1]);
+ ab->data.db[tree_count-1] = DBL_MAX*0.5;
+
+ CV_CALL( err_jk = cvCreateMat( cv_n, tree_count, CV_64F ));
+ err = err_jk->data.db;
+
+ for( j = 0; j < cv_n; j++ )
+ {
+ int tj = 0, tk = 0;
+ for( ; tk < tree_count; tj++ )
+ {
+ double min_alpha = update_tree_rnc(tj, j);
+ if( cut_tree(tj, j, min_alpha) )
+ min_alpha = DBL_MAX;
+
+ for( ; tk < tree_count; tk++ )
+ {
+ if( ab->data.db[tk] > min_alpha )
+ break;
+ err[j*tree_count + tk] = root->tree_error;
+ }
+ }
+ }
+
+ for( ti = 0; ti < tree_count; ti++ )
+ {
+ double sum_err = 0;
+ for( j = 0; j < cv_n; j++ )
+ sum_err += err[j*tree_count + ti];
+ if( ti == 0 || sum_err < min_err )
+ {
+ min_err = sum_err;
+ min_idx = ti;
+ if( use_1se )
+ min_err_se = sqrt( sum_err*(n - sum_err) );
+ }
+ else if( sum_err < min_err + min_err_se )
+ min_idx = ti;
+ }
+ }
+
+ pruned_tree_idx = min_idx;
+ free_prune_data(data->params.truncate_pruned_tree != 0);
+
+ __END__;
+
+ cvReleaseMat( &err_jk );
+ cvReleaseMat( &ab );
+ cvReleaseMat( &temp );
+}
+
+
+double CvDTree::update_tree_rnc( int T, int fold )
+{
+ CvDTreeNode* node = root;
+ double min_alpha = DBL_MAX;
+
+ for(;;)
+ {
+ CvDTreeNode* parent;
+ for(;;)
+ {
+ int t = fold >= 0 ? node->cv_Tn[fold] : node->Tn;
+ if( t <= T || !node->left )
+ {
+ node->complexity = 1;
+ node->tree_risk = node->node_risk;
+ node->tree_error = 0.;
+ if( fold >= 0 )
+ {
+ node->tree_risk = node->cv_node_risk[fold];
+ node->tree_error = node->cv_node_error[fold];
+ }
+ break;
+ }
+ node = node->left;
+ }
+
+ for( parent = node->parent; parent && parent->right == node;
+ node = parent, parent = parent->parent )
+ {
+ parent->complexity += node->complexity;
+ parent->tree_risk += node->tree_risk;
+ parent->tree_error += node->tree_error;
+
+ parent->alpha = ((fold >= 0 ? parent->cv_node_risk[fold] : parent->node_risk)
+ - parent->tree_risk)/(parent->complexity - 1);
+ min_alpha = MIN( min_alpha, parent->alpha );
+ }
+
+ if( !parent )
+ break;
+
+ parent->complexity = node->complexity;
+ parent->tree_risk = node->tree_risk;
+ parent->tree_error = node->tree_error;
+ node = parent->right;
+ }
+
+ return min_alpha;
+}
+
+
+int CvDTree::cut_tree( int T, int fold, double min_alpha )
+{
+ CvDTreeNode* node = root;
+ if( !node->left )
+ return 1;
+
+ for(;;)
+ {
+ CvDTreeNode* parent;
+ for(;;)
+ {
+ int t = fold >= 0 ? node->cv_Tn[fold] : node->Tn;
+ if( t <= T || !node->left )
+ break;
+ if( node->alpha <= min_alpha + FLT_EPSILON )
+ {
+ if( fold >= 0 )
+ node->cv_Tn[fold] = T;
+ else
+ node->Tn = T;
+ if( node == root )
+ return 1;
+ break;
+ }
+ node = node->left;
+ }
+
+ for( parent = node->parent; parent && parent->right == node;
+ node = parent, parent = parent->parent )
+ ;
+
+ if( !parent )
+ break;
+
+ node = parent->right;
+ }
+
+ return 0;
+}
+
+
+void CvDTree::free_prune_data(bool cut_tree)
+{
+ CvDTreeNode* node = root;
+
+ for(;;)
+ {
+ CvDTreeNode* parent;
+ for(;;)
+ {
+ // do not call cvSetRemoveByPtr( cv_heap, node->cv_Tn )
+ // as we will clear the whole cross-validation heap at the end
+ node->cv_Tn = 0;
+ node->cv_node_error = node->cv_node_risk = 0;
+ if( !node->left )
+ break;
+ node = node->left;
+ }
+
+ for( parent = node->parent; parent && parent->right == node;
+ node = parent, parent = parent->parent )
+ {
+ if( cut_tree && parent->Tn <= pruned_tree_idx )
+ {
+ data->free_node( parent->left );
+ data->free_node( parent->right );
+ parent->left = parent->right = 0;
+ }
+ }
+
+ if( !parent )
+ break;
+
+ node = parent->right;
+ }
+
+ if( data->cv_heap )
+ cvClearSet( data->cv_heap );
+}
+
+
+void CvDTree::free_tree()
+{
+ if( root && data && data->shared )
+ {
+ pruned_tree_idx = INT_MIN;
+ free_prune_data(true);
+ data->free_node(root);
+ root = 0;
+ }
+}
+
+CvDTreeNode* CvDTree::predict( const CvMat* _sample,
+ const CvMat* _missing, bool preprocessed_input ) const
+{
+ CvDTreeNode* result = 0;
+ int* catbuf = 0;
+
+ CV_FUNCNAME( "CvDTree::predict" );
+
+ __BEGIN__;
+
+ int i, step, mstep = 0;
+ const float* sample;
+ const uchar* m = 0;
+ CvDTreeNode* node = root;
+ const int* vtype;
+ const int* vidx;
+ const int* cmap;
+ const int* cofs;
+
+ if( !node )
+ CV_ERROR( CV_StsError, "The tree has not been trained yet" );
+
+ if( !CV_IS_MAT(_sample) || CV_MAT_TYPE(_sample->type) != CV_32FC1 ||
+ (_sample->cols != 1 && _sample->rows != 1) ||
+ (_sample->cols + _sample->rows - 1 != data->var_all && !preprocessed_input) ||
+ (_sample->cols + _sample->rows - 1 != data->var_count && preprocessed_input) )
+ CV_ERROR( CV_StsBadArg,
+ "the input sample must be 1d floating-point vector with the same "
+ "number of elements as the total number of variables used for training" );
+
+ sample = _sample->data.fl;
+ step = CV_IS_MAT_CONT(_sample->type) ? 1 : _sample->step/sizeof(sample[0]);
+
+ if( data->cat_count && !preprocessed_input ) // cache for categorical variables
+ {
+ int n = data->cat_count->cols;
+ catbuf = (int*)cvStackAlloc(n*sizeof(catbuf[0]));
+ for( i = 0; i < n; i++ )
+ catbuf[i] = -1;
+ }
+
+ if( _missing )
+ {
+ if( !CV_IS_MAT(_missing) || !CV_IS_MASK_ARR(_missing) ||
+ !CV_ARE_SIZES_EQ(_missing, _sample) )
+ CV_ERROR( CV_StsBadArg,
+ "the missing data mask must be 8-bit vector of the same size as input sample" );
+ m = _missing->data.ptr;
+ mstep = CV_IS_MAT_CONT(_missing->type) ? 1 : _missing->step/sizeof(m[0]);
+ }
+
+ vtype = data->var_type->data.i;
+ vidx = data->var_idx && !preprocessed_input ? data->var_idx->data.i : 0;
+ cmap = data->cat_map ? data->cat_map->data.i : 0;
+ cofs = data->cat_ofs ? data->cat_ofs->data.i : 0;
+
+ while( node->Tn > pruned_tree_idx && node->left )
+ {
+ CvDTreeSplit* split = node->split;
+ int dir = 0;
+ for( ; !dir && split != 0; split = split->next )
+ {
+ int vi = split->var_idx;
+ int ci = vtype[vi];
+ i = vidx ? vidx[vi] : vi;
+ float val = sample[i*step];
+ if( m && m[i*mstep] )
+ continue;
+ if( ci < 0 ) // ordered
+ dir = val <= split->ord.c ? -1 : 1;
+ else // categorical
+ {
+ int c;
+ if( preprocessed_input )
+ c = cvRound(val);
+ else
+ {
+ c = catbuf[ci];
+ if( c < 0 )
+ {
+ int a = c = cofs[ci];
+ int b = (ci+1 >= data->cat_ofs->cols) ? data->cat_map->cols : cofs[ci+1];
+
+ int ival = cvRound(val);
+ if( ival != val )
+ CV_ERROR( CV_StsBadArg,
+ "one of input categorical variable is not an integer" );
+
+ int sh = 0;
+ while( a < b )
+ {
+ sh++;
+ c = (a + b) >> 1;
+ if( ival < cmap[c] )
+ b = c;
+ else if( ival > cmap[c] )
+ a = c+1;
+ else
+ break;
+ }
+
+ if( c < 0 || ival != cmap[c] )
+ continue;
+
+ catbuf[ci] = c -= cofs[ci];
+ }
+ }
+ c = ( (c == 65535) && data->is_buf_16u ) ? -1 : c;
+ dir = CV_DTREE_CAT_DIR(c, split->subset);
+ }
+
+ if( split->inversed )
+ dir = -dir;
+ }
+
+ if( !dir )
+ {
+ double diff = node->right->sample_count - node->left->sample_count;
+ dir = diff < 0 ? -1 : 1;
+ }
+ node = dir < 0 ? node->left : node->right;
+ }
+
+ result = node;
+
+ __END__;
+
+ return result;
+}
+
+
+CvDTreeNode* CvDTree::predict( const Mat& _sample, const Mat& _missing, bool preprocessed_input ) const
+{
+ CvMat sample = _sample, mmask = _missing;
+ return predict(&sample, mmask.data.ptr ? &mmask : 0, preprocessed_input);
+}
+
+
+const CvMat* CvDTree::get_var_importance()
+{
+ if( !var_importance )
+ {
+ CvDTreeNode* node = root;
+ double* importance;
+ if( !node )
+ return 0;
+ var_importance = cvCreateMat( 1, data->var_count, CV_64F );
+ cvZero( var_importance );
+ importance = var_importance->data.db;
+
+ for(;;)
+ {
+ CvDTreeNode* parent;
+ for( ;; node = node->left )
+ {
+ CvDTreeSplit* split = node->split;
+
+ if( !node->left || node->Tn <= pruned_tree_idx )
+ break;
+
+ for( ; split != 0; split = split->next )
+ importance[split->var_idx] += split->quality;
+ }
+
+ for( parent = node->parent; parent && parent->right == node;
+ node = parent, parent = parent->parent )
+ ;
+
+ if( !parent )
+ break;
+
+ node = parent->right;
+ }
+
+ cvNormalize( var_importance, var_importance, 1., 0, CV_L1 );
+ }
+
+ return var_importance;
+}
+
+
+void CvDTree::write_split( CvFileStorage* fs, CvDTreeSplit* split ) const
+{
+ int ci;
+
+ cvStartWriteStruct( fs, 0, CV_NODE_MAP + CV_NODE_FLOW );
+ cvWriteInt( fs, "var", split->var_idx );
+ cvWriteReal( fs, "quality", split->quality );
+
+ ci = data->get_var_type(split->var_idx);
+ if( ci >= 0 ) // split on a categorical var
+ {
+ int i, n = data->cat_count->data.i[ci], to_right = 0, default_dir;
+ for( i = 0; i < n; i++ )
+ to_right += CV_DTREE_CAT_DIR(i,split->subset) > 0;
+
+ // ad-hoc rule when to use inverse categorical split notation
+ // to achieve more compact and clear representation
+ default_dir = to_right <= 1 || to_right <= MIN(3, n/2) || to_right <= n/3 ? -1 : 1;
+
+ cvStartWriteStruct( fs, default_dir*(split->inversed ? -1 : 1) > 0 ?
+ "in" : "not_in", CV_NODE_SEQ+CV_NODE_FLOW );
+
+ for( i = 0; i < n; i++ )
+ {
+ int dir = CV_DTREE_CAT_DIR(i,split->subset);
+ if( dir*default_dir < 0 )
+ cvWriteInt( fs, 0, i );
+ }
+ cvEndWriteStruct( fs );
+ }
+ else
+ cvWriteReal( fs, !split->inversed ? "le" : "gt", split->ord.c );
+
+ cvEndWriteStruct( fs );
+}
+
+
+void CvDTree::write_node( CvFileStorage* fs, CvDTreeNode* node ) const
+{
+ CvDTreeSplit* split;
+
+ cvStartWriteStruct( fs, 0, CV_NODE_MAP );
+
+ cvWriteInt( fs, "depth", node->depth );
+ cvWriteInt( fs, "sample_count", node->sample_count );
+ cvWriteReal( fs, "value", node->value );
+
+ if( data->is_classifier )
+ cvWriteInt( fs, "norm_class_idx", node->class_idx );
+
+ cvWriteInt( fs, "Tn", node->Tn );
+ cvWriteInt( fs, "complexity", node->complexity );
+ cvWriteReal( fs, "alpha", node->alpha );
+ cvWriteReal( fs, "node_risk", node->node_risk );
+ cvWriteReal( fs, "tree_risk", node->tree_risk );
+ cvWriteReal( fs, "tree_error", node->tree_error );
+
+ if( node->left )
+ {
+ cvStartWriteStruct( fs, "splits", CV_NODE_SEQ );
+
+ for( split = node->split; split != 0; split = split->next )
+ write_split( fs, split );
+
+ cvEndWriteStruct( fs );
+ }
+
+ cvEndWriteStruct( fs );
+}
+
+
+void CvDTree::write_tree_nodes( CvFileStorage* fs ) const
+{
+ //CV_FUNCNAME( "CvDTree::write_tree_nodes" );
+
+ __BEGIN__;
+
+ CvDTreeNode* node = root;
+
+ // traverse the tree and save all the nodes in depth-first order
+ for(;;)
+ {
+ CvDTreeNode* parent;
+ for(;;)
+ {
+ write_node( fs, node );
+ if( !node->left )
+ break;
+ node = node->left;
+ }
+
+ for( parent = node->parent; parent && parent->right == node;
+ node = parent, parent = parent->parent )
+ ;
+
+ if( !parent )
+ break;
+
+ node = parent->right;
+ }
+
+ __END__;
+}
+
+
+void CvDTree::write( CvFileStorage* fs, const char* name ) const
+{
+ //CV_FUNCNAME( "CvDTree::write" );
+
+ __BEGIN__;
+
+ cvStartWriteStruct( fs, name, CV_NODE_MAP, CV_TYPE_NAME_ML_TREE );
+
+ //get_var_importance();
+ data->write_params( fs );
+ //if( var_importance )
+ //cvWrite( fs, "var_importance", var_importance );
+ write( fs );
+
+ cvEndWriteStruct( fs );
+
+ __END__;
+}
+
+
+void CvDTree::write( CvFileStorage* fs ) const
+{
+ //CV_FUNCNAME( "CvDTree::write" );
+
+ __BEGIN__;
+
+ cvWriteInt( fs, "best_tree_idx", pruned_tree_idx );
+
+ cvStartWriteStruct( fs, "nodes", CV_NODE_SEQ );
+ write_tree_nodes( fs );
+ cvEndWriteStruct( fs );
+
+ __END__;
+}
+
+
+CvDTreeSplit* CvDTree::read_split( CvFileStorage* fs, CvFileNode* fnode )
+{
+ CvDTreeSplit* split = 0;
+
+ CV_FUNCNAME( "CvDTree::read_split" );
+
+ __BEGIN__;
+
+ int vi, ci;
+
+ if( !fnode || CV_NODE_TYPE(fnode->tag) != CV_NODE_MAP )
+ CV_ERROR( CV_StsParseError, "some of the splits are not stored properly" );
+
+ vi = cvReadIntByName( fs, fnode, "var", -1 );
+ if( (unsigned)vi >= (unsigned)data->var_count )
+ CV_ERROR( CV_StsOutOfRange, "Split variable index is out of range" );
+
+ ci = data->get_var_type(vi);
+ if( ci >= 0 ) // split on categorical var
+ {
+ int i, n = data->cat_count->data.i[ci], inversed = 0, val;
+ CvSeqReader reader;
+ CvFileNode* inseq;
+ split = data->new_split_cat( vi, 0 );
+ inseq = cvGetFileNodeByName( fs, fnode, "in" );
+ if( !inseq )
+ {
+ inseq = cvGetFileNodeByName( fs, fnode, "not_in" );
+ inversed = 1;
+ }
+ if( !inseq ||
+ (CV_NODE_TYPE(inseq->tag) != CV_NODE_SEQ && CV_NODE_TYPE(inseq->tag) != CV_NODE_INT))
+ CV_ERROR( CV_StsParseError,
+ "Either 'in' or 'not_in' tags should be inside a categorical split data" );
+
+ if( CV_NODE_TYPE(inseq->tag) == CV_NODE_INT )
+ {
+ val = inseq->data.i;
+ if( (unsigned)val >= (unsigned)n )
+ CV_ERROR( CV_StsOutOfRange, "some of in/not_in elements are out of range" );
+
+ split->subset[val >> 5] |= 1 << (val & 31);
+ }
+ else
+ {
+ cvStartReadSeq( inseq->data.seq, &reader );
+
+ for( i = 0; i < reader.seq->total; i++ )
+ {
+ CvFileNode* inode = (CvFileNode*)reader.ptr;
+ val = inode->data.i;
+ if( CV_NODE_TYPE(inode->tag) != CV_NODE_INT || (unsigned)val >= (unsigned)n )
+ CV_ERROR( CV_StsOutOfRange, "some of in/not_in elements are out of range" );
+
+ split->subset[val >> 5] |= 1 << (val & 31);
+ CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );
+ }
+ }
+
+ // for categorical splits we do not use inversed splits,
+ // instead we inverse the variable set in the split
+ if( inversed )
+ for( i = 0; i < (n + 31) >> 5; i++ )
+ split->subset[i] ^= -1;
+ }
+ else
+ {
+ CvFileNode* cmp_node;
+ split = data->new_split_ord( vi, 0, 0, 0, 0 );
+
+ cmp_node = cvGetFileNodeByName( fs, fnode, "le" );
+ if( !cmp_node )
+ {
+ cmp_node = cvGetFileNodeByName( fs, fnode, "gt" );
+ split->inversed = 1;
+ }
+
+ split->ord.c = (float)cvReadReal( cmp_node );
+ }
+
+ split->quality = (float)cvReadRealByName( fs, fnode, "quality" );
+
+ __END__;
+
+ return split;
+}
+
+
+CvDTreeNode* CvDTree::read_node( CvFileStorage* fs, CvFileNode* fnode, CvDTreeNode* parent )
+{
+ CvDTreeNode* node = 0;
+
+ CV_FUNCNAME( "CvDTree::read_node" );
+
+ __BEGIN__;
+
+ CvFileNode* splits;
+ int i, depth;
+
+ if( !fnode || CV_NODE_TYPE(fnode->tag) != CV_NODE_MAP )
+ CV_ERROR( CV_StsParseError, "some of the tree elements are not stored properly" );
+
+ CV_CALL( node = data->new_node( parent, 0, 0, 0 ));
+ depth = cvReadIntByName( fs, fnode, "depth", -1 );
+ if( depth != node->depth )
+ CV_ERROR( CV_StsParseError, "incorrect node depth" );
+
+ node->sample_count = cvReadIntByName( fs, fnode, "sample_count" );
+ node->value = cvReadRealByName( fs, fnode, "value" );
+ if( data->is_classifier )
+ node->class_idx = cvReadIntByName( fs, fnode, "norm_class_idx" );
+
+ node->Tn = cvReadIntByName( fs, fnode, "Tn" );
+ node->complexity = cvReadIntByName( fs, fnode, "complexity" );
+ node->alpha = cvReadRealByName( fs, fnode, "alpha" );
+ node->node_risk = cvReadRealByName( fs, fnode, "node_risk" );
+ node->tree_risk = cvReadRealByName( fs, fnode, "tree_risk" );
+ node->tree_error = cvReadRealByName( fs, fnode, "tree_error" );
+
+ splits = cvGetFileNodeByName( fs, fnode, "splits" );
+ if( splits )
+ {
+ CvSeqReader reader;
+ CvDTreeSplit* last_split = 0;
+
+ if( CV_NODE_TYPE(splits->tag) != CV_NODE_SEQ )
+ CV_ERROR( CV_StsParseError, "splits tag must stored as a sequence" );
+
+ cvStartReadSeq( splits->data.seq, &reader );
+ for( i = 0; i < reader.seq->total; i++ )
+ {
+ CvDTreeSplit* split;
+ CV_CALL( split = read_split( fs, (CvFileNode*)reader.ptr ));
+ if( !last_split )
+ node->split = last_split = split;
+ else
+ last_split = last_split->next = split;
+
+ CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );
+ }
+ }
+
+ __END__;
+
+ return node;
+}
+
+
+void CvDTree::read_tree_nodes( CvFileStorage* fs, CvFileNode* fnode )
+{
+ CV_FUNCNAME( "CvDTree::read_tree_nodes" );
+
+ __BEGIN__;
+
+ CvSeqReader reader;
+ CvDTreeNode _root;
+ CvDTreeNode* parent = &_root;
+ int i;
+ parent->left = parent->right = parent->parent = 0;
+
+ cvStartReadSeq( fnode->data.seq, &reader );
+
+ for( i = 0; i < reader.seq->total; i++ )
+ {
+ CvDTreeNode* node;
+
+ CV_CALL( node = read_node( fs, (CvFileNode*)reader.ptr, parent != &_root ? parent : 0 ));
+ if( !parent->left )
+ parent->left = node;
+ else
+ parent->right = node;
+ if( node->split )
+ parent = node;
+ else
+ {
+ while( parent && parent->right )
+ parent = parent->parent;
+ }
+
+ CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );
+ }
+
+ root = _root.left;
+
+ __END__;
+}
+
+
+void CvDTree::read( CvFileStorage* fs, CvFileNode* fnode )
+{
+ CvDTreeTrainData* _data = new CvDTreeTrainData();
+ _data->read_params( fs, fnode );
+
+ read( fs, fnode, _data );
+ get_var_importance();
+}
+
+
+// a special entry point for reading weak decision trees from the tree ensembles
+void CvDTree::read( CvFileStorage* fs, CvFileNode* node, CvDTreeTrainData* _data )
+{
+ CV_FUNCNAME( "CvDTree::read" );
+
+ __BEGIN__;
+
+ CvFileNode* tree_nodes;
+
+ clear();
+ data = _data;
+
+ tree_nodes = cvGetFileNodeByName( fs, node, "nodes" );
+ if( !tree_nodes || CV_NODE_TYPE(tree_nodes->tag) != CV_NODE_SEQ )
+ CV_ERROR( CV_StsParseError, "nodes tag is missing" );
+
+ pruned_tree_idx = cvReadIntByName( fs, node, "best_tree_idx", -1 );
+ read_tree_nodes( fs, tree_nodes );
+
+ __END__;
+}
+
+/* End of file. */
projects / opencv.git / blobdiff