genoprob_computer2.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_GENOPROB_COMPUTER_H_
#define _SPEL_GENOPROB_COMPUTER_H_

#include "geno_matrix.h"
#include "labelled_matrix.h"
#include "data/qtl_chromosome.h"


struct out_of_segment_exception : public std::exception {};


struct geno_prob_computer {
    const std::vector<double>* locus;
    const MatrixXd* LV;
    const geno_matrix* this_gen;
    std::unordered_map<double, MatrixXd> tr_cache;

    std::vector<MatrixXd> TR;
    MatrixXd init_kernel;
    /*MatrixXd redux;*/
    std::map<int, MatrixXd> left, right;
    MatrixXd redux;
    VectorXd inv_stat_dist;
    std::vector<label_type> unique_labels;

    geno_prob_computer()
        : locus(0), LV(0), this_gen(0), TR(), init_kernel(), left(), right(), redux(), inv_stat_dist()
    {}

    geno_prob_computer(const geno_prob_computer& tmp)
        : locus(tmp.locus)
        , LV(tmp.LV)
        , this_gen(tmp.this_gen)
        , left(std::map<int, MatrixXd>(tmp.left))
        , right(std::map<int, MatrixXd>(tmp.right))
        , redux(tmp.redux)
        , inv_stat_dist(tmp.inv_stat_dist)
        , unique_labels(tmp.unique_labels)
    {}

    geno_prob_computer(geno_prob_computer&& tmp)
        : locus(tmp.locus)
        , LV(tmp.LV)
        , this_gen(tmp.this_gen)
        , left(std::forward<std::map<int, MatrixXd>>(tmp.left))
        , right(std::forward<std::map<int, MatrixXd>>(tmp.right))
        , redux(tmp.redux)
        , inv_stat_dist(tmp.inv_stat_dist)
        , unique_labels(std::forward<std::vector<label_type>>(tmp.unique_labels))
    {}

    geno_prob_computer&
        operator = (geno_prob_computer&& tmp)
        {
            locus = tmp.locus;
            LV = tmp.LV;
            this_gen = tmp.this_gen;
            TR.swap(tmp.TR);
            left.swap(tmp.left);
            right.swap(tmp.right);
            redux = tmp.redux;
            inv_stat_dist = tmp.inv_stat_dist;
            unique_labels.swap(tmp.unique_labels);
            return *this;
        }

    geno_prob_computer&
        operator = (const geno_prob_computer& tmp)
        {
            locus = tmp.locus;
            LV = tmp.LV;
            this_gen = tmp.this_gen;
            TR = tmp.TR;
            left = tmp.left;
            right = tmp.right;
            redux = tmp.redux;
            inv_stat_dist = tmp.inv_stat_dist;
            unique_labels = tmp.unique_labels;
            return *this;
        }

    MatrixXd& compute_tr(double dist)
    {
        auto it = tr_cache.find(dist);
        if (it == tr_cache.end()) {
             auto& TR = tr_cache[dist] = this_gen->exp(dist);
//              if ((TR.array() > 1.).any() || (TR.array() < -1.e-13).any()) {
//                  double max = TR.array().abs().maxCoeff();
//                  double min = TR.minCoeff();
//                  MSG_DEBUG("BAD TR @" << dist << "cM min=" << min << " max=" << max << std::endl << TR);
//                  abort();
//              }
            double min = TR.minCoeff();
            if (min < 0) {
                TR.array() += min;
//                 double max = TR.array().abs().maxCoeff();
//                 MSG_DEBUG("BAD TR @" << dist << "cM min=" << min << " max=" << max << std::endl << TR);
//                 abort();
            }
            normalize_by_col(TR);
            return tr_cache[dist];
        }
        return it->second;
    }

    MatrixXd normalize_by_col(const MatrixXd& m)
    {
        MatrixXd ret = MatrixXd::Zero(m.innerSize(), m.outerSize());
        VectorXd sums = m.colwise().sum();
        for (int i = 0; i < m.outerSize(); ++i) {
            if (sums(i) != 0) {
                ret.col(i) = m.col(i) / sums(i);
            }
        }
        return ret;
    }

