search_branch_and_bound.hpp

Go to the documentation of this file.

#ifndef FSTSEARCHBRANCHANDBOUND_H
#define FSTSEARCHBRANCHANDBOUND_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 <cassert>
#include <iostream>
#include <sstream>
#include <vector>
#include <list>
#include "error.hpp"
#include "global.hpp"
#include "stopwatch.hpp"
#include "search.hpp"
#include "result_tracker.hpp"
#include "candidate_evaluator.hpp"
#include "candidate_evaluator_sequential.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 DIMTYPE, class SUBSET, class CRITERION>
class Search_Branch_And_Bound : public Search<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>{
public:
        typedef Search<RETURNTYPE,DIMTYPE,SUBSET,CRITERION> parent; 
        typedef std::shared_ptr<CRITERION> PCriterion; 
        typedef std::shared_ptr<SUBSET> PSubset; 
        typedef Candidate_Evaluator<RETURNTYPE,DIMTYPE,SUBSET,CRITERION> EVALUATORENGINE; 
        typedef std::shared_ptr<EVALUATORENGINE> PEvaluatorEngine; 
        explicit Search_Branch_And_Bound(PEvaluatorEngine engine=PEvaluatorEngine()):Search<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>(),_engine(engine) {
                if(!_engine) _engine.reset(new Candidate_Evaluator_Sequential<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>());
                notify("Search_Branch_And_Bound constructor.");
        }
        virtual ~Search_Branch_And_Bound() {notify("Search_Branch_And_Bound destructor.");}
 
        // NOTE: crit is expected to be already initialize()d before use here
        virtual bool search(const DIMTYPE target_d, RETURNTYPE &result, const PSubset sub, const PCriterion crit, std::ostream& os=std::cout) override; 
 
        virtual std::ostream& print(std::ostream& os) const override {os << "Search_Branch_And_Bound() abstract class"; return os;};
 
protected:
        typedef enum {UNKNOWN, RANDOM, PREDICTED, COMPUTED} NodeType;
        typedef struct {
                DIMTYPE    feature;  
                RETURNTYPE value;    
                NodeType   type;     
                bool       valid;    
        } Node;
        typedef Node* PNode; 
        typedef std::list<Node> ListNode; 
        typedef std::vector<ListNode> VectorListNode; 
        typedef typename ListNode::iterator ListNodeIter; 
        typedef typename std::vector<ListNode>::iterator TreeLevelIter; 
        typedef std::vector<Node> VectorNode; 
 
private:
        // the following structures need to be valid throughout the course of one search
        PCriterion _criterion; 
 
        PEvaluatorEngine _engine; 
 
        VectorListNode _tree; 
 
        ListNodeIter   _nodeiter; 
        DIMTYPE        _availiter; 
 
        VectorNode     _available; 
        DIMTYPE         get_available_size() {DIMTYPE c=0; for(DIMTYPE i=0;i<_available.size();i++) if(_available[i].valid) c++; return c;} 
        void            available_add(const DIMTYPE feature) {assert(_available.size()>feature && !_available[feature].valid); _available[feature].valid=true; _available[feature].type=UNKNOWN;} 
        void            available_remove(const DIMTYPE feature) {assert(_available.size()>feature && _available[feature].valid); _available[feature].valid=false;} 
 
        PSubset        _currentset; 
        void            currentset_add(const DIMTYPE feature) {assert(!_currentset->selected_raw(feature)); _currentset->select_raw(feature);} 
        void            currentset_remove(const DIMTYPE feature) {assert(_currentset->selected_raw(feature)); _currentset->deselect_raw(feature);} 
 
        PSubset        _tempset; 
 
        PSubset        _boundset; 
        RETURNTYPE     _boundvalue; 
        bool           _boundvalid; 
        
        DIMTYPE         k; 
        DIMTYPE         n; 
        DIMTYPE         d; 
        
protected:
        PCriterion const &get_criterion() const {assert(_criterion); return _criterion;} 
        PEvaluatorEngine const &get_engine() const {assert(_engine); return _engine;} 
 
        RETURNTYPE get_bound_value() const {return _boundvalue;} 
        bool       is_bound_valid() const {return _boundvalid;} 
 
        int     get_k() const {return k;} 
        DIMTYPE get_n() const {return n;} 
        DIMTYPE get_d() const {return d;} 
 
        bool getFirstNode(PNode &nod) {assert(k<n-d+1 && _tree.size()==n-d); if(_tree[k].empty()) return false; else {_nodeiter=_tree[k].begin(); nod=&(*_nodeiter); return true;}} 
        bool getNextNode(PNode &nod) {assert(k<n-d+1 && _tree.size()==n-d); if(++_nodeiter==_tree[k].end()) return false; else {nod=&(*_nodeiter); return true;}} 
 
