result_tracker_feature_stats.hpp

Go to the documentation of this file.

#ifndef FSTRESULTTRACKERFEATURESTATS_H
#define FSTRESULTTRACKERFEATURESTATS_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>
#ifdef FST_THREADED
#include <mutex>
#endif
#include <iostream>
#include <sstream>
#include <list>
#include <vector>
#include <algorithm>
#include <cstddef>
#include "error.hpp"
#include "global.hpp"
#include "result_tracker_dupless.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 IDXTYPE, typename DIMTYPE, class SUBSET>
class Result_Tracker_Feature_Stats : public Clonable_Impl<Result_Tracker_Feature_Stats<RETURNTYPE,IDXTYPE,DIMTYPE,SUBSET>,Result_Tracker_Dupless<RETURNTYPE,IDXTYPE,DIMTYPE,SUBSET>,CLONE_INHERIT,CLONE_INHERIT,CLONE_DEEP> {
public:
        typedef Result_Tracker_Dupless<RETURNTYPE,IDXTYPE,DIMTYPE,SUBSET> parent; 
        typedef Clonable_Impl<Result_Tracker_Feature_Stats,parent,CLONE_INHERIT,CLONE_INHERIT,CLONE_DEEP> parent_impl; 
        typedef typename parent::ResultRec ResultRec; 
        typedef ResultRec* PResultRec; 
        typedef std::shared_ptr<SUBSET> PSubset; 
        Result_Tracker_Feature_Stats(const IDXTYPE capacity_limit=0):parent_impl(capacity_limit), _n(0), _pair_itersize(0) {_itersize[0]=0; _itersize[1]=0; _itersize[2]=0; notify("Result_Tracker_Feature_Stats constructor.");}
        Result_Tracker_Feature_Stats(const Result_Tracker_Feature_Stats& rtfs):parent_impl(rtfs), _n(0), _pair_itersize(0) {_itersize[0]=0; _itersize[1]=0; _itersize[2]=0; notify("Result_Tracker_Feature_Stats copy-constructor.");} 
        virtual ~Result_Tracker_Feature_Stats() {notify("Result_Tracker_Feature_Stats destructor.");}

        bool compute_stats(std::ostream& os=std::cout);
        bool print_stats(std::ostream& os=std::cout);

        bool getFirstDAF(RETURNTYPE &value, DIMTYPE &feature, const unsigned int DAFidx=0) const;
        bool getNextDAF(RETURNTYPE &value, DIMTYPE &feature, const unsigned int DAFidx=0) const;
 
        // ---- opt-in pairwise (group) DAF statistics ----
        bool compute_pair_stats(const DIMTYPE top_k, std::ostream& os=std::cout);
        bool print_pair_stats(std::ostream& os=std::cout, const std::size_t top_pairs=0);
        bool getFirstPairDAF(RETURNTYPE &daf0, RETURNTYPE &daf1, bool &daf1_valid, RETURNTYPE &synergy, DIMTYPE &f1, DIMTYPE &f2) const;
        bool getNextPairDAF(RETURNTYPE &daf0, RETURNTYPE &daf1, bool &daf1_valid, RETURNTYPE &synergy, DIMTYPE &f1, DIMTYPE &f2) const;
 
 
 
        virtual std::ostream& print(std::ostream& os) const override {os << "Result_Tracker_Feature_Stats(limit=" << parent::_capacity_limit << ") size " << parent::results.size(); return os;}
 
protected:
        DIMTYPE _n; 
        typedef struct {
                RETURNTYPE freq_is;    
                RETURNTYPE freq_isnot; 
                RETURNTYPE mean_is;    
                RETURNTYPE mean_isnot; 
                RETURNTYPE stdev_is;   
                RETURNTYPE stdev_isnot; 
                RETURNTYPE daf[3];
                bool daf_valid[3]; 
        } FeatureStat;
        std::vector<FeatureStat> _stats; 