    MatrixXd lv(int i)
    {
        return LV->col(i).asDiagonal();
    }

#ifndef WITH_OVERLUMPING
        static
    MatrixXd
        make_collect(const std::set<subset>& lumping_partition, int n_cols)
        {
            /* Creation of L1 */
            MatrixXd C = MatrixXd::Zero(lumping_partition.size(), n_cols);
            int i = 0;
            for (const auto& part: lumping_partition) {
                for (int s: part) {
                    C(i, s) = 1.;
                }
                ++i;
            }
            return C;
        }
#endif

    
    void init(const geno_matrix* gen)
    {
        /*tr_cache.clear();*/
        this_gen = gen;
        /*redux = make_redux(matrix_labels_to_cliques(break_matrix.column_labels));*/
#ifdef WITH_OVERLUMPING
        experimental_lumper el(*this_gen);
        redux = el.make_collect(el.partition_on_labels(), this_gen->inf_mat.cols());
#else
        lumper<MatrixXd, label_type> el(this_gen->inf_mat, this_gen->labels);
        auto pol = el.partition_on_labels();
        redux = make_collect({pol.begin(), pol.end()}, this_gen->inf_mat.cols());
#endif
        inv_stat_dist = this_gen->stat_dist;
        for (int i = 0; i < inv_stat_dist.size(); ++i) {
            inv_stat_dist(i) = 1. / inv_stat_dist(i);
        }
        unique_labels = this_gen->get_unique_labels();
#if 0
        unique_labels.reserve(this_gen->labels.size());
        std::set<label_type> latemp;  //(this_gen->labels.begin(), this_gen->labels.end());
        for (const auto& l: this_gen->labels) {
            if (latemp.find(l) == latemp.end()) {
                latemp.insert(l);
                unique_labels.push_back(l);
            }
        }
        /*unique_labels.assign(latemp.begin(), latemp.end());*/
#endif

        TR.clear();
        for (size_t i = 1; i < locus->size(); ++i) {
            TR.emplace_back(compute_tr((*locus)[i] - (*locus)[i - 1]));
        }
    }

    bool compute_partial_kernels()
    {
        std::map<int, std::string> lstring, rstring;
        left.clear();
        right.clear();

        int l = 0, r = LV->outerSize() - 1;

        MatrixXd kernel;

        left[0] = kernel = (this_gen->stat_dist.array() * LV->col(0).array()).matrix();

        lstring[0] = "LV0.";
        /*MSG_DEBUG("0 " << lstring[0]);*/
        /*MSG_DEBUG(left[0]);*/
        for(++l; l <= r; ++l) {
            MatrixXd old_kernel = kernel;
            kernel = LV->col(l).asDiagonal() * TR[l - 1] * kernel;
            double s = kernel.sum();
            if (s == 0) {
                return false;
            }
            if (s != s) {
                MSG_ERROR("#" << l << " LEFT KERNEL IS NaN", "");
                MSG_DEBUG("LV->col(" << l << ") =" << std::endl << LV->col(l).transpose());
                MSG_DEBUG("TR[" << (l - 1) << "] =" << std::endl << TR[l - 1]);
                MSG_DEBUG("old kernel =" << std::endl << old_kernel);
                abort();
            }
            kernel /= s;
            left[l] = kernel;
            std::stringstream ss;
            ss << "LV" << l << ".TR" << (l-1) << '.' << lstring[l-1];
            lstring[l] = ss.str();
            /*MSG_DEBUG(l << " " << lstring[l]);*/
            /*MSG_DEBUG(left[l]);*/
        }

        right[r+1] = VectorXd::Ones(this_gen->rows()) / this_gen->rows();
        /*right[r+1] = this_gen->stat_dist;*/
        rstring[r+1] = "";
        right[r-1] = right[r] = kernel = (this_gen->stat_dist.array() * LV->col(r).array()).matrix(); //LV->col(r) / LV->col(r).sum();
        { std::stringstream s; s << ".LV" << r; rstring[r-1] = rstring[r] = s.str(); }
        for (--r; r > 0; --r) {
            MatrixXd old_kernel = kernel;
            kernel = LV->col(r).asDiagonal() * TR[r] * kernel;
            double s = kernel.sum();
            if (s == 0) {
                return false;
            }
            if (s != s) {
                MSG_ERROR("#" << r << " RIGHT KERNEL IS NaN", "");
                MSG_DEBUG("LV->col(" << l << ") =" << std::endl << LV->col(l).transpose());
                MSG_DEBUG("TR[" << (l - 1) << "] =" << std::endl << TR[l - 1]);
                MSG_DEBUG("old kernel =" << std::endl << old_kernel);
                abort();
            }
            kernel /= s;
            right[r-1] = kernel;
            std::stringstream ss;
            ss << rstring[r] << ".TR" << r << ".LV" << r;
            rstring[r-1] = ss.str();
        }

        r = LV->outerSize() - 1;
        kernel = init_kernel;
        return true;
    }

