classifier_rf.hpp

Go to the documentation of this file.

#ifndef FSTCLASSIFIERRF_H
#define FSTCLASSIFIERRF_H
 /*
=========================================================================
Copyright and license:
  * FST4 software (with the exception of any bundled or externally linked
    third-party components) is Copyright (C) 2000-present Institute of
    Information Theory and Automation of the Czech Academy of Sciences (UTIA).
  * FST4 is free software: you can redistribute it and/or modify it under
    the terms of the GNU Affero General Public License as published by the
    Free Software Foundation, version 3. The full text of the license is in
    the file LICENSE in the distribution root and at
    https://www.gnu.org/licenses/agpl-3.0.html.
  * FST4 is distributed in the hope that it will be useful, but WITHOUT ANY
    WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public
    License for more details.
  * Closed-source use without the AGPL's copyleft obligations requires a
    commercial license from UTIA; in case of interest write to
    utia@utia.cas.cz. LICENSING.md in the distribution root explains both
    options, the contributor license agreement that every contribution
    requires, and the contacts (technical questions and pull requests:
    somol@utia.cas.cz, with "FST" in the subject line).
========================================================================== */
 
#include <memory>
#include <vector>
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cassert>
#include <iostream>
#include "error.hpp"
#include "global.hpp"
#include "clonable.hpp"
#include "classifier.hpp"
 
/*============== Template parameter type naming conventions ==============
--------- Numeric types: -------------------------------------------------
DATATYPE - data sample values - usually real numbers (but may be integers
          in text processing etc.)
REALTYPE - must be real numbers - for representing intermediate results of
          calculations like mean, covariance etc.
IDXTYPE - index values for enumeration of data samples - (nonnegative) integers,
          extent depends on numbers of samples in data
DIMTYPE - index values for enumeration of features (dimensions), or classes (not
          class sizes) - (nonnegative) integers, usually lower extent than IDXTYPE,
          but be aware of expressions like _classes*_features*_features !
          in linearized representations of feature matrices for all classes
BINTYPE - feature selection marker type - represents ca. <10 different feature
          states (selected, deselected, sel./desel. temporarily 1st nested loop, 2nd...)
RETURNTYPE - criterion value: real value, but may be extended in future to support
          multiple values
--------- Class types: ---------------------------------------------------
SUBSET       - class of class type Subset
CLASSIFIER   - class implementing interface defined in abstract class Classifier
EVALUATOR    - class implementing interface defined in abstract class Sequential_Step
DISTANCE     - class implementing interface defined in abstract class Distance
DATAACCESSOR - class implementing interface defined in abstract class Data_Accessor
INTERVALCONTAINER - class of class type Data_Intervaller
CONTAINER    - STL container of class type Data_Interval
========================================================================== */
 
namespace FST {

template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
class Classifier_Random_Forest : public Clonable_Impl<Classifier_Random_Forest<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>,Classifier<RETURNTYPE,DIMTYPE,SUBSET,DATAACCESSOR> > {
public:
        typedef Classifier<RETURNTYPE,DIMTYPE,SUBSET,DATAACCESSOR> parent; 
        typedef Clonable_Impl<Classifier_Random_Forest,parent> parent_impl; 
        typedef std::shared_ptr<DATAACCESSOR> PDataAccessor;      
        typedef std::shared_ptr<SUBSET> PSubset;                  
        typedef typename DATAACCESSOR::PPattern PPattern;         
        Classifier_Random_Forest() : _tree_count(100), _mtry(0), _max_depth(0), _min_leaf(1), _seed_salt(1),
                _classes(0), _features(0), _subfeatures(0) {notify("Classifier_Random_Forest constructor.");}
        virtual ~Classifier_Random_Forest() {notify("Classifier_Random_Forest destructor.");}
 
        void set_tree_count(const unsigned int t) {assert(t>0); _tree_count=t;} 
        void set_mtry(const DIMTYPE m) {_mtry=m;}                    
        void set_max_depth(const unsigned int dp) {_max_depth=dp;}   
        void set_min_leaf(const IDXTYPE ml) {assert(ml>0); _min_leaf=ml;} 
        void set_seed_salt(const unsigned int s) {_seed_salt=s;}     
 
        void initialize(const PDataAccessor da); 
 
