frontends4.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_FRONTENDS_H_
#define _SPEL_FRONTENDS_H_

#include "cache2.h"
#include "basic_data.h"
#include "model.h"
// #include "../../attic/frontends2.h"
/*#include "model/tests.h"*/
#include <regex>

#include <boost/math/distributions/normal.hpp> // for normal_distribution

// #include "report_frontend.h"


// #include "excel_report.h"

#define SIGNAL_DISPLAY_ONELINER

#define CONSTRAINT_RESIDUAL_EPSILON 1.e-5

  using boost::math::normal; // typedef provides default type is double.
  using boost::math::cdf;
  using boost::math::mean;
  using boost::math::variance;
  using boost::math::quantile;
  using boost::math::complement;


typedef std::pair<const chromosome*, double> selected_locus;

inline bool operator < (const selected_locus& sl1, const selected_locus& sl2) { return sl1.first < sl2.first || (sl1.first == sl2.first && sl1.second < sl2.second); }

inline std::ostream& operator << (std::ostream& os, const selected_locus& sl) { return os << sl.first->name << ':' << sl.second; }


struct chromosome_search_domain {
    const chromosome* chrom;
    std::vector<double> loci;

    chromosome_search_domain(const chromosome* c, const std::vector<double>& l) : chrom(c), loci(l) {}

    struct const_iterator {
        const chromosome_search_domain* this_csd;
        std::vector<double>::const_iterator i;
        const_iterator() : this_csd(NULL), i(__()) {}
        const_iterator(const chromosome_search_domain* t, const std::vector<double>::const_iterator& i_) : this_csd(t), i(i_) {}
        bool operator == (const const_iterator& other) const { return i == other.i; }
        bool operator != (const const_iterator& other) const { return i != other.i; }
        /*bool operator < (const const_iterator& other) const { return i < other.i; }*/
        /*size_t operator - (const const_iterator& other) const { return i - other.i; }*/
        std::pair<const chromosome*, double> operator * () const { return {this_csd->chrom, *i}; }
        const_iterator& operator ++ () { ++i; return *this; }
        const_iterator& operator -- () { --i; return *this; }
        static std::vector<double>::const_iterator __() { static std::vector<double> _; return _.end(); }
    };

    const_iterator begin() const { return {this, loci.begin()}; }
    const_iterator end() const { return {this, loci.end()}; }
    const_iterator cbegin() const { return {this, loci.begin()}; }
    const_iterator cend() const { return {this, loci.end()}; }
};


typedef std::vector<chromosome_search_domain> genome_search_domain;

struct gsd_iterator {
    std::vector<chromosome_search_domain>::const_iterator csd_i, csd_j;
    chromosome_search_domain::const_iterator i, j;

    bool operator == (const gsd_iterator& other) const { return csd_i == other.csd_i && i == other.i; }
    bool operator != (const gsd_iterator& other) const { return csd_i != other.csd_i || i != other.i; }

    gsd_iterator&
        operator ++ ()
        {
            ++i;
            if (i == j) {
                MSG_DEBUG("at end of chromosome!");
                ++csd_i;
                if (csd_i != csd_j) {
                    i = csd_i->begin();
                    j = csd_i->end();
                } else {
                    i = {};
                    j = {};
                }
            }
            return *this;
        }

    std::pair<const chromosome*, double> operator * () const { return *i; }
};

namespace std {
    inline gsd_iterator begin(const genome_search_domain& gsd) { return {gsd.begin(), gsd.end(), gsd.begin()->begin(), gsd.begin()->end()}; }
    inline gsd_iterator end(const genome_search_domain& gsd) { return {gsd.end(), gsd.end(), {}, {}}; }
}


typedef std::vector<selected_locus> locus_set;

typedef std::vector<locus_set> model_descriptor;

std::pair<bool, double>
detect_strongest_qtl(chromosome_value chr, const locus_key& lk,
                     const model& M0, const std::vector<double> pos);

MatrixXd
ftest_along_chromosome(chromosome_value chr, const locus_key& lk,
                       const model& M0, const std::vector<double> pos);



struct signal_display {
#ifdef SIGNAL_DISPLAY_ONELINER
    static const char* tick(double x)
    {
        static const char* ticks[9] = { " ", "\u2581", "\u2582", "\u2583", "\u2584", "\u2585", "\u2586", "\u2587", "\u2588" };
        return ticks[x < 0. ? 0
                            : x >= 1. ? 8
                                      : ((int) floor(x * 9))];
    }

    VectorXd values;
    int imax_;
    bool above_;

    signal_display(const VectorXd& v, int imax, bool above)
        : values(v.innerSize()), imax_(imax), above_(above)
    {
        values = v;
        double vmin = values.minCoeff();
        double vmax = values.maxCoeff();
        if (vmin == vmax) {
            values = (values.array() - vmin).matrix();
        } else {
            values = ((values.array() - vmin) / (vmax - vmin)).matrix();
        }
    }