    std::pair<size_t, size_t> find_lr(double loc)
    {
        size_t l = 0, r = LV->outerSize() - 1;
        if (loc < locus->front() || loc > locus->back()) {
            throw out_of_segment_exception();
        }

        for (; (l + 1) < r && (*locus)[l + 1] <= loc; ++l);
        for (; r > 0 && (*locus)[r - 1] > loc; --r);

        return {l, r};
    }

    VectorXd fast_compute_at(double loc, size_t l, size_t r)
    {
#if 0
        MSG_DEBUG("l=" << l);
        MSG_DEBUG("r=" << r);
        MSG_DEBUG("locus_l=" << (*locus)[l]);
        MSG_DEBUG("locus_r=" << (*locus)[r]);
        MSG_DEBUG("locus=" << loc);
#endif
        MatrixXd& TRL = compute_tr(loc - (*locus)[l]);
        MatrixXd& TRR = compute_tr((*locus)[r] - loc);
#if 0
        MSG_DEBUG(MATRIX_SIZE(TRL));
        MSG_DEBUG(TRL);
        MSG_DEBUG(MATRIX_SIZE(TRR));
        MSG_DEBUG(TRR);
        MSG_DEBUG(MATRIX_SIZE(left[l]));
        MSG_DEBUG(left[l]);
        MSG_DEBUG(MATRIX_SIZE(right[l]));
        MSG_DEBUG(right[l]);
#endif
        auto L = (TRL * left[l]).array();
        auto R = (TRR * right[l]).array();
#if 0
        /*std::cout << loc << std::endl << (L * R).matrix().transpose() << std::endl;*/
        MSG_DEBUG(MATRIX_SIZE(L));
        MSG_DEBUG(L);
        MSG_DEBUG(MATRIX_SIZE(R));
        MSG_DEBUG(R);
        MSG_DEBUG(MATRIX_SIZE(redux));
        MSG_DEBUG(redux);
#endif
        VectorXd vv = redux * (L * R * inv_stat_dist.array()).matrix();
        double s = vv.sum();
        if (s != s) {
            MSG_ERROR("geno_prob_computer spotted a nan!", "");
            MSG_DEBUG(vv);
            abort();
        }
        if (s == 0) {
            return vv;
        }
        return vv / s;
    }

    MatrixXd
        fast_compute_over_segment_raw(const std::vector<double>& steps)
        {
            MatrixXd ret(unique_labels.size(), steps.size());
            MSG_DEBUG(MATRIX_SIZE(ret));
            size_t l = 0, last = steps.size() - 1, r = !!(LV->outerSize() - 1);
            for (size_t i = 0; i < steps.size(); ++i) {
                /* FIXME: optimize the CORRECT initialization of l & r */
                std::tie(l, r) = find_lr(steps[i]);
                ret.col(i) = fast_compute_at(steps[i], l, r);
                for (int k = 0; k < ret.rows(); ++k) {
                    if (ret(k, i) < FLOAT_TOL && ret(k, i) > -FLOAT_TOL) {
                        ret(k, i) = 0;
                    }
                }
                if ((ret.col(i).array() < 0).any()) {
                    MSG_DEBUG("Bad value at column #" << i << std::endl << ret.col(i).transpose());
//                     MSG_DEBUG((*this_gen));
#if 0
                    MSG_DEBUG("steps: " << steps);
                    MSG_DEBUG("LV:" << std::endl << (*LV));
                    MSG_DEBUG(*this_gen);
                    for (const auto& m: TR) {
                        MSG_DEBUG(m);
                    }
                    for (const auto& kv: left) {
                        MSG_DEBUG("left " << kv.first);
                        MSG_DEBUG(kv.second);
                    }
                    for (const auto& kv: right) {
                        MSG_DEBUG("right " << kv.first);
                        MSG_DEBUG(kv.second);
                    }
#endif
                    /*throw 0;*/
                }
                while (steps[i] >= (*locus)[r]) {
                    ++l;
                    if (i < last) {
                        ++r;
                    } else {
                        break;
                    }
                    /*r += i < last;*/
                }
            }
//             MSG_DEBUG("STATE PROBABILITIES" << std::endl << ret);
            return ret;
        }