        virtual bool classify(DIMTYPE &cls, const PPattern &pattern) override;  // majority tree vote; ties resolve to the lowest class index
        virtual bool classify_scores(std::vector<RETURNTYPE> &scores, const PPattern &pattern) override; // per-class tree-vote fractions for ranking metrics (AUC)
        virtual bool train(const PDataAccessor da, const PSubset sub) override; // grows the forest on the designated training part
        // test() is the accuracy loop inherited from the Classifier base
 
        virtual std::ostream& print(std::ostream& os) const override
                {os << "Classifier_Random_Forest(trees=" << _tree_count << ", mtry=" << _mtry << ", max_depth=" << _max_depth << ", min_leaf=" << _min_leaf << ")"; return os;}
protected:
        struct Node {
                DIMTYPE feature;    
                REALTYPE threshold; 
                long left;          
                long right;         
                DIMTYPE leaf_class; 
        };
 
        // configuration
        unsigned int _tree_count; 
        DIMTYPE _mtry;             
        unsigned int _max_depth;   
        IDXTYPE _min_leaf;         
        unsigned int _seed_salt;   
 
        // data shape + trained state
        DIMTYPE _classes;                   
        DIMTYPE _features;                  
        DIMTYPE _subfeatures;               
        std::vector<DIMTYPE> _index;        
        std::vector<IDXTYPE> _votes;        
        std::vector<std::vector<Node> > _forest; 
 
        // training scratch (rebuilt per train(); copied by clones harmlessly)
        std::vector<REALTYPE> _X;           
        std::vector<DIMTYPE> _y;            
 
        // hash-seeded local RNG (xorshift64*) - deliberately NOT the thread_local fst_rand stream
        std::uint64_t _rng_state; 
        std::uint64_t rng_next() {std::uint64_t &s=_rng_state; s^=s>>12; s^=s<<25; s^=s>>27; return s*2685821657736338717ULL;} 
        IDXTYPE rng_below(const IDXTYPE n) {assert(n>0); return (IDXTYPE)(rng_next()%(std::uint64_t)n);} 

