settings.h

/* Spell-QTL  Software suite for the QTL analysis of modern datasets.
 * Copyright (C) 2016,2017  Damien Leroux <damien.leroux@inra.fr>, Sylvain Jasson <sylvain.jasson@inra.fr>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef _SPEL_SETTINGS_H_
#define _SPEL_SETTINGS_H_

#include "input/read_mark.h"
#include "input/read_map.h"
#include "input/read_trait.h"
#include "input/design.h"
#include <exception>
/*#include <Eigen/Core>*/
#include "eigen.h"
#include "proba_config.h"
#include "data/genoprob_computer2.h"
#include "data/qtl_chromosome.h"
#include "input/read_mark.h"
#include "input/read_trait.h"

#include "bayes/output.h"

#include "input/observations.h"
#include <thread>

//#include "../3rd-party/ThreadPool/ThreadPool.h"

#include "basic_file_checks.h"

#include "linear_combination.h"


#include "task_pool.h"

enum ComputationType { NoTest=0, FTest=1, FTestLOD=2, R2=4, Chi2=8, Chi2LOD=16, Mahalanobis=32 };
enum ComputationResults { Test=0, RSS=1, Residuals=2, Coefficients=4, Rank=8 };

constexpr ComputationType operator | (ComputationType a, ComputationType b) { return  (ComputationType) (((int)a) | ((int)b)); }
constexpr ComputationResults operator | (ComputationResults a, ComputationResults b) { return  (ComputationResults) (((int)a) | ((int)b)); }



struct pleiotropy_descr {
    std::string name;
    double tolerance;
    std::vector<std::string> traits;
};


struct marker_observation_spec {
    ObservationDomain domain;
    size_t domain_size;
    std::vector<std::pair<char, std::vector<std::pair<char, char>>>> scores;
    /*std::map<char, VectorXd> compile() const;*/
    const std::vector<std::pair<char, char>>& score(char obs) const
    {
        static std::vector<std::pair<char, char>> empty;
        auto it = std::find_if(scores.begin(), scores.end(), [=] (const std::pair<char, std::vector<std::pair<char, char>>>& e) { return e.first == obs; });
        return it != scores.end() ? it->second : empty;
    }
};


struct population_marker_observation {
    std::string filename;
    std::string format_name;
    size_t process_size;
    marker_data observations;

    std::map<char, VectorXd> observation_vectors;

    MatrixXd raw_observations(const std::vector<char>& observations) const
    {
        MatrixXd ret = MatrixXd::Ones(process_size, observations.size());
        std::string debug_obs(observations.begin(), observations.end());
        /*MSG_DEBUG("raw_observations[start] " << name << " obs=" << debug_obs);*/
        /*MSG_DEBUG(main_process().column_labels);*/
        /*for (auto& kv: obs_vectors) {*/
        /*MSG_DEBUG(" * " << kv.first << "   " << kv.second.transpose());*/
        /*}*/
        for (size_t i = 0; i < observations.size(); ++i) {
            std::map<char, VectorXd>::const_iterator ovi = observation_vectors.find(observations[i]);
            if (ovi != observation_vectors.end()) {
                ret.col(i) = ovi->second;
            } else {
                /* error! */
                MSG_ERROR("OBSERVATION NOT FOUND! obs #" << i << " chr(" << ((int)observations[i]) << ')', "FIX IT");
            }
        }
        /*MSG_DEBUG("raw_observations " << name << " obs=" << debug_obs << std::endl << ret);*/
        /*MSG_DEBUG("raw_observations obs=" << debug_obs << std::endl << ret);*/
        return ret;
    }

    operator marker_data& () { return observations; }
    operator std::string () { return filename; }
};


struct LD_matrices {
    std::vector<double> loci;
    std::vector<MatrixXd> ld;
};


struct bad_settings_exception : public std::exception {
    std::string m;
    bad_settings_exception(const std::string& msgs)
        : m(MESSAGE("One or more errors were found." << std::endl << msgs))
    {}
    const char* what() const throw() { return m.c_str(); }
};


struct trait_metadata_type {
    std::vector<std::string> dim_names;
    MatrixXd P;
};



struct settings_t {
    std::string name;
    std::string notes;

    std::string map_filename;
    std::vector<chromosome> map;
    /*format_specification_t* marker_observation_specs;*/
    /*std::multimap<std::string, qtl_pop_type> populations;*/
    std::vector<std::shared_ptr<qtl_pop_type>> populations;
    std::map<std::string, std::vector<std::shared_ptr<qtl_pop_type>>> pops_by_trait;

