model.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_MODEL_MODEL_H_
#define _SPEL_MODEL_MODEL_H_

#include <stdexcept>

/*#include <Eigen/SVD>*/
/*#include <Eigen/QR>*/
#include "eigen.h"
#include "beta_gamma.h"

#include <cmath>


#include "labelled_matrix.h"
#include "settings.h"
#include "print.h"
#include "data/chromosome.h"

/*#define COMPONENT_EPSILON (1.e-10)*/
#define COMPONENT_EPSILON (active_settings->tolerance)

typedef labelled_matrix<Eigen::Matrix<double, -1, -1>, int, std::vector<char>> model_block_type;

struct model_block_key_struc;
typedef std::shared_ptr<model_block_key_struc> model_block_key;


struct model;

static inline
void f_test(const model& model_current, const model& model_new, int col_num, MatrixXd* pvalue, MatrixXd* lod);

struct key_already_exists : public std::exception {
    std::string m_what;
    key_already_exists(const model&, const model_block_key&);
    const char* what() const throw()
    {
        return m_what.c_str();
    }
};

typedef enum { mbk_Covar, mbk_CI, mbk_POP, mbk_Dominance, mbk_Interaction } key_type;

inline
collection<model_block_type>
disassemble_block_multipop(const model_block_type& block, const collection<const qtl_pop_type*>& all_pops)
{
    VectorXi keep = VectorXi::Zero(block.cols());
    int r0 = 0;
    collection<model_block_type> ret;
    ret.reserve(all_pops.size());
    for (const auto& vpop: all_pops) {
        const qtl_pop_type* pop = *vpop;
        auto sub_block = block.data.middleRows(r0, pop->size());
        size_t ncols = 0;
        for (int col = 0; col < sub_block.cols(); ++col) {
            keep(col) = !sub_block.col(col).isZero(active_settings->tolerance);
            ncols += keep(col);
        }
        /*MSG_DEBUG("disasm block keep = " << keep.transpose());*/
        /*MSG_QUEUE_FLUSH();*/
        ret.emplace_back(model_block_type{});
        auto& pop_block = *ret.back();
        pop_block.column_labels.reserve(ncols);
        pop_block.data.resize(pop->size(), ncols);
        int pcol = 0;
        for (int col = 0; col < sub_block.cols(); ++col) {
            if (keep(col)) {
                pop_block.data.col(pcol++) = sub_block.col(col);
                pop_block.column_labels.push_back(block.column_labels[col]);
            }
        }
        r0 += pop->size();
    }
    return ret;
}

inline
value<model_block_type>
assemble_block_multipop(const collection<model_block_type>& pop_blocks)
{
    std::map<std::vector<char>, int> col_indices;
    for (const auto& pb: pop_blocks) {
        for (const auto& v: pb->column_labels) {
            col_indices[v] = 0;
        }
    }

    int index = 0;
    for (auto& kv: col_indices) {
        kv.second = index++;
    }

    int r0 = 0;

    for (const auto& vpb: pop_blocks) {
        r0 += vpb->rows();
    }

    value<model_block_type> ret = model_block_type{};
    ret->data = MatrixXd::Zero(r0, col_indices.size());
    ret->column_labels.reserve(col_indices.size());
    for (const auto& kv: col_indices) {
        ret->column_labels.push_back(kv.first);
    }

    r0 = 0;

    for (const auto& vpb: pop_blocks) {
        const auto& pb = *vpb;
        for (int i = 0; i < pb.cols(); ++i) {
            int col = col_indices[pb.column_labels[i]];
            ret->data.block(r0, col, pb.rows(), 1) = pb.data.col(i);
        }
        r0 += pb.rows();
    }

    return ret;
}



struct model_block_key_struc {
    key_type type;

    const chromosome* chr;
    locus_key loci;
    model_block_key left, right;

    model_block_key_struc(key_type kt, const chromosome* c, const locus_key& lk, const model_block_key& c1, const model_block_key& c2)
        : type(kt), chr(c), loci(lk), left(c1), right(c2)
    {}

    size_t
        order() const
        {
            switch (type) {
                case mbk_Covar:
                case mbk_CI:
                case mbk_POP: return 1; // loci->depth();
                case mbk_Dominance: return 2; // 2 * loci->depth();
                case mbk_Interaction: return left->order() + right->order();
            };
            return 0;
        }

    static
    model_block_key
    cross_indicator()
    {
        model_block_key void_;
        locus_key void_lk;
        return std::make_shared<model_block_key_struc>(mbk_CI, (const chromosome*) NULL, void_lk, void_, void_);
    }

    static
    model_block_key
    covariables()
    {
        model_block_key void_;
        locus_key void_lk;
        return std::make_shared<model_block_key_struc>(mbk_Covar, (const chromosome*) NULL, void_lk, void_, void_);
    }

    static
        model_block_key
        pop(const chromosome* c, locus_key loci)
        {
            model_block_key void_;
            return std::make_shared<model_block_key_struc>(mbk_POP, c, loci, void_, void_);
        }

    static
        model_block_key
        dominance(const model_block_key& haplo)
        {
            model_block_key void_;
            locus_key void_lk;
            return std::make_shared<model_block_key_struc>(mbk_Dominance, (const chromosome*) NULL, void_lk, haplo, void_);
        }

    static
        model_block_key
        interaction(const model_block_key& l, const model_block_key& r)
        {
            locus_key void_lk;
            if (l < r) {
                return std::make_shared<model_block_key_struc>(mbk_Interaction, (const chromosome*) NULL, void_lk, l, r);
            } else {
                return std::make_shared<model_block_key_struc>(mbk_Interaction, (const chromosome*) NULL, void_lk, r, l);
            }
        }

    bool
    is_valid()
    {
        switch (type) {
            case mbk_Covar:
            case mbk_CI:          return true;
            case mbk_POP:         return !loci->is_empty();
            case mbk_Dominance:   return left->is_valid();
            case mbk_Interaction: return left->is_valid() && right->is_valid();
        };
        return false;  /* make compiler happy. */
    }

