← Back to branch summary

~paparazzi-uav/paparazzi/v5.0-manual

« back to all changes in this revision

Viewing changes to sw/ext/opencv_bebop/opencv/apps/traincascade/old_ml_boost.cpp

  • Committer: Paparazzi buildbot
  • Date: 2016-05-18 15:00:29 UTC
  • Revision ID: felix.ruess+docbot@gmail.com-20160518150029-e8lgzi5kvb4p7un9
Manual import commit 4b8bbb730080dac23cf816b98908dacfabe2a8ec from v5.0 branch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
sw/ext/opencv_bebop/opencv/apps/traincascade/old_ml_boost.cpp
 
1
/*M///////////////////////////////////////////////////////////////////////////////////////
 
2
//
 
3
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
 
4
//
 
5
//  By downloading, copying, installing or using the software you agree to this license.
 
6
//  If you do not agree to this license, do not download, install,
 
7
//  copy or use the software.
 
8
//
 
9
//
 
10
//                        Intel License Agreement
 
11
//
 
12
// Copyright (C) 2000, Intel Corporation, all rights reserved.
 
13
// Third party copyrights are property of their respective owners.
 
14
//
 
15
// Redistribution and use in source and binary forms, with or without modification,
 
16
// are permitted provided that the following conditions are met:
 
17
//
 
18
//   * Redistribution's of source code must retain the above copyright notice,
 
19
//     this list of conditions and the following disclaimer.
 
20
//
 
21
//   * Redistribution's in binary form must reproduce the above copyright notice,
 
22
//     this list of conditions and the following disclaimer in the documentation
 
23
//     and/or other materials provided with the distribution.
 
24
//
 
25
//   * The name of Intel Corporation may not be used to endorse or promote products
 
26
//     derived from this software without specific prior written permission.
 
27
//
 
28
// This software is provided by the copyright holders and contributors "as is" and
 
29
// any express or implied warranties, including, but not limited to, the implied
 
30
// warranties of merchantability and fitness for a particular purpose are disclaimed.
 
31
// In no event shall the Intel Corporation or contributors be liable for any direct,
 
32
// indirect, incidental, special, exemplary, or consequential damages
 
33
// (including, but not limited to, procurement of substitute goods or services;
 
34
// loss of use, data, or profits; or business interruption) however caused
 
35
// and on any theory of liability, whether in contract, strict liability,
 
36
// or tort (including negligence or otherwise) arising in any way out of
 
37
// the use of this software, even if advised of the possibility of such damage.
 
38
//
 
39
//M*/
 
40
 
 
41
#include "old_ml_precomp.hpp"
 
42
 
 
43
static inline double
 
44
log_ratio( double val )
 
45
{
 
46
    const double eps = 1e-5;
 
47
 
 
48
    val = MAX( val, eps );
 
49
    val = MIN( val, 1. - eps );
 
50
    return log( val/(1. - val) );
 
51
}
 
52
 
 
53
 
 
54
CvBoostParams::CvBoostParams()
 
55
{
 
56
    boost_type = CvBoost::REAL;
 
57
    weak_count = 100;
 
58
    weight_trim_rate = 0.95;
 
59
    cv_folds = 0;
 
60
    max_depth = 1;
 
61
}
 
62
 
 
63
 
 
64
CvBoostParams::CvBoostParams( int _boost_type, int _weak_count,
 
65
                                        double _weight_trim_rate, int _max_depth,
 
66
                                        bool _use_surrogates, const float* _priors )
 
67
{
 
68
    boost_type = _boost_type;
 
69
    weak_count = _weak_count;
 
70
    weight_trim_rate = _weight_trim_rate;
 
71
    split_criteria = CvBoost::DEFAULT;
 
72
    cv_folds = 0;
 
73
    max_depth = _max_depth;
 
74
    use_surrogates = _use_surrogates;
 
75
    priors = _priors;
 
76
}
 
77
 
 
78
 
 
79
 
 
80
///////////////////////////////// CvBoostTree ///////////////////////////////////
 
81
 
 
82
CvBoostTree::CvBoostTree()
 
83
{
 
84
    ensemble = 0;
 
85
}
 
86
 
 
87
 
 
88
CvBoostTree::~CvBoostTree()
 
89
{
 
90
    clear();
 
91
}
 
92
 
 
93
 
 
94
void
 
95
CvBoostTree::clear()
 
96
{
 
97
    CvDTree::clear();
 
98
    ensemble = 0;
 
99
}
 
100
 
 
101
 
 
102
bool
 
103
CvBoostTree::train( CvDTreeTrainData* _train_data,
 
104
                    const CvMat* _subsample_idx, CvBoost* _ensemble )
 
105
{
 
106
    clear();
 
107
    ensemble = _ensemble;
 
108
    data = _train_data;
 
109
    data->shared = true;
 
110
    return do_train( _subsample_idx );
 
111
}
 
112
 
 
113
 
 
114
bool
 
115
CvBoostTree::train( const CvMat*, int, const CvMat*, const CvMat*,
 
116
                    const CvMat*, const CvMat*, const CvMat*, CvDTreeParams )
 
117
{
 
118
    assert(0);
 
119
    return false;
 
120
}
 
121
 
 
122
 
 
123
bool
 
124
CvBoostTree::train( CvDTreeTrainData*, const CvMat* )
 
125
{
 
126
    assert(0);
 
127
    return false;
 
128
}
 
129
 
 
130
 
 
131
void
 
132
CvBoostTree::scale( double _scale )
 
133
{
 
134
    CvDTreeNode* node = root;
 
135
 
 
136
    // traverse the tree and scale all the node values
 
137
    for(;;)
 
138
    {
 
139
        CvDTreeNode* parent;
 
140
        for(;;)
 
141
        {
 
142
            node->value *= _scale;
 
143
            if( !node->left )
 
144
                break;
 
145
            node = node->left;
 
146
        }
 
147
 
 
148
        for( parent = node->parent; parent && parent->right == node;
 
149
            node = parent, parent = parent->parent )
 
150
            ;
 
151
 
 
152
        if( !parent )
 
153
            break;
 
154
 
 
155
        node = parent->right;
 
156
    }
 
157
}
 
158
 
 
159
 
 
160
void
 
161
CvBoostTree::try_split_node( CvDTreeNode* node )
 
162
{
 
163
    CvDTree::try_split_node( node );
 
164
 
 
165
    if( !node->left )
 
166
    {
 
167
        // if the node has not been split,
 
168
        // store the responses for the corresponding training samples
 
169
        double* weak_eval = ensemble->get_weak_response()->data.db;
 
170
        cv::AutoBuffer<int> inn_buf(node->sample_count);
 
171
        const int* labels = data->get_cv_labels( node, (int*)inn_buf );
 
172
        int i, count = node->sample_count;
 
173
        double value = node->value;
 
174
 
 
175
        for( i = 0; i < count; i++ )
 
176
            weak_eval[labels[i]] = value;
 
177
    }
 
178
}
 
179
 
 
180
 
 
181
double
 
182
CvBoostTree::calc_node_dir( CvDTreeNode* node )
 
183
{
 
184
    char* dir = (char*)data->direction->data.ptr;
 
185
    const double* weights = ensemble->get_subtree_weights()->data.db;
 
186
    int i, n = node->sample_count, vi = node->split->var_idx;
 
187
    double L, R;
 
188
 
 
189
    assert( !node->split->inversed );
 
190
 
 
191
    if( data->get_var_type(vi) >= 0 ) // split on categorical var
 
192
    {
 
193
        cv::AutoBuffer<int> inn_buf(n);
 
194
        const int* cat_labels = data->get_cat_var_data( node, vi, (int*)inn_buf );
 
195
        const int* subset = node->split->subset;
 
196
        double sum = 0, sum_abs = 0;
 
197
 
 
198
        for( i = 0; i < n; i++ )
 
199
        {
 
200
            int idx = ((cat_labels[i] == 65535) && data->is_buf_16u) ? -1 : cat_labels[i];
 
201
            double w = weights[i];
 
202
            int d = idx >= 0 ? CV_DTREE_CAT_DIR(idx,subset) : 0;
 
203
            sum += d*w; sum_abs += (d & 1)*w;
 
204
            dir[i] = (char)d;
 
205
        }
 
206
 
 
207
        R = (sum_abs + sum) * 0.5;
 
208
        L = (sum_abs - sum) * 0.5;
 
209
    }
 
210
    else // split on ordered var
 
211
    {
 
212
        cv::AutoBuffer<uchar> inn_buf(2*n*sizeof(int)+n*sizeof(float));
 
213
        float* values_buf = (float*)(uchar*)inn_buf;
 
214
        int* sorted_indices_buf = (int*)(values_buf + n);
 
215
        int* sample_indices_buf = sorted_indices_buf + n;
 
216
        const float* values = 0;
 
217
        const int* sorted_indices = 0;
 
218
        data->get_ord_var_data( node, vi, values_buf, sorted_indices_buf, &values, &sorted_indices, sample_indices_buf );
 
219
        int split_point = node->split->ord.split_point;
 
220
        int n1 = node->get_num_valid(vi);
 
221
 
 
222
        assert( 0 <= split_point && split_point < n1-1 );
 
223
        L = R = 0;
 
224
 
 
225
        for( i = 0; i <= split_point; i++ )
 
226
        {
 
227
            int idx = sorted_indices[i];
 
228
            double w = weights[idx];
 
229
            dir[idx] = (char)-1;
 
230
            L += w;
 
231
        }
 
232
 
 
233
        for( ; i < n1; i++ )
 
234
        {
 
235
            int idx = sorted_indices[i];
 
236
            double w = weights[idx];
 
237
            dir[idx] = (char)1;
 
238
            R += w;
 
239
        }
 
240
 
 
241
        for( ; i < n; i++ )
 
242
            dir[sorted_indices[i]] = (char)0;
 
243
    }
 
244
 
 
245
    node->maxlr = MAX( L, R );
 
246
    return node->split->quality/(L + R);
 
247
}
 
248
 
 
249
 
 
250
CvDTreeSplit*
 
251
CvBoostTree::find_split_ord_class( CvDTreeNode* node, int vi, float init_quality,
 
252
                                    CvDTreeSplit* _split, uchar* _ext_buf )
 