    friend std::ostream& operator << (std::ostream& os, const signal_display& sd)
    {
        os << _WHITE << '[';
        for (int i = 0; i < sd.values.innerSize(); ++i) {
            if (i == sd.imax_) {
                os << (sd.above_ ? _GREEN : _RED);
            }
            os << tick(sd.values(i));
            if (i == sd.imax_) {
                os << _WHITE;
            }
        }
        return os << ']' << _NORMAL;
    }
#else
    braille_grid grid;

    signal_display(const chromosome& chr, const std::vector<double>& X, const VectorXd& y, int imax, double threshold)
        : grid(build(chr, X, y, imax, threshold))
    {}

    braille_grid
        build(const chromosome& chr, const std::vector<double>& X, const VectorXd& y, int imax, double threshold)
        {
            std::vector<double> Y(y.data(), y.data() + y.size());
            int padding_left = 0;
            int W = (int) (msg_handler_t::termcols() * .8);
            if (W > 1000) {
                W = 80;
            }
            braille_grid chr_map = chr.pretty_print(W, {}, {}, padding_left, false);

            braille_plot plot(W - padding_left, 5, 0, X.back(), 0, std::max(threshold, Y[imax]));
            plot.plot(X, Y);
            plot.hline(threshold, 1, 1, 0, 255, 0);
            bool above = Y[imax] > threshold;
            plot.vline(X[imax], 1, 0, above ? 0 : 255, above ? 255 : 0, 0);
            return plot.compose_vert(true, chr_map, false);
        }

    friend
        std::ostream&
        operator << (std::ostream& os, const signal_display& sd)
        {
            std::stringstream tmp;
            tmp << sd.grid;
            return os << tmp.str();
        }
#endif
};


struct model_manager;
struct test_result;

void make_report(model_manager& mm, const std::vector<std::string>& cmdline);
void report_computation(model& Mbase, std::string trait_name, chromosome_value chr, std::vector<double>& loci, const computation_along_chromosome& cac, const test_result& tr);
void report_algo_phase(std::string descr);


struct QTL {
    std::string chromosome;
    double locus;
    std::vector<double> LOD_loci;
    std::vector<double> LOD;

    QTL(const std::string& n, double l, const std::vector<double>& x, const MatrixXd& y)
        : chromosome(n), locus(l), LOD_loci(x), LOD(y.data(), y.data() + y.size())
    {
        /*MSG_DEBUG("QTL at " << chromosome << ':' << locus);*/
        /*MSG_DEBUG(y);*/
        /*MSG_DEBUG(MATRIX_SIZE(y));*/
        /*MSG_DEBUG("" << LOD);*/
    }

    static
        double
        interpolate(double x0, double y0, double x1, double y1, double yT)
        {
            double delta_x = x1 - x0;
            double delta_y = y1 - y0;
            return delta_x * (yT - y0) / delta_y + x0;
        }

    std::array<double, 4>
        confidence_interval(const std::string &trait, const std::vector<QTL> &selection, ComputationType lod_test_type);

    std::array<double, 4>
        confidence_interval();
};


#include "report.h"



enum probability_mode { Joint, Single, Skip };



typedef std::pair<double, double> forbidden_interval_type;
typedef std::vector<forbidden_interval_type> forbidden_interval_vector_type;
typedef std::map<chromosome_value, forbidden_interval_vector_type> forbidden_interval_map_type;

inline bool operator < (const forbidden_interval_type& fi1, const forbidden_interval_type& fi2) { return fi1.first < fi2.first; }

struct search_interval_type;
struct search_lg_type;


struct test_result {
    search_lg_type* this_interval;
    const chromosome* chrom;
    double locus;
    double test_value;
    int index;
    bool over_threshold;
    model_block_key pop_block_key;
    value<model_block_type> pop_block;
    model_block_key dom_block_key;
    value<model_block_type> dom_block;

    test_result()
        : chrom(NULL), locus(0), test_value(0), index(0), over_threshold(false), pop_block_key(), pop_block(), dom_block_key(), dom_block()
    {}

    test_result(search_lg_type* ti, const chromosome* c, double l, double tv, int i, bool ot, const model_block_key& mbk, const value<model_block_type>& mb, const model_block_key& mbkd, const value<model_block_type>& mbd)
        : this_interval(ti), chrom(c), locus(l), test_value(tv), index(i), over_threshold(ot), pop_block_key(mbk), pop_block(mb), dom_block_key(mbkd), dom_block(mbd)
    {}

    test_result(const test_result& tr)
        : this_interval(tr.this_interval),
        chrom(tr.chrom), locus(tr.locus), test_value(tr.test_value),
        index(tr.index), over_threshold(tr.over_threshold),
        pop_block_key(tr.pop_block_key), pop_block(tr.pop_block),
        dom_block_key(tr.dom_block_key), dom_block(tr.dom_block)
    {}

    test_result&
        operator = (const test_result& tr)
        {
            this_interval = tr.this_interval;
            chrom = tr.chrom;
            locus = tr.locus;
            test_value = tr.test_value;
            index = tr.index;
            over_threshold = tr.over_threshold;
            pop_block_key = tr.pop_block_key;
            pop_block = tr.pop_block;
            dom_block_key = tr.dom_block_key;
            dom_block = tr.dom_block;
            return *this;
        }