    std::vector<qtl_chromosome> working_set;

    double step;
    int parallel;
    std::string work_directory;
    bool connected;
    bool epistasis;
    bool pleiotropy;
    double pleiotropy_tolerance;
    std::vector<std::string> chromosome_selection;
    std::string cofactor_marker_selection_filename;
    double cofactor_marker_selection_distance;

    std::string epistasis_qtl_selection_filename;

    std::map<std::string, double> qtl_thresholds;
    int n_permutations;
    double qtl_threshold_quantile;
    std::map<std::string, double> cofactor_thresholds;

    std::string skeleton_mode;
    std::vector<std::string> skeleton_markers;
    double skeleton_interval;

    std::string cofactor_algorithm;
    double cofactor_exclusion_window_size;

    std::string qtl_algorithm;

    std::string initial_selection;

#if 0
    enum class detection_method_t : char {
        Undef,
        iQTLm,
        CIM
    };

    bool P_required;
    bool T_required;
    bool E_required;
    bool D_required;
    detection_method_t detection_method;
    double detection_window;
#endif

//	ThreadPool* pool;
	static std::thread::id main_thread;

	double tolerance;
    std::map<const chromosome*, std::vector<double>> estimation_loci;

    /* FIXME LD data is global to the pedigree (to the ancestors), not local to a population */
    std::map<std::string, std::vector<std::string>> ld_parents;
    std::map<std::string, std::map<const chromosome*, LD_matrices>> ld_data;
    std::map<std::string, std::map<const chromosome*, std::map<double, int>>> parent_count_per_pop_per_chr;

    LV_database LV;

    std::map<std::string, std::vector<std::vector<std::shared_ptr<qtl_pop_type>>>> linked_pops;

    std::unordered_map<std::thread::id, std::vector<std::string>> thread_stacks;

    size_t max_order;
    bool cross_indicator_can_interact;

    std::vector<std::string> npoint_gen;

    bool output_npoint;

    std::unique_ptr<pop_data_type> pop_data;
    std::vector<std::pair<std::string, std::string>> datasets;
    std::vector<std::pair<std::string, std::string>> covariables;
    std::vector<pleiotropy_descr> pleiotropic_traits_descr;

    std::vector<std::string> with_covariables;
    std::vector<std::string> with_traits;
    std::vector<std::string> with_lg;

    std::vector<trait> covar_per_pop;

    double lod_support_inner, lod_support_outer;

    std::map<std::string, trait_metadata_type> trait_metadata;

    double Rjoin;
    double Rskip;
    
    std::string report_format_string;
    ComputationResults additional_table_output;

    settings_t()
        : notes()
        , map_filename(), map()
        /*, populations()*/
        , pops_by_trait()
        , working_set()
        , step(1)
        , parallel(0)
        , work_directory(".")
        , connected(false)
        , epistasis(false)
        , pleiotropy(false)
        , pleiotropy_tolerance(.001)
        , chromosome_selection()
        , cofactor_marker_selection_filename()
        , cofactor_marker_selection_distance(10)
        , epistasis_qtl_selection_filename()
        , qtl_thresholds()
        , n_permutations(10000)
        , qtl_threshold_quantile(0.05)
        , cofactor_thresholds()
        , skeleton_mode("auto")
        , skeleton_markers()
        , skeleton_interval(20.)
        , cofactor_algorithm("forward")
        , cofactor_exclusion_window_size(30.)
        , qtl_algorithm("iQTLm")
        /*, detection_method(detection_method_t::Undef)*/
        /*, detection_window(10.)*/
//		, pool(0)
		, tolerance(1.e-10)
        , estimation_loci()
        , ld_parents()
        , ld_data()
        , parent_count_per_pop_per_chr()
        , LV()
        , linked_pops()
        , thread_stacks()
        , max_order(1)
        , cross_indicator_can_interact(false)
        , npoint_gen()
        , output_npoint(false)
        , pop_data()
        , datasets()
        , covariables()
        , pleiotropic_traits_descr()
        , with_covariables()
        , with_traits()
        , with_lg()
        , covar_per_pop()
        , lod_support_inner(1.)
        , lod_support_outer(lod_support_inner * 2)
        , Rjoin(0.)
        , Rskip(1.)
        , report_format_string("text")
//         , report_format_string("excel")
//         , additional_table_output((ComputationResults) 0)
        , additional_table_output(ComputationResults::RSS | ComputationResults::Rank)
    {}

