joint locus block reduction is BROKEN.

Makefile

-include Makefile.conf
-
-
-all: spell-qtl spell-pedigree spell-marker
-
-depend:
-	+cd src && make .depend
-	+cd src/bayes/ && make .depend
-	+cd src/pedigree/ && make .depend
-
-spell-qtl: depend
-	+cd src && $(MAKE)
-
-spell-pedigree: depend
-	+cd src/pedigree && $(MAKE)
-
-spell-marker: depend
-	+cd src/bayes && $(MAKE)
-
-clean:
-	+cd src && make clean
-	+cd src/bayes && make clean
-	+cd src/pedigree && make clean
-
-
-install: spell-qtl spell-pedigree spell-marker
-	cp -v src/spell-qtl $(INSTALL_DIR)/bin
-	cp -v src/bayes/spell-marker $(INSTALL_DIR)/bin
-	cp -v src/pedigree/spell-pedigree $(INSTALL_DIR)/bin
-

cmake_install.cmake

-# Install script for directory: /home/damien/devel/spell-qtl
-
-# Set the install prefix
-if(NOT DEFINED CMAKE_INSTALL_PREFIX)
-  set(CMAKE_INSTALL_PREFIX "/usr/local")
-endif()
-string(REGEX REPLACE "/$" "" CMAKE_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}")
-
-# Set the install configuration name.
-if(NOT DEFINED CMAKE_INSTALL_CONFIG_NAME)
-  if(BUILD_TYPE)
-    string(REGEX REPLACE "^[^A-Za-z0-9_]+" ""
-           CMAKE_INSTALL_CONFIG_NAME "${BUILD_TYPE}")
-  else()
-    set(CMAKE_INSTALL_CONFIG_NAME "Debug")
-  endif()
-  message(STATUS "Install configuration: \"${CMAKE_INSTALL_CONFIG_NAME}\"")
-endif()
-
-# Set the component getting installed.
-if(NOT CMAKE_INSTALL_COMPONENT)
-  if(COMPONENT)
-    message(STATUS "Install component: \"${COMPONENT}\"")
-    set(CMAKE_INSTALL_COMPONENT "${COMPONENT}")
-  else()
-    set(CMAKE_INSTALL_COMPONENT)
-  endif()
-endif()
-
-# Install shared libraries without execute permission?
-if(NOT DEFINED CMAKE_INSTALL_SO_NO_EXE)
-  set(CMAKE_INSTALL_SO_NO_EXE "1")
-endif()
-
-if(CMAKE_INSTALL_COMPONENT)
-  set(CMAKE_INSTALL_MANIFEST "install_manifest_${CMAKE_INSTALL_COMPONENT}.txt")
-else()
-  set(CMAKE_INSTALL_MANIFEST "install_manifest.txt")
-endif()
-
-string(REPLACE ";" "\n" CMAKE_INSTALL_MANIFEST_CONTENT
-       "${CMAKE_INSTALL_MANIFEST_FILES}")
-file(WRITE "/home/damien/devel/spell-qtl/${CMAKE_INSTALL_MANIFEST}"
-     "${CMAKE_INSTALL_MANIFEST_CONTENT}")

include/computations/frontends4.h

@@ -212,7 +212,7 @@ struct model_manager;
 struct test_result;
 
 void make_report(model_manager& mm, const std::vector<std::string>& cmdline);
-void report_computation(std::string trait_name, chromosome_value chr, std::vector<double>& loci, const computation_along_chromosome& cac, const test_result& tr);
+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);
 
 
@@ -453,9 +453,9 @@ struct search_lg_type {
     {
         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]);
+//         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;
         }