    void reset()
    {
        chrom = NULL;
        locus = 0;
        test_value = 0;
        index = 0;
        over_threshold = false;
        /*block_key.selection.clear();*/
        pop_block_key.reset();
        pop_block = value<model_block_type>();
        dom_block_key.reset();
        dom_block = value<model_block_type>();
    }

    friend
        std::ostream& operator << (std::ostream& os, const test_result& tr)
        {
            os << "<result chrom=" << (tr.chrom ? tr.chrom->name : "nil");
            os << " locus=" << tr.locus;
            os << " test=" << tr.test_value;
            os << " at=" << tr.index;
            os << " over?=" << tr.over_threshold;
            if (tr.pop_block_key) {
                os << " block_key=" << tr.pop_block_key;
            } else {
                os << " block_key=none";
            }
            if (tr.dom_block_key) {
                os << " dominance_block_key=" << tr.dom_block_key;
            } else {
                os << " dominance_block_key=none";
            }
            return os << '>';
        }

    bool
        operator < (const test_result& other) const
        {
            return test_value < other.test_value;
        }

    void
        select(model_manager& mm) const;
};


locus_probabilities_type
locus_probabilities(const context_key& ck, const locus_key& lk,
                    /*const MatrixXd& mgo,*/
                    int ind,
                    const std::vector<double>& loci);


struct search_lg_type {
    const chromosome* chrom;
    /* all positions in this interval */
    std::vector<probability_mode> local_modes;
    std::vector<locus_key> local_selections;
    std::map<locus_key, std::vector<double>> loci_by_selection;
    std::vector<double> all_positions, active_loci;
    /* all USED positions in this interval */
    std::vector<double> positions;
    /* positions actually used in the segment test (all_positions \ selection \ forbidden_intervals) */
    std::vector<double> effective_positions;
    /* current selection (subset of base model's selection for this search interval) */
    locus_key selection;
    /* the model blocks along the chromosome currently under study */
    collection<model_block_type> locus_blocks, single_blocks;
    /* dominance matrices per locus per population */
    collection<parental_origin_per_locus_type> dominance_blocks;
    test_result local_max;

    search_lg_type(const chromosome* chr, const std::vector<double>& steps, locus_key sel)
            : chrom(chr)
            , local_modes()
            , local_selections()
            , loci_by_selection()
            , all_positions(steps)
            , active_loci()
            , positions()
            , effective_positions()
            , selection(sel)
            , locus_blocks()
            , dominance_blocks()
            , local_max()
    {
        recompute_modes();
    }

    search_lg_type(const search_lg_type& slg)
            : chrom(slg.chrom)
            , local_modes(slg.local_modes)
            , local_selections(slg.local_selections)
            , loci_by_selection()
            , all_positions(slg.all_positions)
            , active_loci()
            , positions(slg.positions)
            , effective_positions(slg.effective_positions)
            , selection(slg.selection)
            , locus_blocks(slg.locus_blocks)
            , dominance_blocks(slg.dominance_blocks)
            , local_max(slg.local_max)
    {
        recompute_modes();
    }

    void
    compute_blocks(const collection<population_value>& all_pops)
    {
        if (locus_blocks.size()) {
            return;
        }
        for (const auto& vpop: all_pops) {
            context_key ck(new context_key_struc(*vpop, chrom, all_positions));
            for (auto& x: make_collection<Disk>(locus_probabilities,
                                                as_value(ck), as_value(locus_key{}), range<int>(0, (*vpop)->size(), 1),
                                                as_value(ck->loci))) {
                (void)*x;
            }
        }
        locus_blocks
                = compute_parental_origins_multipop(
                all_pops,
                as_value(chrom),
                local_selections,
                all_positions,
                active_loci);
//        MSG_DEBUG("Computing dominance probabilities");
        dominance_blocks
                = compute_dominance_multipop(
                all_pops,
                as_value(chrom),
                local_selections,
                all_positions,
                active_loci);
    }


    size_t
    get_closest_index(double locus)
    {
        size_t i = 0;
        size_t end = all_positions.size() - 1;
//         MSG_DEBUG("get_closest_index(" << locus << ") all_positions.size=" << all_positions.size());
//         MSG_DEBUG("position #" << i << " at " << all_positions[i]);
        for (; i < end && all_positions[i] < locus; ++i);  // MSG_DEBUG("position #" << i << " at " << all_positions[i]);
        if (i == 0 && i == end) {
            return i;
        }
        double l = locus - all_positions[i - 1];
        double r = all_positions[i] - locus;
        return i - (l < r);
    }
    
    double
    match_locus(double locus)
    {
        size_t i = get_closest_index(locus);
        double new_locus = all_positions[i];
        if (new_locus != locus) {
            MSG_WARNING("Locus " << chrom->name << ':' << locus << " adjusted to test position " << new_locus);
        }
        return new_locus;
    }
    