    labelled_matrix<MatrixXd, label_type, double>
        fast_compute_over_segment(const std::vector<double>& steps)
        {
            if (!compute_partial_kernels()) {
                labelled_matrix<MatrixXd, label_type, double>
                    ret(unique_labels, steps);
                ret.data = MatrixXd::Zero(ret.row_labels.size(), ret.column_labels.size());
                return ret;
            }

            labelled_matrix<MatrixXd, label_type, double>
                ret(unique_labels, steps);
            ret.data = fast_compute_over_segment_raw(steps);
            return ret;
        }
};



struct segment_computer_t {
    const geno_matrix* g_this;
    qtl_chromosome chr;
    geno_prob_computer gpc_vec;
    /*std::vector<double> steps;*/
    std::vector<double> loci;
    double noise;

    segment_computer_t() : g_this(0), chr(), gpc_vec(), /*steps(),*/ noise(0) {}

    segment_computer_t(segment_computer_t&& sc)
        : g_this(sc.g_this), chr(sc.chr), gpc_vec(std::move(sc.gpc_vec)), /*steps(std::move(sc.steps)),*/ loci(std::move(sc.loci)), noise(sc.noise)
    {}


    segment_computer_t(const geno_matrix& g, const chromosome* pchr, const locus_key& lk, double /*step*/, double nz)
        : g_this(&g), chr(pchr, lk), gpc_vec(), noise(nz)
    {
        chr.copy_loci(loci);
        /*steps = compute_steps(loci, step);*/
        init();
    }

    segment_computer_t(const geno_matrix& g, const chromosome* pchr, const locus_key& lk, double nz/*, const std::vector<double>& loc_vec*/)
        : g_this(&g), chr(pchr, lk), gpc_vec(), noise(nz)
    {
        chr.copy_loci(loci);
        /*steps = loc_vec;*/
        init();
    }

    bool
        operator == (const segment_computer_t& sc) const { return g_this == sc.g_this /*&& steps == sc.steps*/ && noise == sc.noise; }

    void
        init()
        {
            /*MSG_DEBUG("init !");*/
            gpc_vec.locus = &loci;
            gpc_vec.init(g_this);
        }

    segment_computer_t&
        operator = (segment_computer_t&& sc)
        {
            /*MSG_DEBUG("OPERATOR = :" << __LINE__);*/
            g_this = sc.g_this;
            chr = sc.chr;
            chr.copy_loci(loci);
            /*steps.swap(sc.steps);*/
            /*gpc_vec.swap(sc.gpc_vec);*/
            gpc_vec = std::move(sc.gpc_vec);
            noise = sc.noise;
            init();
            return *this;
        }

    labelled_matrix<MatrixXd, label_type, double>
        compute(const MatrixXd& LV, const std::vector<double>& steps)
        {
            /*DUMP_FILE_LINE();*/
            auto unphased_LV = g_this->unphased_LV();
            auto joint_qtl_iterator = chr.qtl_state_iterator(unphased_LV);
            /*MatrixXd lv = g_this->noisy_LV(LV, noise);*/
            MatrixXd qtl_lv(LV.innerSize(), loci.size());
            /*MSG_DEBUG(MATRIX_SIZE(lv) << std::endl << lv.transpose());*/
            gpc_vec.LV = &qtl_lv;
            joint_qtl_iterator.init_expansion(LV, qtl_lv);
            /*MSG_DEBUG(MATRIX_SIZE(qtl_lv));*/
            /*MSG_DEBUG(qtl_lv.transpose());*/
            auto tmp = compute(steps);
            int n_states = tmp.innerSize();
            std::vector<label_type> row_labels;
            row_labels.reserve(n_states * joint_qtl_iterator.size());
            /*MSG_DEBUG("qtl_iterator.size = " << joint_qtl_iterator.size());*/
            for (int i = 0; i < joint_qtl_iterator.size(); ++i) {
                row_labels.insert(row_labels.end(), tmp.row_labels.begin(), tmp.row_labels.end());
            }
            labelled_matrix<MatrixXd, label_type, double> ret(row_labels, steps);
            /*MSG_DEBUG("ret.data$dim(" << ret.data.innerSize() << ',' << ret.data.outerSize() << ") tmp.data$dim(" << tmp.innerSize() << ',' << tmp.outerSize() << ") n_states=" << n_states);*/
            /*ret.data.topRows(n_states) = compute().data;*/
            ret.data.topRows(n_states) = tmp.data;
            /*ret.emplace_back(compute());*/
            int row_ofs = n_states;
            while (!joint_qtl_iterator.next()) {
                if (joint_qtl_iterator.update(LV, qtl_lv)) {
                    /*ret.emplace_back(compute());*/
                    /*MSG_DEBUG(MATRIX_SIZE(qtl_lv));*/
                    /*MSG_DEBUG(qtl_lv.transpose());*/
                    ret.data.middleRows(row_ofs, n_states) = compute(steps).data;
                } else {
                    ret.data.middleRows(row_ofs, n_states) = MatrixXd::Zero(n_states, ret.data.outerSize());
                }
                row_ofs += n_states;
            }
            return ret;
        }

