basic_data.cc

/* 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/>.
 */

//#include "input.h"
//#include <mutex>
#include "eigen.h"
#include "input.h"
//#include "data/geno_matrix.h"
#include "cache/md5.h"
#include "computations/base.h"
#include "cache2.h"
#include "computations/basic_data.h"
#include "io/output_impl.h"


collection<std::string>
all_traits()
{
    collection<std::string> ret;
    for (auto& kv: active_settings->linked_pops) {
        ret.emplace_back(kv.first);
    }
    return ret;
}


std::vector<std::pair<const chromosome*, double>>
selected_qtls()
{
    std::vector<std::pair<const chromosome*, double>> ret;
    for (auto& qc: active_settings->working_set) {
        for (auto& l: qc.qtl) {
            ret.push_back({qc.chr, l});
        }
    }
    return ret;
}


#if 0
population_value
population_by_name(const std::string& name)
{
    for (auto& kv: active_settings->populations) {
        if (kv.first == name) {
            return &kv.second;
        }
    }
    return NULL;
}
#endif

collection<population_value> all_populations()
{
	collection<population_value> ret;
	ret.reserve(active_settings->populations.size());
	for (auto& p: active_settings->populations) {
		/*ret.push_back({&kv.second});*/
        ret.push_back(p.get());
	}
	return ret;
}

std::vector<double>
selected_qtls_on_chromosome(qtl_chromosome_value chr)
{
    return chr->qtl;
}


#if 0
generation_value
generation_by_name(const std::string& name)
{
    /* FIXME */
    return NULL;
    /*return active_settings->design->generation[name];*/
}


generation_value
qtl_generation(const population_value& pop)
{
    /* FIXME */
    return NULL;
    /*return generation_by_name(pop->qtl_generation_name);*/
}
#else

generation_value
generation_by_name(const population_value& pop, const std::string& /*name*/)
{
    return pop->gen.get();
}

generation_value
qtl_generation(const population_value& pop)
{
    return pop->gen.get();
}

#endif

chromosome_value
chromosome_by_name(const std::string& name)
{
    for (auto& chr: active_settings->map) {
        if (chr.name == name) {
            return &chr;
        }
    }
    return NULL;
}


qtl_chromosome_value
qtl_chromosome_by_name(const std::string& name)
{
    for (auto& chr: active_settings->working_set) {
        if (chr.chr->name == name) {
            return &chr;
        }
    }
    return NULL;
}


collection<std::string>
all_observed_generation_names(population_value /*pop*/)
{
    collection<std::string> ret;
    /* FIXME */
    /*for (auto& kv: pop->marker_observation_filenames) {*/
        /*ret.emplace_back(kv.first);*/
    /*}*/
    return ret;
}


#if 0
/* DEPRECAITIDE */
std::vector<char>
marker_observations(population_value pop, chromosome_value chr, generation_value gen, int ind)
{
    return pop->get_observed_mark(gen->name)
            .observations
            .get_obs(chr->raw.marker_name.begin(), chr->raw.marker_name.end(), ind);
}


collection<double>
test_loci(chromosome_value chr)
{
    std::vector<double> steps = compute_steps(chr->marker_locus, active_settings->step);
    return collection<double>(steps.begin(), steps.end());
}
#endif

value<std::vector<double>>
test_loci(chromosome_value chr)
{
    return compute_steps(chr->condensed.marker_locus, active_settings->step);
}


#if 0
segment_computer_t
genoprob_computer(population_value pop, generation_value gen, qtl_chromosome_value qtl_chr)
{
    /*return gen->segment_computer(qtl_chr, 0, active_settings->step);*/
    segment_computer_t ret(*gen, qtl_chr, active_settings->step, 0.);
    return ret;
}


marker_observation_spec
marker_observation_specifications(generation_value gen)
{
    return active_settings->marker_observation_specs[gen->name];
}
#endif

/*observation_vectors_type*/
/*get_observation_vectors(generation_value gen)*/
/*{*/
        /*return gen->observation_vectors(marker_observation_specifications(gen));*/
/*}*/


/*pedigree_type*/
/*pedigree(population_value pop, generation_value gen, int ind)*/
/*{*/
    /*return pop->get_pedigree(gen, ind);*/
/*}*/