        bool getFirstAvailable(PNode &avail) {assert(_available.size()>0); for(_availiter=0; _availiter<_available.size() && !_available[_availiter].valid; _availiter++) {} if(_availiter<_available.size()) {avail=&(_available[_availiter]); return true;} else return false;} 
        bool getNextAvailable(PNode &avail) {assert(_available.size()>0 && _availiter<_available.size()); _availiter++; while(_availiter<_available.size() && !_available[_availiter].valid) {_availiter++;} if(_availiter<_available.size()) {avail=&(_available[_availiter]); return true;} else return false;} 
        bool getAvailable(DIMTYPE idx, PNode &avail) {assert(idx<_available.size()); if(idx<_available.size()) {avail=&(_available[idx]); return true;} else return false;} 
 
protected:
        PSubset& get_currentset(){return _currentset;} 
        Node&    get_parent_node(){assert(k>0 && _tree.size()>=k && !_tree[k-1].empty()); return _tree[k-1].back();} 
        Node&    get_current_node(){assert(_tree.size()>k && !_tree[k].empty()); return _tree[k].back();} 
        
        std::ostream* shared_os; 
        
        void setup_default_structures(const DIMTYPE d, const DIMTYPE n);
        void update_bound(const RETURNTYPE val, const PSubset &sub);
        void make_tree_level();
        void process_rightmost_path();
 
        std::shared_ptr<StopWatch> swatch; 
protected:
        // the following may/must be overriden in descendant Branch and Bound specialization classes (to implement specific functionality or threading)
        virtual void process_leafs(); 
        virtual void leaf_result(const DIMTYPE feature, const RETURNTYPE value) {} 
        virtual void initialize(const DIMTYPE d, const DIMTYPE n, const PCriterion crit)=0;
        virtual void pre_evaluate_availables()=0;
        virtual void post_process_tree_level()=0;
        virtual bool cut_possible()=0;
};
 
template<class RETURNTYPE, typename DIMTYPE, class SUBSET, class CRITERION>
void Search_Branch_And_Bound<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>::setup_default_structures(const DIMTYPE d, const DIMTYPE n)
{
        assert(d>0);
        assert(d<n);
        _tree.resize(n-d); // logically n-d+1 levels, but the last (leaf) level is not stored - its nodes are evaluated directly by process_leafs()
        for(TreeLevelIter i=_tree.begin();i!=_tree.end();i++) i->clear();
        _available.resize(n); for(DIMTYPE f=0;f<n;f++) {_available[f].feature=f; _available[f].valid=true;}
        _currentset.reset(new SUBSET(n)); _currentset->select_all_raw();
        _tempset.reset(new SUBSET(n));
        _boundset.reset(new SUBSET(n)); _boundvalid=false;
}
 
template<class RETURNTYPE, typename DIMTYPE, class SUBSET, class CRITERION>
void Search_Branch_And_Bound<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>::update_bound(const RETURNTYPE val, const PSubset &sub)
{
        assert(sub);
        assert(_boundset);
        assert(sub->get_n()==_boundset->get_n());
        assert(swatch);
        if(!_boundvalid || val>_boundvalue) {
                _boundvalue=val;
                _boundset->stateless_copy(*sub);
                _boundvalid=true;
                if(parent::output_normal()) this->report_solution(*shared_os, val, *sub, msgbuf() << "NEW BOUND (val="<<val<<"): " << *sub << std::endl << *swatch << std::endl << std::endl << std::flush);
        }
        if(parent::result_tracker_active()) parent::_tracker->add(val,sub);
}       
 
template<class RETURNTYPE, typename DIMTYPE, class SUBSET, class CRITERION>
void Search_Branch_And_Bound<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>::make_tree_level() // make k-th tree level
{
        // based on k and no. of _available features
        assert(k>0);
        assert(k<n-d+1);
        assert(_tree.size()==n-d);
        assert(_tree[k].empty());
        assert(_available.size()==n);
        
        pre_evaluate_availables(); // assign values to candidates to enable sorting into the tree level to be constructed here
        
        const DIMTYPE levelsize=get_available_size()-(n-d-k+1);
        assert(levelsize>0);
        assert(levelsize<=get_available_size());
 
        // sort-in levelsize features according to value (in case of type UNKNOWN feature put at the end of list)
        PNode avail;
        for(bool got=getFirstAvailable(avail);got;got=getNextAvailable(avail))
        {
                DIMTYPE pos=0;
                ListNodeIter it=_tree[k].begin();
                while(it!=_tree[k].end() && it->value<avail->value) {++it; ++pos;}
                if(pos<levelsize) {
                        _tree[k].insert(it,*avail);
                        if(_tree[k].size()>levelsize) _tree[k].pop_back();
                }
        }
        assert(_tree[k].size()==levelsize);
        for(ListNodeIter it=_tree[k].begin();it!=_tree[k].end();it++) {
                available_remove(it->feature);
        }
        
        post_process_tree_level(); // enable to re-evaluate chosen candidate values (e.g., supply missing true _criterion values, etc.)
}
 