    friend
        std::ostream&
        operator << (std::ostream& os, const model_block_key_struc& mbk)
        {
            switch (mbk.type) {
                case mbk_Covar:
                    os << "Covar";
                    break;
                case mbk_CI:
                    os << "Cross";
                    break;
                case mbk_POP:
                    os << mbk.chr->name << ':' << mbk.loci;
                    break;
                case mbk_Dominance:
                    os << "Dom(" << mbk.left << ')';
                    break;
                case mbk_Interaction:
                    os << '(' << mbk.left << ":" << mbk.right << ')';
                    break;
            };
            return os;
        }

    friend
        std::ostream&
        operator << (std::ostream& os, const model_block_key& mbk)
        {
            return os << (*mbk);
        }

    bool
        can_interact_with(const model_block_key& mbk) const
        {
            if (type == mbk_CI || mbk->type == mbk_CI) {
                return active_settings->cross_indicator_can_interact;
            }
            auto v1 = flatten();
            auto v2 = mbk->flatten();
            decltype(v1) inter(std::min(v1.size(), v2.size()));
            auto e = std::set_intersection(v1.begin(), v1.end(), v2.begin(), v2.end(), inter.begin());
            return e == inter.begin();
        }

    bool
        has_chromosome(const chromosome* c) const
        {
            switch (type) {
                case mbk_Covar:
                case mbk_CI:
                    return false;
                case mbk_POP:
                    return chr == c;
                case mbk_Dominance:
                    return left->has_chromosome(c);
                case mbk_Interaction:
                    return left->has_chromosome(c) || right->has_chromosome(c);
            };
            return false;
        }

    bool
        has_locus(const chromosome* c, double l) const
        {
            switch (type) {
                case mbk_Covar:
                case mbk_CI:
                    return false;
                case mbk_POP:
                    return chr == c && loci->has(l);
                case mbk_Dominance:
                    return left->has_locus(c, l);
                case mbk_Interaction:
                    return left->has_locus(c, l) || right->has_locus(c, l);
            };
            return false;
        }

    bool
        has(const model_block_key& mbk) const
        {
            return *this == *mbk
                || (left && left->has(mbk))
                || (right && right->has(mbk));
        }

    bool
        operator < (const model_block_key_struc& other) const
        {
            /* CI < POP < Dom < Inter
             * then sort on haplotypes
             * for interactions, sort on left&right children
             */
            switch (type) {
                case mbk_Covar:
                    return other.type != mbk_Covar;
                case mbk_CI:
                    return other.type != mbk_CI && other.type != mbk_Covar;
                case mbk_POP:
                    switch (other.type) {
                        case mbk_Covar:
                        case mbk_CI:
                            return false;
                        case mbk_POP:
                            return chr < other.chr || (chr == other.chr && loci < other.loci);
                        default:
                            return true;
                    };
                case mbk_Dominance:
                    switch (other.type) {
                        case mbk_Covar:
                        case mbk_CI:
                        case mbk_POP:
                            return false;
                        case mbk_Dominance:
                            return (*left) < (*other.left);
                        case mbk_Interaction:
                            return true;
                    };
                case mbk_Interaction:
                    {
                        auto o1 = order();
                        auto o2 = other.order();
                        if (o1 < o2) {
                            return true;
                        }
                        if (o1 > o2) {
                            return false;
                        }
                    }
                    switch (other.type) {
                        case mbk_Interaction:
                            return (*left) < (*other.left) || ((*left) == (*other.left) && (*right) < (*other.right));
                        default:
                            return false;
                    };
            };
            return false;
        }

    bool remove_position(const chromosome* chr, double locus)
    {
        switch (type) {
            case mbk_Covar:
            case mbk_CI:
                return true;
            case mbk_POP:
                if (has_locus(chr, locus)) {
                    loci = loci - locus;
                }
                return !loci->is_empty();
            case mbk_Dominance:
                return left->remove_position(chr, locus);
            case mbk_Interaction:
                bool a = left->remove_position(chr, locus);
                bool b = right->remove_position(chr, locus);
                return a && b;
        };
        if (!is_valid()) {
            MSG_ERROR("Created an invalid key!", "");
            throw 0;
        }
        return false;
    }

    friend
        inline
        bool
        operator < (const model_block_key& k1, const model_block_key& k2)
        {
            if (!k1) {
                return !!k2;
            }
            if (!k2) {
                return false;
            }
            return (*k1) < (*k2);
        }

    bool
        operator == (const model_block_key_struc& other) const
        {
            if (type != other.type) {
                return false;
            }
            switch (type) {
                case mbk_Interaction:
                    if (!((*right) == (*other.right))) {
                        return false;
                    }
                case mbk_Dominance:
                    if (!((*left) == (*other.left))) {
                        return false;
                    }
                    break;
                case mbk_POP:
                    return chr == other.chr && loci == other.loci;
                case mbk_Covar:
                case mbk_CI:
                    return true;
            };
            return false;
        }

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

    operator std::string () const { std::stringstream ret; ret << (*this); return ret.str(); }

private:
    typedef std::pair<const chromosome*, locus_key> item;

    std::vector<item>
        flatten() const
        {
            std::vector<item> ret;
            flatten(ret);
            std::sort(ret.begin(), ret.end(), [](const item& a, const item& b) { return a.first < b.first || (a.first == b.first && a.second < b.second); });
            return ret;
        }

    void
        flatten(std::vector<item>& ret) const
        {
            switch (type) {
                case mbk_Interaction:
                    right->flatten(ret);
                case mbk_Dominance:
                    left->flatten(ret);
                    break;
                case mbk_POP:
                    ret.emplace_back(chr, loci);
                case mbk_Covar:
                case mbk_CI:
                    ;
            };
        }

};