253
{
 
254
    const float epsilon = FLT_EPSILON*2;
 
255
 
 
256
    const double* weights = ensemble->get_subtree_weights()->data.db;
 
257
    int n = node->sample_count;
 
258
    int n1 = node->get_num_valid(vi);
 
259
 
 
260
    cv::AutoBuffer<uchar> inn_buf;
 
261
    if( !_ext_buf )
 
262
        inn_buf.allocate(n*(3*sizeof(int)+sizeof(float)));
 
263
    uchar* ext_buf = _ext_buf ? _ext_buf : (uchar*)inn_buf;
 
264
    float* values_buf = (float*)ext_buf;
 
265
    int* sorted_indices_buf = (int*)(values_buf + n);
 
266
    int* sample_indices_buf = sorted_indices_buf + n;
 
267
    const float* values = 0;
 
268
    const int* sorted_indices = 0;
 
269
    data->get_ord_var_data( node, vi, values_buf, sorted_indices_buf, &values, &sorted_indices, sample_indices_buf );
 
270
    int* responses_buf = sorted_indices_buf + n;
 
271
    const int* responses = data->get_class_labels( node, responses_buf );
 
272
    const double* rcw0 = weights + n;
 
273
    double lcw[2] = {0,0}, rcw[2];
 
274
    int i, best_i = -1;
 
275
    double best_val = init_quality;
 
276
    int boost_type = ensemble->get_params().boost_type;
 
277
    int split_criteria = ensemble->get_params().split_criteria;
 
278
 
 
279
    rcw[0] = rcw0[0]; rcw[1] = rcw0[1];
 
280
    for( i = n1; i < n; i++ )
 
281
    {
 
282
        int idx = sorted_indices[i];
 
283
        double w = weights[idx];
 
284
        rcw[responses[idx]] -= w;
 
285
    }
 
286
 
 
287
    if( split_criteria != CvBoost::GINI && split_criteria != CvBoost::MISCLASS )
 
288
        split_criteria = boost_type == CvBoost::DISCRETE ? CvBoost::MISCLASS : CvBoost::GINI;
 
289
 
 
290
    if( split_criteria == CvBoost::GINI )
 
291
    {
 
292
        double L = 0, R = rcw[0] + rcw[1];
 
293
        double lsum2 = 0, rsum2 = rcw[0]*rcw[0] + rcw[1]*rcw[1];
 
294
 
 
295
        for( i = 0; i < n1 - 1; i++ )
 
296
        {
 
297
            int idx = sorted_indices[i];
 
298
            double w = weights[idx], w2 = w*w;
 
299
            double lv, rv;
 
300
            idx = responses[idx];
 
301
            L += w; R -= w;
 
302
            lv = lcw[idx]; rv = rcw[idx];
 
303
            lsum2 += 2*lv*w + w2;
 
304
            rsum2 -= 2*rv*w - w2;
 
305
            lcw[idx] = lv + w; rcw[idx] = rv - w;
 
306
 
 
307
            if( values[i] + epsilon < values[i+1] )
 
308
            {
 
309
                double val = (lsum2*R + rsum2*L)/(L*R);
 
310
                if( best_val < val )
 
311
                {
 
312
                    best_val = val;
 
313
                    best_i = i;
 
314
                }
 
315
            }
 
316
        }
 
317
    }
 
318
    else
 
319
    {
 
320
        for( i = 0; i < n1 - 1; i++ )
 
321
        {
 
322
            int idx = sorted_indices[i];
 
323
            double w = weights[idx];
 
324
            idx = responses[idx];
 
325
            lcw[idx] += w;
 
326
            rcw[idx] -= w;
 
327
 
 
328
            if( values[i] + epsilon < values[i+1] )
 
329
            {
 
330
                double val = lcw[0] + rcw[1], val2 = lcw[1] + rcw[0];
 
331
                val = MAX(val, val2);
 
332
                if( best_val < val )
 
333
                {
 
334
                    best_val = val;
 
335
                    best_i = i;
 
336
                }
 
337
            }
 
338
        }
 
339
    }
 
340
 
 
341
    CvDTreeSplit* split = 0;
 
342
    if( best_i >= 0 )
 
343
    {
 
344
        split = _split ? _split : data->new_split_ord( 0, 0.0f, 0, 0, 0.0f );
 
345
        split->var_idx = vi;
 
346
        split->ord.c = (values[best_i] + values[best_i+1])*0.5f;
 
347
        split->ord.split_point = best_i;
 
348
        split->inversed = 0;
 
349
        split->quality = (float)best_val;
 
350
    }
 
351
    return split;
 
352
}
 
353
 
 
354
template<typename T>
 
355
class LessThanPtr
 
356
{
 
357
public:
 
358
    bool operator()(T* a, T* b) const { return *a < *b; }
 
359
};
 
360
 
 
361
CvDTreeSplit*
 
362
CvBoostTree::find_split_cat_class( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split, uchar* _ext_buf )
 
363
{
 
364
    int ci = data->get_var_type(vi);
 
365
    int n = node->sample_count;
 
366
    int mi = data->cat_count->data.i[ci];
 
367
 
 
368
    int base_size = (2*mi+3)*sizeof(double) + mi*sizeof(double*);
 
369
    cv::AutoBuffer<uchar> inn_buf((2*mi+3)*sizeof(double) + mi*sizeof(double*));
 
370
    if( !_ext_buf)
 
371
        inn_buf.allocate( base_size + 2*n*sizeof(int) );
 
372
    uchar* base_buf = (uchar*)inn_buf;
 
373
    uchar* ext_buf = _ext_buf ? _ext_buf : base_buf + base_size;
 
374
 
 
375
    int* cat_labels_buf = (int*)ext_buf;
 
376
    const int* cat_labels = data->get_cat_var_data(node, vi, cat_labels_buf);
 
377
    int* responses_buf = cat_labels_buf + n;
 
378
    const int* responses = data->get_class_labels(node, responses_buf);
 
379
    double lcw[2]={0,0}, rcw[2]={0,0};
 
380
 
 
381
    double* cjk = (double*)cv::alignPtr(base_buf,sizeof(double))+2;
 
382
    const double* weights = ensemble->get_subtree_weights()->data.db;
 
383
    double** dbl_ptr = (double**)(cjk + 2*mi);
 
384
    int i, j, k, idx;
 
385
    double L = 0, R;
 
386
    double best_val = init_quality;
 
387
    int best_subset = -1, subset_i;
 
388
    int boost_type = ensemble->get_params().boost_type;
 
389
    int split_criteria = ensemble->get_params().split_criteria;
 
390
 
 
391
    // init array of counters:
 
392
    // c_{jk} - number of samples that have vi-th input variable = j and response = k.
 
393
    for( j = -1; j < mi; j++ )
 
394
        cjk[j*2] = cjk[j*2+1] = 0;
 
395
 
 
396
    for( i = 0; i < n; i++ )
 
397
    {
 
398
        double w = weights[i];
 
399
        j = ((cat_labels[i] == 65535) && data->is_buf_16u) ? -1 : cat_labels[i];
 
400
        k = responses[i];
 
401
        cjk[j*2 + k] += w;
 
402
    }
 
403
 
 
404
    for( j = 0; j < mi; j++ )
 
405
    {
 
406
        rcw[0] += cjk[j*2];
 
407
        rcw[1] += cjk[j*2+1];
 
408
        dbl_ptr[j] = cjk + j*2 + 1;
 
409
    }
 
410
 
 
411
    R = rcw[0] + rcw[1];
 
412
 
 
413
    if( split_criteria != CvBoost::GINI && split_criteria != CvBoost::MISCLASS )
 
414
        split_criteria = boost_type == CvBoost::DISCRETE ? CvBoost::MISCLASS : CvBoost::GINI;
 
415
 
 
416
    // sort rows of c_jk by increasing c_j,1
 
417
    // (i.e. by the weight of samples in j-th category that belong to class 1)
 
418
    std::sort(dbl_ptr, dbl_ptr + mi, LessThanPtr<double>());
 
419
 
 
420
    for( subset_i = 0; subset_i < mi-1; subset_i++ )
 
421
    {
 
422
        idx = (int)(dbl_ptr[subset_i] - cjk)/2;
 
423
        const double* crow = cjk + idx*2;
 
424
        double w0 = crow[0], w1 = crow[1];
 
425
        double weight = w0 + w1;
 
426
 
 
427
        if( weight < FLT_EPSILON )
 
428
            continue;
 
429
 
 
430
        lcw[0] += w0; rcw[0] -= w0;
 
431
        lcw[1] += w1; rcw[1] -= w1;
 
432
 
 
433
        if( split_criteria == CvBoost::GINI )
 
434
        {
 
435
            double lsum2 = lcw[0]*lcw[0] + lcw[1]*lcw[1];
 
436
            double rsum2 = rcw[0]*rcw[0] + rcw[1]*rcw[1];
 
437
 
 
438
            L += weight;
 
439
            R -= weight;
 
440
 
 
441
            if( L > FLT_EPSILON && R > FLT_EPSILON )
 
442
            {
 
443
                double val = (lsum2*R + rsum2*L)/(L*R);
 
444
                if( best_val < val )
 
445
                {
 
446
                    best_val = val;
 
447
                    best_subset = subset_i;
 
448
                }
 
449
            }
 
450
        }
 
451
        else
 
452
        {
 
453
            double val = lcw[0] + rcw[1];
 
454
            double val2 = lcw[1] + rcw[0];
 
455
 
 
456
            val = MAX(val, val2);
 
457
            if( best_val < val )
 
458
            {
 
459
                best_val = val;
 
460
                best_subset = subset_i;
 
461
            }
 
462
        }
 
463
    }
 
464
 
 
465
    CvDTreeSplit* split = 0;
 
466
    if( best_subset >= 0 )
 
467
    {
 
468
        split = _split ? _split : data->new_split_cat( 0, -1.0f);
 
469
        split->var_idx = vi;
 
470
        split->quality = (float)best_val;
 
471
        memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));
 
472
        for( i = 0; i <= best_subset; i++ )
 
473
        {
 
474
            idx = (int)(dbl_ptr[i] - cjk) >> 1;
 
475
            split->subset[idx >> 5] |= 1 << (idx & 31);
 
476
        }
 
477
    }
 
478
    return split;
 
479
}
 
480
 
 
481
 
 
482
CvDTreeSplit*
 
483
CvBoostTree::find_split_ord_reg( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split, uchar* _ext_buf )
 
484
{
 
485
    const float epsilon = FLT_EPSILON*2;
 
486
    const double* weights = ensemble->get_subtree_weights()->data.db;
 
487
    int n = node->sample_count;
 
488
    int n1 = node->get_num_valid(vi);
 
489
 
 
490
    cv::AutoBuffer<uchar> inn_buf;
 
491
    if( !_ext_buf )
 
492
        inn_buf.allocate(2*n*(sizeof(int)+sizeof(float)));
 
493
    uchar* ext_buf = _ext_buf ? _ext_buf : (uchar*)inn_buf;
 
494
 
 
495
    float* values_buf = (float*)ext_buf;
 
496
    int* indices_buf = (int*)(values_buf + n);
 
497
    int* sample_indices_buf = indices_buf + n;
 
498
    const float* values = 0;
 
499
    const int* indices = 0;
 
500
    data->get_ord_var_data( node, vi, values_buf, indices_buf, &values, &indices, sample_indices_buf );
 
501
    float* responses_buf = (float*)(indices_buf + n);
 
502
    const float* responses = data->get_ord_responses( node, responses_buf, sample_indices_buf );
 
503
 
 
504
    int i, best_i = -1;
 
505
    double L = 0, R = weights[n];
 
506
    double best_val = init_quality, lsum = 0, rsum = node->value*R;
 
507
 
 
508
    // compensate for missing values
 
509
    for( i = n1; i < n; i++ )
 
510
    {
 
511
        int idx = indices[i];
 
512
        double w = weights[idx];
 
513
        rsum -= responses[idx]*w;
 
514
        R -= w;
 
515
    }
 
516
 
 
517
    // find the optimal split
 
518
    for( i = 0; i < n1 - 1; i++ )
 
519
    {
 
520
        int idx = indices[i];
 
521
        double w = weights[idx];
 
522
        double t = responses[idx]*w;
 
523
        L += w; R -= w;
 
524
        lsum += t; rsum -= t;
 
525
 
 
526
        if( values[i] + epsilon < values[i+1] )
 
527
        {
 
528
            double val = (lsum*lsum*R + rsum*rsum*L)/(L*R);
 
529
            if( best_val < val )
 
530
            {
 
531
                best_val = val;
 
532
                best_i = i;
 
533
            }
 
534
        }
 
535
    }
 
536
 
 
537
    CvDTreeSplit* split = 0;
 
538
    if( best_i >= 0 )
 
539
    {
 
540
        split = _split ? _split : data->new_split_ord( 0, 0.0f, 0, 0, 0.0f );
 
541
        split->var_idx = vi;
 
542
        split->ord.c = (values[best_i] + values[best_i+1])*0.5f;
 
543
        split->ord.split_point = best_i;
 
544
        split->inversed = 0;
 
545
        split->quality = (float)best_val;
 
546
    }
 
547
    return split;
 
548
}
 
549
 
 
550
 
 
551
CvDTreeSplit*
 
552
CvBoostTree::find_split_cat_reg( CvDTreeNode* node, int vi, float init_quality, CvDTreeSplit* _split, uchar* _ext_buf )
 