range<int>
individual_range(population_value pop)
{
    /*auto it = pop->families.find(pop->qtl_generation_name);*/
    /*if (it == pop->families.end()) {*/
        /*return {0, 0, 0};*/
    /*}*/
    return {0, (int) pop->size(), 1};
    /*auto it = pop->observed_mark.find(pop->qtl_generation_name);*/
    /*if (it == pop->observed_mark.end()) {*/
        /*return {0, 0, 0};*/
    /*}*/
    /*return {0, (int)it->second.observations.n_obs, 1};*/
}


MatrixXd
population_marker_obs(const context_key& ck, int ind)
{
    /*return {ck->pop, ck->chr, ck->gen, ind};*/
    /*return active_settings->LV(ck->chr->name, ck->pop->qtl_generation_name, ind);*/
    return ck->pop->get_LV(ck->chr->name, ind);
}

#if 0
MatrixXd
population_marker_obs(population_value pop, chromosome_value chr, int i, const pop_mgo_data* pmd)
{
    return {pop, chr, pmd->qtl_gen, i};
}

MatrixXd
population_marker_obs(population_value pop, chromosome_value chr, const pedigree_type& ped, const pop_mgo_data* pmd)
{
    /*generation_value qtl_gen = generation_by_name(pop->qtl_generation_name);*/
    /*collection<std::string> aogn = all_observed_generation_names(pop);*/
    /*collection<generation_value> aog = make_collection(generation_by_name, aogn);*/
    /*collection<observation_vectors_type> aov = make_collection(get_observation_vectors, aog);*/
    MatrixXd mgo(chr->marker_locus.size(), pmd->qtl_gen);
    for (size_t i = 0; i < pmd->aogn.size(); ++i) {
        auto mo = make_value<Mem>(marker_observations, pop, chr, *pmd->aog[i], ped[i]);
        mgo[*pmd->aog[i]]
            = (*pmd->aog[i])->raw_observations(*pmd->aov[i], *mo);
    }
    pmd->qtl_gen->update_locus_vectors(mgo);
    return mgo;
}
#endif

stpom_data
compute_state_to_parental_origin_haplo(const context_key& ck, const locus_key& lk);

value<stpom_data>
get_stpom_data(const context_key& ck, const locus_key& lk)
{
    return make_value<Sync|Mem>(compute_state_to_parental_origin_haplo,
                                value<context_key>{ck}, value<locus_key>{lk});
}


int
get_n_parents(const context_key& ck, const locus_key& lk)
{
    return (int) get_stpom_data(ck, lk)->row_labels.size();
}



MatrixXb
get_contrast_groups(generation_value gen, const locus_key& lk)
{
    /* FIXME! */
#if 1
    context_key ck(new context_key_struc(gen));
    value<stpom_data> sd = get_stpom_data(ck, lk->parent);
    MatrixXb ret(2, sd->haplo1.innerSize());
    MatrixXb collapse = MatrixXb::Ones(sd->haplo1.outerSize(), 1);
    /*MSG_DEBUG((*sd));*/
    /*MSG_DEBUG(MATRIX_SIZE(ret));*/
    /*MSG_DEBUG(MATRIX_SIZE(collapse) << std::endl << collapse);*/
    (void)*sd;
    ret.row(0) = (sd->haplo1 * collapse).transpose();
    ret.row(1) = (sd->haplo2 * collapse).transpose();
    /*MSG_DEBUG("CONTRAST GROUPS " << gen->name << std::endl << ret);*/
#else
    MatrixXb ret;
#endif
    return ret;
}




labelled_matrix<MatrixXd, std::vector<char>, label_type>
compute_state_to_parental_origin(const context_key& ck, const locus_key& lk)
{
    value<stpom_data>
        sd = get_stpom_data(ck, lk);

    labelled_matrix<MatrixXd, std::vector<char>, label_type>
        ret(sd->row_labels, sd->col_labels);

    ret.data = sd->haplo1.cast<double>() + sd->haplo2.cast<double>();
    /*MSG_DEBUG("stfopom " << MATRIX_SIZE(ret));*/
    return ret;
}