struct mbk_comp {
    bool
        operator () (const model_block_key& k1, const model_block_key& k2) const
        {
            return (*k1) < (*k2);
        }
};


typedef std::map<model_block_key, value<model_block_type>, mbk_comp> model_block_collection;

namespace std {
    template <>
        struct hash<model_block_key> {
            size_t operator () (const model_block_key& mbk) const
            {
                std::hash<const chromosome*> hc;
                std::hash<locus_key> hlk;
                size_t accum;
#if 0
                /* WARNING FIXME this MUST NOT be used to hash a function parameter
                 * in a disk-cached task, because a POINTER is HASHED and the order
                 * is not guaranteed to be the same in every run.
                 */
                for (const auto& kv: mbk.selection) {
                    accum = impl::ROTATE<7>(impl::ROTATE<7>(accum * hc(kv.first)) ^ hlk(kv.second));
                }
#endif
                accum = 0xbadc0def;
                switch (mbk->type) {
                    case mbk_Covar:
                    case mbk_CI:
                        return 0;
                    case mbk_POP:
                        accum = 0xdeadbe3f;
                        accum = impl::ROTATE<7>(impl::ROTATE<7>(accum * hc(mbk->chr)) ^ hlk(mbk->loci));
                        break;
                    case mbk_Dominance:
                        accum = impl::ROTATE<7>(accum * operator () (mbk->left));
                        break;
                    case mbk_Interaction:
                        accum = impl::ROTATE<7>(accum * operator () (mbk->left)) ^ impl::ROTATE<13>(accum * operator () (mbk->right));
                };
                return accum;
            }
        };
} // namespace std



static inline
bool
around_zero(double o)
{
    return o < COMPONENT_EPSILON && o > -COMPONENT_EPSILON;
}

static inline
bool
much_smaller_than(double a, double b)
{
    return a < (COMPONENT_EPSILON * b);
}

static inline
void
set_if_much_smaller_than(double& a, double b)
{
    double tmp = COMPONENT_EPSILON * b;
    if (a < tmp) {
        a = tmp;
    }
}


using namespace Eigen;

static inline
MatrixXd concat_right(const std::vector<MatrixXd>& mat_vec)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (auto& m: mat_vec) {
        full_size += m.outerSize();
        /*MSG_DEBUG("preparing concat_right with matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
    }
    ret.resize(mat_vec.front().innerSize(), full_size);
    full_size = 0;
    for (auto& m: mat_vec) {
        /*MSG_DEBUG("concat_right in M(" << ret.innerSize() << ',' << ret.outerSize() << ") at col " << full_size << "matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
        /*ret.block(0, full_size, ret.innerSize(), m->outerSize()) = *m;*/
        ret.middleCols(full_size, m.outerSize()) = m;
        full_size += m.outerSize();
    }
    return ret;
}

static inline
MatrixXd concat_right(const collection<model_block_type>& mat_vec)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (auto m: mat_vec) {
        full_size += m->outerSize();
        /*MSG_DEBUG("preparing concat_right with matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
    }
    ret.resize(mat_vec.front()->innerSize(), full_size);
    full_size = 0;
    for (auto m: mat_vec) {
        /*MSG_DEBUG("concat_right in M(" << ret.innerSize() << ',' << ret.outerSize() << ") at col " << full_size << "matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
        /*ret.block(0, full_size, ret.innerSize(), m->outerSize()) = *m;*/
        ret.middleCols(full_size, m->outerSize()) = m->data;
        full_size += m->outerSize();
    }
    return ret;
}


static inline
MatrixXd concat_right(const model_block_collection& mat_map)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (auto m: mat_map) {
        full_size += m.second->outerSize();
        /*MSG_DEBUG("preparing concat_right with matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
    }
    ret.resize(mat_map.begin()->second->innerSize(), full_size);
    full_size = 0;
    for (auto m: mat_map) {
        /*MSG_DEBUG("concat_right in M(" << ret.innerSize() << ',' << ret.outerSize() << ") at col " << full_size << "matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
        /*ret.block(0, full_size, ret.innerSize(), m->outerSize()) = *m;*/
        ret.middleCols(full_size, m.second->outerSize()) = m.second->data;
        full_size += m.second->outerSize();
    }
    return ret;
}

static inline
MatrixXd concat_down(const std::vector<MatrixXd>& mat_vec)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (auto& m: mat_vec) {
        full_size += m.innerSize();
    }
    ret.resize(full_size, mat_vec.front().outerSize());
    full_size = 0;
    for (auto& m: mat_vec) {
        ret.middleRows(full_size, m.innerSize()) = m;
        full_size += m.innerSize();
    }
    return ret;
}

static inline
MatrixXd concat_down(const std::vector<const MatrixXd*>& mat_vec)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (auto m: mat_vec) {
        full_size += m->innerSize();
    }
    ret.resize(full_size, mat_vec.front()->outerSize());
    full_size = 0;
    for (auto m: mat_vec) {
        /*ret.block(full_size, 0, m->innerSize(), ret.outerSize()) = *m;*/
        ret.middleRows(full_size, m->innerSize()) = *m;
        full_size += m->innerSize();
    }
    return ret;
}


static inline
std::pair<int, MatrixXd>
rank_and_components(const MatrixXd& M)
{
    JacobiSVD<MatrixXd> svd(M, ComputeThinU);

    std::cout << "Singular values " << svd.singularValues().transpose() << std::endl;
    int nzsv = svd.nonzeroSingularValues();

    return {nzsv, svd.matrixU().leftCols(nzsv)};
}


static inline
MatrixXd components(const MatrixXd& M, const MatrixXd& P)
{
    MatrixXd pnorm(P.innerSize(), P.outerSize());
    for (int i = 0; i < P.outerSize(); ++i) {
        pnorm.col(i) = P.col(i).normalized();
    }
    MatrixXd orth = M - pnorm * pnorm.transpose() * M; /* feu ! */
    return rank_and_components(orth).second;
}