    void
    recompute_modes(double force_pos=-1)
    {
        local_modes.clear();
        local_modes.resize(all_positions.size(), Single);
        local_selections.clear();
        local_selections.resize(all_positions.size());
        active_loci.clear();
        active_loci.reserve(all_positions.size());
        if (force_pos != -1) {
            force_pos = match_locus(force_pos);
        }
        for (size_t i = 0; i < all_positions.size(); ++i) {
            double locus = all_positions[i];
            if (!!selection) {
                for (double l: selection) {
                    if (locus != force_pos && locus >= l - active_settings->Rskip && locus <= l + active_settings->Rskip) {
                        local_modes[i] = Skip;
                        local_selections[i].reset();
                        break;
                    } else if (locus >= l - active_settings->Rjoin && locus <= l + active_settings->Rjoin) {
                        local_modes[i] = Joint;
                        local_selections[i] = local_selections[i] + l;
                    }
                }
            }
            if (local_modes[i] != Skip) {
                active_loci.push_back(locus);
            }
            loci_by_selection[local_selections[i]].push_back(locus);
        }
        locus_blocks.clear();
        dominance_blocks.clear();
    }


    void
    test(const collection<population_value>& all_pops, int i0, computation_along_chromosome& cac,
         value<ComputationType> vct, value<ComputationResults> vcr, size_t y_block_cols,
         const value<model>& Mcurrent, const value<model>& Mbase, double threshold)
    {
        if (active_loci.size()) {
            compute_blocks(all_pops);
            compute_along_interval<>(i0, cac, vct, vcr, y_block_cols, Mcurrent, Mbase, local_selections, chrom, active_loci, locus_blocks);
            local_max = find_max(i0, vct, cac, threshold);
        } else {
            local_max.reset();
        }
    }

    std::pair<value<model_block_type>, value<model_block_type>>
    compute_at(const collection<population_value>& all_pops, double position)
    {
        position = match_locus(position);
        recompute_modes(position);
        compute_blocks(all_pops);
        auto it = std::find(active_loci.begin(), active_loci.end(), position);
        if (it != active_loci.end()) {
            auto ofs = it - active_loci.begin();
            return {locus_blocks[ofs], dominance_blocks[ofs]};
        }
        return {};
    }
    
    test_result
    find_max(int i0, value<ComputationType> vct, computation_along_chromosome &cac, double threshold)
    {
        MSG_DEBUG("call to find_max");
        auto ct = *vct;
        auto last_computation =
                (ct == ComputationType::FTest    ? cac.ftest_pvalue
                                                 : ct == ComputationType::FTestLOD ? cac.ftest_lod
                                                                                   : ct == ComputationType::Chi2     ? cac.chi2
                                                                                                                     :/* has to be Chi2LOD */           cac.chi2_lod).middleCols(i0, active_loci.size());
        if (active_loci.size() != (size_t)(last_computation.outerSize())) {
            MSG_ERROR("LOCI INCONSISTENT WITH COMPUTATION RESULT (" << active_loci.size() << " vs " << last_computation.outerSize() << ")", "");
            abort();
        }
        int i_max = -1;
        double max = -1;
        for (int i = 0; i < last_computation.outerSize(); ++i) {
            if (last_computation(0, i) >= max) {
                max = last_computation(0, i);
                i_max = i;
            }
        }

#ifdef SIGNAL_DISPLAY_ONELINER
        signal_display sd(last_computation.transpose(), i_max, max > threshold);
        MSG_DEBUG("[COMPUTATION] " << active_loci.front() << sd << active_loci.back() << " max=" << max << " at " << active_loci[i_max]);
#else
        signal_display sd(*chrom, active_loci, last_computation.transpose(), i_max, threshold);
        MSG_DEBUG("[COMPUTATION] " << active_loci.front() << " ... " << active_loci.back() << " max=" << max << " at " << active_loci[i_max] << std::endl << sd);
#endif

        if (i_max == -1) {
            return {};
        }
        /*model_block_key k = locus_base_key;*/
        model_block_key mbk = model_block_key_struc::pop(chrom, local_selections[i_max] + active_loci[i_max]);

        return {
                this, chrom,
                active_loci[i_max], max, i_max, max > threshold,
                mbk,
                locus_blocks[i_max],
                dominance_blocks[i_max]->cols() ? model_block_key_struc::dominance(mbk) : model_block_key{},
                dominance_blocks[i_max]
        };
    }

    void
    select(const test_result& tr, value<model> M)
    {
        MSG_DEBUG("select test result " << tr << " model " << M->keys());
        assert(chrom == tr.chrom);
        selection = selection + tr.pop_block_key->loci;
        recompute_modes();
        M->add_block(tr.pop_block_key, tr.pop_block);
        if (tr.dom_block_key) {
            model ref = *M;
            ref.compute();
            M->add_block_if_test_is_good(tr.dom_block_key, tr.dom_block, ref);
        }
    }