template<class RETURNTYPE, typename DIMTYPE, class SUBSET, class CRITERION>
void Search_Branch_And_Bound<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>::process_rightmost_path()
{
        assert(_tempset);
        assert(_currentset);
        assert(_available.size()==n);
        assert(get_available_size()>0);
        assert(_criterion);
        _tempset->stateless_copy(*_currentset);
        
        RETURNTYPE val;
        PNode avail;
        
        for(bool got=getFirstAvailable(avail);got;got=getNextAvailable(avail)) _tempset->deselect_raw(avail->feature);
        if(!_criterion->evaluate(val,_tempset)) throw FST::fst_error("Criterion evaluation failure."); 
        update_bound(val,_tempset); //adds to tracker
        if(parent::output_detailed()) this->report(*shared_os, msgbuf() << "rightmost leaf (val="<<val<<"): " << *_tempset << std::endl << std::endl << std::flush);
}
 
template<class RETURNTYPE, typename DIMTYPE, class SUBSET, class CRITERION>
void Search_Branch_And_Bound<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>::process_leafs()
{
        assert(_currentset);
        assert(get_criterion());
        assert(_engine);
        PNode avail;
        // bound updates and progress output fire once per leaf in enumeration order: the
        // sequential engine fires them interleaved with the evaluations, the threaded engine
        // once the batch has completed - either way the order, and hence the result, is the same
        typename EVALUATORENGINE::ResultCallback on_result=[this](const RETURNTYPE value, const PSubset &csub, const DIMTYPE feat) {
                update_bound(value,csub); //adds to tracker
                if(parent::output_detailed()) this->report(*shared_os, msgbuf() << "leaf (val="<<value<<"): " << *csub << std::endl << std::endl << std::flush);
                leaf_result(feat,value); // descendant hook (FBB/BBPP prediction learning) - fires in candidate order under both engine implementations
        };
        _engine->initialize(get_criterion(), get_n(), on_result);
        for(bool got=getFirstAvailable(avail);got;got=getNextAvailable(avail))
        {
                currentset_remove(avail->feature);
                if(!_engine->add_candidate(_currentset,avail->feature)) throw FST::fst_error("Criterion evaluation failure.");
                currentset_add(avail->feature);
        }
        _engine->collect(); // (return value carries no error here: sequential evaluation failures throw at add_candidate above, threaded ones throw in the worker)
}
 