enum class SolverType { QR };

typedef std::vector<std::map<model_block_key, MatrixXd>> constraint_list;

inline
std::ostream& operator << (std::ostream& os, const constraint_list& cl)
{
    for (const auto& map: cl) {
        for (const auto& kv: map) {
            os << kv.first << std::endl << kv.second << std::endl;
        }
        os << "---" << std::endl;
    }
    return os;
}


MatrixXd
contrast_groups(const collection<const qtl_pop_type*>& all_pops, const locus_key& lk);



struct filename_stream {
    struct filtering_streambuf : public std::streambuf {
        std::streambuf* original;
        filtering_streambuf(std::streambuf* _) : original(_) {}
    protected:
        int overflow(int ch) override
        {
            switch (ch) {
                case ',':
                case '}':
                    return ch;
                case '{':
                    ch = '-';
                    break;
                case ' ':
                    ch = '_';
                default:;
            };
            return original->sputc(ch);
        }
    };

    std::ostringstream sstream;
    filtering_streambuf rdbuf;
    std::ostream ostream;

    filename_stream()
        : sstream(), rdbuf(sstream.rdbuf()), ostream(&rdbuf)
    {}

    template <typename T>
        friend
        filename_stream& operator << (filename_stream& fs, T x)
        {
            fs.ostream << x;
            return fs;
        }

    operator
        std::string () const
        {
            return sstream.str();
        }
};


/*#define ULTRA_MEGA_PARANOID*/
static inline
MatrixXd matrix_inverse(const MatrixXd& m)
{
    /*return m.fullPivHouseholderQr().inverse();*/
    /*return m.inverse();*/
    /* FIXME Should NOT use SVD (see note in compute() */
    JacobiSVD<MatrixXd> inverter(m, ComputeFullV);
    auto& V = inverter.matrixV();
    VectorXd inv_sv(inverter.singularValues());
    for (int i = 0; i < inv_sv.innerSize(); ++i) {
        if (!around_zero(inv_sv(i))) {
            inv_sv(i) = 1. / inv_sv(i);
        } else {
            inv_sv(i) = 0.;
        }
    }
#ifdef ULTRA_MEGA_PARANOID
    /* CHECK THAT THIS IS STABLE BY THE BEN ISRAEL SEQUENCE */
    MatrixXd svd_ret = V * inv_sv.asDiagonal() * V.transpose();
    MatrixXd ret = svd_ret;
    MatrixXd test = 2 * ret - ret * m * ret;
    if (!test.isApprox(ret, .00001)) {
        MSG_DEBUG("ret" << std::endl << ret << std::endl << "test" << std::endl << test);
    } else {
        MSG_DEBUG("m^-1 IS GOOD! Yeah man.");
    }
    return ret;
#else
    return V * inv_sv.asDiagonal() * V.transpose();
#endif
}


struct gauss_elimination {
    MatrixXd block;
    std::vector<int> leading;

    gauss_elimination(const MatrixXd& mat)
        : block(mat)
        , leading(mat.rows())
    {
        for (int i = 0; i < block.rows(); ++i) {
            find_leading_coef(i);
        }

        sort_rows(0);
        for (int ref = 0; ref < block.rows(); ++ref) {
            if (is_null(ref)) {
                break;
            }
            anihilate_coef(ref, leading[ref]);
            sort_rows(ref);
        }
    }

    bool is_null(int row) const { return leading[row] == block.cols(); }
    int not_null(int row) const { return (int) (leading[row] != block.cols()); }

    int
        rank() const
        {
            /*MSG_DEBUG("echelon form" << std::endl << block);*/
            int count = 0;
            for (int i = 0; i < block.rows(); ++i) {
                count += not_null(i);
            }
            return count;
        }

    void
        sort_rows(int start)
        {
            std::vector<int> indices(block.rows() - start);
            int i = start;
            std::generate(indices.begin(), indices.end(), [&] () { return i++; });
            std::sort(indices.begin(), indices.end(), [&](int a, int b) { return leading[a] < leading[b]; });
            int sz = (int) indices.size();
            for (i = 0; i < sz; ++i) {
                if (i > indices[i]) {
                    swap_rows(i, indices[i]);
                }
            }
        }

    void
        swap_rows(int r1, int r2)
        {
            VectorXd tmp = block.row(r1);
            block.row(r1) = block.row(r2);
            block.row(r2) = tmp;
            std::swap(leading[r1], leading[r2]);
        }

    void
        anihilate_coef(int ref, int j)
        {
            for (int i = ref + 1; i < block.rows(); ++i) {
                if (leading[i] == j) {
                    anihilate_coef_in_row(ref, i, j);
                }
            }
        }

    void
        anihilate_coef_in_row(int ref, int r, int j)
        {
            block.row(r) *= 1. / block(r, j);
            block.row(r) -= block.row(ref);
            block(r, j) = 0;
            find_leading_coef(r);
            if (not_null(r)) {
                block.row(r) *= 1. / block(r, leading[r]);
            }
        }

    void
        find_leading_coef(int r)
        {
            int j;
            for (j = 0; j < block.cols() && block(r, j) == 0; ++j);
            leading[r] = j;
        }
};


struct model {
    model()
        : m_Y(), m_blocks(), m_X(), m_rank(), m_rss(), m_coefficients(), m_solver_type(), m_computed(false), m_with_constraints(false), m_ghost_constraint(), m_all_pops(), m_ancestor_names(), m_max_order(1), m_threshold(0), m_rank_fixer()

    {}