        typedef struct {
                RETURNTYPE freq;  
                RETURNTYPE mean;  
                RETURNTYPE stdev; 
                bool valid;       
        } SubSizeStat;
        std::vector<SubSizeStat> _dstat; 
 
        typedef std::vector<DIMTYPE> ORDERTYPE; 
        ORDERTYPE _order[3];  
        mutable DIMTYPE _itersize[3]; 
#ifdef FST_THREADED
        mutable std::mutex _itersize_mutex; 
#endif

        typedef struct {
                DIMTYPE f1;              
                DIMTYPE f2;              
                RETURNTYPE freq_both;    
                RETURNTYPE freq_neither; 
                RETURNTYPE daf0;         
                RETURNTYPE daf1;         
                RETURNTYPE synergy;      
                bool valid;              
                bool daf1_valid;         
        } PairStat;
        std::vector<PairStat> _pair_stats;      
        std::vector<DIMTYPE> _pair_features;    
        std::vector<std::size_t> _pair_order;   
        mutable std::size_t _pair_itersize;     
        static std::size_t pair_index(const std::size_t a, const std::size_t b, const std::size_t k) {assert(a<b && b<k); return a*k - a*(a+1)/2 + (b-a-1);}
};
 
 
 
template<class RETURNTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET>
bool Result_Tracker_Feature_Stats<RETURNTYPE, IDXTYPE, DIMTYPE, SUBSET>::compute_stats(std::ostream& os)
{
        if(parent::results.empty()) return false;
        
        typename parent::CONSTRESULTSITER iter=parent::results.begin();
        // get n from first recorded sub (assumes all others have the same n)
        _n=iter->sub->get_n();
        _stats.resize(_n);
        for(DIMTYPE f=0;f<_n;f++) {
                _stats[f].freq_is=_stats[f].freq_isnot=_stats[f].mean_is=_stats[f].mean_isnot=_stats[f].stdev_is=_stats[f].stdev_isnot=0;
                _stats[f].daf_valid[0]=false;
                _stats[f].daf_valid[1]=false;
                _stats[f].daf_valid[2]=false;
        }
        _dstat.resize(_n+1);
        for(DIMTYPE ss=1;ss<=_n;ss++) {
                _dstat[ss].mean=_dstat[ss].stdev=_dstat[ss].freq=0;
                _dstat[ss].valid=false;
        }
 
        if(parent::output_normal()) {std::ostringstream sos; sos << std::endl << "compute_stats() 1/4.." << std::endl << std::flush; syncout::print(os,sos);}
        // gather initial stats
        for(iter=parent::results.begin();iter!=parent::results.end();iter++)
        {
                assert(iter->sub->get_n()==_n);
                DIMTYPE d=iter->sub->get_d();
                assert(d>0 && d<=_n);
                for(DIMTYPE f=0;f<_n;f++) {
                        if(iter->sub->selected_raw(f)) {
                                _stats[f].freq_is++;
                                _stats[f].mean_is+=iter->value;
                        } else {
                                _stats[f].freq_isnot++;
                                _stats[f].mean_isnot+=iter->value;
                        }
                }
                _dstat[d].mean+=iter->value;
                _dstat[d].freq++;
        }
        
        // calculate mean values
        for(DIMTYPE f=0;f<_n;f++)
        {
                if(_stats[f].freq_is>0) _stats[f].mean_is/=_stats[f].freq_is;
                if(_stats[f].freq_isnot>0) _stats[f].mean_isnot/=_stats[f].freq_isnot;
        }
        for(DIMTYPE ss=1;ss<=_n;ss++) if(_dstat[ss].freq>0) {
                _dstat[ss].mean/=_dstat[ss].freq;
        }
        if(parent::output_normal()) {std::ostringstream sos; sos << "compute_stats() 2/4.." << std::endl << std::flush; syncout::print(os,sos);}
                
        // gather more stats
        for(iter=parent::results.begin();iter!=parent::results.end();iter++)
        {
                assert(iter->sub->get_n()==_n);
                DIMTYPE d=iter->sub->get_d();
                assert(d>0 && d<=_n);
                for(DIMTYPE f=0;f<_n;f++) {
                        if(iter->sub->selected_raw(f)) {
                                _stats[f].stdev_is+=(_stats[f].mean_is-iter->value)*(_stats[f].mean_is-iter->value);
                        } else {
                                _stats[f].stdev_isnot+=(_stats[f].mean_isnot-iter->value)*(_stats[f].mean_isnot-iter->value);
                        }
                }
                _dstat[d].stdev+=(_dstat[d].mean-iter->value)*(_dstat[d].mean-iter->value);
        }
        
        // calculate stddev values
        for(DIMTYPE f=0;f<_n;f++)
        {
                if(_stats[f].freq_is>0) _stats[f].stdev_is=sqrt(_stats[f].stdev_is/_stats[f].freq_is);
                // each guard must test the SAME frequency its line divides by: guarding the 'isnot'
                // line with freq_is would leave stdev_isnot unnormalized, or divide it by zero
                if(_stats[f].freq_isnot>0) _stats[f].stdev_isnot=sqrt(_stats[f].stdev_isnot/_stats[f].freq_isnot);
        }
        for(DIMTYPE ss=1;ss<=_n;ss++) if(_dstat[ss].freq>0) {
                _dstat[ss].stdev=sqrt(_dstat[ss].stdev/_dstat[ss].freq);
                _dstat[ss].valid=true;
        }
        
        
        if(parent::output_normal()) {std::ostringstream sos; sos << "compute_stats() 3/4.." << std::endl << std::flush; syncout::print(os,sos);}
        // calculate DAF ranks
        for(DIMTYPE f=0;f<_n;f++)
        {
                // DAF0 needs only both partitions to be non-empty, but it deliberately shares DAF1's
                // validity test: a feature whose criterion values never vary at all carries no evidence,
                // so ranking it alongside features that do would be misleading
                if(_stats[f].freq_is>0 && _stats[f].freq_isnot>0 && (_stats[f].stdev_is>0 || _stats[f].stdev_isnot>0) ) 
                {
                        RETURNTYPE freq=_stats[f].freq_is+_stats[f].freq_isnot;
                        //DAF1
                        _stats[f].daf[1]=(_stats[f].mean_is-_stats[f].mean_isnot)/( (_stats[f].freq_is/freq)*_stats[f].stdev_is + (_stats[f].freq_isnot/freq)*_stats[f].stdev_isnot );
                        _stats[f].daf_valid[1]=true;
                        //DAF0
                        _stats[f].daf[0]=_stats[f].mean_is-_stats[f].mean_isnot;
                        _stats[f].daf_valid[0]=true;
                } 
                else _stats[f].daf_valid[1]=_stats[f].daf_valid[0]=false;
        }
        // DAF2: like DAF0, but every probe's value is first divided by the stdev of all probes of
        // the SAME cardinality, which removes the systematic level difference between small and
        // large subsets. Probes of a cardinality whose stdev is 0 contribute nothing, so a feature
        // can end up with no usable probe on one side - its DAF2 is then 0 (but still marked
        // valid, unlike DAF0/DAF1). mi1/mi2/cnt1/cnt2 are indexed by FEATURE, 0.._n-1.
        std::vector<RETURNTYPE> mi1,mi2; mi1.resize(_n); mi2.resize(_n);
        std::vector<IDXTYPE> cnt1,cnt2; cnt1.resize(_n); cnt2.resize(_n);
        for(DIMTYPE f=0;f<_n;f++) {mi1[f]=mi2[f]=0; cnt1[f]=cnt2[f]=0;}
        for(iter=parent::results.begin();iter!=parent::results.end();iter++)
        {
                assert(iter->sub->get_n()==_n);
                DIMTYPE d=iter->sub->get_d();
                assert(d>0 && d<=_n);
                
                if(_dstat[d].stdev>0) {
                        for(DIMTYPE f=0;f<_n;f++) {
                                if(iter->sub->selected_raw(f)) {mi1[f]+=iter->value/_dstat[d].stdev; cnt1[f]++;}
                                else {mi2[f]+=iter->value/_dstat[d].stdev; cnt2[f]++;}
                        }
                }
        }       
        for(DIMTYPE f=0;f<_n;f++) {
                if(cnt1[f]==0 || cnt2[f]==0) _stats[f].daf[2]=0; else
                {
                        mi1[f]/=cnt1[f];
                        mi2[f]/=cnt2[f];
                        _stats[f].daf[2]=mi1[f]-mi2[f];
                }
                _stats[f].daf_valid[2]=true;
        }
 
        if(parent::output_normal()) {std::ostringstream sos; sos << "compute_stats() 4/4.." << std::endl << std::endl << std::flush; syncout::print(os,sos);}
        
        // order features according to DAF0, DAF1 and DAF2
        for(unsigned int DAFidx=0; DAFidx<3; DAFidx++)
        {
                _order[DAFidx].clear();
                for(DIMTYPE f=0;f<_n;f++)
                {
                        if(_stats[f].daf_valid[DAFidx]) {
                                typename ORDERTYPE::iterator iter=_order[DAFidx].begin();
                                while(iter!=_order[DAFidx].end() && _stats[f].daf[DAFidx]<_stats[*iter].daf[DAFidx]) iter++;
                                _order[DAFidx].insert(iter,f);
                        }
                }
                if(_order[DAFidx].size()<_n && parent::output_normal()) {std::ostringstream sos; sos << "WARNING: DAF"<<DAFidx<<" rank could not be computed for all features due to insufficient number of evaluated probe subsets." << std::endl << std::flush; syncout::print(os,sos);}
        }
        if(parent::output_normal()) {
                RETURNTYPE minfreq=_stats[0].freq_is; if(_stats[0].freq_isnot<_stats[0].freq_is) minfreq=_stats[0].freq_isnot;
                for(DIMTYPE f=1;f<_n;f++) {
                        if(_stats[f].freq_isnot<minfreq) minfreq=_stats[f].freq_isnot;
                        if(_stats[f].freq_is<minfreq) minfreq=_stats[f].freq_is;
                }
                {
                        std::ostringstream sos; 
                        sos << "Lowest feature evaluation frequency over all probes is "<<minfreq<<"." << std::endl;
                        if(minfreq<100) sos << "WARNING: To ensure reliable DAF estimates the minimum feature evaluation frequency should be at least in the order of hundreds. More probes may be needed." << std::endl;
                        sos << std::endl << std::flush;
                        syncout::print(os,sos);
                }
        }
 
        return true;
}
 