@@ -604,7 +604,8 @@ struct search_lg_type {
     void
     deselect(double position, const collection<population_value>& all_pops, value<model> M)
     {
-        position = match_locus(position);
+        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) {
@@ -617,13 +618,18 @@ struct search_lg_type {
 
         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 {
-                reduce(all_pops, position, it->second);
+                MSG_DEBUG("Need to reduce block (remaining key is " << key << ')');
+                reduce(all_pops, i_locus, it->second);
             }
         }
-        selection = selection - position;
+        for (auto& ls: local_selections) {
+            ls = ls - position;
+        }
+//         local_selections[i_locus] = local_selections[i_locus] - position;
         recompute_modes();
         /*deselect(position);*/
     }
@@ -671,29 +677,32 @@ struct search_lg_type {
     }
 
     void
-    reduce(const collection<population_value>& all_pops, double position, value<model_block_type>& vblock)
+    reduce(const collection<population_value>& all_pops, size_t i_locus, value<model_block_type>& vblock)
     {
-        position = match_locus(position);
+        double position = all_positions[i_locus];
+        auto local_sel = local_selections[i_locus] + position;
         /* FIXME what about dominance?? */
-        locus_key lk2 = selection - position;
+        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, selection->parent, *vblock, all_pops);
+        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 = selection->reduce(active_settings->parent_count_per_pop_per_chr
+            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->parent)->row_labels;
+            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->parent),
+                                                    as_value(lk2),
                                                     all_popl,
                                                     all_pops,
                                                     true)[0];
@@ -1093,7 +1102,7 @@ struct model_manager {
                         last_computation[chrom] = NULL;
                 };
 
-                report_computation(trait_name, chrom, si.active_loci, cac[chrom], si.local_max);
+                report_computation(*Mbase, trait_name, chrom, si.active_loci, cac[chrom], si.local_max);
             }
             return cac;
         }

include/computations/report.h

@@ -147,7 +147,7 @@ struct analysis_report : public analysis_report_computations {
     }
 
     analysis_report&
-    report_algo_selection(std::string chrom, double score, double locus, double threshold);
+    report_algo_selection(model& Mbase, std::string chrom, double score, double locus, double threshold);
 
     analysis_report&
     report_full_map(region_of_interest_type& R, point_of_interest_type& P)
@@ -192,7 +192,7 @@ struct analysis_report : public analysis_report_computations {
     }
     
     analysis_report&
-    report_computation(std::string chrom, std::string trait_name, std::vector<double>& loci, const computation_along_chromosome& cac, const test_result& tr);
+    report_computation(model& Mbase, std::string chrom, std::string trait_name, std::vector<double>& loci, const computation_along_chromosome& cac, const test_result& tr);
     analysis_report&
     report_model(const model& Mcurrent);
     analysis_report&

include/model/model.h

@@ -138,14 +138,15 @@ assemble_block_multipop(const collection<model_block_type>& pop_blocks)
 
 
 struct model_block_key_struc {
+    int n_col;
     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)
+    model_block_key_struc(int nc, key_type kt, const chromosome* c, const locus_key& lk, const model_block_key& c1, const model_block_key& c2)
+        : n_col(nc), type(kt), chr(c), loci(lk), left(c1), right(c2)
     {}
 
     size_t
@@ -167,7 +168,7 @@ struct model_block_key_struc {
     {
         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_);
+        return std::make_shared<model_block_key_struc>(0, mbk_CI, (const chromosome*) NULL, void_lk, void_, void_);
     }
 
     static
@@ -176,7 +177,7 @@ struct model_block_key_struc {
     {
         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_);
+        return std::make_shared<model_block_key_struc>(0, mbk_Covar, (const chromosome*) NULL, void_lk, void_, void_);
     }
 
     static
@@ -184,7 +185,7 @@ struct model_block_key_struc {
         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_);
+            return std::make_shared<model_block_key_struc>(0, mbk_POP, c, loci, void_, void_);
         }
 
     static
@@ -193,7 +194,7 @@ struct model_block_key_struc {
         {
             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_);
+            return std::make_shared<model_block_key_struc>(0, mbk_Dominance, (const chromosome*) NULL, void_lk, haplo, void_);
         }
 
     static