553
{
 
554
    const double* weights = ensemble->get_subtree_weights()->data.db;
 
555
    int ci = data->get_var_type(vi);
 
556
    int n = node->sample_count;
 
557
    int mi = data->cat_count->data.i[ci];
 
558
    int base_size = (2*mi+3)*sizeof(double) + mi*sizeof(double*);
 
559
    cv::AutoBuffer<uchar> inn_buf(base_size);
 
560
    if( !_ext_buf )
 
561
        inn_buf.allocate(base_size + n*(2*sizeof(int) + sizeof(float)));
 
562
    uchar* base_buf = (uchar*)inn_buf;
 
563
    uchar* ext_buf = _ext_buf ? _ext_buf : base_buf + base_size;
 
564
 
 
565
    int* cat_labels_buf = (int*)ext_buf;
 
566
    const int* cat_labels = data->get_cat_var_data(node, vi, cat_labels_buf);
 
567
    float* responses_buf = (float*)(cat_labels_buf + n);
 
568
    int* sample_indices_buf = (int*)(responses_buf + n);
 
569
    const float* responses = data->get_ord_responses(node, responses_buf, sample_indices_buf);
 
570
 
 
571
    double* sum = (double*)cv::alignPtr(base_buf,sizeof(double)) + 1;
 
572
    double* counts = sum + mi + 1;
 
573
    double** sum_ptr = (double**)(counts + mi);
 
574
    double L = 0, R = 0, best_val = init_quality, lsum = 0, rsum = 0;
 
575
    int i, best_subset = -1, subset_i;
 
576
 
 
577
    for( i = -1; i < mi; i++ )
 
578
        sum[i] = counts[i] = 0;
 
579
 
 
580
    // calculate sum response and weight of each category of the input var
 
581
    for( i = 0; i < n; i++ )
 
582
    {
 
583
        int idx = ((cat_labels[i] == 65535) && data->is_buf_16u) ? -1 : cat_labels[i];
 
584
        double w = weights[i];
 
585
        double s = sum[idx] + responses[i]*w;
 
586
        double nc = counts[idx] + w;
 
587
        sum[idx] = s;
 
588
        counts[idx] = nc;
 
589
    }
 
590
 
 
591
    // calculate average response in each category
 
592
    for( i = 0; i < mi; i++ )
 
593
    {
 
594
        R += counts[i];
 
595
        rsum += sum[i];
 
596
        sum[i] = fabs(counts[i]) > DBL_EPSILON ? sum[i]/counts[i] : 0;
 
597
        sum_ptr[i] = sum + i;
 
598
    }
 
599
 
 
600
    std::sort(sum_ptr, sum_ptr + mi, LessThanPtr<double>());
 
601
 
 
602
    // revert back to unnormalized sums
 
603
    // (there should be a very little loss in accuracy)
 
604
    for( i = 0; i < mi; i++ )
 
605
        sum[i] *= counts[i];
 
606
 
 
607
    for( subset_i = 0; subset_i < mi-1; subset_i++ )
 
608
    {
 
609
        int idx = (int)(sum_ptr[subset_i] - sum);
 
610
        double ni = counts[idx];
 
611
 
 
612
        if( ni > FLT_EPSILON )
 
613
        {
 
614
            double s = sum[idx];
 
615
            lsum += s; L += ni;
 
616
            rsum -= s; R -= ni;
 
617
 
 
618
            if( L > FLT_EPSILON && R > FLT_EPSILON )
 
619
            {
 
620
                double val = (lsum*lsum*R + rsum*rsum*L)/(L*R);
 
621
                if( best_val < val )
 
622
                {
 
623
                    best_val = val;
 
624
                    best_subset = subset_i;
 
625
                }
 
626
            }
 
627
        }
 
628
    }
 
629
 
 
630
    CvDTreeSplit* split = 0;
 
631
    if( best_subset >= 0 )
 
632
    {
 
633
        split = _split ? _split : data->new_split_cat( 0, -1.0f);
 
634
        split->var_idx = vi;
 
635
        split->quality = (float)best_val;
 
636
        memset( split->subset, 0, (data->max_c_count + 31)/32 * sizeof(int));
 
637
        for( i = 0; i <= best_subset; i++ )
 
638
        {
 
639
            int idx = (int)(sum_ptr[i] - sum);
 
640
            split->subset[idx >> 5] |= 1 << (idx & 31);
 
641
        }
 
642
    }
 
643
    return split;
 
644
}
 
645
 
 
646
 
 
647
CvDTreeSplit*
 
648
CvBoostTree::find_surrogate_split_ord( CvDTreeNode* node, int vi, uchar* _ext_buf )
 
649
{
 
650
    const float epsilon = FLT_EPSILON*2;
 
651
    int n = node->sample_count;
 
652
    cv::AutoBuffer<uchar> inn_buf;
 
653
    if( !_ext_buf )
 
654
        inn_buf.allocate(n*(2*sizeof(int)+sizeof(float)));
 
655
    uchar* ext_buf = _ext_buf ? _ext_buf : (uchar*)inn_buf;
 
656
    float* values_buf = (float*)ext_buf;
 
657
    int* indices_buf = (int*)(values_buf + n);
 
658
    int* sample_indices_buf = indices_buf + n;
 
659
    const float* values = 0;
 
660
    const int* indices = 0;
 
661
    data->get_ord_var_data( node, vi, values_buf, indices_buf, &values, &indices, sample_indices_buf );
 
662
 
 
663
    const double* weights = ensemble->get_subtree_weights()->data.db;
 
664
    const char* dir = (char*)data->direction->data.ptr;
 
665
    int n1 = node->get_num_valid(vi);
 
666
    // LL - number of samples that both the primary and the surrogate splits send to the left
 
667
    // LR - ... primary split sends to the left and the surrogate split sends to the right
 
668
    // RL - ... primary split sends to the right and the surrogate split sends to the left
 
669
    // RR - ... both send to the right
 
670
    int i, best_i = -1, best_inversed = 0;
 
671
    double best_val;
 
672
    double LL = 0, RL = 0, LR, RR;
 
673
    double worst_val = node->maxlr;
 
674
    double sum = 0, sum_abs = 0;
 
675
    best_val = worst_val;
 
676
 
 
677
    for( i = 0; i < n1; i++ )
 
678
    {
 
679
        int idx = indices[i];
 
680
        double w = weights[idx];
 
681
        int d = dir[idx];
 
682
        sum += d*w; sum_abs += (d & 1)*w;
 
683
    }
 
684
 
 
685
    // sum_abs = R + L; sum = R - L
 
686
    RR = (sum_abs + sum)*0.5;
 
687
    LR = (sum_abs - sum)*0.5;
 
688
 
 
689
    // initially all the samples are sent to the right by the surrogate split,
 
690
    // LR of them are sent to the left by primary split, and RR - to the right.
 
691
    // now iteratively compute LL, LR, RL and RR for every possible surrogate split value.
 
692
    for( i = 0; i < n1 - 1; i++ )
 
693
    {
 
694
        int idx = indices[i];
 
695
        double w = weights[idx];
 
696
        int d = dir[idx];
 
697
 
 
698
        if( d < 0 )
 
699
        {
 
700
            LL += w; LR -= w;
 
701
            if( LL + RR > best_val && values[i] + epsilon < values[i+1] )
 
702
            {
 
703
                best_val = LL + RR;
 
704
                best_i = i; best_inversed = 0;
 
705
            }
 
706
        }
 
707
        else if( d > 0 )
 
708
        {
 
709
            RL += w; RR -= w;
 
710
            if( RL + LR > best_val && values[i] + epsilon < values[i+1] )
 
711
            {
 
712
                best_val = RL + LR;
 
713
                best_i = i; best_inversed = 1;
 
714
            }
 
715
        }
 
716
    }
 
717
 
 
718
    return best_i >= 0 && best_val > node->maxlr ? data->new_split_ord( vi,
 
719
        (values[best_i] + values[best_i+1])*0.5f, best_i,
 
720
        best_inversed, (float)best_val ) : 0;
 
721
}
 
722
 
 
723
 
 
724
CvDTreeSplit*
 
725
CvBoostTree::find_surrogate_split_cat( CvDTreeNode* node, int vi, uchar* _ext_buf )
 
726
{
 
727
    const char* dir = (char*)data->direction->data.ptr;
 
728
    const double* weights = ensemble->get_subtree_weights()->data.db;
 
729
    int n = node->sample_count;
 
730
    int i, mi = data->cat_count->data.i[data->get_var_type(vi)];
 
731
 
 
732
    int base_size = (2*mi+3)*sizeof(double);
 
733
    cv::AutoBuffer<uchar> inn_buf(base_size);
 
734
    if( !_ext_buf )
 
735
        inn_buf.allocate(base_size + n*sizeof(int));
 
736
    uchar* ext_buf = _ext_buf ? _ext_buf : (uchar*)inn_buf;
 
737
    int* cat_labels_buf = (int*)ext_buf;
 
738
    const int* cat_labels = data->get_cat_var_data(node, vi, cat_labels_buf);
 
739
 
 
740
    // LL - number of samples that both the primary and the surrogate splits send to the left
 
741
    // LR - ... primary split sends to the left and the surrogate split sends to the right
 
742
    // RL - ... primary split sends to the right and the surrogate split sends to the left
 
743
    // RR - ... both send to the right
 
744
    CvDTreeSplit* split = data->new_split_cat( vi, 0 );
 
745
    double best_val = 0;
 
746
    double* lc = (double*)cv::alignPtr(cat_labels_buf + n, sizeof(double)) + 1;
 
747
    double* rc = lc + mi + 1;
 
748
 
 
749
    for( i = -1; i < mi; i++ )
 
750
        lc[i] = rc[i] = 0;
 
751
 
 
752
    // 1. for each category calculate the weight of samples
 
753
    // sent to the left (lc) and to the right (rc) by the primary split
 
754
    for( i = 0; i < n; i++ )
 
755
    {
 
756
        int idx = ((cat_labels[i] == 65535) && data->is_buf_16u) ? -1 : cat_labels[i];
 
757
        double w = weights[i];
 
758
        int d = dir[i];
 
759
        double sum = lc[idx] + d*w;
 
760
        double sum_abs = rc[idx] + (d & 1)*w;
 
761
        lc[idx] = sum; rc[idx] = sum_abs;
 
762
    }
 
763
 
 
764
    for( i = 0; i < mi; i++ )
 
765
    {
 
766
        double sum = lc[i];
 
767
        double sum_abs = rc[i];
 
768
        lc[i] = (sum_abs - sum) * 0.5;
 
769
        rc[i] = (sum_abs + sum) * 0.5;
 
770
    }
 
771
 
 
772
    // 2. now form the split.
 
773
    // in each category send all the samples to the same direction as majority
 
774
    for( i = 0; i < mi; i++ )
 
775
    {
 
776
        double lval = lc[i], rval = rc[i];
 
777
        if( lval > rval )
 
778
        {
 
779
            split->subset[i >> 5] |= 1 << (i & 31);
 
780
            best_val += lval;
 
781
        }
 
782
        else
 
783
            best_val += rval;
 
784
    }
 
785
 
 
786
    split->quality = (float)best_val;
 
787
    if( split->quality <= node->maxlr )
 
788
        cvSetRemoveByPtr( data->split_heap, split ), split = 0;
 
789
 
 
790
    return split;
 
791
}
 
792
 
 
793
 
 
794
void
 
795
CvBoostTree::calc_node_value( CvDTreeNode* node )
 