template<class RETURNTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET>
bool Result_Tracker_Feature_Stats<RETURNTYPE, IDXTYPE, DIMTYPE, SUBSET>::print_stats(std::ostream& os)
{
        std::ostringstream sos;
        if(_n<1) {sos << "Result_Tracker_Feature_Stats:: Nothing to print out." << std::endl; syncout::print(os,sos); return false;}
        {
                sos << "Probe subset stats for each (valid) cardinality:"<< std::endl;
                for(DIMTYPE ss=1;ss<=_n;ss++) 
                        if(_dstat[ss].valid) sos << " d: "<< ss << "  mean=" << _dstat[ss].mean << ", stddev=" << _dstat[ss].stdev << ", freq=" << _dstat[ss].freq << std::endl;
                sos << std::endl << "Feature stats over all probe subsets:"<< std::endl;
                for(DIMTYPE f=0;f<_n;f++) {
                        sos << " F:"<<f<< " daf0 = "; if(_stats[f].daf_valid[0]) sos << _stats[f].daf[0]; else sos << "n/a"; 
                        sos << ", daf1 = "; if(_stats[f].daf_valid[1]) sos << _stats[f].daf[1]; else sos << "n/a";
                        sos << ", daf2 = "; if(_stats[f].daf_valid[2]) sos << _stats[f].daf[2]; else sos << "n/a";
                        sos<<",  freq "<<_stats[f].freq_is<<"|"<<_stats[f].freq_isnot<<",  mean ";
                        if(_stats[f].freq_is>0) sos<<_stats[f].mean_is; else sos<< "n/a";
                        sos <<"|";
                        if(_stats[f].freq_isnot>0) sos<<_stats[f].mean_isnot; else sos<< "n/a";
                        sos<<", stdev ";
                        if(_stats[f].freq_is>0) sos<<_stats[f].stdev_is; else sos<<"n/a";
                        sos<<"|";
                        if(_stats[f].freq_isnot>0) sos<<_stats[f].stdev_isnot; else sos << "n/a";
                        sos<<std::endl;
                }
                sos<<std::endl;
        }
        syncout::print(os,sos);
        return true;
}
 