    labelled_matrix<MatrixXd, label_type, double>
        compute(const std::vector<double>& steps)
        {
            return gpc_vec.fast_compute_over_segment(steps);
        }

    friend
        inline
        std::ostream&
        operator << (std::ostream& os, const segment_computer_t& sc)
        {
            return os << '<' << sc.g_this->name
                /*<< " n_steps=" << sc.steps.size()*/
                << " chr=" << sc.chr.chr->name
                << " noise=" << sc.noise << '>';
        }

};


struct optim_segment_computer_t {
    const geno_matrix* g_this;
    qtl_chromosome chr;
    std::vector<geno_prob_computer> gpc_vec;
    std::vector<bool> gpc_not_zero;
    std::vector<double> loci;
    double noise;
    std::vector<label_type> row_labels;
    std::vector<MatrixXd> lv_vec;
    size_t n_states;

    optim_segment_computer_t()
        : g_this(0), chr(), gpc_vec(), gpc_not_zero(), loci(), noise(), row_labels(), lv_vec(), n_states(0) {}

    optim_segment_computer_t(optim_segment_computer_t&& other)
        : g_this(other.g_this)
        , chr(other.chr)
        , gpc_vec(std::move(other.gpc_vec))
        , gpc_not_zero(std::move(other.gpc_not_zero))
        , loci(std::move(other.loci))
        , noise(other.noise)
        , row_labels(std::move(other.row_labels))
        , lv_vec(std::move(other.lv_vec))
        , n_states(other.n_states)
    {
        for (size_t i = 0; i < gpc_vec.size(); ++i) {
            gpc_vec[i].locus = &loci;
            gpc_vec[i].LV = &lv_vec[i];
        }
    }

    optim_segment_computer_t(const optim_segment_computer_t& other)
        : g_this(other.g_this)
        , chr(other.chr)
        , gpc_vec(other.gpc_vec)
        , gpc_not_zero(other.gpc_not_zero)
        , loci(other.loci)
        , noise(other.noise)
        , row_labels(other.row_labels)
        , lv_vec(other.lv_vec)
        , n_states(other.n_states)
    {
        for (size_t i = 0; i < gpc_vec.size(); ++i) {
            gpc_vec[i].locus = &loci;
            gpc_vec[i].LV = &lv_vec[i];
        }
        /*MSG_DEBUG("optim_segment_computer_t copy " << gpc_vec.size() << " parts; labels " << row_labels);*/
    }

    optim_segment_computer_t&
        operator = (const optim_segment_computer_t& other)
        {
            g_this = other.g_this;
            chr = other.chr;
            gpc_vec = other.gpc_vec;
            gpc_not_zero = other.gpc_not_zero;
            loci = other.loci;
            noise = other.noise;
            row_labels = other.row_labels;
            lv_vec = other.lv_vec;
            n_states = other.n_states;
            for (size_t i = 0; i < gpc_vec.size(); ++i) {
                gpc_vec[i].locus = &loci;
                gpc_vec[i].LV = &lv_vec[i];
            }
            /*MSG_DEBUG("optim_segment_computer_t = " << gpc_vec.size() << " parts; labels " << row_labels);*/
            return *this;
        }