796
{
 
797
    int i, n = node->sample_count;
 
798
    const double* weights = ensemble->get_weights()->data.db;
 
799
    cv::AutoBuffer<uchar> inn_buf(n*(sizeof(int) + ( data->is_classifier ? sizeof(int) : sizeof(int) + sizeof(float))));
 
800
    int* labels_buf = (int*)(uchar*)inn_buf;
 
801
    const int* labels = data->get_cv_labels(node, labels_buf);
 
802
    double* subtree_weights = ensemble->get_subtree_weights()->data.db;
 
803
    double rcw[2] = {0,0};
 
804
    int boost_type = ensemble->get_params().boost_type;
 
805
 
 
806
    if( data->is_classifier )
 
807
    {
 
808
        int* _responses_buf = labels_buf + n;
 
809
        const int* _responses = data->get_class_labels(node, _responses_buf);
 
810
        int m = data->get_num_classes();
 
811
        int* cls_count = data->counts->data.i;
 
812
        for( int k = 0; k < m; k++ )
 
813
            cls_count[k] = 0;
 
814
 
 
815
        for( i = 0; i < n; i++ )
 
816
        {
 
817
            int idx = labels[i];
 
818
            double w = weights[idx];
 
819
            int r = _responses[i];
 
820
            rcw[r] += w;
 
821
            cls_count[r]++;
 
822
            subtree_weights[i] = w;
 
823
        }
 
824
 
 
825
        node->class_idx = rcw[1] > rcw[0];
 
826
 
 
827
        if( boost_type == CvBoost::DISCRETE )
 
828
        {
 
829
            // ignore cat_map for responses, and use {-1,1},
 
830
            // as the whole ensemble response is computes as sign(sum_i(weak_response_i)
 
831
            node->value = node->class_idx*2 - 1;
 
832
        }
 
833
        else
 
834
        {
 
835
            double p = rcw[1]/(rcw[0] + rcw[1]);
 
836
            assert( boost_type == CvBoost::REAL );
 
837
 
 
838
            // store log-ratio of the probability
 
839
            node->value = 0.5*log_ratio(p);
 
840
        }
 
841
    }
 
842
    else
 
843
    {
 
844
        // in case of regression tree:
 
845
        //  * node value is 1/n*sum_i(Y_i), where Y_i is i-th response,
 
846
        //    n is the number of samples in the node.
 
847
        //  * node risk is the sum of squared errors: sum_i((Y_i - <node_value>)^2)
 
848
        double sum = 0, sum2 = 0, iw;
 
849
        float* values_buf = (float*)(labels_buf + n);
 
850
        int* sample_indices_buf = (int*)(values_buf + n);
 
851
        const float* values = data->get_ord_responses(node, values_buf, sample_indices_buf);
 
852
 
 
853
        for( i = 0; i < n; i++ )
 
854
        {
 
855
            int idx = labels[i];
 
856
            double w = weights[idx]/*priors[values[i] > 0]*/;
 
857
            double t = values[i];
 
858
            rcw[0] += w;
 
859
            subtree_weights[i] = w;
 
860
            sum += t*w;
 
861
            sum2 += t*t*w;
 
862
        }
 
863
 
 
864
        iw = 1./rcw[0];
 
865
        node->value = sum*iw;
 
866
        node->node_risk = sum2 - (sum*iw)*sum;
 
867
 
 
868
        // renormalize the risk, as in try_split_node the unweighted formula
 
869
        // sqrt(risk)/n is used, rather than sqrt(risk)/sum(weights_i)
 
870
        node->node_risk *= n*iw*n*iw;
 
871
    }
 
872
 
 
873
    // store summary weights
 
874
    subtree_weights[n] = rcw[0];
 
875
    subtree_weights[n+1] = rcw[1];
 
876
}
 
877
 
 
878
 
 
879
void CvBoostTree::read( CvFileStorage* fs, CvFileNode* fnode, CvBoost* _ensemble, CvDTreeTrainData* _data )
 
880
{
 
881
    CvDTree::read( fs, fnode, _data );
 
882
    ensemble = _ensemble;
 
883
}
 
884
 
 
885
void CvBoostTree::read( CvFileStorage*, CvFileNode* )
 
886
{
 
887
    assert(0);
 
888
}
 
889
 
 
890
void CvBoostTree::read( CvFileStorage* _fs, CvFileNode* _node,
 
891
                        CvDTreeTrainData* _data )
 
892
{
 
893
    CvDTree::read( _fs, _node, _data );
 
894
}
 
895
 
 
896
 
 
897
/////////////////////////////////// CvBoost /////////////////////////////////////
 
898
 
 
899
CvBoost::CvBoost()
 
900
{
 
901
    data = 0;
 
902
    weak = 0;
 
903
    default_model_name = "my_boost_tree";
 
904
 
 
905
    active_vars = active_vars_abs = orig_response = sum_response = weak_eval =
 
906
        subsample_mask = weights = subtree_weights = 0;
 
907
    have_active_cat_vars = have_subsample = false;
 
908
 
 
909
    clear();
 
910
}
 
911
 
 
912
 
 
913
void CvBoost::prune( CvSlice slice )
 
914
{
 
915
    if( weak && weak->total > 0 )
 
916
    {
 
917
        CvSeqReader reader;
 
918
        int i, count = cvSliceLength( slice, weak );
 
919
 
 
920
        cvStartReadSeq( weak, &reader );
 
921
        cvSetSeqReaderPos( &reader, slice.start_index );
 
922
 
 
923
        for( i = 0; i < count; i++ )
 
924
        {
 
925
            CvBoostTree* w;
 
926
            CV_READ_SEQ_ELEM( w, reader );
 
927
            delete w;
 
928
        }
 
929
 
 
930
        cvSeqRemoveSlice( weak, slice );
 
931
    }
 
932
}
 
933
 
 
934
 
 
935
void CvBoost::clear()
 
936
{
 
937
    if( weak )
 
938
    {
 
939
        prune( CV_WHOLE_SEQ );
 
940
        cvReleaseMemStorage( &weak->storage );
 
941
    }
 
942
    if( data )
 
943
        delete data;
 
944
    weak = 0;
 
945
    data = 0;
 
946
    cvReleaseMat( &active_vars );
 
947
    cvReleaseMat( &active_vars_abs );
 
948
    cvReleaseMat( &orig_response );
 
949
    cvReleaseMat( &sum_response );
 
950
    cvReleaseMat( &weak_eval );
 
951
    cvReleaseMat( &subsample_mask );
 
952
    cvReleaseMat( &weights );
 
953
    cvReleaseMat( &subtree_weights );
 
954
 
 
955
    have_subsample = false;
 
956
}
 
957
 
 
958
 
 
959
CvBoost::~CvBoost()
 
960
{
 
961
    clear();
 
962
}
 
963
 
 
964
 
 
965
CvBoost::CvBoost( const CvMat* _train_data, int _tflag,
 
966
                  const CvMat* _responses, const CvMat* _var_idx,
 
967
                  const CvMat* _sample_idx, const CvMat* _var_type,
 
968
                  const CvMat* _missing_mask, CvBoostParams _params )
 
969
{
 
970
    weak = 0;
 
971
    data = 0;
 
972
    default_model_name = "my_boost_tree";
 
973
 
 
974
    active_vars = active_vars_abs = orig_response = sum_response = weak_eval =
 
975
        subsample_mask = weights = subtree_weights = 0;
 
976
 
 
977
    train( _train_data, _tflag, _responses, _var_idx, _sample_idx,
 
978
           _var_type, _missing_mask, _params );
 
979
}
 
980
 
 
981
 
 
982
bool
 
983
CvBoost::set_params( const CvBoostParams& _params )
 
984
{
 
985
    bool ok = false;
 
986
 
 
987
    CV_FUNCNAME( "CvBoost::set_params" );
 
988
 
 
989
    __BEGIN__;
 
990
 
 
991
    params = _params;
 
992
    if( params.boost_type != DISCRETE && params.boost_type != REAL &&
 
993
        params.boost_type != LOGIT && params.boost_type != GENTLE )
 
994
        CV_ERROR( CV_StsBadArg, "Unknown/unsupported boosting type" );
 
995
 
 
996
    params.weak_count = MAX( params.weak_count, 1 );
 
997
    params.weight_trim_rate = MAX( params.weight_trim_rate, 0. );
 
998
    params.weight_trim_rate = MIN( params.weight_trim_rate, 1. );
 
999
    if( params.weight_trim_rate < FLT_EPSILON )
 
1000
        params.weight_trim_rate = 1.f;
 
1001
 
 
1002
    if( params.boost_type == DISCRETE &&
 
1003
        params.split_criteria != GINI && params.split_criteria != MISCLASS )
 
1004
        params.split_criteria = MISCLASS;
 
1005
    if( params.boost_type == REAL &&
 
1006
        params.split_criteria != GINI && params.split_criteria != MISCLASS )
 
1007
        params.split_criteria = GINI;
 
1008
    if( (params.boost_type == LOGIT || params.boost_type == GENTLE) &&
 
1009
        params.split_criteria != SQERR )
 
1010
        params.split_criteria = SQERR;
 
1011
 
 
1012
    ok = true;
 
1013
 
 
1014
    __END__;
 
1015
 
 
1016
    return ok;
 
1017
}
 
1018
 
 
1019
 
 
1020
bool
 
1021
CvBoost::train( const CvMat* _train_data, int _tflag,
 
1022
              const CvMat* _responses, const CvMat* _var_idx,
 
1023
              const CvMat* _sample_idx, const CvMat* _var_type,
 
1024
              const CvMat* _missing_mask,
 
1025
              CvBoostParams _params, bool _update )
 
1026
{
 
1027
    bool ok = false;
 
1028
    CvMemStorage* storage = 0;
 
1029
 
 
1030
    CV_FUNCNAME( "CvBoost::train" );
 
1031
 
 
1032
    __BEGIN__;
 
1033
 
 
1034
    int i;
 
1035
 
 
1036
    set_params( _params );
 
1037
 
 
1038
    cvReleaseMat( &active_vars );
 
1039
    cvReleaseMat( &active_vars_abs );
 
1040
 
 
1041
    if( !_update || !data )
 
1042
    {
 
1043
        clear();
 
1044
        data = new CvDTreeTrainData( _train_data, _tflag, _responses, _var_idx,
 
1045
            _sample_idx, _var_type, _missing_mask, _params, true, true );
 
1046
 
 
1047
        if( data->get_num_classes() != 2 )
 
1048
            CV_ERROR( CV_StsNotImplemented,
 
1049
            "Boosted trees can only be used for 2-class classification." );
 
1050
        CV_CALL( storage = cvCreateMemStorage() );
 
1051
        weak = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvBoostTree*), storage );
 
1052
        storage = 0;
 
1053
    }
 
1054
    else
 
1055
    {
 
1056
        data->set_data( _train_data, _tflag, _responses, _var_idx,
 
1057
            _sample_idx, _var_type, _missing_mask, _params, true, true, true );
 
1058
    }
 
1059
 
 
1060
    if ( (_params.boost_type == LOGIT) || (_params.boost_type == GENTLE) )
 
1061
        data->do_responses_copy();
 
1062
 
 
1063
    update_weights( 0 );
 
1064
 
 
1065
    for( i = 0; i < params.weak_count; i++ )
 
1066
    {
 
1067
        CvBoostTree* tree = new CvBoostTree;
 
1068
        if( !tree->train( data, subsample_mask, this ) )
 
1069
        {
 
1070
            delete tree;
 
1071
            break;
 
1072
        }
 
1073
        //cvCheckArr( get_weak_response());
 
1074
        cvSeqPush( weak, &tree );
 
1075
        update_weights( tree );
 
1076
        trim_weights();
 
1077
        if( cvCountNonZero(subsample_mask) == 0 )
 
1078
            break;
 
1079
    }
 
1080
 
 
1081
    if(weak->total > 0)
 
1082
    {
 
1083
        get_active_vars(); // recompute active_vars* maps and condensed_idx's in the splits.
 
1084
        data->is_classifier = true;
 
1085
        data->free_train_data();
 
1086
        ok = true;
 
1087
    }
 
1088
    else
 
1089
        clear();
 
1090
 
 
1091
    __END__;
 
1092
 
 
1093
    return ok;
 
1094
}
 
1095
 
 
1096
bool CvBoost::train( CvMLData* _data,
 
1097
             CvBoostParams _params,
 
1098
             bool update )
 