    model(const value<MatrixXd>& y, double threshold, const collection<const qtl_pop_type*>& pops, const std::map<char, std::string>& anam, size_t mo=1, SolverType st=SolverType::QR)
        : m_Y(y)
        , m_blocks(), m_X()
		, m_rank(), m_rss(), m_coefficients(), m_residuals()
		, m_solver_type(st)
        , m_computed(false)
        , m_with_constraints(false)
        , m_ghost_constraint()
        , m_all_pops(pops)
        , m_ancestor_names(anam)
        , m_max_order(mo)
        , m_threshold(threshold)
        , m_rank_fixer()
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y.innerSize() << ',' << y.outerSize() << ')'); }*/
    {}

    model(const value<MatrixXd>& y, double threshold, const collection<const qtl_pop_type*>& pops, size_t mo=1, SolverType st=SolverType::QR)
        : m_Y(y)
        , m_blocks(), m_X()
		, m_rank(), m_rss(), m_coefficients(), m_residuals()
		, m_solver_type(st)
        , m_computed(false)
        , m_with_constraints(false)
        , m_ghost_constraint()
        , m_all_pops(pops)
        , m_ancestor_names((*pops.front())->ancestor_names)
        , m_max_order(mo)
        , m_threshold(threshold)
        , m_rank_fixer()
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y.innerSize() << ',' << y.outerSize() << ')'); }*/
    {}

    model(const model& mo)
        : m_Y(mo.m_Y)
        , m_blocks(), m_X(mo.m_X)
		, m_rank(mo.m_rank), m_rss(mo.m_rss), m_coefficients(mo.m_coefficients)
        , m_residuals(mo.m_residuals)
		, m_solver_type(mo.m_solver_type)
        , m_computed(mo.m_computed)
        , m_with_constraints(mo.m_with_constraints)
        , m_ghost_constraint()
        , m_all_pops(mo.m_all_pops)
        , m_ancestor_names(mo.m_ancestor_names)
        , m_max_order(mo.m_max_order)
        , m_threshold(mo.m_threshold)
        , m_rank_fixer()
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')'); }*/
    {
        for (const auto& kv: mo.m_blocks) {
            model_block_key mbk = std::make_shared<model_block_key_struc>(*kv.first);
            m_blocks.emplace(mbk, kv.second);
        }
    }

    model(const value<MatrixXd>& y, const model& mo)
        : m_Y(y)
        , m_blocks(), m_X(mo.m_X)
		, m_rank(), m_rss(), m_coefficients()
        , m_residuals()
		, m_solver_type(mo.m_solver_type)
        , m_computed(false)
        , m_with_constraints(mo.m_with_constraints)
        , m_ghost_constraint()
        , m_all_pops(mo.m_all_pops)
        , m_ancestor_names(mo.m_ancestor_names)
        , m_max_order(mo.m_max_order)
        , m_threshold(mo.m_threshold)
        , m_rank_fixer()
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')'); }*/
    {
        for (const auto& kv: mo.m_blocks) {
            model_block_key mbk = std::make_shared<model_block_key_struc>(*kv.first);
            m_blocks.emplace(mbk, kv.second);
        }
    }

# if 0
    model&
        operator = (const model& mo)
        = delete;
#else
    model&
        operator = (const model& mo)
        {
            m_Y = mo.m_Y;
            m_blocks = mo.m_blocks;
            m_computed = mo.m_computed;
            m_with_constraints = mo.m_with_constraints;
            m_ghost_constraint = mo.m_ghost_constraint;
            m_X = mo.m_X;
            m_rss = mo.m_rss;
            m_coefficients = mo.m_coefficients;
            m_residuals = mo.m_residuals;
            m_rank = mo.m_rank;
            m_all_pops = mo.m_all_pops;
            m_ancestor_names = mo.m_ancestor_names;
            m_max_order = mo.m_max_order;
            m_threshold = mo.m_threshold;
            return *this;
        }

    model&
        operator = (model&& mo)
        {
            m_Y = mo.m_Y;
            m_blocks.clear();
            m_blocks.swap(mo.m_blocks);
            m_computed = mo.m_computed;
            m_with_constraints = mo.m_with_constraints;
            m_ghost_constraint = mo.m_ghost_constraint;
            m_X = mo.m_X;
            m_rss = mo.m_rss;
            m_coefficients = mo.m_coefficients;
            m_residuals = mo.m_residuals;
            m_rank = mo.m_rank;
            m_ancestor_names.swap(mo.m_ancestor_names);
            m_max_order = mo.m_max_order;
            m_threshold = mo.m_threshold;
            return *this;
        }