template<class RETURNTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET>
bool Result_Tracker_Feature_Stats<RETURNTYPE, IDXTYPE, DIMTYPE, SUBSET>::getFirstDAF(RETURNTYPE &value, DIMTYPE &feature, const unsigned int DAFidx) const
{
        assert(DAFidx<3);
#ifdef FST_THREADED
        std::unique_lock<std::mutex> lock(_itersize_mutex);
#endif
        _itersize[DAFidx]=0;
        if(!_order[DAFidx].empty())     {
                feature=_order[DAFidx][0];
                assert(feature<_n);
                value=_stats[feature].daf[DAFidx];
                _itersize[DAFidx]=1;
                return true;
        } else return false;
}
 
template<class RETURNTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET>
bool Result_Tracker_Feature_Stats<RETURNTYPE, IDXTYPE, DIMTYPE, SUBSET>::getNextDAF(RETURNTYPE &value, DIMTYPE &feature, const unsigned int DAFidx) const
{
        assert(DAFidx<3);
#ifdef FST_THREADED
        std::unique_lock<std::mutex> lock(_itersize_mutex);
#endif
        if(0<_itersize[DAFidx] && _itersize[DAFidx]<_order[DAFidx].size()) {
                feature=_order[DAFidx][_itersize[DAFidx]];
                assert(feature<_n);
                value=_stats[feature].daf[DAFidx];
                ++_itersize[DAFidx];
                return true;
        }
        return false;
}
 