1099
{
 
1100
    bool result = false;
 
1101
 
 
1102
    CV_FUNCNAME( "CvBoost::train" );
 
1103
 
 
1104
    __BEGIN__;
 
1105
 
 
1106
    const CvMat* values = _data->get_values();
 
1107
    const CvMat* response = _data->get_responses();
 
1108
    const CvMat* missing = _data->get_missing();
 
1109
    const CvMat* var_types = _data->get_var_types();
 
1110
    const CvMat* train_sidx = _data->get_train_sample_idx();
 
1111
    const CvMat* var_idx = _data->get_var_idx();
 
1112
 
 
1113
    CV_CALL( result = train( values, CV_ROW_SAMPLE, response, var_idx,
 
1114
        train_sidx, var_types, missing, _params, update ) );
 
1115
 
 
1116
    __END__;
 
1117
 
 
1118
    return result;
 
1119
}
 
1120
 
 
1121
void CvBoost::initialize_weights(double (&p)[2])
 
1122
{
 
1123
    p[0] = 1.;
 
1124
    p[1] = 1.;
 
1125
}
 
1126
 
 
1127
void
 
1128
CvBoost::update_weights( CvBoostTree* tree )
 
1129
{
 
1130
    CV_FUNCNAME( "CvBoost::update_weights" );
 
1131
 
 
1132
    __BEGIN__;
 
1133
 
 
1134
    int i, n = data->sample_count;
 
1135
    double sumw = 0.;
 
1136
    int step = 0;
 
1137
    float* fdata = 0;
 
1138
    int *sample_idx_buf;
 
1139
    const int* sample_idx = 0;
 
1140
    cv::AutoBuffer<uchar> inn_buf;
 
1141
    size_t _buf_size = (params.boost_type == LOGIT) || (params.boost_type == GENTLE) ? (size_t)(data->sample_count)*sizeof(int) : 0;
 
1142
    if( !tree )
 
1143
        _buf_size += n*sizeof(int);
 
1144
    else
 
1145
    {
 
1146
        if( have_subsample )
 
1147
            _buf_size += data->get_length_subbuf()*(sizeof(float)+sizeof(uchar));
 
1148
    }
 
1149
    inn_buf.allocate(_buf_size);
 
1150
    uchar* cur_buf_pos = (uchar*)inn_buf;
 
1151
 
 
1152
    if ( (params.boost_type == LOGIT) || (params.boost_type == GENTLE) )
 
1153
    {
 
1154
        step = CV_IS_MAT_CONT(data->responses_copy->type) ?
 
1155
            1 : data->responses_copy->step / CV_ELEM_SIZE(data->responses_copy->type);
 
1156
        fdata = data->responses_copy->data.fl;
 
1157
        sample_idx_buf = (int*)cur_buf_pos;
 
1158
        cur_buf_pos = (uchar*)(sample_idx_buf + data->sample_count);
 
1159
        sample_idx = data->get_sample_indices( data->data_root, sample_idx_buf );
 
1160
    }
 
1161
    CvMat* dtree_data_buf = data->buf;
 
1162
    size_t length_buf_row = data->get_length_subbuf();
 
1163
    if( !tree ) // before training the first tree, initialize weights and other parameters
 
1164
    {
 
1165
        int* class_labels_buf = (int*)cur_buf_pos;
 
1166
        cur_buf_pos = (uchar*)(class_labels_buf + n);
 
1167
        const int* class_labels = data->get_class_labels(data->data_root, class_labels_buf);
 
1168
        // in case of logitboost and gentle adaboost each weak tree is a regression tree,
 
1169
        // so we need to convert class labels to floating-point values
 
1170
 
 
1171
        double w0 = 1./ n;
 
1172
        double p[2] = { 1., 1. };
 
1173
        initialize_weights(p);
 
1174
 
 
1175
        cvReleaseMat( &orig_response );
 
1176
        cvReleaseMat( &sum_response );
 
1177
        cvReleaseMat( &weak_eval );
 
1178
        cvReleaseMat( &subsample_mask );
 
1179
        cvReleaseMat( &weights );
 
1180
        cvReleaseMat( &subtree_weights );
 
1181
 
 
1182
        CV_CALL( orig_response = cvCreateMat( 1, n, CV_32S ));
 
1183
        CV_CALL( weak_eval = cvCreateMat( 1, n, CV_64F ));
 
1184
        CV_CALL( subsample_mask = cvCreateMat( 1, n, CV_8U ));
 
1185
        CV_CALL( weights = cvCreateMat( 1, n, CV_64F ));
 
1186
        CV_CALL( subtree_weights = cvCreateMat( 1, n + 2, CV_64F ));
 
1187
 
 
1188
        if( data->have_priors )
 
1189
        {
 
1190
            // compute weight scale for each class from their prior probabilities
 
1191
            int c1 = 0;
 
1192
            for( i = 0; i < n; i++ )
 
1193
                c1 += class_labels[i];
 
1194
            p[0] = data->priors->data.db[0]*(c1 < n ? 1./(n - c1) : 0.);
 
1195
            p[1] = data->priors->data.db[1]*(c1 > 0 ? 1./c1 : 0.);
 
1196
            p[0] /= p[0] + p[1];
 
1197
            p[1] = 1. - p[0];
 
1198
        }
 
1199
 
 
1200
        if (data->is_buf_16u)
 
1201
        {
 
1202
            unsigned short* labels = (unsigned short*)(dtree_data_buf->data.s + data->data_root->buf_idx*length_buf_row +
 
1203
                data->data_root->offset + (size_t)(data->work_var_count-1)*data->sample_count);
 
1204
            for( i = 0; i < n; i++ )
 
1205
            {
 
1206
                // save original categorical responses {0,1}, convert them to {-1,1}
 
1207
                orig_response->data.i[i] = class_labels[i]*2 - 1;
 
1208
                // make all the samples active at start.
 
1209
                // later, in trim_weights() deactivate/reactive again some, if need
 
1210
                subsample_mask->data.ptr[i] = (uchar)1;
 
1211
                // make all the initial weights the same.
 
1212
                weights->data.db[i] = w0*p[class_labels[i]];
 
1213
                // set the labels to find (from within weak tree learning proc)
 
1214
                // the particular sample weight, and where to store the response.
 
1215
                labels[i] = (unsigned short)i;
 
1216
            }
 
1217
        }
 
1218
        else
 
1219
        {
 
1220
            int* labels = dtree_data_buf->data.i + data->data_root->buf_idx*length_buf_row +
 
1221
                data->data_root->offset + (size_t)(data->work_var_count-1)*data->sample_count;
 
1222
 
 
1223
            for( i = 0; i < n; i++ )
 
1224
            {
 
1225
                // save original categorical responses {0,1}, convert them to {-1,1}
 
1226
                orig_response->data.i[i] = class_labels[i]*2 - 1;
 
1227
                // make all the samples active at start.
 
1228
                // later, in trim_weights() deactivate/reactive again some, if need
 
1229
                subsample_mask->data.ptr[i] = (uchar)1;
 
1230
                // make all the initial weights the same.
 
1231
                weights->data.db[i] = w0*p[class_labels[i]];
 
1232
                // set the labels to find (from within weak tree learning proc)
 
1233
                // the particular sample weight, and where to store the response.
 
1234
                labels[i] = i;
 
1235
            }
 
1236
        }
 
1237
 
 
1238
        if( params.boost_type == LOGIT )
 
1239
        {
 
1240
            CV_CALL( sum_response = cvCreateMat( 1, n, CV_64F ));
 
1241
 
 
1242
            for( i = 0; i < n; i++ )
 
1243
            {
 
1244
                sum_response->data.db[i] = 0;
 
1245
                fdata[sample_idx[i]*step] = orig_response->data.i[i] > 0 ? 2.f : -2.f;
 
1246
            }
 
1247
 
 
1248
            // in case of logitboost each weak tree is a regression tree.
 
1249
            // the target function values are recalculated for each of the trees
 
1250
            data->is_classifier = false;
 
1251
        }
 
1252
        else if( params.boost_type == GENTLE )
 
1253
        {
 
1254
            for( i = 0; i < n; i++ )
 
1255
                fdata[sample_idx[i]*step] = (float)orig_response->data.i[i];
 
1256
 
 
1257
            data->is_classifier = false;
 
1258
        }
 
1259
    }
 
1260
    else
 
1261
    {
 
1262
        // at this moment, for all the samples that participated in the training of the most
 
1263
        // recent weak classifier we know the responses. For other samples we need to compute them
 
1264
        if( have_subsample )
 
1265
        {
 
1266
            float* values = (float*)cur_buf_pos;
 
1267
            cur_buf_pos = (uchar*)(values + data->get_length_subbuf());
 
1268
            uchar* missing = cur_buf_pos;
 
1269
            cur_buf_pos = missing + data->get_length_subbuf() * (size_t)CV_ELEM_SIZE(data->buf->type);
 
1270
 
 
1271
            CvMat _sample, _mask;
 
1272
 
 
1273
            // invert the subsample mask
 
1274
            cvXorS( subsample_mask, cvScalar(1.), subsample_mask );
 
1275
            data->get_vectors( subsample_mask, values, missing, 0 );
 
1276
 
 
1277
            _sample = cvMat( 1, data->var_count, CV_32F );
 
1278
            _mask = cvMat( 1, data->var_count, CV_8U );
 
1279
 
 
1280
            // run tree through all the non-processed samples
 
1281
            for( i = 0; i < n; i++ )
 
1282
                if( subsample_mask->data.ptr[i] )
 
1283
                {
 
1284
                    _sample.data.fl = values;
 
1285
                    _mask.data.ptr = missing;
 
1286
                    values += _sample.cols;
 
1287
                    missing += _mask.cols;
 
1288
                    weak_eval->data.db[i] = tree->predict( &_sample, &_mask, true )->value;
 
1289
                }
 
1290
        }
 
1291
 
 
1292
        // now update weights and other parameters for each type of boosting
 
1293
        if( params.boost_type == DISCRETE )
 
1294
        {
 
1295
            // Discrete AdaBoost:
 
1296
            //   weak_eval[i] (=f(x_i)) is in {-1,1}
 
1297
            //   err = sum(w_i*(f(x_i) != y_i))/sum(w_i)
 
1298
            //   C = log((1-err)/err)
 
1299
            //   w_i *= exp(C*(f(x_i) != y_i))
 
1300
 
 
1301
            double C, err = 0.;
 
1302
            double scale[] = { 1., 0. };
 
1303
 
 
1304
            for( i = 0; i < n; i++ )
 
1305
            {
 
1306
                double w = weights->data.db[i];
 
1307
                sumw += w;
 
1308
                err += w*(weak_eval->data.db[i] != orig_response->data.i[i]);
 
1309
            }
 
1310
 
 
1311
            if( sumw != 0 )
 
1312
                err /= sumw;
 
1313
            C = err = -log_ratio( err );
 
1314
            scale[1] = exp(err);
 
1315
 
 
1316
            sumw = 0;
 
1317
            for( i = 0; i < n; i++ )
 
1318
            {
 
1319
                double w = weights->data.db[i]*
 
1320
                    scale[weak_eval->data.db[i] != orig_response->data.i[i]];
 
1321
                sumw += w;
 
1322
                weights->data.db[i] = w;
 
1323
            }
 
1324
 
 
1325
            tree->scale( C );
 
1326
        }
 
1327
        else if( params.boost_type == REAL )
 
1328
        {
 
1329
            // Real AdaBoost:
 
1330
            //   weak_eval[i] = f(x_i) = 0.5*log(p(x_i)/(1-p(x_i))), p(x_i)=P(y=1|x_i)
 
1331
            //   w_i *= exp(-y_i*f(x_i))
 
1332
 
 
1333
            for( i = 0; i < n; i++ )
 
1334
                weak_eval->data.db[i] *= -orig_response->data.i[i];
 
1335
 
 
1336
            cvExp( weak_eval, weak_eval );
 
1337
 
 
1338
            for( i = 0; i < n; i++ )
 
1339
            {
 
1340
                double w = weights->data.db[i]*weak_eval->data.db[i];
 
1341
                sumw += w;
 
1342
                weights->data.db[i] = w;
 
1343
            }
 
1344
        }
 
1345
        else if( params.boost_type == LOGIT )
 
1346
        {
 
1347
            // LogitBoost:
 
1348
            //   weak_eval[i] = f(x_i) in [-z_max,z_max]
 
1349
            //   sum_response = F(x_i).
 
1350
            //   F(x_i) += 0.5*f(x_i)
 
1351
            //   p(x_i) = exp(F(x_i))/(exp(F(x_i)) + exp(-F(x_i))=1/(1+exp(-2*F(x_i)))
 
1352
            //   reuse weak_eval: weak_eval[i] <- p(x_i)
 
1353
            //   w_i = p(x_i)*1(1 - p(x_i))
 
1354
            //   z_i = ((y_i+1)/2 - p(x_i))/(p(x_i)*(1 - p(x_i)))
 
1355
            //   store z_i to the data->data_root as the new target responses
 
1356
 
 
1357
            const double lb_weight_thresh = FLT_EPSILON;
 
1358
            const double lb_z_max = 10.;
 
1359
            /*float* responses_buf = data->get_resp_float_buf();
 
1360
            const float* responses = 0;
 
1361
            data->get_ord_responses(data->data_root, responses_buf, &responses);*/
 
1362
 
 
1363
            /*if( weak->total == 7 )
 
1364
                putchar('*');*/
 
1365
 
 
1366
            for( i = 0; i < n; i++ )
 
1367
            {
 
1368
                double s = sum_response->data.db[i] + 0.5*weak_eval->data.db[i];
 
1369
                sum_response->data.db[i] = s;
 
1370
                weak_eval->data.db[i] = -2*s;
 
1371
            }
 
1372
 
 
1373
            cvExp( weak_eval, weak_eval );
 
1374
 
 
1375
            for( i = 0; i < n; i++ )
 
1376
            {
 
1377
                double p = 1./(1. + weak_eval->data.db[i]);
 
1378
                double w = p*(1 - p), z;
 
1379
                w = MAX( w, lb_weight_thresh );
 
1380
                weights->data.db[i] = w;
 
1381
                sumw += w;
 
1382
                if( orig_response->data.i[i] > 0 )
 
1383
                {
 
1384
                    z = 1./p;
 
1385
                    fdata[sample_idx[i]*step] = (float)MIN(z, lb_z_max);
 
1386
                }
 
1387
                else
 
1388
                {
 
1389
                    z = 1./(1-p);
 
1390
                    fdata[sample_idx[i]*step] = (float)-MIN(z, lb_z_max);
 
1391
                }
 
1392
            }
 
1393
        }
 
1394
        else
 
1395
        {
 
1396
            // Gentle AdaBoost:
 
1397
            //   weak_eval[i] = f(x_i) in [-1,1]
 
1398
            //   w_i *= exp(-y_i*f(x_i))
 
1399
            assert( params.boost_type == GENTLE );
 
1400
 
 
1401
            for( i = 0; i < n; i++ )
 
1402
                weak_eval->data.db[i] *= -orig_response->data.i[i];
 
1403
 
 
1404
            cvExp( weak_eval, weak_eval );
 
1405
 
 
1406
            for( i = 0; i < n; i++ )
 
1407
            {
 
1408
                double w = weights->data.db[i] * weak_eval->data.db[i];
 
1409
                weights->data.db[i] = w;
 
1410
                sumw += w;
 
1411
            }
 
1412
        }
 
1413
    }
 
1414
 
 
1415
    // renormalize weights
 
1416
    if( sumw > FLT_EPSILON )
 
1417
    {
 
1418
        sumw = 1./sumw;
 
1419
        for( i = 0; i < n; ++i )
 
1420
            weights->data.db[i] *= sumw;
 
1421
    }
 
1422
 
 
1423
    __END__;
 
1424
}
 
