Bayesable.h

Go to the documentation of this file.

#pragma once 

#include <utility>
#include <iomanip>
#include <signal.h>

#include "FleetArgs.h"
#include "Errors.h"
#include "Datum.h"
#include "IO.h"
#include "Statistics/FleetStatistics.h"
#include "Miscellaneous.h"
#include "Statistics/StreamingSum.h"

// if we define this, then we won't use breakouts
//#define NO_BREAKOUT 1

extern std::atomic<bool> CTRL_C;

//cache::lru_cache<size_t, std::tuple<double,double,double>> posterior_cache(10);

template<typename _datum_t, typename _data_t=std::vector<_datum_t>>
class Bayesable {
public:
    
    // We'll define a vector of pairs of inputs and outputs
    // this may be used externally to define what data is
    typedef _datum_t datum_t;
    typedef _data_t  data_t; 

    // log values for all
    double prior; 
    double likelihood;
    double posterior; // Posterior always stores at temperature 1
    
    // some born info/statistics
    uintmax_t born; // what count were you born at?
    size_t born_chain_idx; // if we came from a ChainPool, what chain index were we born on?
    
    Bayesable() : prior(NaN), likelihood(NaN), posterior(NaN), born(++FleetStatistics::hypothesis_births), born_chain_idx(0) {  }
    
    // Stuff for subclasses to implement: 
    virtual size_t hash() const=0;
    virtual std::string string(std::string prefix="") const = 0; 
    virtual double compute_prior() = 0;     
    
    virtual double compute_single_likelihood(const datum_t& datum)  {
        throw NotImplementedError();
    }
    
    virtual void clear_bayes() {
        // necessary for inserting into big collections not by prior
        prior = 0.0;
        likelihood = 0.0;
        posterior = 0.0;
    }
    
    virtual double compute_likelihood(const data_t& data, const double breakout=-infinity) {
    
        // include this in case a subclass overrides to make it non-iterable -- then it must define its own likelihood
        if constexpr (is_iterable_v<data_t>) { 
            // defaultly a sum over datums in data (e.g. assuming independence)
            likelihood = 0.0;
            for(const auto& d : data) {
                
                
                auto sll = compute_single_likelihood(d);
                
                likelihood += sll;
                
                if(likelihood == -infinity or std::isnan(likelihood)) break; // no need to continue
                
                // This is a breakout in case our ll is too low
                if(FleetArgs::LIKELIHOOD_BREAKOUT){
                    assert((sll <= 0 or breakout==-infinity) && "*** Cannot use breakout if likelihoods are positive");
                    if(likelihood < breakout) {
                        return likelihood = -infinity; // should not matter what value, but let's make it -infinity
                    }
                }
                
                // add a break here 
                if(CTRL_C) break; 
            }
            
            return likelihood;      
        }
        else {
            // should never execute this because it must be defined
            throw NotImplementedError("*** If you use a non-iterable data_t, then you must define compute_likelihood on your own."); 
        }
    }
    
    
    virtual double compute_tempered_likelihood(const data_t& data, int ladder_rank, const double breakout=-infinity) {
        return compute_likelihood(data,breakout); 
    }

    virtual double compute_posterior(const data_t& data, const std::pair<double,double> breakoutpair=std::make_pair(-infinity,1.0)) {
        
        ++FleetStatistics::posterior_calls; // just keep track of how many calls 
        
        // Always compute a prior
        prior = compute_prior();
        
        // if the prior is -inf, then we don't have to compute likelihood
        if(prior == -infinity) {
            likelihood = NaN;
            posterior = -infinity;
        }
        else {      
            // NOTE: here we *subtract* off prior so the breakout in compute_likelihood is ONLY for the likelihood
            auto [b,t] = breakoutpair; // see MCMCChain.h to understand this and hte next line
            likelihood = compute_likelihood(data, (b-prior)*t);
            posterior = prior + likelihood; 
        }
        
        return posterior;
    }
    
    virtual double at_temperature(double t) const {
        // temperature here applies to the likelihood only
        return prior + likelihood/t;
    }
    
    // Defaultly we search by posterior 
    auto operator<=>(const Bayesable<datum_t,data_t>& other) const {
        if(posterior != other.posterior) {
            return fp_ordering(posterior, other.posterior); 
        }
        else if(prior != other.prior) {
            return fp_ordering(prior,other.prior); 
        }
        else {
            // we could be equal, but we only want to say that if we are. 
            // otherwise this is decided by the hash function 
            // so we must ensure that the hash function is only equal for equal values
            return this->hash() <=> other.hash(); 
        }
    }
    
    virtual void show(std::string prefix="") {
        if(prefix != "") {
            print(prefix, this->posterior, this->prior, this->likelihood, QQ(this->string())); 
        }
        else {
            print(this->posterior, this->prior, this->likelihood, QQ(this->string())); 
        }
    }
};



template<typename _datum_t, typename _data_t>
std::ostream& operator<<(std::ostream& o, Bayesable<_datum_t,_data_t>& x) {
    o << x.string();
    return o;
}


// little helper functions 
template<typename HYP>
std::function get_posterior = [](const HYP& h) {return h.posterior; };

template<typename HYP>
std::function get_prior = [](const HYP& h) {return h.prior; };

template<typename HYP>
std::function get_likelihood = [](const HYP& h) {return h.likelihood; };


// Interface for std::fmt
//#include "fmt/format.h"
//
//template<typename _datum_t, typename _data_t>
//struct fmt::formatter<Bayesable<_datum_t,_data_t>&> {
//  constexpr auto parse(format_parse_context& ctx) { return ctx.begin(); }
//
//  template <typename FormatContext>
//  auto format(const Bayesable<_datum_t,_data_t>& h, FormatContext& ctx) {
//    return format_to(ctx.out(), "{}", h.string());
//  }
//};