template<class RETURNTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET>
bool Result_Tracker_Feature_Stats<RETURNTYPE, IDXTYPE, DIMTYPE, SUBSET>::compute_pair_stats(const DIMTYPE top_k, std::ostream& os)
{
        _pair_stats.clear(); _pair_features.clear(); _pair_order.clear(); _pair_itersize=0;
        if(parent::results.empty() || _n<1 || _order[0].empty()) return false; // compute_stats() must have run first
        if(top_k<2) return false;
        const std::size_t K = std::min<std::size_t>(top_k, _order[0].size());
        if(K<2) return false;
 
        // the top-K DAF0-ranked features define the pair set; local index <-> global feature
        _pair_features.resize(K);
        for(std::size_t k=0;k<K;k++) _pair_features[k]=_order[0][k];
 
        // accumulation pass: only the both-present cell needs per-probe work; the
        // one-only and both-absent cells derive by inclusion-exclusion below
        const std::size_t P=K*(K-1)/2;
        std::vector<RETURNTYPE> both_freq(P,0), both_sum(P,0), both_sumsq(P,0);
        std::vector<RETURNTYPE> f_freq(K,0), f_sum(K,0), f_sumsq(K,0);
        RETURNTYPE t_freq=0, t_sum=0, t_sumsq=0;
        std::vector<std::size_t> present; present.reserve(K);
        for(typename parent::CONSTRESULTSITER iter=parent::results.begin();iter!=parent::results.end();iter++)
        {
                assert(iter->sub->get_n()==_n);
                const RETURNTYPE v=iter->value;
                t_freq+=1; t_sum+=v; t_sumsq+=v*v;
                present.clear();
                for(std::size_t k=0;k<K;k++) if(iter->sub->selected_raw(_pair_features[k])) present.push_back(k);
                for(std::size_t i=0;i<present.size();i++) {
                        const std::size_t a=present[i];
                        f_freq[a]+=1; f_sum[a]+=v; f_sumsq[a]+=v*v;
                        for(std::size_t j=i+1;j<present.size();j++) {
                                const std::size_t ij=pair_index(a,present[j],K);
                                both_freq[ij]+=1; both_sum[ij]+=v; both_sumsq[ij]+=v*v;
                        }
                }
        }
 
        // derive the 2x2 cells and the pair indicators
        _pair_stats.resize(P);
        std::size_t valid_count=0;
        for(std::size_t a=0;a<K;a++) for(std::size_t b=a+1;b<K;b++) {
                const std::size_t ij=pair_index(a,b,K);
                PairStat &ps=_pair_stats[ij];
                ps.f1=_pair_features[a]; ps.f2=_pair_features[b];
                const RETURNTYPE fb=both_freq[ij], sb=both_sum[ij], qb=both_sumsq[ij];
                const RETURNTYPE f1o=f_freq[a]-fb, s1o=f_sum[a]-sb;       // f1 present, f2 absent
                const RETURNTYPE f2o=f_freq[b]-fb, s2o=f_sum[b]-sb;       // f2 present, f1 absent
                const RETURNTYPE fn=t_freq-f_freq[a]-f_freq[b]+fb;        // neither present
                const RETURNTYPE sn=t_sum-f_sum[a]-f_sum[b]+sb;
                const RETURNTYPE qn=t_sumsq-f_sumsq[a]-f_sumsq[b]+qb;
                ps.freq_both=fb; ps.freq_neither=fn;
                ps.daf0=ps.daf1=ps.synergy=0;
                ps.valid=(fb>0 && fn>0 && f1o>0 && f2o>0);
                ps.daf1_valid=false;
                if(ps.valid) {
                        const RETURNTYPE mb=sb/fb, mn=sn/fn, m1o=s1o/f1o, m2o=s2o/f2o;
                        ps.daf0=mb-mn;
                        ps.synergy=mb+mn-m1o-m2o;
                        RETURNTYPE vb=qb/fb-mb*mb; if(vb<0) vb=0;             // guard FP negatives
                        RETURNTYPE vn=qn/fn-mn*mn; if(vn<0) vn=0;
                        const RETURNTYPE denom=(fb/(fb+fn))*sqrt(vb) + (fn/(fb+fn))*sqrt(vn);
                        if(denom>0) {ps.daf1=ps.daf0/denom; ps.daf1_valid=true;}
                        ++valid_count;
                }
        }
 
        // order valid pairs by decreasing pair-DAF0 (deterministic tie-break on feature indices)
        _pair_order.reserve(valid_count);
        for(std::size_t ij=0;ij<P;ij++) if(_pair_stats[ij].valid) _pair_order.push_back(ij);
        std::sort(_pair_order.begin(),_pair_order.end(),[this](const std::size_t x, const std::size_t y)->bool {
                if(_pair_stats[x].daf0!=_pair_stats[y].daf0) return _pair_stats[x].daf0>_pair_stats[y].daf0;
                if(_pair_stats[x].f1!=_pair_stats[y].f1) return _pair_stats[x].f1<_pair_stats[y].f1;
                return _pair_stats[x].f2<_pair_stats[y].f2;
        });
 
        if(parent::output_normal()) {
                std::ostringstream sos;
                sos << "compute_pair_stats(): pairwise stats over the top " << K << " DAF0-ranked features (" << P << " pairs, " << valid_count << " valid)." << std::endl;
                if(valid_count<P) sos << "WARNING: " << (P-valid_count) << " pairs have an empty presence cell (insufficient probes); they carry no pair-DAF value." << std::endl;
                sos << std::flush;
                syncout::print(os,sos);
        }
        return valid_count>0;
}
 