    void
    deselect(double position, const collection<population_value>& all_pops, value<model> M)
    {
        size_t i_locus = get_closest_index(position);
        position = all_positions[i_locus];
        auto i = M->m_blocks.begin(), j = M->m_blocks.end();
        std::vector<decltype(M->m_blocks.begin())> to_remove;
        for (; i != j; ++i) {
            auto& key = i->first;
//             auto& value = i->second;
            if (key->has_locus(chrom, position)) {
                to_remove.emplace_back(i);
            }
        }

        for (auto it: to_remove) {
            auto& key = it->first;
            MSG_DEBUG("Deselecting " << position << " in block " << key);
            if (!key->remove_position(chrom, position)) {
                M->m_blocks.erase(it);
            } else {
                MSG_DEBUG("Need to reduce block (remaining key is " << key << ')');
                reduce(all_pops, i_locus, it->second);
            }
        }
        for (auto& ls: local_selections) {
            ls = ls - position;
        }
//         local_selections[i_locus] = local_selections[i_locus] - position;
        recompute_modes();
        /*deselect(position);*/
    }

    std::pair<model_block_key, model_block_key>
    select(double pos)
    {
        pos = match_locus(pos);
        MSG_DEBUG("select position " << pos);
        size_t i_select = std::find(all_positions.begin(), all_positions.end(), pos) - all_positions.begin();
        std::pair<model_block_key, model_block_key> ret;
        if (local_modes[i_select] == Joint) {
            ret.second = model_block_key_struc::pop(chrom, selection);
        }
        selection = selection + pos;
        if (local_modes[i_select] == Single) {
            ret.first = model_block_key_struc::pop(chrom, locus_key{} + pos);
        } else {
            ret.first = model_block_key_struc::pop(chrom, selection);
        }
        locus_blocks.clear();
        recompute_modes();
        return ret;
    }

    std::pair<model_block_key, model_block_key>
    deselect(double pos)
    {
        pos = match_locus(pos);
        MSG_DEBUG("deselect position " << pos);
        size_t i_select = std::find(all_positions.begin(), all_positions.end(), pos) - all_positions.begin();
        std::pair<model_block_key, model_block_key> ret;
        if (local_modes[i_select] == Joint) {
            ret.second = model_block_key_struc::pop(chrom, selection);
        } else {
            ret.second = model_block_key_struc::pop(chrom, locus_key{} + pos);
        }
        selection = selection - pos;
        if (local_modes[i_select] == Joint && selection) {
            ret.first = model_block_key_struc::pop(chrom, selection);
        }
        locus_blocks.clear();
        recompute_modes();
        return ret;
    }

    void
    reduce(const collection<population_value>& all_pops, size_t i_locus, value<model_block_type>& vblock)
    {
        double position = all_positions[i_locus];
        auto local_sel = local_selections[i_locus] + position;
        /* FIXME what about dominance?? */
        locus_key lk2 = local_selections[i_locus];
        MSG_DEBUG("CALLING REDUCE :: Have position=" << position << " local_sel=" << local_sel << " lk2=" << lk2 << "dim(block)=(" << vblock->rows() << ", " << vblock->cols() << ')');
        auto pop_it = all_pops.begin();
        auto pb = disassemble_parental_origins_multipop(chrom, local_sel->parent, *vblock, all_pops);
        std::vector<collection<parental_origin_per_locus_type>> all_popl;
        all_popl.reserve(pb.size());
        for (auto& vmat: pb) {
            const qtl_pop_type* pop = **pop_it++;
            context_key ck(new context_key_struc(pop, chrom, std::vector<double>()));
            MatrixXd red = local_sel->reduce(active_settings->parent_count_per_pop_per_chr
                                                     .find(pop->qtl_generation_name)->second.find(chrom)->second,
                                             position);
            MatrixXd data = vmat->data * red;
            vmat->data = data;
            vmat->column_labels = get_stpom_data(ck, lk2)->row_labels;
            all_popl.emplace_back();
            all_popl.back().push_back(vmat);
            MSG_DEBUG("Disassembled (" << vmat->rows() << ", " << vmat->cols() << ')');
        }
        vblock = assemble_parental_origins_multipop(as_value(chrom),
                                                    as_value(lk2),
                                                    all_popl,
                                                    all_pops,
                                                    true)[0];
    }

    void
    clear()
    {
        selection = locus_key{};
        recompute_modes();
    }

    friend
    std::ostream&
    operator << (std::ostream& os, const search_lg_type& sl)
    {
        os << '{' << sl.selection << ' ';
        for (size_t i = 0; i < sl.all_positions.size(); ++i) {
            double locus = sl.all_positions[i];
            probability_mode mode = sl.local_modes[i];
            bool selected = !!sl.selection && sl.selection->has(locus);
            if (i) {
                os << ' ';
            }
            if (selected) {
                os << '[';
            }
            os << locus << ':' << (mode == Skip ? "Skip" : mode == Joint ? "Joint" : "Single");
            if (selected) {
                os << ']';
            }
        }
        return os << '}';
    }
};


typedef std::map<chromosome_value, std::vector<search_interval_type>> search_interval_map_type;
typedef std::map<chromosome_value, search_lg_type> search_lg_map_type;