        void grow_tree(std::vector<Node> &tree, std::vector<IDXTYPE> &samples, const DIMTYPE mtry);
        DIMTYPE tree_vote(const std::vector<Node> &tree, const PPattern &pattern) const;
};
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
void Classifier_Random_Forest<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::initialize(const PDataAccessor da)
{
        notify("Classifier_Random_Forest::initialize().");
        assert(da);
        assert(da->getNoOfFeatures()>0);
        if(da->getNoOfClasses()<2) throw fst_error("Classifier_Random_Forest requires at least 2 classes");
        _classes=da->getNoOfClasses();
        _features=da->getNoOfFeatures();
        _index.resize(_features);
        _subfeatures=0; // no model yet
        _forest.clear();
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_Random_Forest<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::train(const PDataAccessor da, const PSubset sub)
{
        notify("Classifier_Random_Forest::train().");
        assert(da);
        assert(sub);
        assert(_classes==da->getNoOfClasses());
        assert(_features==da->getNoOfFeatures());
        assert(sub->get_n_raw()==da->getNoOfFeatures());
 
        // buffer the feature subset index
        _subfeatures=0;
        for(DIMTYPE fsel : sub->selected_features()) {assert(_subfeatures<_features); _index[_subfeatures++]=fsel;}
        assert(_subfeatures==sub->get_d_raw());
        assert(_subfeatures>0);
        const std::size_t d=_subfeatures;
 
        // gather raw selected training values and labels (CART is scale-invariant — no standardization)
        const DIMTYPE stride=da->getNoOfFeatures();
        _X.clear(); _y.clear();
        typename DATAACCESSOR::PPattern p;
        IDXTYPE s,i;
        for(DIMTYPE c=0;c<_classes;c++)
        {
                da->setClass(c);
                for(const auto &blk : da->blocks(TRAIN)) for(i=0,p=blk.data,s=blk.size;i<s;i++)
                {
                        const auto *pat=&p[(std::size_t)i*stride];
                        for(std::size_t f=0;f<d;f++) _X.push_back((REALTYPE)pat[_index[f]]);
                        _y.push_back(c);
                }
        }
        const IDXTYPE N=(IDXTYPE)_y.size();
        if(N==0) return false;
 
        // seed the local RNG from (subset content, `--seed` salt, split index): FNV-1a 64-bit
        std::uint64_t h=14695981039346656037ULL;
        const std::uint64_t prime=1099511628211ULL;
        auto mix=[&h](std::uint64_t v) {for(int b=0;b<8;b++) {h^=(v&0xffu); h*=prime; v>>=8;}};
        mix((std::uint64_t)_seed_salt);
        mix((std::uint64_t)da->getSplitIndex());
        for(std::size_t f=0;f<d;f++) mix((std::uint64_t)_index[f]);
        if(h==0) h=0x9E3779B97F4A7C15ULL; // xorshift state must be nonzero
        _rng_state=h;
 
        // per-node feature-subsample size
        DIMTYPE mtry=_mtry;
        if(mtry==0) mtry=(DIMTYPE)std::floor(std::sqrt((double)d));
        if(mtry<1) mtry=1;
        if(mtry>(DIMTYPE)d) mtry=(DIMTYPE)d;
 
        // grow the forest on bootstrap samples (sequential, fixed order — deterministic)
        _forest.assign(_tree_count,std::vector<Node>());
        std::vector<IDXTYPE> samples(N);
        for(unsigned int t=0;t<_tree_count;t++)
        {
                for(IDXTYPE n=0;n<N;n++) samples[n]=rng_below(N);
                grow_tree(_forest[t],samples,mtry);
        }
        return true;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
void Classifier_Random_Forest<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::grow_tree(std::vector<Node> &tree, std::vector<IDXTYPE> &samples, const DIMTYPE mtry)
{
        const std::size_t d=_subfeatures;
        const DIMTYPE K=_classes;
        struct Work {long node; IDXTYPE lo, hi; unsigned int depth;}; // samples[lo..hi) belong to the node
        std::vector<Work> stack;
        std::vector<IDXTYPE> cnt(K), cntL(K);
        std::vector<DIMTYPE> feats(d);
        std::vector<IDXTYPE> sortbuf;
 
        tree.clear();
        tree.push_back(Node());
        stack.push_back(Work{0,0,(IDXTYPE)samples.size(),0});
        while(!stack.empty())
        {
                const Work w=stack.back(); stack.pop_back();
                const IDXTYPE n=w.hi-w.lo;
                assert(n>0);
 
                // class histogram; majority with lowest-index tie-break
                std::fill(cnt.begin(),cnt.end(),(IDXTYPE)0);
                for(IDXTYPE ii=w.lo;ii<w.hi;ii++) cnt[_y[samples[ii]]]++;
                DIMTYPE major=0;
                for(DIMTYPE c=1;c<K;c++) if(cnt[c]>cnt[major]) major=c;
                const bool pure=(cnt[major]==n);
 
                // n<2*_min_leaf is a shortcut only: no split of n samples can give both children at
                // least _min_leaf. The per-split-point check further below enforces _min_leaf exactly.
                bool make_leaf=pure || n<2*_min_leaf || (_max_depth>0 && w.depth>=_max_depth);
                DIMTYPE best_f=0; REALTYPE best_thr=0.0, best_gain=(REALTYPE)-1.0;
                if(!make_leaf)
                {
                        // draw mtry distinct candidate features (partial Fisher-Yates, fixed order)
                        for(std::size_t f=0;f<d;f++) feats[f]=(DIMTYPE)f;
                        for(DIMTYPE k=0;k<mtry;k++) {const IDXTYPE j=k+rng_below((IDXTYPE)(d-k)); std::swap(feats[k],feats[j]);}
 
                        for(DIMTYPE k=0;k<mtry;k++)
                        {
                                const DIMTYPE f=feats[k];
                                sortbuf.assign(samples.begin()+w.lo,samples.begin()+w.hi);
                                // total order (value, sample index) — deterministic under equal values
                                std::sort(sortbuf.begin(),sortbuf.end(),[this,d,f](const IDXTYPE a, const IDXTYPE b)
                                        {const REALTYPE va=_X[(std::size_t)a*d+f], vb=_X[(std::size_t)b*d+f]; return va<vb || (va==vb && a<b);});
                                // sweep split points between distinct consecutive values; with sL,sR the sums of
                                // squared class counts on each side, maximizing sL/nL + sR/nR is exactly the same
                                // ordering as minimizing the weighted Gini impurity (nL*giniL + nR*giniR), and
                                // avoids the divisions the impurity form would need
                                std::fill(cntL.begin(),cntL.end(),(IDXTYPE)0);
                                for(IDXTYPE ii=0;ii+1<n;ii++)
                                {
                                        cntL[_y[sortbuf[ii]]]++;
                                        const REALTYPE v0=_X[(std::size_t)sortbuf[ii]*d+f], v1=_X[(std::size_t)sortbuf[ii+1]*d+f];
                                        if(!(v0<v1)) continue; // tie group — no split point here
                                        const IDXTYPE nL=ii+1, nR=n-nL;
                                        if(nL<_min_leaf || nR<_min_leaf) continue;
                                        REALTYPE sL=0.0, sR=0.0;
                                        for(DIMTYPE c=0;c<K;c++) {
                                                const REALTYPE l=(REALTYPE)cntL[c], r=(REALTYPE)(cnt[c]-cntL[c]);
                                                sL+=l*l; sR+=r*r;
                                        }
                                        const REALTYPE gain=sL/(REALTYPE)nL+sR/(REALTYPE)nR;
                                        if(gain>best_gain) { // strict: equal gain keeps the first candidate encountered
                                                REALTYPE thr=(v0+v1)/(REALTYPE)2.0;
                                                // FP-midpoint guard: for adjacent representable values the midpoint can round
                                                // up to v1, which would put the v1 samples on the left as well and empty the
                                                // right child; falling back to v0 keeps the partition as the sweep assumed
                                                if(!(thr<v1)) thr=v0;
                                                best_gain=gain; best_f=f; best_thr=thr;
                                        }
                                }
                        }
                        if(best_gain<(REALTYPE)0.0) make_leaf=true; // no valid split point on any candidate (e.g. all-tied values)
                }
 
                if(make_leaf)
                {
                        tree[w.node].left=-1; tree[w.node].right=-1;
                        tree[w.node].feature=0; tree[w.node].threshold=0.0;
                        tree[w.node].leaf_class=major;
                        continue;
                }
 
                // partition the slice (stable — deterministic order for the children)
                auto mid=std::stable_partition(samples.begin()+w.lo,samples.begin()+w.hi,
                        [this,d,best_f,best_thr](const IDXTYPE a){return _X[(std::size_t)a*d+best_f]<=best_thr;});
                const IDXTYPE cut=(IDXTYPE)(mid-samples.begin());
                assert(cut>w.lo && cut<w.hi); // guaranteed by the split search (both children hold >=_min_leaf>=1 samples);
                                              // were it ever violated with asserts disabled, an empty child slice would
                                              // trip the assert(n>0) at the top of the next iteration
 
                const long left=(long)tree.size(); tree.push_back(Node());
                const long right=(long)tree.size(); tree.push_back(Node());
                tree[w.node].feature=best_f; tree[w.node].threshold=best_thr;
                tree[w.node].left=left; tree[w.node].right=right; tree[w.node].leaf_class=0;
                stack.push_back(Work{right,cut,w.hi,w.depth+1});
                stack.push_back(Work{left,w.lo,cut,w.depth+1});
        }
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
DIMTYPE Classifier_Random_Forest<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::tree_vote(const std::vector<Node> &tree, const PPattern &pattern) const
{
        long node=0;
        while(tree[node].left>=0)
        {
                const Node &nd=tree[node];
                node=((REALTYPE)pattern[_index[nd.feature]]<=nd.threshold)?nd.left:nd.right;
        }
        return tree[node].leaf_class;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_Random_Forest<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::classify(DIMTYPE &cls, const PPattern &pattern)
{
        assert(_subfeatures>0 && !_forest.empty()); // train() must have succeeded
        _votes.assign(_classes,(IDXTYPE)0); // member scratch: avoids a per-pattern allocation
        for(const auto &tree : _forest) _votes[tree_vote(tree,pattern)]++;
        cls=0;
        for(DIMTYPE c=1;c<_classes;c++) if(_votes[c]>_votes[cls]) cls=c; // ties resolve to the lowest class index
        return true;
}
 
template<class RETURNTYPE, typename DATATYPE, typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET, class DATAACCESSOR>
bool Classifier_Random_Forest<RETURNTYPE,DATATYPE,REALTYPE,IDXTYPE,DIMTYPE,SUBSET,DATAACCESSOR>::classify_scores(std::vector<RETURNTYPE> &scores, const PPattern &pattern)
{
        assert(_subfeatures>0 && !_forest.empty());
        scores.assign(_classes,(RETURNTYPE)0.0);
        for(const auto &tree : _forest) scores[tree_vote(tree,pattern)]+=1.0;
        for(DIMTYPE c=0;c<_classes;c++) scores[c]/=(RETURNTYPE)_forest.size();
        return true;
}
 
 
} // namespace
#endif // FSTCLASSIFIERRF_H ///:~