    std::vector<std::pair<const chromosome*, double>>
        selected_qtl() const
        {
            std::vector<std::pair<const chromosome*, double>> ret;
            size_t n = 0;
            for (auto& qc: working_set) {
                n += qc.qtl.size();
            }
            ret.reserve(n);
            for (auto& qc: working_set) {
                for (auto& q: qc.qtl) {
                    ret.emplace_back(qc.chr, q);
                }
            }
            return ret;
        }

    std::vector<double>
        selected_qtl(const std::string& name) const
        {
            for (auto& qc: working_set) {
                if (qc.chr->name == name) {
                    return qc.qtl;
                }
            }
            return {};
        }

    ~settings_t()
    {
//        if (pool) { delete pool; }
        /*if (marker_observation_specs) { delete marker_observation_specs; }*/
    }

    const chromosome*
        find_chromosome(const std::string& name) const
        {
            for (auto& c: map) {
                if (c.name == name) {
                    return &c;
                }
            }
            return NULL;
        }

    static settings_t* from_args(int argc, const char** argv);

    bool sanity_check() const;

    void finalize();

//	template <class F, class... Args>
//    auto enqueue(int _Sync, F&& f, Args&&... args)
//        -> std::future<typename std::result_of<F(Args...)>::type>
//		{
//            return TaskPool::enqueue(std::forward<F>(f),
//                                     std::forward<Args>(args)...);
//
////			if (_Sync || parallel < 2 || std::this_thread::get_id() != main_thread) {
////                /*MSG_DEBUG(std::this_thread::get_id()  << " RUN SYNC");*/
////				return std::async(std::launch::deferred,
////								  std::forward<F>(f),
////								  std::forward<Args>(args)...);
////			}
////			return pool->enqueue(std::forward<F>(f),
////								 std::forward<Args>(args)...);
//		}

    std::string set_title(const std::string& t)
    {
        MSG_INFO(GREEN << "Current task: " << t);
        return {};
//        std::string ret;
//        if (pool) {
//            ret = pool->get_title();
//        }
//        if (pool && msg_handler_t::color()) {
//            pool->set_title(t);
//        }
//        return ret;
    }

    void set_parallel(const std::string& s)
    {
        if (s == "auto") {
            parallel = std::thread::hardware_concurrency();
        } else {
            std::stringstream ss(s);
            ss >> parallel;
        }
        TaskPool::set_slot_count(parallel);
    }

    double get_threshold(const std::string& trait, size_t y_block_size);
};


struct format_specification_t {
    std::string filename;
    ObservationDomain domain;
    std::map<std::string, marker_observation_spec> map;

    void operator () (std::map<std::string, marker_observation_spec>& settings)
    {
        settings.insert(map.begin(), map.end());
    }
};


std::ostream& operator << (std::ostream&, const settings_t&);

namespace read_data {
    settings_t* read_settings(std::istream& is);
}


extern settings_t* active_settings;


struct context_key_struc;
typedef std::shared_ptr<context_key_struc> context_key;

struct context_key_struc {
    const qtl_pop_type* pop;
    const geno_matrix* gen;
    const chromosome* chr;
    std::vector<double> loci;
    std::unordered_map<double, size_t> locus_indices;

    void _init_locus_indices()
    {
        size_t sz = loci.size();
        for (size_t i = 0; i < sz; ++i) {
            locus_indices[loci[i]] = i;
        }
    }

    context_key_struc(const qtl_pop_type* p, const chromosome* c, const std::vector<double>& l)
        : pop(p)
        , gen(p->gen.get())
        , chr(c)
        , loci(l)
    {
        _init_locus_indices();
    }

    context_key_struc(const geno_matrix* g)
        : pop(NULL)
        , gen(g)
        , chr(NULL)
        , loci()
    {
    }

    context_key_struc(const qtl_pop_type* p, const geno_matrix* g, const chromosome* c)
        : pop(p)
        , gen(g)
        , chr(c)
        , loci(compute_steps(chr->condensed.marker_locus, active_settings->step))
    {
        _init_locus_indices();
    }

#if 1
    context_key_struc(const qtl_pop_type* p, const chromosome* c)
        : pop(p)
        , gen(p->gen.get())
        , chr(c)
        , loci(compute_steps(chr->condensed.marker_locus, active_settings->step))
    {
        _init_locus_indices();
    }
#endif