inline
search_lg_map_type
skeleton_search_intervals()
{
    search_lg_map_type ret;

    if (active_settings->skeleton_mode == "auto") {
        for (const chromosome& chr: active_settings->map) {
            std::vector<double> pos;
            double accept = -1.;
            for (double l: chr.condensed.marker_locus) {
                if (l > accept) {
                    pos.push_back(l);
                    accept = l + active_settings->skeleton_interval;
                }
            }
            if (pos.size() > 0) {
                ret.emplace(&chr, search_lg_type{&chr, pos, locus_key{}});
            }
        }
    } else if (active_settings->skeleton_mode == "manual") {
        std::map<chromosome_value, std::vector<double>> pos;
        for (const auto& name: active_settings->skeleton_markers) {
            bool found = false;
            for (const chromosome& chr: active_settings->map) {
                auto li = chr.condensed.marker_locus.begin();
                auto lj = chr.condensed.marker_locus.end();
                auto ni = chr.condensed.marker_name.begin();
                for (; li != lj; ++li, ++ni) {
                    if (*ni == name) {
                        pos[&chr].push_back(*li);
                        found = true;
                        break;
                    }
                }
                if (found) {
                    break;
                }
            }
        }
        for (const auto& kv: pos) {
            ret.emplace(kv.first, search_lg_type{kv.first, kv.second, {}});
        }
    }

    return ret;
}




inline
search_lg_map_type
full_search_intervals(const std::map<chromosome_value, locus_key>& selection)
{
    search_lg_map_type ret;
    /* TODO: add definition of intervals to command line: chromosome, mode, beginning, end */

    for (const chromosome& chr: active_settings->map) {
        MSG_DEBUG("full_search_intervals on " << chr.name);
        auto si = selection.find(&chr);
        locus_key sel;
        if (si != selection.end()) {
            sel = si->second;
        }
        ret.emplace(&chr, search_lg_type{&chr, active_settings->estimation_loci[&chr], sel});
        MSG_DEBUG("Current interval for chromosome " << chr.name << ": " << ret.find(&chr)->second);
    }

    return ret;
}


void
read_locus_list(std::string& s, settings_t* target);


struct model_manager {

    /* name of the studied single_trait */
    std::string trait_name;
    double threshold;
    /* populations used in this model */
    collection<population_value> all_pops;
    /* thy reference model */
    value<model> vMcurrent;
    size_t y_block_cols;
    /* thy maximum order for interaction blocks */
    size_t max_order;

    /* all steps, split by haplotypic intervals, inside which joint probabilities must be computed */
    search_lg_map_type search_intervals;
    /* the current selection split by haplotypic interval */
    std::map<chromosome_value, std::vector<locus_key>> selection;
    /* a collection of intervals to not search in, if the user wishes so */
    forbidden_interval_map_type forbidden_intervals;



    /* structure to cache the results of the last computation along the chromosome */
    std::map<chromosome_value, computation_along_chromosome> cac;

    std::map<chromosome_value, std::vector<double>> testpos;

    /* additiona, colpops.front()->ancestor_namesl data for output and postprocessing */
    std::map<std::string, std::map<double, std::pair<double, double>>> qtl_confidence_interval;

    std::map<chromosome_value, MatrixXd*> last_computation;
    std::map<chromosome_value, std::vector<double>> last_test_positions;
    std::map<chromosome_value, test_result> last_best;

    analysis_report* reporter;

/* Construction
 */
    model_manager(const std::string& trait, const collection<population_value>& colpops,
                   const value<MatrixXd>& y,
                   size_t maximum_order = 1,
                   /*ComputationType ct = ComputationType::FTest,*/
                   SolverType st=SolverType::QR)
        : trait_name(trait)
        , all_pops(colpops)
        , vMcurrent(model{y, active_settings->get_threshold(trait, y->cols()), colpops, (*colpops.front())->ancestor_names, maximum_order, st})
        , y_block_cols(vMcurrent->Y().cols())
        , max_order(maximum_order)
        , search_intervals()
        , cac()
        , testpos()
        , qtl_confidence_interval()
        , reporter(NULL)
    {
        vMcurrent->add_block(model_block_key_struc::cross_indicator(), cross_indicator(trait));
        auto covars = make_covariables_block(trait, active_settings->with_covariables);
        MSG_DEBUG("(model_manager) covariables " << active_settings->with_covariables << " covars.cols=" << covars->data.cols());
        if (covars->data.cols() > 0) {
            vMcurrent->add_block(model_block_key_struc::covariables(), covars);
        }
        vMcurrent->compute();
        /*MSG_DEBUG("INITIAL MODEL SIZE: " << MATRIX_SIZE(vMcurrent->X()));*/
        /*MSG_QUEUE_FLUSH();*/
        active_settings->cross_indicator_can_interact = colpops.size() > 1;
    }

    void
        compute_model()
        {
            vMcurrent->compute();
        }