1425
 
 
1426
 
 
1427
void
 
1428
CvBoost::trim_weights()
 
1429
{
 
1430
    //CV_FUNCNAME( "CvBoost::trim_weights" );
 
1431
 
 
1432
    __BEGIN__;
 
1433
 
 
1434
    int i, count = data->sample_count, nz_count = 0;
 
1435
    double sum, threshold;
 
1436
 
 
1437
    if( params.weight_trim_rate <= 0. || params.weight_trim_rate >= 1. )
 
1438
        EXIT;
 
1439
 
 
1440
    // use weak_eval as temporary buffer for sorted weights
 
1441
    cvCopy( weights, weak_eval );
 
1442
 
 
1443
    std::sort(weak_eval->data.db, weak_eval->data.db + count);
 
1444
 
 
1445
    // as weight trimming occurs immediately after updating the weights,
 
1446
    // where they are renormalized, we assume that the weight sum = 1.
 
1447
    sum = 1. - params.weight_trim_rate;
 
1448
 
 
1449
    for( i = 0; i < count; i++ )
 
1450
    {
 
1451
        double w = weak_eval->data.db[i];
 
1452
        if( sum <= 0 )
 
1453
            break;
 
1454
        sum -= w;
 
1455
    }
 
1456
 
 
1457
    threshold = i < count ? weak_eval->data.db[i] : DBL_MAX;
 
1458
 
 
1459
    for( i = 0; i < count; i++ )
 
1460
    {
 
1461
        double w = weights->data.db[i];
 
1462
        int f = w >= threshold;
 
1463
        subsample_mask->data.ptr[i] = (uchar)f;
 
1464
        nz_count += f;
 
1465
    }
 
1466
 
 
1467
    have_subsample = nz_count < count;
 
1468
 
 
1469
    __END__;
 
1470
}
 
1471
 
 
1472
 
 
1473
const CvMat*
 
1474
CvBoost::get_active_vars( bool absolute_idx )
 
1475
{
 
1476
    CvMat* mask = 0;
 
1477
    CvMat* inv_map = 0;
 
1478
    CvMat* result = 0;
 
1479
 
 
1480
    CV_FUNCNAME( "CvBoost::get_active_vars" );
 
1481
 
 
1482
    __BEGIN__;
 
1483
 
 
1484
    if( !weak )
 
1485
        CV_ERROR( CV_StsError, "The boosted tree ensemble has not been trained yet" );
 
1486
 
 
1487
    if( !active_vars || !active_vars_abs )
 
1488
    {
 
1489
        CvSeqReader reader;
 
1490
        int i, j, nactive_vars;
 
1491
        CvBoostTree* wtree;
 
1492
        const CvDTreeNode* node;
 
1493
 
 
1494
        assert(!active_vars && !active_vars_abs);
 
1495
        mask = cvCreateMat( 1, data->var_count, CV_8U );
 
1496
        inv_map = cvCreateMat( 1, data->var_count, CV_32S );
 
1497
        cvZero( mask );
 
1498
        cvSet( inv_map, cvScalar(-1) );
 
1499
 
 
1500
        // first pass: compute the mask of used variables
 
1501
        cvStartReadSeq( weak, &reader );
 
1502
        for( i = 0; i < weak->total; i++ )
 
1503
        {
 
1504
            CV_READ_SEQ_ELEM(wtree, reader);
 
1505
 
 
1506
            node = wtree->get_root();
 
1507
            assert( node != 0 );
 
1508
            for(;;)
 
1509
            {
 
1510
                const CvDTreeNode* parent;
 
1511
                for(;;)
 
1512
                {
 
1513
                    CvDTreeSplit* split = node->split;
 
1514
                    for( ; split != 0; split = split->next )
 
1515
                        mask->data.ptr[split->var_idx] = 1;
 
1516
                    if( !node->left )
 
1517
                        break;
 
1518
                    node = node->left;
 
1519
                }
 
1520
 
 
1521
                for( parent = node->parent; parent && parent->right == node;
 
1522
                    node = parent, parent = parent->parent )
 
1523
                    ;
 
1524
 
 
1525
                if( !parent )
 
1526
                    break;
 
1527
 
 
1528
                node = parent->right;
 
1529
            }
 
1530
        }
 
1531
 
 
1532
        nactive_vars = cvCountNonZero(mask);
 
1533
 
 
1534
        //if ( nactive_vars > 0 )
 
1535
        {
 
1536
            active_vars = cvCreateMat( 1, nactive_vars, CV_32S );
 
1537
            active_vars_abs = cvCreateMat( 1, nactive_vars, CV_32S );
 
1538
 
 
1539
            have_active_cat_vars = false;
 
1540
 
 
1541
            for( i = j = 0; i < data->var_count; i++ )
 
1542
            {
 
1543
                if( mask->data.ptr[i] )
 
1544
                {
 
1545
                    active_vars->data.i[j] = i;
 
1546
                    active_vars_abs->data.i[j] = data->var_idx ? data->var_idx->data.i[i] : i;
 
1547
                    inv_map->data.i[i] = j;
 
1548
                    if( data->var_type->data.i[i] >= 0 )
 
1549
                        have_active_cat_vars = true;
 
1550
                    j++;
 
1551
                }
 
1552
            }
 
1553
 
 
1554
 
 
1555
            // second pass: now compute the condensed indices
 
1556
            cvStartReadSeq( weak, &reader );
 
1557
            for( i = 0; i < weak->total; i++ )
 
1558
            {
 
1559
                CV_READ_SEQ_ELEM(wtree, reader);
 
1560
                node = wtree->get_root();
 
1561
                for(;;)
 
1562
                {
 
1563
                    const CvDTreeNode* parent;
 
1564
                    for(;;)
 
1565
                    {
 
1566
                        CvDTreeSplit* split = node->split;
 
1567
                        for( ; split != 0; split = split->next )
 
1568
                        {
 
1569
                            split->condensed_idx = inv_map->data.i[split->var_idx];
 
1570
                            assert( split->condensed_idx >= 0 );
 
1571
                        }
 
1572
 
 
1573
                        if( !node->left )
 
1574
                            break;
 
1575
                        node = node->left;
 
1576
                    }
 
1577
 
 
1578
                    for( parent = node->parent; parent && parent->right == node;
 
1579
                        node = parent, parent = parent->parent )
 
1580
                        ;
 
1581
 
 
1582
                    if( !parent )
 
1583
                        break;
 
1584
 
 
1585
                    node = parent->right;
 
1586
                }
 
1587
            }
 
1588
        }
 
1589
    }
 
1590
 
 
1591
    result = absolute_idx ? active_vars_abs : active_vars;
 
1592
 
 
1593
    __END__;
 
1594
 
 
1595
    cvReleaseMat( &mask );
 
1596
    cvReleaseMat( &inv_map );
 
1597
 
 
1598
    return result;
 
1599
}
 
1600
 
 
1601
 
 
1602
float
 
1603
CvBoost::predict( const CvMat* _sample, const CvMat* _missing,
 
1604
                  CvMat* weak_responses, CvSlice slice,
 
1605
                  bool raw_mode, bool return_sum ) const
 