    bool operator == (const context_key_struc& ck) const
    {
        return pop == ck.pop && chr == ck.chr;
    }

    bool operator != (const context_key_struc& ck) const { return !(*this == ck); }

    friend
        bool operator == (const context_key& k1, const context_key& k2) { return *k1 == *k2; }

    friend
        bool operator != (const context_key& k1, const context_key& k2) { return *k1 != *k2; }

    friend
        std::ostream&
        operator << (std::ostream& os, const context_key& ck)
        {
            return os << "<context pop="
                      << (ck->pop ? ck->pop->name : "nil")
                      << " gen=" << (ck->gen ? ck->gen->name : "nil")
                      << " chr=" << (ck->chr ? ck->chr->name : "nil")
                      << ", " << ck->loci.size() << " loci>";
        }
};

namespace std {
    template <>
        struct hash<context_key> {
            size_t operator () (const context_key& ck) const {
                hash<string> h;
                size_t hn = 0;
                size_t hc = 0;
                size_t hs = 0;
                size_t hv = 0;
                if (ck->pop) {
                    hn = h(ck->pop->name);
                    hs = hash<size_t>()(ck->pop->size());
                    hv = hash<std::vector<std::string>>()(ck->pop->covariables.dim_names);
                }
                if (ck->chr) {
                    hc = h(ck->chr->name);
                }
                /*return (((hs * 0xdeadbe35) ^ hn) * 0xdeadbe35) ^ hc;*/
                return hv ^ hs ^ hn ^ hc;
            }
        };
}


int
get_n_parents(const context_key& ck, const locus_key& lk);


inline
std::vector<trait>
assemble_traits(const std::vector<trait>& src, const std::vector<std::string>& with_traits, const std::vector<pleiotropy_descr>& pleio) {
    std::vector<trait> ret(src.size() + pleio.size());
    MSG_DEBUG("with_traits.empty() ? " << with_traits.empty());
    auto it = with_traits.empty()
              ? std::copy(src.begin(), src.end(), ret.begin())
              : std::copy_if(src.begin(), src.end(), ret.begin(),
                             [&] (const trait& t) { return std::find(with_traits.begin(), with_traits.end(), t.name) != with_traits.end(); });
    ret.resize(it - ret.begin());
    for (const auto& d: pleio) {
        ret.emplace_back(src, d.traits, d.name, d.tolerance);
    }
    MSG_DEBUG("assembled traits " << ret);
    return ret;
}