template<class RETURNTYPE, typename DIMTYPE, class SUBSET, class CRITERION>
bool Search_Branch_And_Bound<RETURNTYPE,DIMTYPE,SUBSET,CRITERION>::search(const DIMTYPE target_d, RETURNTYPE &result, const PSubset sub, const PCriterion crit, std::ostream& os) // returns found subset + _criterion value
{
        swatch.reset(new StopWatch);
        notify("Search_Branch_And_Bound::search().");
        assert(sub);
        assert(crit);
        n=sub->get_n();
        assert(target_d>0 && target_d<n);
        d=target_d;
        assert(d>0 && d<=n);
        shared_os=&os;
 
        if(parent::result_tracker_active()) parent::_tracker->set_output_detail(parent::get_output_detail());
                
        if(parent::output_normal()) {
                msgbuf sos;
                sos << "---------------------------------------" << std::endl;
                sos << "Starting " << *this << std::endl;
                sos << "with Criterion: " << *crit << std::endl;
                sos << "Target subset size set to: " << target_d << std::endl;
                sos << std::endl << std::flush;
                this->report(os,sos);
        }
        // user-supplied criteria default to is_monotonic()==false and get this warning
        // automatically; the criteria shipped with FST4 that are monotonic declare the
        // property, so they stay silent
        if(parent::output_normal() && !crit->is_monotonic()) {
                this->report(os, msgbuf() << "WARNING: the supplied criterion does not declare monotonicity w.r.t. subset cardinality" << std::endl
                   << "(Criterion::is_monotonic()==false). Branch & Bound optimality guarantees do NOT apply -" << std::endl
                   << "the result may be suboptimal. Override is_monotonic() to true only if removing a feature" << std::endl
                   << "provably never increases the criterion value; otherwise prefer a non-optimal search." << std::endl << std::endl << std::flush);
        }
 
        if(n==d) {
                sub->select_all_raw();
                if(!crit->evaluate(result,sub)) throw FST::fst_error("Criterion evaluation failure."); 
                return true;
        }
        
        setup_default_structures(d,n);
        
        // initialize root
        Node root; // Node::feature and Node::valid are not used in root
        _tempset->select_all_raw(); 
        if(!crit->evaluate(root.value,_tempset)) throw FST::fst_error("Criterion evaluation failure."); 
        root.type=COMPUTED;
        _tree[0].push_back(root);
                
        initialize(d,n,crit); // to enable overloaded functionality
        _criterion=crit;
                
        // Deepest level at which a CUT message was last printed. Cuts near the leaves are both
        // frequent and uninteresting, so a CUT is reported only when it happens at that level or
        // higher up the tree; this throttles the message stream without affecting the search.
        DIMTYPE lastcut_k=n-d;
 
        DIMTYPE feature;
        // TRAVERSAL INVARIANT, maintained by every block below. At any point between iterations:
        //  - the features removed on the path from the root to level k are exactly the features
        //    deselected in _currentset, one per level 1..k-1, namely _tree[j].back().feature;
        //  - a feature is in exactly one of three places: selected in _currentset, sorted into some
        //    level of _tree, or valid in _available. Descending one level moves _tree[k].back().feature
        //    out of _currentset; backtracking one level moves it back into _currentset AND back into
        //    _available, and pops it off the level, so the level's next node becomes the current one.
        //  - a level that becomes empty means its subtree is fully explored, hence the "pop and go up"
        //    loops keep climbing while _tree[k] is empty. Reaching k==0 ends the search.
        k=1; // start with the first level below root
        do{
                // traverse tree branch down
                while(k>0 && k<n-d) {
 
                        make_tree_level(); // make k-th level, choose |_available|-(n-d-k+1) nodes from _available and use them to make next level
                        process_rightmost_path();
                        
                        // recursively cut all sub-branches that can be cut
                        while(k>0 && cut_possible())
                        {
                                assert(!_tree[k].empty());
                                if(k<=lastcut_k) { // report only cuts at or above the last reported level (see lastcut_k)
                                        if(parent::output_normal()) this->report(os, msgbuf() << "CUT (k="<<k<<", val="<<_tree[k].back().value<<", bound="<<get_bound_value()<<", feature=" << _tree[k].back().feature<<")  " << std::endl << std::endl);
                                        lastcut_k=k;
                                }
                                feature=_tree[k].back().feature; assert(feature<n);
                                available_add(feature);
                                _tree[k].pop_back();
                                // the cut emptied this level: climb until a level with an untried node is found,
                                // restoring each passed level's feature into _currentset and _available on the way
                                while(k>0 && _tree[k].empty()) {
                                        --k;
                                        if(k>0) {
                                                assert(!_tree[k].empty());
                                                feature=_tree[k].back().feature; assert(feature<n);
                                                currentset_add(feature);
                                                available_add(feature);
                                                _tree[k].pop_back();
                                        }
                                } 
                        }
                        // no more cutting possible, proceed down to next level (or finish if k==0)
                        if(k>0) {
                                feature=_tree[k].back().feature; assert(feature<n);
                                currentset_remove(feature);
                                ++k;
                        }
                }
                // evaluate all leafs that are descendants of the just visited node
                if(k==n-d) {
                        process_leafs(); //depends on _currentset, _available, not on k
                        // the leaf level has no stored nodes, so there is nothing left to try here:
                        // backtrack up to the closest ancestor that still has an untried node, undoing
                        // each level's feature removal as we climb
                        do {
                                --k;
                                if(k>0) {
                                        assert(!_tree[k].empty());
                                        feature=_tree[k].back().feature; assert(feature<n);
                                        currentset_add(feature);
                                        available_add(feature);
                                        _tree[k].pop_back();
                                }
                        } while(k>0 && _tree[k].empty());
                }
                // proceed down through the next rightmost subbranch to the left from the last one
                if(k>0) {
                        assert(!_tree[k].empty());
                        feature=_tree[k].back().feature; assert(feature<n);
                        currentset_remove(feature);
                        ++k;
                }
                
        } while(k>0);
 
        assert(_boundset);
        assert(sub->get_n()==_boundset->get_n());
        sub->stateless_copy(*_boundset);
        result=_boundvalue;
        if(parent::output_normal()) {
                msgbuf sos; 
                sos << "......................................." << std::endl;
                sos << "Branch and Bound search finished. " << *swatch << std::endl;
                sos << "Search result: "<< std::endl << *sub << std::endl << "Criterion value: " << result << std::endl << std::endl << std::flush;
                this->report(os,sos);
        }
        swatch.reset();
        return true;
}
 
} // namespace
#endif // FSTSEARCHBRANCHANDBOUND_H ///:~