1606
{
 
1607
    float value = -FLT_MAX;
 
1608
 
 
1609
    CvSeqReader reader;
 
1610
    double sum = 0;
 
1611
    int wstep = 0;
 
1612
    const float* sample_data;
 
1613
 
 
1614
    if( !weak )
 
1615
        CV_Error( CV_StsError, "The boosted tree ensemble has not been trained yet" );
 
1616
 
 
1617
    if( !CV_IS_MAT(_sample) || CV_MAT_TYPE(_sample->type) != CV_32FC1 ||
 
1618
        (_sample->cols != 1 && _sample->rows != 1) ||
 
1619
        (_sample->cols + _sample->rows - 1 != data->var_all && !raw_mode) ||
 
1620
        (active_vars && _sample->cols + _sample->rows - 1 != active_vars->cols && raw_mode) )
 
1621
            CV_Error( CV_StsBadArg,
 
1622
        "the input sample must be 1d floating-point vector with the same "
 
1623
        "number of elements as the total number of variables or "
 
1624
        "as the number of variables used for training" );
 
1625
 
 
1626
    if( _missing )
 
1627
    {
 
1628
        if( !CV_IS_MAT(_missing) || !CV_IS_MASK_ARR(_missing) ||
 
1629
            !CV_ARE_SIZES_EQ(_missing, _sample) )
 
1630
            CV_Error( CV_StsBadArg,
 
1631
            "the missing data mask must be 8-bit vector of the same size as input sample" );
 
1632
    }
 
1633
 
 
1634
    int i, weak_count = cvSliceLength( slice, weak );
 
1635
    if( weak_count >= weak->total )
 
1636
    {
 
1637
        weak_count = weak->total;
 
1638
        slice.start_index = 0;
 
1639
    }
 
1640
 
 
1641
    if( weak_responses )
 
1642
    {
 
1643
        if( !CV_IS_MAT(weak_responses) ||
 
1644
            CV_MAT_TYPE(weak_responses->type) != CV_32FC1 ||
 
1645
            (weak_responses->cols != 1 && weak_responses->rows != 1) ||
 
1646
            weak_responses->cols + weak_responses->rows - 1 != weak_count )
 
1647
            CV_Error( CV_StsBadArg,
 
1648
            "The output matrix of weak classifier responses must be valid "
 
1649
            "floating-point vector of the same number of components as the length of input slice" );
 
1650
        wstep = CV_IS_MAT_CONT(weak_responses->type) ? 1 : weak_responses->step/sizeof(float);
 
1651
    }
 
1652
 
 
1653
    int var_count = active_vars->cols;
 
1654
    const int* vtype = data->var_type->data.i;
 
1655
    const int* cmap = data->cat_map->data.i;
 
1656
    const int* cofs = data->cat_ofs->data.i;
 
1657
 
 
1658
    cv::Mat sample = cv::cvarrToMat(_sample);
 
1659
    cv::Mat missing;
 
1660
    if(!_missing)
 
1661
        missing = cv::cvarrToMat(_missing);
 
1662
 
 
1663
    // if need, preprocess the input vector
 
1664
    if( !raw_mode )
 
1665
    {
 
1666
        int sstep, mstep = 0;
 
1667
        const float* src_sample;
 
1668
        const uchar* src_mask = 0;
 
1669
        float* dst_sample;
 
1670
        uchar* dst_mask;
 
1671
        const int* vidx = active_vars->data.i;
 
1672
        const int* vidx_abs = active_vars_abs->data.i;
 
1673
        bool have_mask = _missing != 0;
 
1674
 
 
1675
        sample = cv::Mat(1, var_count, CV_32FC1);
 
1676
        missing = cv::Mat(1, var_count, CV_8UC1);
 
1677
 
 
1678
        dst_sample = sample.ptr<float>();
 
1679
        dst_mask = missing.ptr<uchar>();
 
1680
 
 
1681
        src_sample = _sample->data.fl;
 
1682
        sstep = CV_IS_MAT_CONT(_sample->type) ? 1 : _sample->step/sizeof(src_sample[0]);
 
1683
 
 
1684
        if( _missing )
 
1685
        {
 
1686
            src_mask = _missing->data.ptr;
 
1687
            mstep = CV_IS_MAT_CONT(_missing->type) ? 1 : _missing->step;
 
1688
        }
 
1689
 
 
1690
        for( i = 0; i < var_count; i++ )
 
1691
        {
 
1692
            int idx = vidx[i], idx_abs = vidx_abs[i];
 
1693
            float val = src_sample[idx_abs*sstep];
 
1694
            int ci = vtype[idx];
 
1695
            uchar m = src_mask ? src_mask[idx_abs*mstep] : (uchar)0;
 
1696
 
 
1697
            if( ci >= 0 )
 
1698
            {
 
1699
                int a = cofs[ci], b = (ci+1 >= data->cat_ofs->cols) ? data->cat_map->cols : cofs[ci+1],
 
1700
                    c = a;
 
1701
                int ival = cvRound(val);
 
1702
                if ( (ival != val) && (!m) )
 
1703
                    CV_Error( CV_StsBadArg,
 
1704
                        "one of input categorical variable is not an integer" );
 
1705
 
 
1706
                while( a < b )
 
1707
                {
 
1708
                    c = (a + b) >> 1;
 
1709
                    if( ival < cmap[c] )
 
1710
                        b = c;
 
1711
                    else if( ival > cmap[c] )
 
1712
                        a = c+1;
 
1713
                    else
 
1714
                        break;
 
1715
                }
 
1716
 
 
1717
                if( c < 0 || ival != cmap[c] )
 
1718
                {
 
1719
                    m = 1;
 
1720
                    have_mask = true;
 
1721
                }
 
1722
                else
 
1723
                {
 
1724
                    val = (float)(c - cofs[ci]);
 
1725
                }
 
1726
            }
 
1727
 
 
1728
            dst_sample[i] = val;
 
1729
            dst_mask[i] = m;
 
1730
        }
 
1731
 
 
1732
        if( !have_mask )
 
1733
            missing.release();
 
1734
    }
 
1735
    else
 
1736
    {
 
1737
        if( !CV_IS_MAT_CONT(_sample->type & (_missing ? _missing->type : -1)) )
 
1738
            CV_Error( CV_StsBadArg, "In raw mode the input vectors must be continuous" );
 
1739
    }
 
1740
 
 
1741
    cvStartReadSeq( weak, &reader );
 
1742
    cvSetSeqReaderPos( &reader, slice.start_index );
 
1743
 
 
1744
    sample_data = sample.ptr<float>();
 
1745
 
 
1746
    if( !have_active_cat_vars && missing.empty() && !weak_responses )
 
1747
    {
 
1748
        for( i = 0; i < weak_count; i++ )
 
1749
        {
 
1750
            CvBoostTree* wtree;
 
1751
            const CvDTreeNode* node;
 
1752
            CV_READ_SEQ_ELEM( wtree, reader );
 
1753
 
 
1754
            node = wtree->get_root();
 
1755
            while( node->left )
 
1756
            {
 
1757
                CvDTreeSplit* split = node->split;
 
1758
                int vi = split->condensed_idx;
 
1759
                float val = sample_data[vi];
 
1760
                int dir = val <= split->ord.c ? -1 : 1;
 
1761
                if( split->inversed )
 
1762
                    dir = -dir;
 
1763
                node = dir < 0 ? node->left : node->right;
 
1764
            }
 
1765
            sum += node->value;
 
1766
        }
 
1767
    }
 
1768
    else
 
1769
    {
 
1770
        const int* avars = active_vars->data.i;
 
1771
        const uchar* m = !missing.empty() ? missing.ptr<uchar>() : 0;
 
1772
 
 
1773
        // full-featured version
 
1774
        for( i = 0; i < weak_count; i++ )
 
1775
        {
 
1776
            CvBoostTree* wtree;
 
1777
            const CvDTreeNode* node;
 
1778
            CV_READ_SEQ_ELEM( wtree, reader );
 
1779
 
 
1780
            node = wtree->get_root();
 
1781
            while( node->left )
 
1782
            {
 
1783
                const CvDTreeSplit* split = node->split;
 
1784
                int dir = 0;
 
1785
                for( ; !dir && split != 0; split = split->next )
 
1786
                {
 
1787
                    int vi = split->condensed_idx;
 
1788
                    int ci = vtype[avars[vi]];
 
1789
                    float val = sample_data[vi];
 
1790
                    if( m && m[vi] )
 
1791
                        continue;
 
1792
                    if( ci < 0 ) // ordered
 
1793
                        dir = val <= split->ord.c ? -1 : 1;
 
1794
                    else // categorical
 
1795
                    {
 
1796
                        int c = cvRound(val);
 
1797
                        dir = CV_DTREE_CAT_DIR(c, split->subset);
 
1798
                    }
 
1799
                    if( split->inversed )
 
1800
                        dir = -dir;
 
1801
                }
 
1802
 
 
1803
                if( !dir )
 
1804
                {
 
1805
                    int diff = node->right->sample_count - node->left->sample_count;
 
1806
                    dir = diff < 0 ? -1 : 1;
 
1807
                }
 
1808
                node = dir < 0 ? node->left : node->right;
 
1809
            }
 
1810
            if( weak_responses )
 
1811
                weak_responses->data.fl[i*wstep] = (float)node->value;
 
1812
            sum += node->value;
 
1813
        }
 
1814
    }
 
1815
 
 
1816
    if( return_sum )
 
1817
        value = (float)sum;
 
1818
    else
 
1819
    {
 
1820
        int cls_idx = sum >= 0;
 
1821
        if( raw_mode )
 
1822
            value = (float)cls_idx;
 
1823
        else
 
1824
            value = (float)cmap[cofs[vtype[data->var_count]] + cls_idx];
 
1825
    }
 
1826
 
 
1827
    return value;
 
1828
}
 
1829
 
 
1830
float CvBoost::calc_error( CvMLData* _data, int type, std::vector<float> *resp )
 
1831
{
 
1832
    float err = 0;
 
1833
    const CvMat* values = _data->get_values();
 
1834
    const CvMat* response = _data->get_responses();
 
1835
    const CvMat* missing = _data->get_missing();
 
1836
    const CvMat* sample_idx = (type == CV_TEST_ERROR) ? _data->get_test_sample_idx() : _data->get_train_sample_idx();
 
1837
    const CvMat* var_types = _data->get_var_types();
 
1838
    int* sidx = sample_idx ? sample_idx->data.i : 0;
 
1839
    int r_step = CV_IS_MAT_CONT(response->type) ?
 
1840
                1 : response->step / CV_ELEM_SIZE(response->type);
 
1841
    bool is_classifier = var_types->data.ptr[var_types->cols-1] == CV_VAR_CATEGORICAL;
 
1842
    int sample_count = sample_idx ? sample_idx->cols : 0;
 
1843
    sample_count = (type == CV_TRAIN_ERROR && sample_count == 0) ? values->rows : sample_count;
 
1844
    float* pred_resp = 0;
 
1845
    if( resp && (sample_count > 0) )
 
1846
    {
 
1847
        resp->resize( sample_count );
 
1848
        pred_resp = &((*resp)[0]);
 
1849
    }
 
1850
    if ( is_classifier )
 
1851
    {
 
1852
        for( int i = 0; i < sample_count; i++ )
 
1853
        {
 
1854
            CvMat sample, miss;
 
1855
            int si = sidx ? sidx[i] : i;
 
1856
            cvGetRow( values, &sample, si );
 
1857
            if( missing )
 
1858
                cvGetRow( missing, &miss, si );
 
1859
            float r = (float)predict( &sample, missing ? &miss : 0 );
 
1860
            if( pred_resp )
 
1861
                pred_resp[i] = r;
 
1862
            int d = fabs((double)r - response->data.fl[si*r_step]) <= FLT_EPSILON ? 0 : 1;
 
1863
            err += d;
 
1864
        }
 
1865
        err = sample_count ? err / (float)sample_count * 100 : -FLT_MAX;
 
1866
    }
 
1867
    else
 
1868
    {
 
1869
        for( int i = 0; i < sample_count; i++ )
 
1870
        {
 
1871
            CvMat sample, miss;
 
1872
            int si = sidx ? sidx[i] : i;
 
1873
            cvGetRow( values, &sample, si );
 
1874
            if( missing )
 
1875
                cvGetRow( missing, &miss, si );
 
1876
            float r = (float)predict( &sample, missing ? &miss : 0 );
 
1877
            if( pred_resp )
 
1878
                pred_resp[i] = r;
 
1879
            float d = r - response->data.fl[si*r_step];
 
1880
            err += d*d;
 
1881
        }
 
1882
        err = sample_count ? err / (float)sample_count : -FLT_MAX;
 
1883
    }
 
1884
    return err;
 
1885
}
 
1886
 
 
1887
void CvBoost::write_params( CvFileStorage* fs ) const
 
