geese-bones.hpp

Go to the documentation of this file.

#ifndef GEESE_BONES_HPP
#define GEESE_BONES_HPP 1
 
// #include <vector>
// #include <algorithm>
// #include <random>
// #include <stdexcept>
 
template<typename Ta, typename Tb>
inline std::vector< Ta > vector_caster(const std::vector< Tb > & x) {
 
    std::vector< Ta > ans;
    ans.reserve(x.size());
 
    for (auto i = x.begin(); i != x.end(); ++i)
        ans.push_back(static_cast< Ta >(*i));
 
    return ans;
    
}
 
#define INITIALIZED() if (!this->initialized) \
    throw std::logic_error("The model has not been initialized yet.");
 
// The same need to be locked
RULE_FUNCTION(rule_empty_free) {
 
    return Array(i, j) == 9u;
    
}
 
 
 
// Hasher
 
inline std::vector< double > keygen_full(
    const PhyloArray & array,
    const PhyloCounterData * d
    ) {
 
    // Baseline data: nrows and columns
    std::vector< double > dat = {
        static_cast<double>(array.nrow()) * 100000 +
         static_cast<double>(array.ncol()),
         // state of the parent
         1000000.0, 
         // type of the parent
         array.D_ptr()->duplication ? 1.0 : 0.0, 
         // Annotations with zeros
         0.0 
    };
 
    // State of the parent
    double pow10 = 1.0;
    for (bool i : array.D_ptr()->states) {
        dat[1u] += (i ? 1.0 : 0.0) * pow10;
        pow10 *= 10.0;
    }
 
    // // Annotations with zeros
    // pow10 = 1.0;
    // for (const auto & cell: array.get_data()) {
    //     dat[3u] += (cell == 9u ? 2.0 : static_cast<double>(cell)) * pow10;
    //     pow10 *= 10.0;
    // }
 
    return dat;
    
}
 
inline bool vec_diff(
    const std::vector< size_t > & s,
    const std::vector< size_t > & a
) {
 
    for (size_t i = 0u; i < a.size(); ++i)
        if ((a[i] != 9u) && (a[i] != s[i]))
            return true;
 
    return false;
}
 
class Flock;
 
class Geese {
    friend Flock;
private:
 
    std::mt19937 *                     rengine = nullptr;
    PhyloModel * model   = nullptr;
    std::vector< std::vector< bool > > states;
    size_t n_zeros       = 0u; 
    size_t n_ones        = 0u; 
    size_t n_dupl_events = 0u; 
    size_t n_spec_events = 0u; 
 
public:
 
    // Data
    size_t                       nfunctions;
    std::map< size_t, Node >     nodes;
    
    barry::MapVec_type< size_t > map_to_state_id;
    std::vector< std::vector< std::vector< size_t > > > pset_loc;    
 
    // Tree-traversal sequence
    std::vector< size_t > sequence;
    std::vector< size_t > reduced_sequence;  
 
    // Admin-related objects
    bool initialized     = false;
    bool delete_rengine  = false;
    bool delete_support  = false;
 
    // Information about the type of event
    
    /***
     * @name Information about the type of event
     * @details
     * The type of event is stored in the `etype` member. The possible values
     * are `etype_default`, `etype_speciation`, `etype_duplication`, and
     * `etype_either`.
     * 
    */
    static const size_t etype_default     = 1ul;
    static const size_t etype_speciation  = 0ul;
    static const size_t etype_duplication = 1ul;
    static const size_t etype_either      = 2ul;
 
    Geese();
 
    Geese(
        std::vector< std::vector<size_t> > & annotations,
        std::vector< size_t > &              geneid,
        std::vector< int > &                 parent,
        std::vector< bool > &                duplication
        );
 
    // Copy constructor
    Geese(const Geese & model_, bool copy_data = true);
    
    // Constructor move
    Geese(Geese && x) noexcept;
 
    // Copy assignment
    Geese & operator=(const Geese & model_) = delete;
 
    // // Move assignment
    Geese & operator=(Geese && model_) noexcept = delete;
 
 
    ~Geese();
 
    void init(size_t bar_width = BARRY_PROGRESS_BAR_WIDTH);
 
    void inherit_support(const Geese & model_, bool delete_support_ = false);
 
    // Node * operator()(size_t & nodeid);
    void calc_sequence(Node * n = nullptr);
    void calc_reduced_sequence();
 
    double likelihood(
        const std::vector< double > & par,
        bool as_log = false,
        bool use_reduced_sequence = true,
        size_t ncores = 1u,
        bool no_update_pset_probs = false
        );
 
    double likelihood_exhaust(const std::vector< double > & par);
 
    std::vector< double > get_probabilities() const;
 
    void set_seed(const size_t & s);
    std::vector< std::vector< size_t > > simulate(
        const std::vector< double > & par
        );
 
    size_t nfuns() const noexcept;             
    size_t nnodes() const noexcept;            
    size_t nleafs() const noexcept;            
    size_t nterms() const;                     
    size_t support_size() const noexcept;      
    std::vector< size_t > nannotations() const noexcept;      
    std::vector< std::string > colnames() const;     
    size_t parse_polytomies(
        bool verb = true,
        std::vector< size_t > * dist = nullptr
        ) const noexcept;  
 
 
    std::vector< std::vector<double> > observed_counts();
    void print_observed_counts();
 
    void print() const;
    void print_nodes() const;
 
 
    std::vector< std::vector< double > > predict(
        const std::vector< double > & par,
        std::vector< std::vector< double > > * res_prob = nullptr,
        bool leave_one_out        = false,
        bool only_annotated       = false,
        bool use_reduced_sequence = true
        );
    
    std::vector< std::vector<double> > predict_backend(
        const std::vector< double > & par,
        bool use_reduced_sequence,
        const std::vector< size_t > & preorder
        );
 
    std::vector< std::vector< double > > predict_exhaust_backend(
        const std::vector< double > & par,
        const std::vector< size_t > & preorder
        );
 
    std::vector< std::vector< double > > predict_exhaust(
        const std::vector< double > & par
        );
 
    std::vector< std::vector< double > > predict_sim(
        const std::vector< double > & par,
        bool only_annotated       = false,
        size_t nsims        = 10000u
        );
 
    void init_node(Node & n);
    void update_annotations(
        size_t nodeid,
        std::vector< size_t > newann
    );
 
    std::mt19937 *  get_rengine();
    PhyloCounters * get_counters();
    PhyloModel *    get_model();
    PhyloSupport *  get_support_fun();
    
    std::vector< std::vector< bool > > get_states() const;  
    std::vector< size_t > get_annotated_nodes() const; 
    std::vector< size_t > get_annotations() const; 
 
};
 
#endif