Commit 60f0b923 authored by Damien Leroux's avatar Damien Leroux
Browse files

Late commit with fixes and updates. See the diff…

parent bd47a1a9
...@@ -28,6 +28,7 @@ ...@@ -28,6 +28,7 @@
#include <iostream> #include <iostream>
#include "file.h" #include "file.h"
#include <cstring> #include <cstring>
#include <type_traits>
#include "eigen.h" #include "eigen.h"
#include "error.h" #include "error.h"
...@@ -956,9 +957,11 @@ struct EM_helpers { ...@@ -956,9 +957,11 @@ struct EM_helpers {
struct EM_map { struct EM_map {
std::vector<std::string> marker_names; std::vector<std::string> marker_names;
std::vector<double> distances; std::vector<double> distances;
std::vector<double> r;
double likelihood; double likelihood;
int n_iterations; int n_iterations;
bool converged; bool converged;
double delta;
}; };
...@@ -1079,7 +1082,7 @@ struct gamete_LV_database : public EM_helpers { ...@@ -1079,7 +1082,7 @@ struct gamete_LV_database : public EM_helpers {
gam.get(*name_i++, tmp); gam.get(*name_i++, tmp);
// accum = tmp / tmp.sum(); // accum = tmp / tmp.sum();
// // accum = tmp / (tmp.array() > 0).template cast<double>().sum(); // // accum = tmp / (tmp.array() > 0).template cast<double>().sum();
accum = tmp * .25; accum = tmp/* * .25*/;
for (const auto& t: TR) { for (const auto& t: TR) {
gam.get(*name_i++, tmp); gam.get(*name_i++, tmp);
MSG_DEBUG("accum = " << accum.transpose() << " LV = " << tmp.transpose()); MSG_DEBUG("accum = " << accum.transpose() << " LV = " << tmp.transpose());
...@@ -1315,6 +1318,67 @@ struct gamete_LV_database : public EM_helpers { ...@@ -1315,6 +1318,67 @@ struct gamete_LV_database : public EM_helpers {
MSG_DEBUG(em_state_per_mark_per_ind); MSG_DEBUG(em_state_per_mark_per_ind);
} }
double
EM_forward_likelihood()
{
scoped_indent _("[EM forward likelihood] ");
double accum = 0;
for (size_t i = 0; i < em_chain_is_single.size(); ++i) {
if (em_chain_is_single[i]) {
size_t m = 0;
Eigen::Vector2d v2 = Eigen::Vector2d::Constant(1);
// forward.
for (; m < em_distances.size(); ++m) {
v2 = em_TR2[m] * (v2.array() * em_obs[m][i].unary.array()).matrix();
}
v2 = (v2.array() * em_obs[m][i].unary.array()).matrix();
accum += log(v2.sum());
} else {
size_t m = 0;
Eigen::Vector4d v4 = Eigen::Vector4d::Constant(1);
// forward.
for (; m < em_distances.size(); ++m) {
v4 = em_TR4[m] * (v4.array() * em_obs[m][i].binary.array()).matrix();
}
v4 = (v4.array() * em_obs[m][i].binary.array()).matrix();
accum += log(v4.sum());
}
}
return accum;
}
template <typename M>
double
L1(const M& m)
{
double d = m.sum();
return d < DBL_EPSILON ? 1. / d : 1.;
}
double
twoMarkerLikelihoodAtInf()
{
return twoMarkerLikelihood(0, .5);
auto o1 = em_obs[0].begin();
auto done = em_obs[0].end();
auto o2 = em_obs[1].begin();
auto sing = em_chain_is_single.begin();
double accum = 0;
for (; o1 != done; ++sing/*, ++m1, ++m2*/, ++o1, ++o2) {
if (*sing) {
// accum += log(o1->unary.transpose() * L1(o1->unary) * Eigen::Matrix2d::Constant(.5) * o2->unary);
// MSG_INFO("sing " << o1->unary.transpose() << " ; " << o2->unary.transpose());
accum += log(.5 * o1->unary.sum() * o2->unary.sum());
} else {
// accum += log(o1->binary.transpose() * L1(o1->binary) * Eigen::Matrix4d::Constant(.25) * o2->binary);
// MSG_INFO("doub " << o1->binary.transpose() << " ; " << o2->binary.transpose());
accum += log(.25 * o1->binary.sum() * o2->binary.sum());
}
}
return accum;
}
double double
twoMarkerLikelihood(size_t i, double r) twoMarkerLikelihood(size_t i, double r)
{ {
...@@ -1326,10 +1390,19 @@ struct gamete_LV_database : public EM_helpers { ...@@ -1326,10 +1390,19 @@ struct gamete_LV_database : public EM_helpers {
double accum = 0; double accum = 0;
for (; sing != done; ++sing, ++m1, ++m2, ++o1, ++o2) { for (; sing != done; ++sing, ++m1, ++m2, ++o1, ++o2) {
if (*sing) { if (*sing) {
accum += log(mult_and_norm(m1->unary, o1->unary).transpose() * get_TR_unary(r, true) * mult_and_norm(m2->unary, o2->unary)); // accum += log(mult_and_norm(m1->unary, o1->unary).transpose() * get_TR_unary(r, true) * mult_and_norm(m2->unary, o2->unary));
accum += log((2 * m1->unary.array() * o1->unary.array()).matrix().transpose()
* get_TR_unary(r, true)
* (2 * m2->unary.array() * o2->unary.array()).matrix());
// MSG_INFO("o1 " << o1->unary.transpose());
// MSG_INFO("o2 " << o2->unary.transpose());
} else { } else {
Eigen::Vector4d mm1 = mult_and_norm(m1->binary, o1->binary); // Eigen::Vector4d mm1 = mult_and_norm(m1->binary, o1->binary);
Eigen::Vector4d mm2 = mult_and_norm(m2->binary, o2->binary); // Eigen::Vector4d mm2 = mult_and_norm(m2->binary, o2->binary);
// MSG_INFO("o1 " << o1->binary.transpose());
// MSG_INFO("o2 " << o2->binary.transpose());
Eigen::Vector4d mm1 = 4 * (m1->binary.array() * o1->binary.array()).matrix();
Eigen::Vector4d mm2 = 4 * (m2->binary.array() * o2->binary.array()).matrix();
/*MSG_DEBUG("mm1 " << mm1.transpose());*/ /*MSG_DEBUG("mm1 " << mm1.transpose());*/
/*MSG_DEBUG("mm2 " << mm2.transpose());*/ /*MSG_DEBUG("mm2 " << mm2.transpose());*/
accum += log(mm1.transpose() * get_TR_binary(r, true) * mm2); accum += log(mm1.transpose() * get_TR_binary(r, true) * mm2);
...@@ -1337,7 +1410,7 @@ struct gamete_LV_database : public EM_helpers { ...@@ -1337,7 +1410,7 @@ struct gamete_LV_database : public EM_helpers {
} }
return accum; return accum;
} }
double double
compute_2pt_dist(size_t i) compute_2pt_dist(size_t i)
{ {
...@@ -1370,7 +1443,7 @@ struct gamete_LV_database : public EM_helpers { ...@@ -1370,7 +1443,7 @@ struct gamete_LV_database : public EM_helpers {
return ret; return ret;
//*/ //*/
} }
void void
EM_update_distances() EM_update_distances()
{ {
...@@ -1396,7 +1469,7 @@ struct gamete_LV_database : public EM_helpers { ...@@ -1396,7 +1469,7 @@ struct gamete_LV_database : public EM_helpers {
EM_map EM_map
EM(const std::vector<std::string>& marker_names, double convergence_threshold = 1.e-12, int max_iterations = 100) EM(const std::vector<std::string>& marker_names, bool dist_if_true_else_r = true, double convergence_threshold = 1.e-12, int max_iterations = 100)
{ {
double delta; double delta;
EM_map ret; EM_map ret;
...@@ -1411,12 +1484,25 @@ struct gamete_LV_database : public EM_helpers { ...@@ -1411,12 +1484,25 @@ struct gamete_LV_database : public EM_helpers {
delta = EM_delta(); delta = EM_delta();
++ret.n_iterations; ++ret.n_iterations;
} while (delta > convergence_threshold && ret.n_iterations < max_iterations); } while (delta > convergence_threshold && ret.n_iterations < max_iterations);
ret.r = em_distances;
ret.distances = em_distances; ret.distances = em_distances;
for (double& d: ret.distances) { if (dist_if_true_else_r) {
d = -50. * log(1 - 2. * d); for (double& d: ret.distances) {
d = -50. * log(1 - 2. * d);
}
} }
ret.likelihood = map_likelihood(marker_names, em_distances); // ret.likelihood = map_likelihood(marker_names, em_distances);
ret.likelihood = EM_forward_likelihood();
// auto sing = em_chain_is_single.begin();
// for (const auto& v: em_forward.back()) {
// if (*sing++) {
// ret.likelihood += log(v.unary.sum());
// } else {
// ret.likelihood += log(v.binary.sum());
// }
// }
ret.converged = delta < convergence_threshold; ret.converged = delta < convergence_threshold;
ret.delta = delta;
return ret; return ret;
} }
......
...@@ -571,7 +571,7 @@ void message_queue::run() ...@@ -571,7 +571,7 @@ void message_queue::run()
#ifdef SPELL_UNSAFE_OUTPUT #ifdef SPELL_UNSAFE_OUTPUT
#define MSG_ERROR(_msg_expr_, _workaround_expr_) do { std::cerr << _msg_expr_ << std::endl; } while (0) #define MSG_ERROR(_msg_expr_, _workaround_expr_) do { std::cerr << msg_handler_t::e() << _msg_expr_ << msg_handler_t::n() << std::endl; } while (0)
#define MSG_WARNING(_msg_expr_) do { std::cout << _msg_expr_ << std::endl; } while (0) #define MSG_WARNING(_msg_expr_) do { std::cout << _msg_expr_ << std::endl; } while (0)
#define MSG_INFO(_msg_expr_) do { std::cout << _msg_expr_ << std::endl; } while (0) #define MSG_INFO(_msg_expr_) do { std::cout << _msg_expr_ << std::endl; } while (0)
#define MSG_QUEUE_FLUSH() #define MSG_QUEUE_FLUSH()
...@@ -706,12 +706,18 @@ struct scoped_indent { ...@@ -706,12 +706,18 @@ struct scoped_indent {
inline bool msg_handler_t::state_t::check(bool fatal) inline bool msg_handler_t::state_t::check(bool fatal)
{ {
msg_handler_t::scoped_lock_type _(msg_handler_t::mutex); msg_handler_t::scoped_lock_type _(msg_handler_t::mutex);
// MSG_INFO("count = " << count);
if (count > 0) { if (count > 0) {
cerr() << info() << "[MSG] " << count << " error" cerr() << info() << "[MSG] " << count << " error"
<< (count > 1 ? "s were" : " was") << (count > 1 ? "s were" : " was")
<< " reported. Suggestions to fix this:" << std::endl; << " reported.";
for (auto& w: workarounds) { if (workarounds.size()) {
cerr() << info() << " - " << w << normal() << std::endl; cerr() << " Suggestions to fix this:" << std::endl;
for (auto& w: workarounds) {
cerr() << info() << " - " << w << normal() << std::endl;
}
} else {
cerr() << normal() << std::endl;
} }
if (fatal) { if (fatal) {
cout() << normal() <<"At least one fatal error encountered. Aborting process." << std::endl; cout() << normal() <<"At least one fatal error encountered. Aborting process." << std::endl;
......
...@@ -150,7 +150,8 @@ struct marker_obs_formats { ...@@ -150,7 +150,8 @@ struct marker_obs_formats {
const nlohmann::json& formats = get_internal_formats_json(); const nlohmann::json& formats = get_internal_formats_json();
auto it = formats.find(key); auto it = formats.find(key);
if (it == formats.end()) { if (it == formats.end()) {
MSG_ERROR("Marker observation format is not defined: " << key, SPELL_STRING("Use one of the predefined formats " << build_all_internal_format_names() << " or add a -mos <format-filename> argument to the commandline.")); MSG_ERROR("Marker observation format is not defined: " << key << ". Use one of the predefined formats " << build_all_internal_format_names() << " or add a -mos <format-filename> argument to the commandline.", "");
// MSG_ERROR("Marker observation format is not defined: " << key, SPELL_STRING("Use one of the predefined formats " << build_all_internal_format_names() << " or add a -mos <format-filename> argument to the commandline."));
return {}; return {};
} }
return it.value(); return it.value();
......
...@@ -570,14 +570,6 @@ bn_settings_t* bn_settings_t::from_args(int argc, const char** argv) ...@@ -570,14 +570,6 @@ bn_settings_t* bn_settings_t::from_args(int argc, const char** argv)
if (ret->output_mode == 0) { if (ret->output_mode == 0) {
ret->output_mode = bn_settings_t::OutputPopData; ret->output_mode = bn_settings_t::OutputPopData;
} }
if (ret->output_generations.size() == 0) {
for (const auto& kv: ret->pedigree.generation_names) {
ret->output_generations.emplace_back(kv.second);
}
}
std::sort(ret->output_generations.begin(), ret->output_generations.end());
auto it = std::unique(ret->output_generations.begin(), ret->output_generations.end());
ret->output_generations.resize(it - ret->output_generations.begin());
} }
return ret; return ret;
......
...@@ -194,7 +194,7 @@ job_registry = { ...@@ -194,7 +194,7 @@ job_registry = {
/*MSG_DEBUG("Pedigree" << std::endl << settings->pedigree);*/ /*MSG_DEBUG("Pedigree" << std::endl << settings->pedigree);*/
/*MSG_DEBUG("COMPUTING FACTOR GRAPH FOR " << unique_n_alleles[n] << " ALLELES.");*/ /*MSG_DEBUG("COMPUTING FACTOR GRAPH FOR " << unique_n_alleles[n] << " ALLELES.");*/
/*factor_graph fg(settings->pedigree, unique_n_alleles[n], settings->noise);*/ /*factor_graph fg(settings->pedigree, unique_n_alleles[n], settings->noise);*/
auto fg = factor_graph_type::from_pedigree(settings->pedigree.items, unique_n_alleles[n], settings->output_mode & bn_settings_t::OutputGametes, settings->input_generations, /*settings->output_generations*/ {}); auto fg = factor_graph_type::from_pedigree(settings->pedigree.items, unique_n_alleles[n], settings->output_mode & bn_settings_t::OutputGametes, settings->input_generations, settings->output_generations);
/*MSG_DEBUG("COMPUTED FACTOR GRAPH" << std::endl << fg);*/ /*MSG_DEBUG("COMPUTED FACTOR GRAPH" << std::endl << fg);*/
/*fg.finalize();*/ /*fg.finalize();*/
/*fg.save(settings->job_filename("factor-graph", unique_n_alleles[n]));*/ /*fg.save(settings->job_filename("factor-graph", unique_n_alleles[n]));*/
......
...@@ -42,7 +42,9 @@ int SPELL_BAYES_MAIN(int argc, const char** argv) ...@@ -42,7 +42,9 @@ int SPELL_BAYES_MAIN(int argc, const char** argv)
bn_settings_t* settings = NULL; bn_settings_t* settings = NULL;
try { try {
settings = bn_settings_t::from_args(argc, argv); settings = bn_settings_t::from_args(argc, argv);
msg_handler_t::check(true); if (msg_handler_t::check(false)) {
return -1;
}
} catch (file::error& fe) { } catch (file::error& fe) {
MSG_ERROR("An error happened while reading input files. " << fe.what(), ""); MSG_ERROR("An error happened while reading input files. " << fe.what(), "");
return -1; return -1;
...@@ -62,6 +64,25 @@ int SPELL_BAYES_MAIN(int argc, const char** argv) ...@@ -62,6 +64,25 @@ int SPELL_BAYES_MAIN(int argc, const char** argv)
MSG_DEBUG("loaded from " << settings->pedigree_filename); MSG_DEBUG("loaded from " << settings->pedigree_filename);
MSG_DEBUG("" << settings->pedigree); MSG_DEBUG("" << settings->pedigree);
MSG_DEBUG_DEDENT; MSG_DEBUG_DEDENT;
MSG_INFO("Output gen: " << settings->output_generations);
MSG_INFO("settings->output_generations.size = " << settings->output_generations.size());
MSG_INFO("settings->output_mode & OutputGametes = " << (settings->output_mode & bn_settings_t::OutputGametes));
if (settings->output_generations.size() == 0 || (settings->output_mode & bn_settings_t::OutputGametes)) {
settings->output_generations.clear();
settings->output_generations.reserve(settings->pedigree.generations.size());
for (const auto& gen: settings->pedigree.generations) {
if (gen) {
settings->output_generations.emplace_back(gen->name);
}
}
}
MSG_INFO("Output gen: " << settings->output_generations);
std::sort(settings->output_generations.begin(), settings->output_generations.end());
auto it = std::unique(settings->output_generations.begin(), settings->output_generations.end());
settings->output_generations.resize(it - settings->output_generations.begin());
for (const auto& kv: settings->observed_mark) { for (const auto& kv: settings->observed_mark) {
settings->marker_observation_specs[kv.second.format_name] settings->marker_observation_specs[kv.second.format_name]
......
...@@ -35,7 +35,7 @@ int SPELL_PEDIGREE_MAIN(int argc, const char** argv) ...@@ -35,7 +35,7 @@ int SPELL_PEDIGREE_MAIN(int argc, const char** argv)
print_usage(); print_usage();
} else { } else {
pedigree_type ped = read_pedigree(settings->pedigree_filename, settings->csv_sep, settings->max_states); pedigree_type ped = read_pedigree(settings->pedigree_filename, settings->csv_sep, settings->max_states);
if (!msg_handler_t::check(false)) { if (msg_handler_t::check(false)) {
return -1; return -1;
} }
std::string ped_output; std::string ped_output;
......
Package: spellmaptools Package: spellmaptools
Version: 1.0 Version: 0.3
Date: 2018-03-05 Date: 2018-03-05
Title: SpellMapTools Title: SpellMapTools
Authors@R: c(person("Damien", "Leroux", role = c("aut", "cre"), email = "damien.leroux@inra.fr"), Authors@R: c(person("Damien", "Leroux", role = c("aut", "cre"), email = "damien.leroux@inra.fr"),
......
...@@ -15,7 +15,7 @@ copy_sources_1: ...@@ -15,7 +15,7 @@ copy_sources_1:
@./cp.py . $(DOT_SRC) Makevars package/src @./cp.py . $(DOT_SRC) Makevars package/src
@./cp.py /usr/include/RWrap '*.h' package/src/include/RWrap @./cp.py /usr/include/RWrap '*.h' package/src/include/RWrap
@./cp.py ../../3rd-party json.hpp package/3rd-party @./cp.py ../../3rd-party json.hpp package/3rd-party
@./cp.py . DESCRIPTION package @./cp.py . .Rbuildignore DESCRIPTION package
copy_sources_2: .FORCE copy_sources_2: .FORCE
@./cp.py ~/include/eigen3/ '*' '*/*' '*/*/*' '*/*/*/*' package/src/include/eigen3 @./cp.py ~/include/eigen3/ '*' '*/*' '*/*/*' '*/*/*/*' package/src/include/eigen3
...@@ -49,7 +49,7 @@ package: generate ...@@ -49,7 +49,7 @@ package: generate
install: package install: package
+R CMD INSTALL $(TARGET)_1.0.tar.gz +R CMD INSTALL $(TARGET)_`grep Version DESCRIPTION | sed 's/.* //'`.tar.gz
clean: clean:
......
...@@ -41,11 +41,12 @@ class Spell2PtMatrix; ...@@ -41,11 +41,12 @@ class Spell2PtMatrix;
class SpellMapTools { class SpellMapTools {
private: private:
gamete_LV_database gamete_LV; gamete_LV_database gamete_LV;
double m_convergence_threshold = 1.e-12; double m_convergence_threshold = 1.e-6;
int m_max_iterations = 100; int m_max_iterations = 100;
bool m_active = false;
public: public:
/* observations_specs must have three columns: generation, format, filename */ /* observations_specs must have three columns: generation, format, filename */
SpellMapTools(std::string ped_filename, Rwrap::DataFrame observation_specs); SpellMapTools(std::string ped_filename, Rwrap::DataFrame observation_specs, int mt);
Rwrap::List Rwrap::List
SEM(marker_vec order); SEM(marker_vec order);
Rwrap::List Rwrap::List
...@@ -54,12 +55,14 @@ public: ...@@ -54,12 +55,14 @@ public:
Flips(marker_vec order, int window_size); Flips(marker_vec order, int window_size);
Spell2PtMatrix* Spell2PtMatrix*
LOD2pt(marker_vec rows, marker_vec cols); LOD2pt(marker_vec rows, marker_vec cols);
Spell2PtMatrix*
R2pt(marker_vec rows, marker_vec cols);
double double
map_likelihood(const marker_vec& order, const std::vector<double>& distances) map_likelihood(const marker_vec& order, const std::vector<double>& distances)
{ return gamete_LV.map_likelihood(order, distances); } { return gamete_LV.map_likelihood(order, distances); }
EM_map EM_map
EM(const marker_vec& order) EM(const marker_vec& order, bool dist_if_true_else_r=true)
{ return gamete_LV.EM(order, m_convergence_threshold, m_max_iterations); } { return gamete_LV.EM(order, dist_if_true_else_r, m_convergence_threshold, m_max_iterations); }
// void init_2pt_tr_at_inf() { gamete_LV.init_2pt_tr_at_inf(); } // void init_2pt_tr_at_inf() { gamete_LV.init_2pt_tr_at_inf(); }
double twoMarkerLikelihood() { return gamete_LV.twoMarkerLikelihoodAtInf(); } double twoMarkerLikelihood() { return gamete_LV.twoMarkerLikelihoodAtInf(); }
...@@ -75,7 +78,17 @@ public: ...@@ -75,7 +78,17 @@ public:
if (t) { std::swap(m_convergence_threshold, t); return t; } if (t) { std::swap(m_convergence_threshold, t); return t; }
else { return m_convergence_threshold; } else { return m_convergence_threshold; }
} }
bool active_check()
{
if (!m_active) {
MSG_ERROR("This instance is not active. There was a problem during initialization.", "");
return true;
}
return false;
}
bool is_active() { return m_active; }
}; };
...@@ -86,9 +99,11 @@ private: ...@@ -86,9 +99,11 @@ private:
std::vector<std::string> rownames, colnames; std::vector<std::string> rownames, colnames;
std::map<std::string, int> row_indices, col_indices; std::map<std::string, int> row_indices, col_indices;
std::map<std::pair<std::string, std::string>, double> cache; std::map<std::pair<std::string, std::string>, double> cache;
std::map<std::pair<std::string, std::string>, double> r_cache;
bool lod_if_true_else_r;
public: public:
Spell2PtMatrix(SpellMapTools* smt, std::vector<std::string> r, std::vector<std::string> c) Spell2PtMatrix(SpellMapTools* smt, std::vector<std::string> r, std::vector<std::string> c)
: instance(smt), rownames(r), colnames(c), row_indices(), col_indices(), cache() : instance(smt), rownames(r), colnames(c), row_indices(), col_indices(), cache(), r_cache(), lod_if_true_else_r(true)
{ {
for (const auto& m: rownames) { for (const auto& m: rownames) {
row_indices.emplace(m, row_indices.size()); row_indices.emplace(m, row_indices.size());
...@@ -98,6 +113,12 @@ public: ...@@ -98,6 +113,12 @@ public:
} }
} }
void rate() { lod_if_true_else_r = false; }
void lod() { lod_if_true_else_r = true; }
int is_rate() { return (int) !lod_if_true_else_r; }
int is_lod() { return (int) lod_if_true_else_r; }
marker_vec getRownames() { return rownames; } marker_vec getRownames() { return rownames; }
marker_vec getColnames() { return colnames; } marker_vec getColnames() { return colnames; }
...@@ -107,14 +128,15 @@ public: ...@@ -107,14 +128,15 @@ public:
return {(int) rownames.size(), (int) colnames.size()}; return {(int) rownames.size(), (int) colnames.size()};
} }
double get2PtValue(int i, int j) void
get_value_impl(int i, int j, double& lod, double& r)
{ {
if (i < 1 || j < 1 || i > rownames.size() || j > colnames.size()) { if (i < 1 || j < 1 || i > rownames.size() || j > colnames.size()) {
return NA_REAL; lod = r = NA_REAL;
} }
--i; --j; --i; --j;
if (rownames[i] == colnames[j]) { if (rownames[i] == colnames[j]) {
return 0; lod = r = 0;
} }
std::pair<std::string, std::string> key = {rownames[i], colnames[j]}; std::pair<std::string, std::string> key = {rownames[i], colnames[j]};
if (key.first > key.second) { if (key.first > key.second) {
...@@ -125,19 +147,31 @@ public: ...@@ -125,19 +147,31 @@ public:
// MSG_INFO("Key " << key.first << ':' << key.second << " not in cache!"); // MSG_INFO("Key " << key.first << ':' << key.second << " not in cache!");
static const std::vector<double> at_worlds_end = {.5}; static const std::vector<double> at_worlds_end = {.5};
marker_vec minimap = {rownames[i], colnames[j]}; marker_vec minimap = {rownames[i], colnames[j]};
auto em = instance->EM(minimap); auto em = instance->EM(minimap, false);
double best = em.likelihood; double best = em.likelihood;
r = em.distances[0];
r_cache.emplace(key, r);
// instance->init_2pt_tr_at_inf(); // instance->init_2pt_tr_at_inf();
// double inf_ref = instance->map_likelihood(minimap, at_worlds_end); // double inf_ref = instance->map_likelihood(minimap, at_worlds_end);
double inf_ref = instance->twoMarkerLikelihood(); double inf_ref = instance->twoMarkerLikelihood();
double value = best - inf_ref; lod = best - inf_ref;
// MSG_INFO("LOD(" << minimap[0] << " … " << minimap[1] << " => d=" << em.distances[0] << " best=" << best << " inf_ref=" << inf_ref << " value=" << value); // MSG_INFO("LOD(" << minimap[0] << " … " << minimap[1] << " => d=" << em.distances[0] << " best=" << best << " inf_ref=" << inf_ref << " value=" << lod);
cache.emplace(key, value); cache.emplace(key, lod);
return value;
} else { } else {
// MSG_INFO("Key " << key.first << ':' << key.second << " in cache! " << it->first.first << ':' << it->first.second); // MSG_INFO("Key " << key.first << ':' << key.second << " in cache! " << it->first.first << ':' << it->first.second);
return it->second; lod = it->second;
r = r_cache[key];
}
}
double get2PtValue(int i, int j)
{
if (!instance->is_active()) {
return NA_REAL;
} }
double lod, r;
get_value_impl(i, j, lod, r);
return lod_if_true_else_r ? lod : r;
} }
}; };
...@@ -181,7 +215,7 @@ extern "C" { ...@@ -181,7 +215,7 @@ extern "C" {
SpellMapTools::SpellMapTools(std::string ped_filename, Rwrap::DataFrame observation_specs) SpellMapTools::SpellMapTools(std::string ped_filename, Rwrap::DataFrame observation_specs, int mt)
: gamete_LV() : gamete_LV()