1st version of SMT/R. Had to rename MESSAGE into SPELL_STRING because of conflict.

include/bayes/cli.h

@@ -133,7 +133,7 @@ struct bn_settings_t {
             /*if (!unconst_this) {*/
                 /*MSG_ERROR("Couldn't unconst the pointer to settings.", "Sacrifice a goat and hope for the best.");*/
             /*}*/
-            std::string ret = MESSAGE(work_directory << '/' << job_name << '.' << job << ".data");
+            std::string ret = SPELL_STRING(work_directory << '/' << job_name << '.' << job << ".data");
             std::lock_guard<std::mutex> lock(/*unconst_this->*/job_filenames_mutex);
             /*unconst_this->*/job_filenames.insert(ret);
             return ret;

include/bayes/factor_var2.h

@@ -444,7 +444,7 @@ struct bayesian_network {
 
     size_t new_variable(size_t n_par, size_t n_al, bool observed=false)
     {
-        m_variable_names.emplace_back(MESSAGE(m_variables.size()));
+        m_variable_names.emplace_back(SPELL_STRING(m_variables.size()));
         size_t ret = m_variables.size();
         m_variables.emplace_back(ret, n_par, n_al, observed);
         return ret;
@@ -641,12 +641,12 @@ struct bayesian_network {
     std::string debug_fac(size_t f) const
     {
         const auto& F = m_factors[f];
-        std::string msg1 = MESSAGE(F.type() << '[' << debug_var(F.var1()) << ", " << debug_var(F.var2()));
+        std::string msg1 = SPELL_STRING(F.type() << '[' << debug_var(F.var1()) << ", " << debug_var(F.var2()));
         std::string msg2;
         if (F.var3() != novar) {
-            msg2 = MESSAGE(", " << debug_var(F.var3()));
+            msg2 = SPELL_STRING(", " << debug_var(F.var3()));
         }
-        return MESSAGE(msg1 << msg2 << ']');
+        return SPELL_STRING(msg1 << msg2 << ']');
     }
 
     std::string debug_message(MessageType t, size_t i, bool full=false, const double_buffer_type& msg_buf_table={}) const
@@ -655,17 +655,17 @@ struct bayesian_network {
         switch (t) {
             case MT_F2V:
                 std::tie(f, v) = m_corresp_F2V[i];
-                return MESSAGE(debug_fac(f) << " --> " << debug_var(v));
+                return SPELL_STRING(debug_fac(f) << " --> " << debug_var(v));
             case MT_V2F:
                 std::tie(v, f) = m_corresp_V2F[i];
                 if (full) {
                     std::vector<std::string> inc;
                     for (const auto& b: m_var_messages[i].m_incoming) {
-                        inc.emplace_back(MESSAGE(std::endl << msg_buf_table.back(b).transpose()));
+                        inc.emplace_back(SPELL_STRING(std::endl << msg_buf_table.back(b).transpose()));
                     }
-                    return MESSAGE(debug_var(v) << " --> " << debug_fac(f) << inc);
+                    return SPELL_STRING(debug_var(v) << " --> " << debug_fac(f) << inc);
                 } else {
-                    return MESSAGE(debug_var(v) << " --> " << debug_fac(f));
+                    return SPELL_STRING(debug_var(v) << " --> " << debug_fac(f));
                 }
         };
         return "gcc is dumb";
@@ -1023,7 +1023,7 @@ struct bayesian_network {
             /*return m_var_messages[m_evidence_table[v]].m_observation;*/
             auto it = m_bn->m_evidence_table.find(v);
             if (it == m_bn->m_evidence_table.end()) {
-                throw std::runtime_error(MESSAGE("Variable " << m_bn->debug_var(v) << " is not observed"));
+                throw std::runtime_error(SPELL_STRING("Variable " << m_bn->debug_var(v) << " is not observed"));
             }
             return m_observation[it->second];
         }

include/bayes/factor_var3.h

@@ -196,7 +196,7 @@ struct bn_message_type {
                 kv.second = idx++;
             }
             for (const auto& kv: m_map) {
-                std::string str = MESSAGE(kv.second);
+                std::string str = SPELL_STRING(kv.second);
                 if (str.size() > width) {
                     width = str.size();
                 }
@@ -208,7 +208,7 @@ struct bn_message_type {
             for (const auto& kv: m_map) {
                 bn_label_type l = kv.first.keys.front().state;
                 size_t p = kv.first.keys.front().parent;
-                std::string str = MESSAGE(kv.second);
+                std::string str = SPELL_STRING(kv.second);
                 if (str.size() > width) {
                     width = str.size();
                 }
@@ -222,12 +222,12 @@ struct bn_message_type {
 
             header("");
             for (const auto& kv: var_index) {
-                cell(MESSAGE('#' << kv.first));
+                cell(SPELL_STRING('#' << kv.first));
             }
             out << std::endl;
             auto li = label_index.begin();
             for (size_t row = 0; row < label_index.size(); ++row) {
-                header(MESSAGE("" << li->first));
+                header(SPELL_STRING("" << li->first));
                 for (size_t col = 0; col < var_index.size(); ++col) {
                     cell(values[col][row]);
                 }
@@ -1530,7 +1530,7 @@ struct factor_graph {
             for (const auto& descr_i: message_map) {
                 MSG_DEBUG("setting descr " << descr_i.second << " with " << descr_i.first);
                 MSG_QUEUE_FLUSH();
-                m_message_headers[descr_i.second] = MESSAGE("" << descr_i.first);
+                m_message_headers[descr_i.second] = SPELL_STRING("" << descr_i.first);
             }
 
             for (const auto& upd: m_updaters) {

include/bayes/factor_var4.h

@@ -220,7 +220,7 @@ struct bn_message_type {
                 kv.second = idx++;
             }
             for (const auto& kv: m_map) {
-                std::string str = MESSAGE(kv.second);
+                std::string str = SPELL_STRING(kv.second);
                 if (str.size() > width) {
                     width = str.size();
                 }
@@ -232,7 +232,7 @@ struct bn_message_type {
             for (const auto& kv: m_map) {
                 bn_label_type l = kv.first.keys.front().state;
                 size_t p = kv.first.keys.front().parent;
-                std::string str = MESSAGE(kv.second);
+                std::string str = SPELL_STRING(kv.second);
                 if (str.size() > width) {
                     width = str.size();
                 }
@@ -246,12 +246,12 @@ struct bn_message_type {
 
             header("");
             for (const auto& kv: var_index) {
-                cell(MESSAGE('#' << kv.first));
+                cell(SPELL_STRING('#' << kv.first));
             }
             out << std::endl;
             auto li = label_index.begin();
             for (size_t row = 0; row < label_index.size(); ++row) {
-                header(MESSAGE("" << li->first));
+                header(SPELL_STRING("" << li->first));
                 for (size_t col = 0; col < var_index.size(); ++col) {
                     cell(values[col][row]);
                 }
@@ -2244,7 +2244,7 @@ struct factor_graph {
             for (const auto& iptr: m_interfaces) {
                 const auto& fvec = iptr->factors();
                 for (const auto& fptr: fvec) {
-                    m_message_headers[message_map[voidpair{iptr.get(), fptr.get()}]] = MESSAGE("i{" << iptr->variables() << "} -> F{" << fptr->variables() << '}');
+                    m_message_headers[message_map[voidpair{iptr.get(), fptr.get()}]] = SPELL_STRING("i{" << iptr->variables() << "} -> F{" << fptr->variables() << '}');
                     /*MSG_DEBUG("MESSAGE " << (m_message_headers[message_map[voidpair{iptr.get(), fptr.get()}]]));*/
                     m_updaters.emplace_back(std::static_pointer_cast<bn_message_computer_type>(iptr), message_map[voidpair{iptr.get(), fptr.get()}], iptr->variables());
                     updater_by_message_index[m_updaters.back().output_message()] = m_updaters.size() - 1;
@@ -2258,7 +2258,7 @@ struct factor_graph {
             for (const auto& fptr: m_factors) {
                 const auto& ivec = fptr->interfaces();
                 for (const auto& iptr: ivec) {
-                    m_message_headers[message_map[voidpair{fptr.get(), iptr.get()}]] = MESSAGE("F{" << fptr->variables() << "} -> i{" << iptr->variables() << '}');
+                    m_message_headers[message_map[voidpair{fptr.get(), iptr.get()}]] = SPELL_STRING("F{" << fptr->variables() << "} -> i{" << iptr->variables() << '}');
                     /*MSG_DEBUG("MESSAGE " << (m_message_headers[message_map[voidpair{fptr.get(), iptr.get()}]]));*/
                     m_updaters.emplace_back(std::static_pointer_cast<bn_message_computer_type>(fptr), message_map[voidpair{fptr.get(), iptr.get()}], iptr->variables());
                     updater_by_message_index[m_updaters.back().output_message()] = m_updaters.size() - 1;

include/bayes/generalized_product.h

@@ -492,7 +492,7 @@ struct joint_variable_product_type {
             static std::vector<size_t> ret;
             ret.resize(domains.size());
             compute_coordinates(cursor, ret);
-            return MESSAGE(WHITE << ret << NORMAL);
+            return SPELL_STRING(WHITE << ret << NORMAL);
         }
 
     inline
@@ -1232,7 +1232,7 @@ inline
 genotype_comb_type
 generate_factor(const std::vector<int>& parents, bool dh, int spawnling, const genotype_comb_type& parent_domain, std::map<std::vector<int>, genotype_comb_type>& domains, bool keep_gametes)
 {
-    scoped_indent _(MESSAGE("[generate_factor " << (dh ? "dh " : "") << spawnling << " from parents " << parents << " domain " << parent_domain << "] "));
+    scoped_indent _(SPELL_STRING("[generate_factor " << (dh ? "dh " : "") << spawnling << " from parents " << parents << " domain " << parent_domain << "] "));
     std::vector<bn_label_type> label_g;
     std::vector<int> sorted_parents(parents);
     if (parents.size() == 2) {

include/bayes/graphnode.h

@@ -1734,7 +1734,7 @@ struct graph_type {
         to_image(const std::string& prefix, const std::string& format) const
         {
             auto A = active_nodes();
-            std::string dotfilename = MESSAGE("/tmp/prout" << prefix); //std::tempnam(nullptr);
+            std::string dotfilename = SPELL_STRING("/tmp/prout" << prefix); //std::tempnam(nullptr);
             std::vector<std::string> sub_filenames;
             {
                 std::ofstream dotfile(dotfilename);
@@ -1758,7 +1758,7 @@ struct graph_type {
                     if (is_interface(n)) {
                         dotfile << "shape=diamond, label=\"" << n  << " [" << inner_nodes[n] << "] " << variables_of(n) << "\"";
                     } else if (is_aggregate(n)) {
-                        std::string filename = subgraphs[n]->to_image(MESSAGE("_sub_" << prefix << '_' << n), format);
+                        std::string filename = subgraphs[n]->to_image(SPELL_STRING("_sub_" << prefix << '_' << n), format);
                         dotfile << "shape=box,labeljust=\"l\", labelloc=\"t\", label=\"" << n << " [" << inner_nodes[n] << "]\",image=\"" << filename << "\"";
                         sub_filenames.push_back(filename);
                     } else {
@@ -1773,8 +1773,8 @@ struct graph_type {
                 }
                 dotfile << '}' << std::endl;
             }
-            std::string ret = MESSAGE(prefix << '.' << format);
-            std::string cmd = MESSAGE("dot -T" << format << ' ' << dotfilename << " -o " << ret);
+            std::string ret = SPELL_STRING(prefix << '.' << format);
+            std::string cmd = SPELL_STRING("dot -T" << format << ' ' << dotfilename << " -o " << ret);
             if (system(cmd.c_str())) {
                 MSG_ERROR("Dot failed " << strerror(errno), "");
             }
@@ -2088,7 +2088,7 @@ struct graph_type {
     void
         compute_node_domain(node_index_type n)
         {
-            scoped_indent _(MESSAGE("[compute_node_domain #" << n << "] "));
+            scoped_indent _(SPELL_STRING("[compute_node_domain #" << n << "] "));
             auto inputs = nei_in(n);
             auto vv = variables_of(n);
             if (is_interface(n)) {
@@ -2109,7 +2109,7 @@ struct graph_type {
             } else if (is_aggregate(n)) {
                 subgraphs[n]->compute_domains_and_factors();
             } else /* factor */ {
-                scoped_indent _(MESSAGE("[compute_factor {" << rules[n] << "} => " << variables[n] << "] "));
+                scoped_indent _(SPELL_STRING("[compute_factor {" << rules[n] << "} => " << variables[n] << "] "));
                 MSG_DEBUG(" ### ### COMPUTING FACTOR {" << rules[n] << "} => " << variables[n] << " ### ###");
                 MSG_QUEUE_FLUSH();
                 joint_variable_product_type jvp;
@@ -2141,7 +2141,7 @@ struct graph_type {
                 genotype_comb_type jpar_dom;
                 MSG_QUEUE_FLUSH();
                 if (stack.size() > 1) {
-                    MSG_DEBUG("Domains " << (parent ? MESSAGE(parent << "->" << index_in_parent) : std::string("top-level")) << ' ' << domains);
+                    MSG_DEBUG("Domains " << (parent ? SPELL_STRING(parent << "->" << index_in_parent) : std::string("top-level")) << ' ' << domains);
                     jpar_dom = compute_product(stack.begin(), stack.end(), varset, domains);
                 } else {
                     jpar_dom = stack.front();
@@ -2166,7 +2166,7 @@ struct graph_type {
     void
         compute_domains_and_factors()
         {
-            scoped_indent _(MESSAGE("[compute_domains_and_factors " << (parent ? MESSAGE(parent << "->" << index_in_parent) : std::string("top-level")) << "] "));
+            scoped_indent _(SPELL_STRING("[compute_domains_and_factors " << (parent ? SPELL_STRING(parent << "->" << index_in_parent) : std::string("top-level")) << "] "));
             /*MSG_DEBUG("n_alleles = " << n_alleles);*/
             if (parent) {
                 domains = parent->domains;
@@ -2354,7 +2354,7 @@ struct graph_type {
                 }
             }
             /*tmp->finalize();*/
-            tmp->to_image(MESSAGE(pfx << '_' << this << '_' << std::setw(4) << std::setfill('0') << x), "png");
+            tmp->to_image(SPELL_STRING(pfx << '_' << this << '_' << std::setw(4) << std::setfill('0') << x), "png");
         }
 
     void
@@ -2527,7 +2527,7 @@ struct graph_type {
             : tables(g->node_domains.size()), evidence(g->node_domains.size()), sub_instances(g->node_domains.size()), parent(par)
         {
             var_vec varz; for (node_index_type n: g->active_nodes()) { varz = varz + g->variables_of(n); }
-            scoped_indent _(MESSAGE("[instance " << varz << "] "));
+            scoped_indent _(SPELL_STRING("[instance " << varz << "] "));
             size_t mi = 0;
             for (node_index_type n: g->active_nodes()) {
                 if (g->is_interface(n) || g->is_factor(n)) {

include/bayes/graphnode2.h

@@ -528,7 +528,7 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
         to_image(const std::string& prefix, const std::string& format) const
         {
             auto A = nodes();
-            std::string dotfilename = MESSAGE("/tmp/prout" << prefix); //std::tempnam(nullptr);
+            std::string dotfilename = SPELL_STRING("/tmp/prout" << prefix); //std::tempnam(nullptr);
             std::vector<std::string> sub_filenames;
             {
                 std::ofstream dotfile(dotfilename);
@@ -552,7 +552,7 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
                     if (node_is_interface(n)) {
                         dotfile << "shape=diamond, label=\"" << n  << " [" << inner_nodes(n) << "]   " << variables_of(n) << "\"";
                     } else if (node_is_subgraph(n)) {
-                        std::string filename = subgraph(n)->to_image(MESSAGE("_sub_" << prefix << '_' << n), format);
+                        std::string filename = subgraph(n)->to_image(SPELL_STRING("_sub_" << prefix << '_' << n), format);
                         dotfile << "shape=box,labeljust=\"l\", labelloc=\"t\", label=\"" << n << " [" << inner_nodes(n) << "]\",image=\"" << filename << "\"";
                         sub_filenames.push_back(filename);
                     } else {
@@ -567,8 +567,8 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
                 }
                 dotfile << '}' << std::endl;
             }
-            std::string ret = MESSAGE(prefix << '.' << format);
-            std::string cmd = MESSAGE("dot -T" << format << ' ' << dotfilename << " -o " << ret);
+            std::string ret = SPELL_STRING(prefix << '.' << format);
+            std::string cmd = SPELL_STRING("dot -T" << format << ' ' << dotfilename << " -o " << ret);
             if (system(cmd.c_str())) {
                 MSG_ERROR("Dot failed " << strerror(errno), "");
             }
@@ -603,7 +603,7 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
     void
         finalize()
         {
-            scoped_indent _(MESSAGE("[finalize] "));
+            scoped_indent _(SPELL_STRING("[finalize] "));
 //            to_image("pre-optim", "png");
             MSG_DEBUG("DOMAINS");
             compute_domains_and_factors();
@@ -1025,7 +1025,7 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
             } else if (node_is_subgraph(n)) {
                 subgraph(n)->compute_domains_and_factors();
             } else /* factor */ {
-                scoped_indent _(MESSAGE("[compute_factor {" << rule_of(n) << "} => " << own_variable_of(n) << "] "));
+                scoped_indent _(SPELL_STRING("[compute_factor {" << rule_of(n) << "} => " << own_variable_of(n) << "] "));
                 MSG_DEBUG(" ### ### COMPUTING FACTOR {" << rule_of(n) << "} => " << own_variable_of(n) << " ### ###");
                 MSG_QUEUE_FLUSH();
                 joint_variable_product_type jvp;
@@ -1059,7 +1059,7 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
                 genotype_comb_type jpar_dom;
                 /*MSG_QUEUE_FLUSH();*/
                 if (stack.size() > 1) {
-                    MSG_DEBUG("Domains " << (parent() ? MESSAGE(parent() << "->" << index_in_parent()) : std::string("top-level")) << ' ' << domains);
+                    MSG_DEBUG("Domains " << (parent() ? SPELL_STRING(parent() << "->" << index_in_parent()) : std::string("top-level")) << ' ' << domains);
                     jpar_dom = compute_product(stack.begin(), stack.end(), varset, domains);
                 } else {
                     jpar_dom = stack.front();
@@ -1513,7 +1513,7 @@ struct instance_type {
     message_type
         compute_message(const message_compute_operation_type& op, const std::map<var_vec, genotype_comb_type>& domains)
         {
-            scoped_indent _(MESSAGE("[compute message " << op.emitter << " % " << op.output << " -> " << op.receiver << "] "));
+            scoped_indent _(SPELL_STRING("[compute message " << op.emitter << " % " << op.output << " -> " << op.receiver << "] "));
             MSG_DEBUG("" << op);
             if (sub_instances[op.emitter]) {
                 return sub_instances[op.emitter]->compute(op.receiver, domains);
@@ -1560,7 +1560,7 @@ struct instance_type {
             }
             for (size_t o: variant.operations) {
                 const auto& op = operations[o];
-                scoped_indent _(MESSAGE("[" << get_message_index(op.emitter, op.receiver) << "] "));
+                scoped_indent _(SPELL_STRING("[" << get_message_index(op.emitter, op.receiver) << "] "));
                 messages[get_message_index(op.emitter, op.receiver)] = compute_message(op, domains);
                 MSG_DEBUG("MESSAGES");
                 for (const auto& kv: message_index) {

include/bayes/inet.h

@@ -69,7 +69,7 @@ struct inet_server {
     std::string
     address() const
     {
-        return MESSAGE(hostname << ':' << port);
+        return SPELL_STRING(hostname << ':' << port);
     }
 
     inet_server(uint16_t start_port, uint16_t end_port)
@@ -243,7 +243,8 @@ struct inet_client {
         if (connect(sock,  (struct sockaddr*) &name, sizeof(name)) < 0) {
             MSG_ERROR("Couldn't connect to server: " << strerror(errno), "");
         }
-        write(sock, s.c_str(), s.size());
+        int r = write(sock, s.c_str(), s.size());
+        (void) r;
         close(sock);
     }
 };

include/bayes/nd_matrix.h

@@ -71,7 +71,7 @@ struct nd_matrix_base {
     {
 #ifndef NDEBUG
         if (idx.size() != dim.size()) {
-            throw std::runtime_error(MESSAGE("Wrong number of dimensions (" << idx.size() << " != " << dim.size() << ')'));
+            throw std::runtime_error(SPELL_STRING("Wrong number of dimensions (" << idx.size() << " != " << dim.size() << ')'));
         }
 #endif
         size_t accum = 0;
@@ -79,7 +79,7 @@ struct nd_matrix_base {
         for (size_t i = 0; i < n_dim; ++i) {
 #ifndef NDEBUG
             if (idx[i] >= dim[i]) {
-                throw std::runtime_error(MESSAGE("Index #" << i << " greater than dimension (" << idx[i] << " > " << dim[i] << ')'));
+                throw std::runtime_error(SPELL_STRING("Index #" << i << " greater than dimension (" << idx[i] << " > " << dim[i] << ')'));
             }
 #endif
             accum += idx[i] * strides[i];
@@ -1014,20 +1014,20 @@ struct hadamard_product_type : public nd_matrix_base, public printable<hadamard_
     {
 #ifndef NDEBUG
         if (m1.dimensions().size() < m2.dimensions().size()) {
-            throw std::runtime_error(MESSAGE("Left matrix must have a greater dimensionality than right matrix"));
+            throw std::runtime_error(SPELL_STRING("Left matrix must have a greater dimensionality than right matrix"));
         }
         if (m2.dimensions().size() != common_dim.size()) {
-            throw std::runtime_error(MESSAGE("All dimensions of right matrix must be mapped to dimensions of the left matrix"));
+            throw std::runtime_error(SPELL_STRING("All dimensions of right matrix must be mapped to dimensions of the left matrix"));
         }
         for (const auto& dim: m_common_dim) {
             if (dim.first > m1.dimensions().size()) {
-                throw std::runtime_error(MESSAGE("Invalid dimension number for left matrix (" << dim.first << " > " << m1.dimensions().size() << ")"));
+                throw std::runtime_error(SPELL_STRING("Invalid dimension number for left matrix (" << dim.first << " > " << m1.dimensions().size() << ")"));
             }
             if (dim.second > m2.dimensions().size()) {
-                throw std::runtime_error(MESSAGE("Invalid dimension number for right matrix (" << dim.second << " > " << m2.dimensions().size() << ")"));
+                throw std::runtime_error(SPELL_STRING("Invalid dimension number for right matrix (" << dim.second << " > " << m2.dimensions().size() << ")"));
             }
             if (m1.dimension(dim.first) != m2.dimension(dim.second)) {
-                throw std::runtime_error(MESSAGE("Incompatible dimensions #" << dim.first << "[" << m1.dimension(dim.first) << "] and #" << dim.second << "[" << m1.dimension(dim.second) << "]"));
+                throw std::runtime_error(SPELL_STRING("Incompatible dimensions #" << dim.first << "[" << m1.dimension(dim.first) << "] and #" << dim.second << "[" << m1.dimension(dim.second) << "]"));
             }
         }
 #endif
@@ -1366,7 +1366,7 @@ void print_mat(std::ostream& os, const MAT_TYPE& mat)
     size_t width = 0;
     auto ac = mat.accessor();
     for (size_t i = 0; i < mat.size(); ++i) {
-        size_t tmp = MESSAGE(ac.get(i)).size();
+        size_t tmp = SPELL_STRING(ac.get(i)).size();
         if (tmp > width) {
             width = tmp;
         }

include/bayes/output.h

@@ -953,6 +953,14 @@ struct EM_helpers {
 };
 
 
+struct EM_map {
+    std::vector<std::string> marker_names;
+    std::vector<double> distances;
+    double likelihood;
+    int n_iterations;
+    bool converged;
+};
+
 
 struct gamete_LV_database : public EM_helpers {
     /* IND.ped / MARK => 1 or 2 gametes (.second will be empty when cross type is DH) */
@@ -1016,6 +1024,15 @@ struct gamete_LV_database : public EM_helpers {
     void get_TR(double dist, Eigen::Matrix2d& result) { result = get_TR_unary(dist); }
     void get_TR(double dist, Eigen::Matrix4d& result) { result = get_TR_binary(dist); }
 
+    void
+    init_2pt_tr_at_inf()
+    {
+        tr2.clear();
+        tr4.clear();
+        tr2.push_back(MatrixXd::Constant(2, 2, .5));
+        tr4.push_back(MatrixXd::Constant(4, 4, .25));
+    }
+    
     void
     add_gametes(const std::string& mark, int ind, const VectorXd& lv, bool dh)
     {
@@ -1320,7 +1337,7 @@ struct gamete_LV_database : public EM_helpers {
         }
         return accum;
     }
-
+    
     double
     compute_2pt_dist(size_t i)
     {
@@ -1376,15 +1393,8 @@ struct gamete_LV_database : public EM_helpers {
         MSG_DEBUG("distances old " << em_old_distances << " new " << em_distances << " => " << accum);
         return accum;
     }
-    
-    
-    struct EM_map {
-        std::vector<std::string> marker_names;
-        std::vector<double> distances;
-        double likelihood;
-        int n_iterations;
-    };
-    
+
+
     EM_map
     EM(const std::vector<std::string>& marker_names, double convergence_threshold = 1.e-12, int max_iterations = 100)
     {
@@ -1405,10 +1415,8 @@ struct gamete_LV_database : public EM_helpers {
         for (double& d: ret.distances) {
             d = -50. * log(1 - 2. * d);
         }
-        ret.likelihood = /* ? */ -1;
-        if (ret.n_iterations == max_iterations && delta > convergence_threshold) {
-            ret.n_iterations = -1;
-        }
+        ret.likelihood = map_likelihood(marker_names, em_distances);
+        ret.converged = delta < convergence_threshold;
         return ret;
     }
     

include/cache2.h

@@ -125,7 +125,7 @@ template <typename Ret, typename... Args>
                         std::stringstream ss;
                         ss << "INVOKING " << func_name << "(" << std::endl;
                         std::vector<std::string> avec;
-                        do_with_arg_pack(avec.push_back(MESSAGE(args)));
+                        do_with_arg_pack(avec.push_back(SPELL_STRING(args)));
                         for (size_t i = 0; i < avec.size() - 1; ++i) {
                             ss << avec[i] << ", ";
                         }
@@ -201,7 +201,7 @@ template <typename Ret, typename... Args>
                         std::stringstream ss;
                         ss << "INVOKING " << func_name << "(" << std::endl;
                         std::vector<std::string> avec;
-                        do_with_arg_pack(avec.push_back(MESSAGE(args)));
+                        do_with_arg_pack(avec.push_back(SPELL_STRING(args)));
                         for (size_t i = 0; i < avec.size() - 1; ++i) {
                             ss << avec[i] << ", ";
                         }
@@ -399,7 +399,7 @@ template <typename Ret, typename... Args>
                 std::stringstream ss;
                 ss << m_name << '(';
                 std::vector<std::string> avec;
-                do_with_arg_pack(avec.push_back(MESSAGE(args)));
+                do_with_arg_pack(avec.push_back(SPELL_STRING(args)));
                 for (size_t i = 0; i < avec.size() - 1; ++i) {
                     ss << avec[i] << ", ";
                 }
@@ -425,7 +425,7 @@ template <typename Ret, typename... Args>
 
         Ret operator () (Args... x)
         {
-            static disk_hashtable<Ret(Args...)> dht(MESSAGE(active_settings->work_directory << '/' << active_settings->name << ".cache/" << m_name));
+            static disk_hashtable<Ret(Args...)> dht(SPELL_STRING(active_settings->work_directory << '/' << active_settings->name << ".cache/" << m_name));
             return dht.get_or_compute(m_func, std::forward<Args>(x)...);
         }
 

include/commandline.h

@@ -119,7 +119,7 @@ struct default_value_display {
     {
         C s;
         std::stringstream ss;
-        ss << "Default: " << WHITE << MESSAGE(s.*x) << NORMAL;
+        ss << "Default: " << WHITE << SPELL_STRING(s.*x) << NORMAL;
         str = ss.str();
     }
 

include/computations/excel_report.h

@@ -48,7 +48,7 @@ namespace excel_report {
 
 
 inline excel::stream& operator << (excel::stream& s, const std::vector<char>& vc) { std::string str(vc.begin(), vc.end()); return s << str; }
-inline excel::stream& operator << (excel::stream& s, const model_block_key& mbk) { return s << MESSAGE(mbk); }
+inline excel::stream& operator << (excel::stream& s, const model_block_key& mbk) { return s << SPELL_STRING(mbk); }
 
 namespace excel {
 template <typename RS, typename RF, typename CS, typename CF, typename M>

include/computations/report.h

@@ -162,7 +162,7 @@ struct analysis_report : public analysis_report_computations {
     report_trait(std::string name, const MatrixXd& values)
     {
         m_stream.mode(ReportMode::Trait);
-        select_book("report").select_sheet(MESSAGE("trait_" << name));
+        select_book("report").select_sheet(SPELL_STRING("trait_" << name));
         m_stream.matrix(values);
         return *this;
     }

include/computations/report_output.h

@@ -101,11 +101,11 @@ qtl_descr_to_mat(const std::vector<qtl_description>& descr)
     int r = 0;
     for (const auto& d: descr) {
         mat(r, 0) = d.chromosome;
-        mat(r, 1) = MESSAGE(d.position);
-        mat(r, 2) = MESSAGE(d.inner_interval_start);
-        mat(r, 3) = MESSAGE(d.inner_interval_end);
-        mat(r, 4) = MESSAGE(d.outer_interval_start);
-        mat(r, 5) = MESSAGE(d.outer_interval_end);
+        mat(r, 1) = SPELL_STRING(d.position);
+        mat(r, 2) = SPELL_STRING(d.inner_interval_start);
+        mat(r, 3) = SPELL_STRING(d.inner_interval_end);
+        mat(r, 4) = SPELL_STRING(d.outer_interval_start);
+        mat(r, 5) = SPELL_STRING(d.outer_interval_end);
         ++r;
     }
 //         MSG_DEBUG("QTL REPORT" << std::endl << mat);
@@ -184,28 +184,28 @@ inline
 std::string
 make_report_path(const std::string& trait_name, const std::string& type, const std::string& book)
 {
-    return MESSAGE(active_settings->work_directory << '/' << active_settings->name << ".report/" << type << '/' << trait_name << '/' << book);
+    return SPELL_STRING(active_settings->work_directory << '/' << active_settings->name << ".report/" << type << '/' << trait_name << '/' << book);
 }
 
 inline
 std::string
 make_report_path(const std::string& trait_name, const std::string& type, const std::string& book, const std::string& sheet, const std::string& ext)
 {
-    return MESSAGE(make_report_path(trait_name, type, book) << '/' << sheet << '.' << ext);
+    return SPELL_STRING(make_report_path(trait_name, type, book) << '/' << sheet << '.' << ext);
 }
 
 inline
 std::string
 make_report_path(const std::string& trait_name, const std::string& type, const std::string& book, const std::string& sheet, const std::string& title, const std::string& ext)
 {
-    return MESSAGE(make_report_path(trait_name, type, book) << '/' << sheet << '/' << title << '.' << ext);
+    return SPELL_STRING(make_report_path(trait_name, type, book) << '/' << sheet << '/' << title << '.' << ext);
 }