    optim_segment_computer_t(const geno_matrix& g, const chromosome* pchr, const locus_key& lk, double nz, const MatrixXd& LV)
        : g_this(&g), chr(pchr, lk), gpc_vec(), gpc_not_zero(), loci(), noise(nz), row_labels(), lv_vec(), n_states(0)
    {
        auto unphased_LV = g_this->unphased_LV();
        chr.copy_loci(loci);
        MatrixXd qtl_lv(LV.innerSize(), loci.size());
        auto joint_qtl_iterator = chr.qtl_state_iterator(unphased_LV);
        bool lv_ok = joint_qtl_iterator.init_expansion(LV, qtl_lv);
        gpc_vec.reserve(joint_qtl_iterator.size());
        lv_vec.reserve(joint_qtl_iterator.size());
        {
            geno_prob_computer tmp;
            tmp.locus = &loci;
            tmp.init(g_this);
            n_states = tmp.unique_labels.size();
            row_labels.reserve(n_states * joint_qtl_iterator.size());
            for (int i = 0; i < joint_qtl_iterator.size(); ++i) {
                row_labels.insert(row_labels.end(), tmp.unique_labels.begin(), tmp.unique_labels.end());
            }
        }

        auto
            add_gpc = [&] (bool lv_ok)
            {
                if (lv_ok) {
                    lv_vec.push_back(qtl_lv);
                    gpc_vec.emplace_back();
                    gpc_vec.back().locus = &loci;
                    gpc_vec.back().LV = &lv_vec.back();
                    gpc_vec.back().init(g_this);
                    if (gpc_vec.back().compute_partial_kernels()) {
                        gpc_not_zero.push_back(true);
                        /*MSG_DEBUG("GPC::LV" << std::endl << lv_vec.back());*/
                    } else {
                        gpc_not_zero.push_back(false);
                        gpc_vec.pop_back();
                        gpc_vec.emplace_back();
                        lv_vec.pop_back();
                        /*MSG_DEBUG("GPC::LV none. Update was OK but partial kernel computation failed.");*/
                    }
                } else {
                    gpc_not_zero.push_back(false);
                    gpc_vec.emplace_back();
                    lv_vec.emplace_back();
                    /*MSG_DEBUG("GPC::LV none. Update failed.");*/
                }
            };

        add_gpc(lv_ok);

        for (; !joint_qtl_iterator.next(); ) {
            add_gpc(joint_qtl_iterator.update(LV, qtl_lv));
        }

        /*MSG_DEBUG("optim_segment_computer_t NEW " << gpc_vec.size() << " parts; labels " << row_labels);*/

        /*MSG_DEBUG("NEW GPC row_labels " << row_labels);*/
#if 0
        for (; !joint_qtl_iterator.next();) {
            lv_vec.push_back(qtl_lv);
            gpc_vec.emplace_back();
            gpc_vec.back().locus = &loci;
            gpc_vec.back().LV = &lv_vec.back();
            gpc_vec.back().init(g_this);
            if (joint_qtl_iterator.update(LV, lv_vec.back()) && (lv_vec.push_back(qtl_lv), gpc_vec.back().compute_partial_kernels())) {
                gpc_not_zero.push_back(true);
            } else {
                gpc_not_zero.push_back(false);
                gpc_vec.pop_back();
                gpc_vec.emplace_back();
                lv_vec.pop_back();
                lv_vec.emplace_back();
            }
            MSG_DEBUG("GPC::LV" << std::endl << lv_vec.back());
        }
#endif
    }

    labelled_matrix<MatrixXd, label_type, double>
        compute(const std::vector<double>& steps)
        {
            labelled_matrix<MatrixXd, label_type, double> ret(row_labels, steps);
            /*MSG_DEBUG("segment_computer::compute => " << MATRIX_SIZE(ret.data));*/
            auto nzi = gpc_not_zero.begin(), nzj = gpc_not_zero.end();
            auto gpci = gpc_vec.begin();
            int row_ofs = 0;
            for (; nzi != nzj; ++nzi, ++gpci, row_ofs += n_states) {
                if (*nzi) {
                    ret.data.middleRows(row_ofs, n_states) = gpci->fast_compute_over_segment_raw(steps);
                } else {
                    ret.data.middleRows(row_ofs, n_states) = MatrixXd::Zero(n_states, steps.size());
                }
            }
            return ret;
        }
};

#endif