template<class RETURNTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET>
bool Result_Tracker_Feature_Stats<RETURNTYPE, IDXTYPE, DIMTYPE, SUBSET>::print_pair_stats(std::ostream& os, const std::size_t top_pairs)
{
        std::ostringstream sos;
        if(_pair_order.empty()) {sos << "Result_Tracker_Feature_Stats:: No pairwise stats to print out." << std::endl; syncout::print(os,sos); return false;}
        const std::size_t count = (top_pairs==0 || top_pairs>_pair_order.size()) ? _pair_order.size() : top_pairs;
        sos << "Pairwise DAF stats over the top " << _pair_features.size() << " DAF0-ranked features (showing " << count << " of " << _pair_order.size() << " valid pairs, best pair-DAF0 first):" << std::endl;
        for(std::size_t i=0;i<count;i++) {
                const PairStat &ps=_pair_stats[_pair_order[i]];
                sos << " P:("<<ps.f1<<","<<ps.f2<<") pairdaf0 = " << ps.daf0;
                sos << ", pairdaf1 = "; if(ps.daf1_valid) sos << ps.daf1; else sos << "n/a";
                sos << ", synergy = " << ps.synergy;
                sos << ",  freq both|neither " << ps.freq_both << "|" << ps.freq_neither << std::endl;
        }
        sos << std::endl;
        syncout::print(os,sos);
        return true;
}
 