labelled_matrix<MatrixXd, std::vector<char>, label_type>
compute_state_to_dominance(const context_key& ck, const locus_key& lk)
{
    auto stdom = compute_state_to_parental_origin(ck, lk);
    std::map<std::vector<char>, VectorXd> dominance_vectors;
    int nrows = stdom.rows();
    for (int r1 = 0; r1 < nrows; ++r1) {
        for (int r2 = r1 + 1; r2 < nrows; ++r2) {
            VectorXd prod = (stdom.data.row(r1).array() * stdom.data.row(r2).array()).matrix();
            if (prod.sum() > 0) {
                std::vector<char> label;
                label.push_back('(');
                label.insert(label.end(), stdom.row_labels[r1].begin(), stdom.row_labels[r1].end());
                label.push_back('/');
                label.insert(label.end(), stdom.row_labels[r2].begin(), stdom.row_labels[r2].end());
                label.push_back(')');
                dominance_vectors[label] = prod;
            }
        }
    }

    stdom.row_labels.clear();
    stdom.row_labels.reserve(dominance_vectors.size());
    stdom.data.resize(dominance_vectors.size(), stdom.data.cols());
    for (const auto& lv: dominance_vectors) {
        stdom.data.row(stdom.row_labels.size()) = lv.second;
        stdom.row_labels.emplace_back(lv.first);
    }

    /*MSG_DEBUG("Computed stdom matrix" << std::endl << stdom);*/

    return stdom;
}


stpom_data
compute_state_to_parental_origin_haplo(const context_key& ck, const locus_key& lk)
{
    /*const impl::generation_rs* gen = ck->gen;*/
    /* FIXME! */
    /*const impl::generation_rs* gen = NULL;*/
    const geno_matrix* gen = ck->pop ? ck->pop->gen.get() : ck->gen;
    stpom_data sd;
    /*MSG_DEBUG("Computing stfopom order=" << lk->depth() << " gen@" << gen);*/
    /*MSG_DEBUG_INDENT;*/
    /*MSG_DEBUG((*gen));*/
    if (!(lk && lk->locus != locus_key_struc::no_locus)) {
        std::set<char> parents;
        const std::vector<label_type>& ap_labels = gen->get_unique_labels();
        for (const label_type& ap: ap_labels) {
            parents.insert(ap.first());
            parents.insert(ap.second());
        }
        std::vector<std::vector<char>> parent_labels;
        std::map<char, int> parent_index;
        parent_labels.reserve(parents.size());
        for (char c: parents) {
            parent_labels.push_back({c});
            size_t sz = parent_index.size();
            parent_index[c] = sz;
        }
        sd.row_labels.assign(parent_labels.begin(), parent_labels.end());
        sd.col_labels.assign(ap_labels.begin(), ap_labels.end());
        sd.haplo1 = MatrixXb::Zero(parent_labels.size(), ap_labels.size());
        sd.haplo2 = MatrixXb::Zero(parent_labels.size(), ap_labels.size());
        for (size_t i = 0; i < ap_labels.size(); ++i) {
            sd.haplo1(parent_index[ap_labels[i].first()], i) = 1;
            sd.haplo2(parent_index[ap_labels[i].second()], i) = 1;
        }
        /*MSG_DEBUG("prout");*/
        /*MSG_DEBUG(sd);*/
    } else {
        value<stpom_data> first = get_stpom_data(ck, NULL);
        value<stpom_data> pred = get_stpom_data(ck, lk->parent);

        sd.col_labels.reserve(gen->size()
                              * pred->col_labels.size());
        sd.row_labels.reserve(pred->row_labels.size()
                              * first->row_labels.size());

        for (auto& pred_par_lab: pred->row_labels) {
            for (auto& row: first->row_labels) {
                sd.row_labels.push_back(pred_par_lab);
                sd.row_labels.back().push_back(row.front());
            }
        }

        size_t sz = gen->get_unique_labels().size();
        for (size_t i = 0; i < sz; ++i) {
            sd.col_labels.insert(sd.col_labels.end(), pred->col_labels.begin(), pred->col_labels.end());
        }

        /* NB : c'est d'une puissance... */

        sd.haplo1 = kroneckerProduct(pred->haplo1, first->haplo1);
        sd.haplo2 = kroneckerProduct(pred->haplo2, first->haplo2);
    }

    /*MSG_DEBUG("stfopom " << lk);*/
    /*MSG_DEBUG(sd);*/

    /*MSG_DEBUG_DEDENT;*/
    /*MSG_QUEUE_FLUSH();*/

    return sd;
}