1888
{
 
1889
    const char* boost_type_str =
 
1890
        params.boost_type == DISCRETE ? "DiscreteAdaboost" :
 
1891
        params.boost_type == REAL ? "RealAdaboost" :
 
1892
        params.boost_type == LOGIT ? "LogitBoost" :
 
1893
        params.boost_type == GENTLE ? "GentleAdaboost" : 0;
 
1894
 
 
1895
    const char* split_crit_str =
 
1896
        params.split_criteria == DEFAULT ? "Default" :
 
1897
        params.split_criteria == GINI ? "Gini" :
 
1898
        params.boost_type == MISCLASS ? "Misclassification" :
 
1899
        params.boost_type == SQERR ? "SquaredErr" : 0;
 
1900
 
 
1901
    if( boost_type_str )
 
1902
        cvWriteString( fs, "boosting_type", boost_type_str );
 
1903
    else
 
1904
        cvWriteInt( fs, "boosting_type", params.boost_type );
 
1905
 
 
1906
    if( split_crit_str )
 
1907
        cvWriteString( fs, "splitting_criteria", split_crit_str );
 
1908
    else
 
1909
        cvWriteInt( fs, "splitting_criteria", params.split_criteria );
 
1910
 
 
1911
    cvWriteInt( fs, "ntrees", weak->total );
 
1912
    cvWriteReal( fs, "weight_trimming_rate", params.weight_trim_rate );
 
1913
 
 
1914
    data->write_params( fs );
 
1915
}
 
1916
 
 
1917
 
 
1918
void CvBoost::read_params( CvFileStorage* fs, CvFileNode* fnode )
 
1919
{
 
1920
    CV_FUNCNAME( "CvBoost::read_params" );
 
1921
 
 
1922
    __BEGIN__;
 
1923
 
 
1924
    CvFileNode* temp;
 
1925
 
 
1926
    if( !fnode || !CV_NODE_IS_MAP(fnode->tag) )
 
1927
        return;
 
1928
 
 
1929
    data = new CvDTreeTrainData();
 
1930
    CV_CALL( data->read_params(fs, fnode));
 
1931
    data->shared = true;
 
1932
 
 
1933
    params.max_depth = data->params.max_depth;
 
1934
    params.min_sample_count = data->params.min_sample_count;
 
1935
    params.max_categories = data->params.max_categories;
 
1936
    params.priors = data->params.priors;
 
1937
    params.regression_accuracy = data->params.regression_accuracy;
 
1938
    params.use_surrogates = data->params.use_surrogates;
 
1939
 
 
1940
    temp = cvGetFileNodeByName( fs, fnode, "boosting_type" );
 
1941
    if( !temp )
 
1942
        return;
 
1943
 
 
1944
    if( temp && CV_NODE_IS_STRING(temp->tag) )
 
1945
    {
 
1946
        const char* boost_type_str = cvReadString( temp, "" );
 
1947
        params.boost_type = strcmp( boost_type_str, "DiscreteAdaboost" ) == 0 ? DISCRETE :
 
1948
                            strcmp( boost_type_str, "RealAdaboost" ) == 0 ? REAL :
 
1949
                            strcmp( boost_type_str, "LogitBoost" ) == 0 ? LOGIT :
 
1950
                            strcmp( boost_type_str, "GentleAdaboost" ) == 0 ? GENTLE : -1;
 
1951
    }
 
1952
    else
 
1953
        params.boost_type = cvReadInt( temp, -1 );
 
1954
 
 
1955
    if( params.boost_type < DISCRETE || params.boost_type > GENTLE )
 
1956
        CV_ERROR( CV_StsBadArg, "Unknown boosting type" );
 
1957
 
 
1958
    temp = cvGetFileNodeByName( fs, fnode, "splitting_criteria" );
 
1959
    if( temp && CV_NODE_IS_STRING(temp->tag) )
 
1960
    {
 
1961
        const char* split_crit_str = cvReadString( temp, "" );
 
1962
        params.split_criteria = strcmp( split_crit_str, "Default" ) == 0 ? DEFAULT :
 
1963
                                strcmp( split_crit_str, "Gini" ) == 0 ? GINI :
 
1964
                                strcmp( split_crit_str, "Misclassification" ) == 0 ? MISCLASS :
 
1965
                                strcmp( split_crit_str, "SquaredErr" ) == 0 ? SQERR : -1;
 
1966
    }
 
1967
    else
 
1968
        params.split_criteria = cvReadInt( temp, -1 );
 
1969
 
 
1970
    if( params.split_criteria < DEFAULT || params.boost_type > SQERR )
 
1971
        CV_ERROR( CV_StsBadArg, "Unknown boosting type" );
 
1972
 
 
1973
    params.weak_count = cvReadIntByName( fs, fnode, "ntrees" );
 
1974
    params.weight_trim_rate = cvReadRealByName( fs, fnode, "weight_trimming_rate", 0. );
 
1975
 
 
1976
    __END__;
 
1977
}
 
1978
 
 
1979
 
 
1980
 
 
1981
void
 
1982
CvBoost::read( CvFileStorage* fs, CvFileNode* node )
 
1983
{
 
1984
    CV_FUNCNAME( "CvBoost::read" );
 
1985
 
 
1986
    __BEGIN__;
 
1987
 
 
1988
    CvSeqReader reader;
 
1989
    CvFileNode* trees_fnode;
 
1990
    CvMemStorage* storage;
 
1991
    int i, ntrees;
 
1992
 
 
1993
    clear();
 
1994
    read_params( fs, node );
 
1995
 
 
1996
    if( !data )
 
1997
        EXIT;
 
1998
 
 
1999
    trees_fnode = cvGetFileNodeByName( fs, node, "trees" );
 
2000
    if( !trees_fnode || !CV_NODE_IS_SEQ(trees_fnode->tag) )
 
2001
        CV_ERROR( CV_StsParseError, "<trees> tag is missing" );
 
2002
 
 
2003
    cvStartReadSeq( trees_fnode->data.seq, &reader );
 
2004
    ntrees = trees_fnode->data.seq->total;
 
2005
 
 
2006
    if( ntrees != params.weak_count )
 
2007
        CV_ERROR( CV_StsUnmatchedSizes,
 
2008
        "The number of trees stored does not match <ntrees> tag value" );
 
2009
 
 
2010
    CV_CALL( storage = cvCreateMemStorage() );
 
2011
    weak = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvBoostTree*), storage );
 
2012
 
 
2013
    for( i = 0; i < ntrees; i++ )
 
2014
    {
 
2015
        CvBoostTree* tree = new CvBoostTree();
 
2016
        CV_CALL(tree->read( fs, (CvFileNode*)reader.ptr, this, data ));
 
2017
        CV_NEXT_SEQ_ELEM( reader.seq->elem_size, reader );
 
2018
        cvSeqPush( weak, &tree );
 
2019
    }
 
2020
    get_active_vars();
 
2021
 
 
2022
    __END__;
 
2023
}
 
2024
 
 
2025
 
 
2026
void
 
2027
CvBoost::write( CvFileStorage* fs, const char* name ) const
 
2028
{
 
2029
    CV_FUNCNAME( "CvBoost::write" );
 
2030
 
 
2031
    __BEGIN__;
 
2032
 
 
2033
    CvSeqReader reader;
 
2034
    int i;
 
2035
 
 
2036
    cvStartWriteStruct( fs, name, CV_NODE_MAP, CV_TYPE_NAME_ML_BOOSTING );
 
2037
 
 
2038
    if( !weak )
 
2039
        CV_ERROR( CV_StsBadArg, "The classifier has not been trained yet" );
 
2040
 
 
2041
    write_params( fs );
 
2042
    cvStartWriteStruct( fs, "trees", CV_NODE_SEQ );
 
2043
 
 
2044
    cvStartReadSeq( weak, &reader );
 
2045
 
 
2046
    for( i = 0; i < weak->total; i++ )
 
2047
    {
 
2048
        CvBoostTree* tree;
 
2049
        CV_READ_SEQ_ELEM( tree, reader );
 
2050
        cvStartWriteStruct( fs, 0, CV_NODE_MAP );
 
2051
        tree->write( fs );
 
2052
        cvEndWriteStruct( fs );
 
2053
    }
 
2054
 
 
2055
    cvEndWriteStruct( fs );
 
2056
    cvEndWriteStruct( fs );
 
2057
 
 
2058
    __END__;
 
2059
}
 
2060
 
 
2061
 
 
2062
CvMat*
 
2063
CvBoost::get_weights()
 
2064
{
 
2065
    return weights;
 
2066
}
 
2067
 
 
2068
 
 
2069
CvMat*
 
2070
CvBoost::get_subtree_weights()
 
2071
{
 
2072
    return subtree_weights;
 
2073
}
 
2074
 
 
2075
 
 
2076
CvMat*
 
2077
CvBoost::get_weak_response()
 
2078
{
 
2079
    return weak_eval;
 
2080
}
 
2081
 
 
2082
 
 
2083
const CvBoostParams&
 
2084
CvBoost::get_params() const
 
2085
{
 
2086
    return params;
 
2087
}
 
2088
 
 
2089
CvSeq* CvBoost::get_weak_predictors()
 
2090
{
 
2091
    return weak;
 
2092
}
 
2093
 
 
2094
const CvDTreeTrainData* CvBoost::get_data() const
 
2095
{
 
2096
    return data;
 
2097
}
 
2098
 
 
2099
using namespace cv;
 
2100
 
 
2101
CvBoost::CvBoost( const Mat& _train_data, int _tflag,
 
2102
               const Mat& _responses, const Mat& _var_idx,
 
2103
               const Mat& _sample_idx, const Mat& _var_type,
 
2104
               const Mat& _missing_mask,
 
2105
               CvBoostParams _params )
 
2106
{
 
2107
    weak = 0;
 
2108
    data = 0;
 
2109
    default_model_name = "my_boost_tree";
 
2110
    active_vars = active_vars_abs = orig_response = sum_response = weak_eval =
 
2111
        subsample_mask = weights = subtree_weights = 0;
 
2112
 
 
2113
    train( _train_data, _tflag, _responses, _var_idx, _sample_idx,
 
2114
          _var_type, _missing_mask, _params );
 
2115
}
 
2116
 
 
2117
 
 
2118
bool
 
2119
CvBoost::train( const Mat& _train_data, int _tflag,
 
2120
               const Mat& _responses, const Mat& _var_idx,
 
2121
               const Mat& _sample_idx, const Mat& _var_type,
 
2122
               const Mat& _missing_mask,
 
2123
               CvBoostParams _params, bool _update )
 
2124
{
 
2125
    train_data_hdr = _train_data;
 
2126
    train_data_mat = _train_data;
 
2127
    responses_hdr = _responses;
 
2128
    responses_mat = _responses;
 
2129
 
 
2130
    CvMat vidx = _var_idx, sidx = _sample_idx, vtype = _var_type, mmask = _missing_mask;
 
2131
 
 
2132
    return train(&train_data_hdr, _tflag, &responses_hdr, vidx.data.ptr ? &vidx : 0,
 
2133
          sidx.data.ptr ? &sidx : 0, vtype.data.ptr ? &vtype : 0,
 
2134
          mmask.data.ptr ? &mmask : 0, _params, _update);
 
2135
}
 
2136
 
 
2137
float
 
2138
CvBoost::predict( const Mat& _sample, const Mat& _missing,
 
2139
                  const Range& slice, bool raw_mode, bool return_sum ) const
 
2140
{
 
2141
    CvMat sample = _sample, mmask = _missing;
 
2142
    /*if( weak_responses )
 
2143
    {
 
2144
        int weak_count = cvSliceLength( slice, weak );
 
2145
        if( weak_count >= weak->total )
 
2146
        {
 
2147
            weak_count = weak->total;
 
2148
            slice.start_index = 0;
 
2149
        }
 
2150
 
 
2151
        if( !(weak_responses->data && weak_responses->type() == CV_32FC1 &&
 
2152
              (weak_responses->cols == 1 || weak_responses->rows == 1) &&
 
2153
              weak_responses->cols + weak_responses->rows - 1 == weak_count) )
 
2154
            weak_responses->create(weak_count, 1, CV_32FC1);
 
2155
        pwr = &(wr = *weak_responses);
 
2156
    }*/
 
2157
    return predict(&sample, _missing.empty() ? 0 : &mmask, 0,
 
2158
                   slice == Range::all() ? CV_WHOLE_SEQ : cvSlice(slice.start, slice.end),
 
2159
                   raw_mode, return_sum);
 
2160
}
 
2161
 
 
2162
/* End of file. */