#endif

    struct rank_fixer_type {
        bool m_enabled = false;
        model* m_model = NULL;
        typedef std::map<model_block_key, std::vector<int>> col_order_map_type;
        col_order_map_type m_col_order, m_ctf;
        col_order_map_type::reverse_iterator m_current_block;
        std::vector<int>::iterator m_current_col;
        
        rank_fixer_type() : m_enabled(false), m_model(NULL), m_col_order(), m_ctf() {}
        
        std::vector<int> cols_to_fix(const model_block_key& mbk) { return m_ctf[mbk]; }
        bool enabled() const { return m_enabled; }
        void disable() { m_enabled = false; }
        void enable(model& this_model)
        {
            m_enabled = true;
            m_model = &this_model;
            m_ctf.clear();
            for (const auto& kv: m_model->m_blocks) {
                auto& mb = kv.second;
                auto mass = mb->data.colwise().sum();
                std::vector<int> col_order(mb->cols());
                std::iota(col_order.begin(), col_order.end(), 0);
                std::sort(col_order.begin(), col_order.end(), [&] (int a, int b) { return mass(a) < mass(b); });
                std::swap(m_col_order[kv.first], col_order);
            }
            MSG_DEBUG("Have col_order: " << m_col_order);
            m_current_block = m_col_order.rbegin();
            m_current_col = m_current_block->second.begin();
            if (m_current_block != m_col_order.rend()) {
                m_ctf[m_current_block->first].push_back(*m_current_col);
            }
        }
        
        bool next(bool validate_previous_column)
        {
            MSG_DEBUG("[rank_fixer] next (did" << (validate_previous_column ? "" : "n't") << " keep previously constrained column)");
            if (!validate_previous_column) {
                m_ctf[m_current_block->first].pop_back();
            }
            if (m_current_block == m_col_order.rend()) {
                MSG_DEBUG("[rank_fixer] at end.");
                return false;
            }
            if (m_current_col == m_current_block->second.end()) {
                MSG_DEBUG("[rank_fixer] next block.");
                ++m_current_block;
                if (m_current_block != m_col_order.rend()) {
                    m_current_col = m_current_block->second.begin();
                }
            } else {
                MSG_DEBUG("[rank_fixer] next column.");
                ++m_current_col;
                if (m_current_col == m_current_block->second.end()) {
                    MSG_DEBUG("[rank_fixer] next block.");
                    ++m_current_block;
                    if (m_current_block != m_col_order.rend()) {
                        m_current_col = m_current_block->second.begin();
                    }
                }
            }
            if (m_current_block != m_col_order.rend()) {
                MSG_DEBUG("[rank_fixer] adding " << m_current_block->first << '#' << (*m_current_col));
                m_ctf[m_current_block->first].push_back(*m_current_col);
                return true;
            }
            return false;
        }
    };

    constraint_list
    compute_constraint(const model_block_key& mbk, const model_block_type& mb)
    {
        auto ret = compute_constraint_for_block_type(mbk, mb);
//         MSG_DEBUG("compute_constraint");
        if (m_rank_fixer.enabled()) {
//             MSG_DEBUG("constraint count before column zeroing: " << ret.size());
            std::vector<int> cols_to_fix = m_rank_fixer.cols_to_fix(mbk);
            for (int c: cols_to_fix) {
                MatrixXd colzero = MatrixXd::Zero(1, mb.cols());
                colzero(0, c) = 1;
                ret.emplace_back(std::map<model_block_key, MatrixXd>{{mbk, colzero}});
            }
//             MSG_DEBUG("constraint count after column zeroing: " << ret.size());
        }
        return ret;
    }

    constraint_list
        compute_constraint_for_block_type(const model_block_key& mbk, const model_block_type& mb)
        {
            constraint_list ret;
            if (mbk->type == mbk_CI) {
                return ret;
            /*} else if (mbk.selection.size() == 1) {*/
            } else if (mbk->type == mbk_POP) {
//                size_t n_qtl = mb.column_labels.front().size();
//                if (n_qtl == 1) {
                if (1) {
                    /* only one chromosome: do the contrast groups trick */
                    ret = {{{mbk, contrast_groups(m_all_pops, mbk->loci)}}};
//                } else {
//                    MSG_DEBUG("Creating constraints for " << mbk);
//                    std::vector<std::set<char>> uniq_letters_per_qtl(mb.column_labels.front().size());
//                    /* need the epistasis magic here */
//                    for (const auto& vec: mb.column_labels) {
//                        auto i = uniq_letters_per_qtl.begin();
//                        for (char c: vec) { i->insert(c); ++i; }
//                    }
//                    std::vector<size_t> letter_counts(uniq_letters_per_qtl.size());
//                    for (size_t i = 0; i < letter_counts.size(); ++i) { letter_counts[i] = uniq_letters_per_qtl[i].size(); }
//                    size_t start = 1;
//                    size_t finish = (1 << letter_counts.size()) - 1;
//                    for (size_t variant = start; variant < finish; ++variant) {
//                        MatrixXd constraint = MatrixXd::Identity(1, 1);
//                        for (size_t i = 0; i < letter_counts.size(); ++i) {
//                            bool flat = (variant >> i) & 1;
//                            MatrixXd tmp;
//                            if (flat) {
//                                tmp = kroneckerProduct(constraint, MatrixXd::Ones(1, letter_counts[i]));
//                            } else {
//                                tmp = kroneckerProduct(constraint, MatrixXd::Identity(letter_counts[i], letter_counts[i]));
//                            }
//                            constraint = tmp;
//                        }
//                        ret.emplace_back();
//                        MSG_DEBUG("Created constraint for " << mbk);
//                        MSG_DEBUG(constraint);
//                        ret.back().insert({{mbk, constraint}});
//                    }
                }
            } else if (mbk->type == mbk_Interaction) {
                /*MSG_DEBUG("Computing constraint for interaction " << mbk);*/
                /*MSG_QUEUE_FLUSH();*/
                auto b1 = m_blocks[mbk->left];
                auto b2 = m_blocks[mbk->right];
                int cols1 = b1->data.cols();
                int cols2 = b2->data.cols();
                constraint_list
                    C1 = compute_constraint(mbk->left, *b1),
                    C2 = compute_constraint(mbk->right, *b2);
                if (C1.size() == 0) {
                    C1 = {{{{mbk->left, MatrixXd::Ones(1, b1->cols())}}}};
                }
                if (C2.size() == 0) {
                    C2 = {{{{mbk->left, MatrixXd::Ones(1, b2->cols())}}}};
                }
                /*MSG_DEBUG("Computing interaction constraints from" << std::endl*/
                        /*<< "C1" << std::endl << C1*/
                        /*<< "C2" << std::endl << C2*/
                        /*);*/
                MatrixXd constraint;
                MatrixXd ghost_columns;
                if (mb.cols() < b1->cols() * b2->cols()) {
                    ghost_columns = MatrixXd::Zero(b1->cols() * b2->cols(), mb.cols());
                    int i = 0;
                    int j = 0;
                    /*MSG_DEBUG("Computing ghost columns");*/
                    /*MSG_DEBUG("L1 " << b1->column_labels);*/
                    /*MSG_DEBUG("L2 " << b2->column_labels);*/
                    for (const auto& l1: b1->column_labels) {
                        for (const auto& l2: b2->column_labels) {
                            auto il = make_interaction_label(l1, l2);
                            /*MSG_DEBUG("-- compare " << il << " and " << mb.column_labels[j]);*/
                            if (il == mb.column_labels[j]) {
                                ghost_columns(i, j++) = 1;
                            }
                            i++;
                        }
                    }
                    /*MSG_DEBUG("ghost_columns redux" << std::endl << ghost_columns);*/
                }
                for (const auto& map1: C1) {
                    for (const auto& c1: map1) {
                        constraint = kroneckerProduct(c1.second, MatrixXd::Identity(cols2, cols2));
                        if (ghost_columns.size()) {
                            constraint = constraint * ghost_columns;
                        }
                        ret.emplace_back();
                        ret.back().insert({{mbk, constraint}});
                    }
                }
                for (const auto& map2: C2) {
                    for (const auto& c2: map2) {
                        constraint = kroneckerProduct(MatrixXd::Identity(cols1, cols1), c2.second);
                        if (ghost_columns.size()) {
                            constraint = constraint * ghost_columns;
                        }
                        ret.emplace_back();
                        ret.back().insert({{mbk, constraint}});
                    }
                }
            }
            if (mbk == m_ghost_constraint.first) {
                int n_par = m_blocks[mbk]->outerSize();
                ret.emplace_back();
                ret.back().insert({{m_ghost_constraint.first, MatrixXd::Identity(n_par, n_par)},
                        {m_ghost_constraint.second, -1. * MatrixXd::Identity(n_par, n_par)}});
            }
            /*MSG_DEBUG("Constraints for block " << mbk);*/
            /*MSG_DEBUG("" << ret);*/
            return ret;
        }


    /*bool operator == (const model& m) const*/
    /*{*/
        /*return m_Y == m.m_Y && m_blocks == m.m_blocks && m_rss == m.m_rss;*/
    /*}*/

    void
        add_block(const model_block_key& key, const value<model_block_type>& x)
        {
            /*MSG_DEBUG("add_block(" << key << ", <some block>)");*/
            if (!x) {
                MSG_ERROR("Attempting to insert a null block in the model.", "");
                abort();
            }
            if (!key->is_valid()) {
                MSG_ERROR("Attempting to insert an invalid block key in the model.", "");
                abort();
            }
            if (m_blocks.find(key) != m_blocks.end()) {
                /*MSG_DEBUG("Key already exists");*/
                return;
                /*throw key_already_exists(*this, key);*/
            }
            if (!x) {
                MSG_ERROR("Attempting to insert a null block in the model.", "");
                abort();
            }
            m_computed = false;
            m_blocks.insert({key, x});
            /*MSG_DEBUG("selection now: " << keys());*/
            /*m_blocks.push_back(x);*/
            /*m_keys.push_back(key);*/
        }

    bool
        add_block_if_test_is_good(const model_block_key& key, const value<model_block_type>& x, const model& ref)
        {
            MatrixXd pvalue(m_Y->cols(), 1);
            model tmp(*this);
            tmp.add_block(key, x);
            tmp.compute();
            f_test(ref, tmp, 0, &pvalue, NULL);
            if (pvalue(0, 0) < m_threshold) {
                return false;
            }
            add_block(key, x);
            MSG_DEBUG("New block " << key << " has an f-test pvalue of " << pvalue(0, 0) << " >= " << m_threshold);
            return true;
        }

    void
        add_blocks(const model_block_collection& mbc)
        {
            m_computed = false;
            for (auto kv: mbc) {
                if (!kv.second) {
                    MSG_ERROR("Attempting to insert a null block in the model.", "");
                    abort();
                }
                if (!kv.first->is_valid()) {
                    MSG_ERROR("Attempting to insert an invalid block key in the model.", "");
                    abort();
                }
            }
            m_blocks.insert(mbc.begin(), mbc.end());
        }

    void
        remove_block(const value<model_block_type>& x)
        {
            auto i = m_blocks.begin();
            auto j = m_blocks.end();
            while (i != j && i->second != x) ++i;
            if (i != j) {
                m_computed = false;
                std::vector<model_block_collection::iterator> to_remove;
                /*auto mbk = i->first;*/
                /*remove_dependencies(mbk);*/
                for (auto b: reverse_dependencies(to_remove, i->first)) {
                    m_blocks.erase(b);
                }
            }
            /*auto it = m_blocks.erase(std::find(m_blocks.begin(), m_blocks.end(), x));*/
            /*m_keys.erase(m_keys.begin() + (it - m_blocks.begin()));*/
        }

    /* FIXME need dependency graph for block keys */
    void
        remove_block(const model_block_key& x)
        {
            std::vector<model_block_collection::iterator> to_remove;
            m_computed = false;
            /*remove_dependencies(x);*/
            for (auto b: reverse_dependencies(to_remove, x)) {
                m_blocks.erase(b);
            }
            /*m_blocks.erase(x);*/
            /*auto it = m_keys.erase(std::find(m_keys.begin(), m_keys.end(), x));*/
            /*m_blocks.erase(m_blocks.begin() + (it - m_keys.begin()));*/
        }

    const std::vector<model_block_collection::iterator>&
        reverse_dependencies(std::vector<model_block_collection::iterator>& to_remove, const model_block_key& x)
        {
            /*std::vector<model_block_key> to_remove;*/
            for (auto i = m_blocks.begin(), j = m_blocks.end(); i != j; ++i) {
                if (/*(i->first->type == mbk_Dominance || i->first->type == mbk_Interaction) &&*/ i->first->has(x)) {
                    to_remove.push_back(i);
                }
            }
            /*for (const auto& mbk: to_remove) {*/
                /*remove_block(mbk);*/
            /*}*/
            return to_remove;
        }

    void
        remove_dependencies(const model_block_key& mbk)
        {
#if 0
            auto i = m_blocks.begin(), j = m_blocks.end();
            for (; i != j;) {
                if (i->first->contains(mbk)) {
                    auto tmp = i;
                    ++tmp;
                    m_blocks.erase(i);
                    i = tmp;
                } else {
                    ++i;
                }
            }
#else
            switch (mbk->type) {
                case mbk_Interaction:
                    remove_dependencies(mbk->right);
                case mbk_Dominance:
                    remove_dependencies(mbk->left);
                case mbk_POP:
                case mbk_Covar:
                case mbk_CI:
                    ;
            };
            auto it = m_blocks.find(mbk);
            if (it != m_blocks.end()) {
                m_blocks.erase(it);
            }
#endif
        }

    void enable_constraints()
    {
        if (!m_with_constraints) {
            m_with_constraints = true;
            m_computed = false;
        }
    }

    void disable_constraints()
    {
        if (m_with_constraints) {
            m_with_constraints = false;
            m_computed = false;
        }
    }

    void ghost_constraint(const chromosome* chr, double l1, double l2)
    {
        /*model_block_key k1, k2;*/
        /*k1 += std::make_pair(chr, l1);*/
        /*k2 += std::make_pair(chr, l2);*/
        locus_key lk, k1, k2;
        k1 = lk + l1;
        k2 = lk + l2;
        m_ghost_constraint.first = model_block_key_struc::pop(chr, k1);
        m_ghost_constraint.second = model_block_key_struc::pop(chr, k2);
    }

    SolverType solver_type() const
    {
        return m_solver_type;
    }

    void compute()
    {
        compute_impl();
        if (m_with_constraints && rank() < m_X->cols()) {
            MSG_WARNING("Insufficient rank in model matrix. Attempting to fix this automatically.");
            m_rank_fixer.enable(*this);
            bool keep;
//             MSG_DEBUG(X());
//             MSG_DEBUG("rank: " << rank());
            do {
                auto old_rank = rank();
                m_computed = false;
                compute_impl();
//                 MSG_DEBUG(X());
//                 MSG_DEBUG("rank: " << rank());
                keep = old_rank < rank();
            } while (rank() < m_X->cols() && m_rank_fixer.next(keep));
            m_rank_fixer.disable();
        }
    }
    
	void compute_impl()
	{
//        MSG_DEBUG('[' << std::this_thread::get_id() << "] Recomputing model with selection " << keys());
        if (m_computed) {
            return;
        }

        {
            auto i = m_blocks.begin(), j = m_blocks.end();
            for (; i != j;) {
                if (i->first->type == mbk_POP && i->first->loci->depth() == 0) {
                    auto tmp = i;
                    ++i;
                    m_blocks.erase(tmp);
                } else {
                    ++i;
                }
            }
        }
//        MSG_DEBUG('[' << std::this_thread::get_id() << "] Filtered selection " << keys());

		m_computed = true;
		m_X = new immediate_value<MatrixXd>(MatrixXd());
		m_coefficients = new immediate_value<MatrixXd>(MatrixXd());
		m_residuals = new immediate_value<MatrixXd>(MatrixXd());
		m_rss = new immediate_value<VectorXd>(VectorXd());

        if (m_with_constraints) {
            constraint_list constraints;
            size_t n_rows;
            size_t n_pop_rows = 0;
            size_t n_cols = 0;
            /* establish constraints */
            for (const auto& b: m_blocks) {
                auto tmp = compute_constraint(b.first, *b.second);
                constraints.insert(constraints.end(), tmp.begin(), tmp.end());
                n_cols += b.second->outerSize();
                n_pop_rows = b.second->innerSize();
            }
            /* compute total size */
            n_rows = n_pop_rows;
            for (const auto& cmap: constraints) {
                int block_rows = 0;
                for (const auto& kv: cmap) {
                    if (kv.second.innerSize() > block_rows) {
                        block_rows = kv.second.innerSize();
                    }
                }
                n_rows += block_rows;
            }
            /**m_X = MatrixXd(n_rows, n_cols);*/
            *m_X = MatrixXd::Zero(n_rows, n_cols);
            /* now fill the matrix */
            n_rows = 0;
            n_cols = 0;
            for (const auto& b: m_blocks) {
                const model_block_type& mb = *b.second;
                m_X->block(n_rows, n_cols, mb.innerSize(), mb.outerSize()) = mb.data;
//                MSG_DEBUG("DEBUG X" << std::endl << (*m_X));
//                MSG_QUEUE_FLUSH();
                n_cols += mb.outerSize();
            }
            n_rows = n_pop_rows;
//            MSG_DEBUG("n_rows=" << n_rows << " n_cols=" << n_cols);
//            MSG_DEBUG("DEBUG X" << std::endl << (*m_X));
            if (0) {
                /*os << "<model Y(" << m.m_Y->innerSize() << ',' << m.m_Y->outerSize() << "), " << m.m_blocks.size() << " blocks>";*/
                MatrixXd big = MatrixXd::Zero(1, X().cols());
                model_print::matrix_with_sections<std::string, void, std::string, std::vector<char>> mws(big);
                /*model_print::matrix_with_sections<std::string, void, model_block_key, void> mws(m.X());*/
                set_column_sections(mws);
                MSG_DEBUG("DEBUG X" << std::endl << mws);
            }
            for (const auto& cmap: constraints) {
                n_cols = 0;
                size_t row_incr = 0;
//                MSG_DEBUG("row=" << n_rows << " col=" << n_cols);
//                MSG_QUEUE_FLUSH();
                for (const auto& b: m_blocks) {
                    auto it = cmap.find(b.first);
                    if (it != cmap.end()) {
                        const MatrixXd& constraint = it->second;
//                        MSG_DEBUG('[' << std::this_thread::get_id() << "] Adding (" << constraint.innerSize() << ',' << constraint.outerSize() << ") constraint at " << n_rows << "," << n_cols << " within (" << m_X->rows() << ',' << m_X->cols() << ')');
//                        MSG_DEBUG(constraint);