inline
void
settings_t::finalize()
{
    {
        std::string log_file = MESSAGE(active_settings->work_directory << '/' << active_settings->name << ".spell-qtl.log");
        ensure_directories_exist(active_settings->work_directory);
        msg_handler_t::set_log_file(log_file);
        std::chrono::system_clock::time_point now = std::chrono::system_clock::now();
        std::time_t now_c = std::chrono::system_clock::to_time_t(now);
        (*msg_handler_t::instance().log_file)
                << "========================================================" << std::endl
                << "Log started on " <<  std::put_time(std::localtime(&now_c), "%F %T") << std::endl
                << "========================================================" << std::endl;
    }

    std::string pop_data_filename = work_directory + "/" + name + ".cache/" + name + ".spell-marker.data";
    if (!check_file(pop_data_filename, false, false, false)) {
        MSG_ERROR("The spell-marker data file was not found at " << pop_data_filename << ".", "Please run spell-pedigree and spell-marker prior to running spell-qtl.");
    }

    if (with_lg.size()) {
        std::vector<chromosome> with_map;
        for (const auto& chr: map) {
            if (std::find(with_lg.begin(), with_lg.end(), chr.name) != with_lg.end()) {
                with_map.push_back(chr);
            }
        }
        map.swap(with_map);
    }

    pop_data.reset(pop_data_type::load_from(pop_data_filename));
    pop_data->assemble_chromosomes(map);

    auto all_variables = with_traits + with_covariables;
    std::vector<std::string> in_pleio;
    for (const auto& p: pleiotropic_traits_descr) {
        in_pleio = in_pleio + p.traits;
    }
    all_variables = all_variables + in_pleio;
    auto removed_traits = with_traits % in_pleio;
    auto removed_covar = with_covariables % in_pleio;
    with_traits = with_traits - in_pleio;
    with_covariables = with_covariables - in_pleio;
    if (removed_covar.size()) {
        MSG_ERROR("Some covariables are already used in pleiotropic traits: " << removed_covar, "Specify covariable names that are different from the trait names.");
        return;
    }
    if (removed_traits.size()) {
        MSG_ERROR("Some single traits are already used in pleiotropic traits: " << removed_traits, "Specify single trait names that are different from the trait names used in pleiotropic traits.");
        return;
    }
    std::set<std::string> actual_single_traits;
    for (const auto& gt: datasets) {
        ifile ifs(gt.second);
        auto all_traits = read_data::read_trait(ifs, all_variables);
        auto observed_traits = assemble_traits(all_traits, with_traits, pleiotropic_traits_descr);
        for (const auto& t: observed_traits) {
            if (t.dim_names.size() == 1) {
                actual_single_traits.insert(t.name);
            }
        }
        pop_data->set_qtl_generation(gt.first);
        MSG_DEBUG("with_covariables " << with_covariables);
        MSG_DEBUG("Creating covariables " << " from traits " << with_covariables << " in list " << observed_traits);
        covar_per_pop.emplace_back(all_traits, with_covariables, "Covariables", 0);
        MSG_DEBUG("COVAR struct " << covar_per_pop.back());
        pop_data->extract_subpops(populations, gt.second, observed_traits, linked_pops, covar_per_pop.back());
    }
    for (const auto& n_p: linked_pops) {
        for (const auto& vp: n_p.second) {
            for (const auto& p: vp) {
                MSG_DEBUG("trait " << n_p.first << " pop " << p->name << ' ' << MATRIX_SIZE(p->gen->inf_mat) << " labels " << p->gen->labels);
            }
        }
    }
    for (const auto& pleio: pleiotropic_traits_descr) {
        auto& pops = linked_pops[pleio.name];
        MatrixXd xtx = MatrixXd::Zero(pleio.traits.size(), pleio.traits.size());
        for (const auto& vp: pops) {
            for (const auto& p: vp) {
                if (!p->observed_traits.front().raw.size()) {
                    MSG_ERROR("Invalid dataset", "");
                } else {
                    xtx += p->observed_traits.front().XtX();
                }
            }
        }
        auto P = trait::do_PCA(xtx, pleio.tolerance);
        for (const auto& vp: pops) {
            for (const auto& p: vp) {
                p->observed_traits.front().setP(P);
            }
        }
        trait_metadata[pleio.name].P = P;
        trait_metadata[pleio.name].dim_names = pleio.traits;
    }
    for (const auto& t: actual_single_traits) {
        trait_metadata[t].P = MatrixXd::Ones(1, 1);
        trait_metadata[t].dim_names = {t};
    }

    for (const auto& chr: map) {
        estimation_loci[&chr] = compute_steps(chr.condensed.marker_locus, step);
    }

    for (const auto& pptr: populations) {
        auto& parent_count_per_chr = parent_count_per_pop_per_chr[pptr->qtl_generation_name];
        auto ldi = ld_data.find(pptr->qtl_generation_name);
        if (ldi != ld_data.end()) {
            for (const auto& chr: map) {
                auto& parent_count = parent_count_per_chr[&chr];
                auto& ld = (*ldi).second[&chr];
                const auto& loci = estimation_loci[&chr];
                for (size_t i = 0; i < loci.size(); ++i) {
                    if (loci[i] != ld.loci[i]) {
                        MSG_ERROR("Discrepancy between LD locus (" << ld.loci[i] << ") and estimation locus (" << loci[i] << ") on chromosome " << chr.name, "Ensure the LD data is consistent with the map and step size you provided.");
                    }
                    parent_count[ld.loci[i]] = ld.ld[i].cols();
                }
            }
        } else {
            if (!pptr->gen) {
                MSG_ERROR("PROUT", "");
                continue;
            }
            context_key ck(new context_key_struc(pptr->gen.get()));
            locus_key lk;
            size_t n_par = get_n_parents(ck, lk);
            for (const auto& chr: map) {
                auto& parent_count = parent_count_per_chr[&chr];
                const auto& loci = estimation_loci[&chr];
                for (size_t i = 0; i < loci.size(); ++i) {
                    parent_count[loci[i]] = (int) n_par;
                }
            }
        }
    }
}


#endif