    void
        cleanup_intervals(forbidden_interval_vector_type& intervals)
        {
            std::sort(intervals.begin(), intervals.end());
            forbidden_interval_vector_type output;
            output.reserve(intervals.size());
            output.push_back(intervals.front());
            for (size_t i = 1; i < intervals.size(); ++i) {
                if (output.back().second > intervals[i].first) {
                    output.back().second = std::max(output.back().second, intervals[i].second);
                } else {
                    output.push_back(intervals[i]);
                }
            }
            intervals.swap(output);
        }

    void
        add_forbidden_interval(chromosome_value chr, double start, double end)
        {
            auto& i = forbidden_intervals[chr];
            i.emplace_back(start, end);
            cleanup_intervals(i);
        }

    /* FIXME this can't remove a sub-interval or split intervals. need previous/next position info for each locus to do so. */
    void remove_forbidden_interval(chromosome_value chr, double start, double end)
        {
            auto& i = forbidden_intervals[chr];
            std::pair<double, double> p {start, end};
            auto it = std::find(i.begin(), i.end(), p);
            if (it != i.end()) {
                i.erase(it);
            }
        }

    void
        select(search_lg_type& si, double pos)
        {
            auto blocks = si.compute_at(all_pops, pos);
            auto ins_rm_key = si.select(pos);
            si.compute_blocks(all_pops);
            if (ins_rm_key.second) {
                vMcurrent->remove_block(ins_rm_key.second);
            }
            vMcurrent->add_block_with_interactions(ins_rm_key.first, blocks.first);
            if (blocks.second) {
                model ref(*vMcurrent);
                ref.compute();
                vMcurrent->add_block_with_interactions(model_block_key_struc::dominance(ins_rm_key.first), blocks.second, &ref);
            }
            debug_selection(SPELL_STRING("select(" << si.chrom->name << '[' << si.active_loci.front() << '-' << si.active_loci.back() << "]:" << pos << ')'));
        }

    void
        select(chromosome_value chr, double pos)
        {
            auto it = search_intervals.find(chr);
            if (it != search_intervals.end()) {
                select(it->second, pos);
            }
            debug_selection(SPELL_STRING("select(" << chr->name << ':' << pos << ')'));
        }

    void
        deselect(chromosome_value chr, double pos)
        {
            auto it = search_intervals.find(chr);
            if (it != search_intervals.end()) {
                it->second.deselect(pos, all_pops, vMcurrent);
            }
            it->second.compute_blocks(all_pops);
            debug_selection(SPELL_STRING("deselect(" << chr->name << ':' << pos << ')'));
        }

    void
        clear_selection()
        {
            for (auto& chr_vsi: search_intervals) {
                vMcurrent->remove_blocks_with_chromosome(chr_vsi.first);
                chr_vsi.second.clear();
            }
            debug_selection("clear_selection");
        }

    void
    debug_selection(const std::string& pfx)
    {
        MSG_DEBUG("OPERATION " << pfx << " Current selection " << get_selection() << " current keys " << vMcurrent->keys());
    }

    std::vector<std::string>
        split(const std::string& str, const std::string& regex)
        {
            std::regex re(regex);
            return {std::sregex_token_iterator(str.begin(), str.end(), re, -1), std::sregex_token_iterator()};
        }

    void
        set_selection(const std::string& selstr)
        {
            for (const auto& qtl: split(selstr, "\\s+")) {
                auto chr_loc = split(qtl, ":");
                double d;
                std::stringstream(chr_loc[1]) >> d;
                select(active_settings->find_chromosome(chr_loc[0]), d);
            }
        }

    value<model_block_type>
        create_interaction_block(const value<model_block_type>& b1, const value<model_block_type>& b2)
        {
            value<model_block_type> new_block = model_block_type{};
            auto& labels = new_block->column_labels;
            auto& mat = new_block->data;
            labels.reserve(b1->column_labels.size() * b2->column_labels.size());
            for (const auto& v1: b1->column_labels) {
                for (const auto& v2: b2->column_labels) {
                    /*MSG_DEBUG("creating interaction label from [" << v1 << "] and [" << v2 << ']');*/
                    labels.emplace_back();
                    labels.back().reserve(3 + v1.size() + v2.size());
                    labels.back().push_back('(');
                    labels.back().insert(labels.back().end(), v1.begin(), v1.end());
                    labels.back().push_back('/');
                    labels.back().insert(labels.back().end(), v2.begin(), v2.end());
                    labels.back().push_back(')');
                }
            }
            mat.resize(b1->data.rows(), b1->data.cols() * b2->data.cols());
            for (int i = 0; i < mat.rows(); ++i) {
                mat.row(i) = kroneckerProduct(b1->data.row(i), b2->data.row(i));
            }
            return new_block;
        }

    void
        add_block_with_interactions(model_block_key mbk, const value<model_block_type>& block)
        {
            vMcurrent->add_block_with_interactions(mbk, block);
        }