template<class RETURNTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET>
bool Result_Tracker_Feature_Stats<RETURNTYPE, IDXTYPE, DIMTYPE, SUBSET>::getFirstPairDAF(RETURNTYPE &daf0, RETURNTYPE &daf1, bool &daf1_valid, RETURNTYPE &synergy, DIMTYPE &f1, DIMTYPE &f2) const
{
#ifdef FST_THREADED
        std::unique_lock<std::mutex> lock(_itersize_mutex);
#endif
        _pair_itersize=0;
        if(!_pair_order.empty()) {
                const PairStat &ps=_pair_stats[_pair_order[0]];
                daf0=ps.daf0; daf1=ps.daf1; daf1_valid=ps.daf1_valid; synergy=ps.synergy; f1=ps.f1; f2=ps.f2;
                _pair_itersize=1;
                return true;
        } else return false;
}
 
template<class RETURNTYPE, typename IDXTYPE, typename DIMTYPE, class SUBSET>
bool Result_Tracker_Feature_Stats<RETURNTYPE, IDXTYPE, DIMTYPE, SUBSET>::getNextPairDAF(RETURNTYPE &daf0, RETURNTYPE &daf1, bool &daf1_valid, RETURNTYPE &synergy, DIMTYPE &f1, DIMTYPE &f2) const
{
#ifdef FST_THREADED
        std::unique_lock<std::mutex> lock(_itersize_mutex);
#endif
        if(0<_pair_itersize && _pair_itersize<_pair_order.size()) {
                const PairStat &ps=_pair_stats[_pair_order[_pair_itersize]];
                daf0=ps.daf0; daf1=ps.daf1; daf1_valid=ps.daf1_valid; synergy=ps.synergy; f1=ps.f1; f2=ps.f2;
                ++_pair_itersize;
                return true;
        }
        return false;
}
 
 
} // namespace
#endif // FSTRESULTTRACKERFEATURESTATS_H ///:~