@@ -202,9 +203,9 @@ struct model_block_key_struc {
         {
             locus_key void_lk;
             if (l < r) {
-                return std::make_shared<model_block_key_struc>(mbk_Interaction, (const chromosome*) NULL, void_lk, l, r);
+                return std::make_shared<model_block_key_struc>(1, 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);
+                return std::make_shared<model_block_key_struc>(1, mbk_Interaction, (const chromosome*) NULL, void_lk, r, l);
             }
         }
 
@@ -425,6 +426,59 @@ struct model_block_key_struc {
             return false;
         }
 
+        
+#if 0
+    /* helper to create the reduction matrix when removing a locus from a complex model block. */
+    Eigen::MatrixXd reduce(std::function<int(double)> n_ancestors_at_locus, double loc_remove) const
+    {
+        MatrixXd ret;
+
+        int n_geno;
+
+        switch (type) {
+            case mbk_CI:
+            case mbk_Covar:
+                return {};  // should never happen.
+            case mbk_Dominance:
+                
+            case mbk_Interaction:
+            case mbk_POP:
+                if (locus == no_locus) {
+                    return MatrixXd::Ones(1, 1);
+                }
+
+                n_geno = n_ancestors_at_locus(locus);
+                
+                if (loc_remove == locus) {
+                    /*MSG_DEBUG("this locus");*/
+                    ret = MatrixXd::Ones(n_geno, 1);
+                } else {
+                    /*MSG_DEBUG("not this locus");*/
+                    ret = VectorXd::Ones(n_geno).asDiagonal();
+                }
+                if (parent) {
+                    return kroneckerProduct(
+        #ifdef LK_REDUCE_R2L
+                                ret
+                                ,
+        #endif
+                                parent->reduce(n_ancestors_at_locus, loc_remove)
+        #ifndef LK_REDUCE_R2L
+                                ,
+                                ret
+        #endif
+                                );
+                } else {
+                    return ret;
+                }
+        };
+        
+    }
+    
+#endif
+
+
+        
     friend
         inline
         bool
@@ -1311,17 +1365,17 @@ struct model {
         }
     }
 
-    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);
-    }
+//     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
     {

include/settings.h

@@ -286,7 +286,7 @@ struct settings_t {
         , covar_per_pop()
         , lod_support_inner(1.)
         , lod_support_outer(lod_support_inner * 2)
-        , Rjoin(10.)
+        , Rjoin(0.)
         , Rskip(1.)
         , report_format_string("text")
 //         , report_format_string("excel")

src/computations/model.cc

@@ -795,7 +795,7 @@ disassemble_parental_origins_multipop(
 		{
 			return pop.observed_traits.front().values.size();
 		};
-
+        
     std::function<std::vector<std::vector<char>>(const qtl_pop_type&)>
 		get_labels = [=](const qtl_pop_type& pop)
 		{

src/computations/report.cc

@@ -261,9 +261,9 @@ analysis_report_computations::contrast_matrices(const std::vector<contrast_descr
 
 
 analysis_report&
-analysis_report::report_computation(std::string chrom, std::string trait_name, std::vector<double>& loci, const computation_along_chromosome& cac, const test_result& tr)
+analysis_report::report_computation(model& Mbase, std::string chrom, std::string trait_name, std::vector<double>& loci, const computation_along_chromosome& cac, const test_result& tr)
 {
-    report_algo_selection(chrom, tr.test_value, tr.locus, mm.threshold);
+    report_algo_selection(Mbase, chrom, tr.test_value, tr.locus, mm.threshold);
     std::vector<std::string> colnames = {"locus"};
     MatrixXd values = MatrixXd::Zero(loci.size(),
                                      1
@@ -307,19 +307,19 @@ analysis_report::init_selection_report()
 
 
 analysis_report&
-analysis_report::report_algo_selection(std::string chrom, double score, double locus, double threshold)
+analysis_report::report_algo_selection(model& Mbase, std::string chrom, double score, double locus, double threshold)
 {
     init_selection_report();
-    mm.vMcurrent->compute();
+    Mbase.compute();
     ++m_selection_index;
     m_stream.mode(ReportMode::AlgoSel);
     select_book("algo_trace").select_sheet("Computations");
     m_stream
-    .selection(m_selection_index, chrom, MESSAGE(mm.vMcurrent->keys()), score, locus, threshold);
+    .selection(m_selection_index, chrom, MESSAGE(Mbase.keys()), score, locus, threshold);
 //         if (m_what & AR::Model) {
         select_book("algo_trace").select_sheet(MESSAGE(m_selection_index << " M0"));
         m_stream
-        .matrix("", mm.vMcurrent->printable_X());
+        .matrix("", Mbase.printable_X());
 //         }
     select_previous_book();
     return *this;
@@ -476,9 +476,9 @@ void make_report(model_manager& mm, const std::vector<std::string>& cmdline)
 }
 
 
-void report_computation(std::string trait_name, chromosome_value chr, std::vector<double>& loci, const computation_along_chromosome& cac, const test_result& tr)
+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)
 {
-    ar_instance().report_computation(chr->name, trait_name, loci, cac, tr);
+    ar_instance().report_computation(Mbase, chr->name, trait_name, loci, cac, tr);
 }