    const std::map<chromosome_value, computation_along_chromosome>&
        custom_test(ComputationType ct, ComputationResults cr, value<model>& Mbase,
                    const forbidden_interval_map_type& forbidden, bool clear_chromosome=false, chromosome_value chr=NULL)
        {
            std::pair<model_block_collection, std::vector<locus_key>> chrom_sel_backup;
            value<ComputationType> vct = ct;

//             if (reporter && reporter->output_rss) { cr = cr | RSS; }
//             if (reporter && reporter->output_rank) { cr = cr | Rank; }

            value<ComputationResults> vcr = cr | active_settings->additional_table_output;

            cac.clear();
            static forbidden_interval_vector_type no_intervals;
            last_best.clear();
            last_computation.clear();
            testpos.clear();
            last_test_positions.clear();
            for (auto& chr_intervals: search_intervals) {
                if (chr != NULL && chr != chr_intervals.first) {
                    continue;
                }
                chromosome_value chrom = chr_intervals.first;

                if (clear_chromosome) {
                    chrom_sel_backup = remove_chromosome_selection(chrom);
                }

                Mbase->compute();

                auto& si = chr_intervals.second;

//                MSG_DEBUG("active_loci.size=" << si.active_loci.size() << " :: " << si.active_loci);

                init_computation(cac[chrom], ct, cr | active_settings->additional_table_output, si.active_loci.size(), Mbase->Y().rows(), Mbase->Y().cols());

//                auto i0i = i0_vec.begin();

                si.test(all_pops, 0, cac[chrom], vct, vcr, y_block_cols, vMcurrent, Mbase, threshold);
                testpos[chrom].insert(testpos[chrom].end(), si.active_loci.begin(), si.active_loci.end());
                if (last_best[chrom] < si.local_max) {
                    last_best[chrom] = si.local_max;
                }
                last_test_positions[chrom] = si.active_loci;
//                MSG_DEBUG("Last Test Positions on chromosome " << chrom->name << ": " << last_test_positions[chrom]);

                if (clear_chromosome) {
                    restore_chromosome_selection(chrom, chrom_sel_backup);
                }

                switch(ct) {
                    case ComputationType::FTest:
                        last_computation[chrom] = &cac[chrom].ftest_pvalue;
                        break;
                    case ComputationType::FTestLOD:
                        last_computation[chrom] = &cac[chrom].ftest_lod;
                        break;
                    case ComputationType::Chi2:
                        last_computation[chrom] = &cac[chrom].chi2;
                        break;
                    case ComputationType::Chi2LOD:
                        last_computation[chrom] = &cac[chrom].chi2_lod;
                        break;
                    default:
                        last_computation[chrom] = NULL;
                };

                report_computation(*Mbase, trait_name, chrom, si.active_loci, cac[chrom], si.local_max);
            }
            return cac;
        }

    std::vector<model_block_key>
        keys() const
        {
            std::vector<model_block_key> ret;
            return ret;
        }

    locus_key
        get_selection(const chromosome* chr) const
        {
            return search_intervals.find(chr)->second.selection;
        }

    std::pair<model_block_collection, std::vector<locus_key>>
        remove_chromosome_selection(const chromosome* chr)
        {
            auto& si = search_intervals.find(chr)->second;
            return {vMcurrent->extract_blocks_with_chromosome(chr), {si.selection}};
        }

    void
        restore_chromosome_selection(chromosome_value chrom, const std::pair<model_block_collection, std::vector<locus_key>>& chrom_sel)
        {
            MSG_DEBUG("before restore_chromosome_selection " << vMcurrent->keys());
            auto& si = search_intervals.find(chrom)->second;
            si.selection = chrom_sel.second.front();
            si.recompute_modes();
            vMcurrent->add_blocks(chrom_sel.first);
            MSG_DEBUG("after restore_chromosome_selection " << vMcurrent->keys());
        }

    std::map<chromosome_value, locus_key>
        get_selection() const
        {
            std::map<chromosome_value, locus_key> ret;
            for (const auto& chr_vsi: search_intervals) {
                chromosome_value chr = chr_vsi.first;
//                for (const auto& si: chr_vsi.second) {
                    if (chr_vsi.second.selection->depth() > 0) {
                        ret[chr] = /*ret[chr] +*/ chr_vsi.second.selection;
                    }
//                }
            }
            return ret;
        }

    void
    set_selection(const std::map<chromosome_value, locus_key>& sel)
    {
        for (const auto& kv: sel) {
            chromosome_value chr = kv.first;
            for (double l: kv.second) {
                select(chr, l);
            }
        }
    }

    friend
        std::ostream&
        operator << (std::ostream& os, const std::map<chromosome_value, locus_key>& sel)
        {
            if (sel.size() == 0) {
                return os << "(empty selection)";
            }
            static std::string empty, sep(", ");
            const std::string* _ = &empty;
            for (const auto& kv: sel) {
                for (double l: kv.second) {
                    os << (*_) << kv.first << ':' << l;
                    _ = &sep;
                }
            }
            return os;
        }

    std::vector<QTL>
        QTLs()
        {
            std::vector<QTL> ret;

            auto sel = get_selection();
//             MSG_DEBUG("Computing LOD values for QTLs in selec