Commit f8749dd1 authored by Damien Leroux's avatar Damien Leroux
Browse files

Compiles with clang. executes 10% faster with clang.

parent 143174a8
...@@ -20,7 +20,7 @@ project(spell_qtl) ...@@ -20,7 +20,7 @@ project(spell_qtl)
#LIST(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/CMake") #LIST(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/CMake")
#INCLUDE(pandocology) #INCLUDE(pandocology)
find_package(Boost 1.55.0 REQUIRED) find_package(Boost 1.56.0 REQUIRED)
#set(CMAKE_CONFIGURATION_TYPES Debug Release CACHE INTERNAL FORCE) #set(CMAKE_CONFIGURATION_TYPES Debug Release CACHE INTERNAL FORCE)
set(CMAKE_VERBOSE_MAKEFILE ON) set(CMAKE_VERBOSE_MAKEFILE ON)
......
...@@ -126,7 +126,7 @@ struct state_index_type { ...@@ -126,7 +126,7 @@ struct state_index_type {
state_index_type state_index_type
bad() bad()
{ {
return {{~0ULL, ~0ULL, ~0ULL, ~0ULL}}; return {{{~0ULL, ~0ULL, ~0ULL, ~0ULL}}};
} }
friend friend
......
...@@ -31,6 +31,7 @@ ...@@ -31,6 +31,7 @@
#include "labelled_matrix.h" #include "labelled_matrix.h"
#if 0
namespace detail { namespace detail {
template<typename T> template<typename T>
struct has_const_iterator struct has_const_iterator
...@@ -73,7 +74,7 @@ namespace detail { ...@@ -73,7 +74,7 @@ namespace detail {
&& has_begin_end<T>::beg_value && has_begin_end<T>::beg_value
&& has_begin_end<T>::end_value> {}; && has_begin_end<T>::end_value> {};
} }
#endif
#ifdef USE_OPENSSL_MD5 #ifdef USE_OPENSSL_MD5
#include <openssl/evp.h> #include <openssl/evp.h>
......
...@@ -29,6 +29,8 @@ ...@@ -29,6 +29,8 @@
#define SIGNAL_DISPLAY_ONELINER #define SIGNAL_DISPLAY_ONELINER
#define CONSTRAINT_RESIDUAL_EPSILON 1.e-5
using boost::math::normal; // typedef provides default type is double. using boost::math::normal; // typedef provides default type is double.
using boost::math::cdf; using boost::math::cdf;
using boost::math::mean; using boost::math::mean;
...@@ -2380,6 +2382,22 @@ void analysis_report::report_final_model(model_manager& mm) ...@@ -2380,6 +2382,22 @@ void analysis_report::report_final_model(model_manager& mm)
mm.vMcurrent->compute(); mm.vMcurrent->compute();
/* Checks */
{
auto X = mm.vMcurrent->X();
if (mm.vMcurrent->rank() < X.cols()) {
MSG_ERROR("Insufficient rank for model matrix.", "Report to the developers");
}
int n_constraints = X.rows();
for (const auto& p: mm.all_pops) {
n_constraints -= (*p)->size();
}
auto res = mm.vMcurrent->residuals().bottomRows(n_constraints);
if (res.lpNorm<0>() > CONSTRAINT_RESIDUAL_EPSILON) {
MSG_ERROR("Inconsistent contraints for model matrix.", "Report to the developers");
}
}
MSG_DEBUG(__FILE__ << ':' << __LINE__ << " keys " << mm.vMcurrent->keys()); MSG_DEBUG(__FILE__ << ':' << __LINE__ << " keys " << mm.vMcurrent->keys());
if (mm.vMcurrent->n_obs() <= mm.vMcurrent->dof() + 1) { if (mm.vMcurrent->n_obs() <= mm.vMcurrent->dof() + 1) {
......
...@@ -122,13 +122,38 @@ struct geno_prob_computer { ...@@ -122,13 +122,38 @@ struct geno_prob_computer {
return LV->col(i).asDiagonal(); 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) void init(const geno_matrix* gen)
{ {
/*tr_cache.clear();*/ /*tr_cache.clear();*/
this_gen = gen; this_gen = gen;
/*redux = make_redux(matrix_labels_to_cliques(break_matrix.column_labels));*/ /*redux = make_redux(matrix_labels_to_cliques(break_matrix.column_labels));*/
#ifdef WITH_OVERLUMPING
experimental_lumper el(*this_gen); experimental_lumper el(*this_gen);
redux = el.make_collect(el.partition_on_labels(), this_gen->inf_mat.cols()); 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; inv_stat_dist = this_gen->stat_dist;
for (int i = 0; i < inv_stat_dist.size(); ++i) { for (int i = 0; i < inv_stat_dist.size(); ++i) {
inv_stat_dist(i) = 1. / inv_stat_dist(i); inv_stat_dist(i) = 1. / inv_stat_dist(i);
......
...@@ -30,6 +30,10 @@ ...@@ -30,6 +30,10 @@
struct locus_key_struc; struct locus_key_struc;
typedef std::shared_ptr<const locus_key_struc> locus_key; typedef std::shared_ptr<const locus_key_struc> locus_key;
inline bool _notnull(const void* p) { return !!p; }
inline bool _null(const void* p) { return !p; }
template <class C> bool _notnull(const C* ptr) { return _notnull(static_cast<const void*>(ptr)); }
struct locus_key_struc { struct locus_key_struc {
struct no_locus_exception { struct no_locus_exception {
const char* what() const { return "Bad locus"; } const char* what() const { return "Bad locus"; }
...@@ -64,7 +68,7 @@ struct locus_key_struc { ...@@ -64,7 +68,7 @@ struct locus_key_struc {
bool is_empty() const bool is_empty() const
{ {
return !this || locus == no_locus; return _null(this) || locus == no_locus;
} }
size_t locus_order(double loc) const size_t locus_order(double loc) const
...@@ -81,7 +85,7 @@ struct locus_key_struc { ...@@ -81,7 +85,7 @@ struct locus_key_struc {
size_t depth() const size_t depth() const
{ {
if (!this || locus == no_locus) { if (_null(this) || locus == no_locus) {
return 0; return 0;
} }
if (parent) { if (parent) {
...@@ -92,7 +96,7 @@ struct locus_key_struc { ...@@ -92,7 +96,7 @@ struct locus_key_struc {
bool has(double loc) const bool has(double loc) const
{ {
return !!this && (loc == locus || (parent && parent->has(loc))); return _notnull(this) && (loc == locus || (parent && parent->has(loc)));
} }
void _fill_vec(std::vector<double>& vec) const void _fill_vec(std::vector<double>& vec) const
......
...@@ -85,12 +85,9 @@ struct builder { ...@@ -85,12 +85,9 @@ struct builder {
builder(builder&& b) : _s(std::move(b._s)), out(b.out) {} builder(builder&& b) : _s(std::move(b._s)), out(b.out) {}
~builder() { auto synchronize_outputs = lock(); out << _s.str(); out.flush(); } ~builder() { auto synchronize_outputs = lock(); out << _s.str(); out.flush(); }
builder<S>& operator << (const char* x) { std::move(_s) << x; return *this; } builder<S>& operator << (const char* x);
template <typename X> builder<S>& operator << (X&& x) { _s << x; return *this; } template <typename X> builder<S>& operator << (X&& x);
// template <typename X> builder<S>& operator << (const X& x) { std::move(_s) << x; return *this; } builder<S>& operator << (_Stream&(*x)(_Stream&));
// template <typename X> unless_integral<X, builder<S>&> operator << (X&& x) { std::move(_s) << x; return *this; }
// template <typename X> if_integral<X, builder<S>&> operator << (X x) { std::move(_s) << x; return *this; }
builder<S>& operator << (_Stream&(*x)(_Stream&)) { _s << x; return *this; }
}; };
...@@ -102,12 +99,9 @@ struct endpoint { ...@@ -102,12 +99,9 @@ struct endpoint {
endpoint(_Stream& _) : out(_) {} endpoint(_Stream& _) : out(_) {}
builder<S> operator << (const char* x) { builder<S> b(out); b._s << x; return b; } builder<S> operator << (const char* x);
template <typename X> builder<S> operator << (X&& x) { builder<S> b(out); b._s << x; return b; } template <typename X> builder<S> operator << (X&& x);
// template <typename X> builder<S> operator << (const X& x) { builder<S> b(out); b._s << ((const X&) x); return b; } builder<S> operator << (_Stream&(*x)(_Stream&));
// template <typename X> unless_integral<X, builder<S>> operator << (const X&& x) { builder<S> b(out); b._s << x; return b; }
// template <typename X> if_integral<X, builder<S>> operator << (X x) { builder<S> b(out); b._s << x; return b; }
builder<S> operator << (_Stream&(*x)(_Stream&)) { builder<S> b(out); b._s << x; return b; }
}; };
......
...@@ -109,35 +109,14 @@ struct file { ...@@ -109,35 +109,14 @@ struct file {
template <typename ARG> template <typename ARG>
friend friend
file& operator << (file& f, ARG&& arg) file& operator << (file& f, ARG&& arg);
{
if (f.m_impl) {
f.m_impl << arg;
f.check_state();
}
return f;
}
template <typename ARG> template <typename ARG>
friend friend
file& operator >> (file& f, ARG&& arg) file& operator >> (file& f, ARG&& arg);
{
if (f.m_impl) {
f.check_state();
f.m_impl >> arg;
}
return f;
}
friend friend
file& operator << (file& f, std::ostream& (&func) (std::ostream&)) file& operator << (file& f, std::ostream& (&func) (std::ostream&));
{
if (f.m_impl) {
f.m_impl << func;
f.check_state();
}
return f;
}
file& file&
write(const char* s, std::streamsize n) write(const char* s, std::streamsize n)
......
...@@ -30,7 +30,7 @@ ...@@ -30,7 +30,7 @@
#undef WITH_OVERLUMPING #undef WITH_OVERLUMPING
#define WITH_OVERLUMPING // #define WITH_OVERLUMPING
/*#define eq_dbl(_x_, _y_) (fabs((_x_) - (_y_)) < FLOAT_TOL)*/ /*#define eq_dbl(_x_, _y_) (fabs((_x_) - (_y_)) < FLOAT_TOL)*/
...@@ -834,7 +834,7 @@ geno_matrix recall_matrix(label_type::letter_type n) ...@@ -834,7 +834,7 @@ geno_matrix recall_matrix(label_type::letter_type n)
#if 0 #ifndef WITH_OVERLUMPING
inline inline
geno_matrix lump_using_partition_weighted(const geno_matrix& m, const std::set<subset>& lumping_partition) geno_matrix lump_using_partition_weighted(const geno_matrix& m, const std::set<subset>& lumping_partition)
{ {
...@@ -1197,6 +1197,7 @@ struct experimental_lumper { ...@@ -1197,6 +1197,7 @@ struct experimental_lumper {
auto CP = make_collect(PP, M.cols()); auto CP = make_collect(PP, M.cols());
auto DP = make_dispatch(CP); auto DP = make_dispatch(CP);
auto M_hat = lump_using_partition_weighted(M, CP, DP); auto M_hat = lump_using_partition_weighted(M, CP, DP);
#ifdef WITH_OVERLUMPING
if (M_hat.cols() > max_size) { if (M_hat.cols() > max_size) {
experimental_lumper perfect(M_hat); experimental_lumper perfect(M_hat);
MSG_DEBUG_INDENT_EXPR("[Lumped geno_matrix] "); MSG_DEBUG_INDENT_EXPR("[Lumped geno_matrix] ");
...@@ -1207,8 +1208,11 @@ struct experimental_lumper { ...@@ -1207,8 +1208,11 @@ struct experimental_lumper {
auto D = perfect.make_dispatch(C); auto D = perfect.make_dispatch(C);
return lump_using_partition_weighted(perfect.M, C, D); return lump_using_partition_weighted(perfect.M, C, D);
} else { } else {
#endif
return M_hat; return M_hat;
#ifdef WITH_OVERLUMPING
} }
#endif
} }
bool insert_split(std::set<subset>& P, const subset& part, const subset& split) bool insert_split(std::set<subset>& P, const subset& part, const subset& split)
......
...@@ -686,7 +686,7 @@ express_conditional_probability(const combination_type<PARENT_TYPE, STATE_TYPE>& ...@@ -686,7 +686,7 @@ express_conditional_probability(const combination_type<PARENT_TYPE, STATE_TYPE>&
/*std::map<key_list, double> prob_expressed = extract_probability_of(table, expressed);*/ /*std::map<key_list, double> prob_expressed = extract_probability_of(table, expressed);*/
std::map<key_list, double> prob_condition = extract_probability_of(table, condition); std::map<key_list, double> prob_condition = extract_probability_of(table, condition);
for (auto& kv: prob_condition) { kv.second = 1. / kv.second; } for (auto& kv: prob_condition) { kv.second = 1. / kv.second; }
double n = 0; // double n = 0;
/*for (const auto& kv: prob_expressed) { n += kv.second; }*/ /*for (const auto& kv: prob_expressed) { n += kv.second; }*/
/*n = 1. / n;*/ /*n = 1. / n;*/
/*for (auto& kv: prob_expressed) { kv.second *= n; }*/ /*for (auto& kv: prob_expressed) { kv.second *= n; }*/
......
...@@ -482,8 +482,10 @@ struct section_renderer { ...@@ -482,8 +482,10 @@ struct section_renderer {
} /* namespace model_print */ } /* namespace model_print */
using model_print::matrix_with_sections;
template <typename RS, typename RF, typename CS, typename CF, typename MT> template <typename RS, typename RF, typename CS, typename CF, typename MT>
std::ostream& operator << (std::ostream& os, const model_print::matrix_with_sections<RS, RF, CS, CF, MT>& mws) std::ostream& operator << (std::ostream& os, const matrix_with_sections<RS, RF, CS, CF, MT>& mws)
{ {
model_print::section_renderer<RS, RF, CS, CF, MT> sr(mws); model_print::section_renderer<RS, RF, CS, CF, MT> sr(mws);
return os << sr; return os << sr;
......
...@@ -583,6 +583,7 @@ struct pedigree_type { ...@@ -583,6 +583,7 @@ struct pedigree_type {
individual_index_type dh(const std::string& name, individual_index_type p1) { return fill_db(name, spawn(name, {p1})); } individual_index_type dh(const std::string& name, individual_index_type p1) { return fill_db(name, spawn(name, {p1})); }
individual_index_type ancestor(const std::string& name) { return fill_db(name, spawn(name, {})); } individual_index_type ancestor(const std::string& name) { return fill_db(name, spawn(name, {})); }
#if 0
void propagate_symmetries(int n, geno_matrix& gen) void propagate_symmetries(int n, geno_matrix& gen)
{ {
MSG_DEBUG_INDENT_EXPR("[propagate symmetries #" << n << "] "); MSG_DEBUG_INDENT_EXPR("[propagate symmetries #" << n << "] ");
...@@ -630,6 +631,7 @@ struct pedigree_type { ...@@ -630,6 +631,7 @@ struct pedigree_type {
MSG_DEBUG_DEDENT; MSG_DEBUG_DEDENT;
MSG_DEBUG_DEDENT; MSG_DEBUG_DEDENT;
} }
#endif
void compute_generation(const std::string& generation_name, int n) void compute_generation(const std::string& generation_name, int n)
{ {
...@@ -730,7 +732,7 @@ struct pedigree_type { ...@@ -730,7 +732,7 @@ struct pedigree_type {
/*if (tree[n].is_crossing()) {*/ /*if (tree[n].is_crossing()) {*/
/*propagate_symmetries(n, *generations.back());*/ /*propagate_symmetries(n, *generations.back());*/
/*} else if (tree[n].is_ancestor()) {*/ /*} else if (tree[n].is_ancestor()) {*/
generations.back()->symmetries = symmetry_group_type(generations.back()->labels); // generations.back()->symmetries = symmetry_group_type(generations.back()->labels);
/*}*/ /*}*/
/**node_generations[n] = lump(new_gen);*/ /**node_generations[n] = lump(new_gen);*/
if (cached_gen) { if (cached_gen) {
...@@ -996,7 +998,7 @@ struct pedigree_type { ...@@ -996,7 +998,7 @@ struct pedigree_type {
}; };
typedef vector_iterator<label_type> label_iterator; typedef vector_iterator<label_type> label_iterator;
typedef vector_iterator<symmetry_table_type> symmetry_iterator; // typedef vector_iterator<symmetry_table_type> symmetry_iterator;
template <typename VALUE_TYPE> template <typename VALUE_TYPE>
VALUE_TYPE eval_one(size_t node, const std::vector<bool>& recompute, const std::vector<vector_iterator<VALUE_TYPE>>& iterators, VALUE_TYPE eval_one(size_t node, const std::vector<bool>& recompute, const std::vector<vector_iterator<VALUE_TYPE>>& iterators,
...@@ -1080,7 +1082,7 @@ struct pedigree_type { ...@@ -1080,7 +1082,7 @@ struct pedigree_type {
static bn_label_type reentrant_bn_label(size_t, const bn_label_type& l) { return l; } static bn_label_type reentrant_bn_label(size_t, const bn_label_type& l) { return l; }
static gencomb_type reentrant_LC(size_t, const gencomb_type&) { return {1.}; } static gencomb_type reentrant_LC(size_t, const gencomb_type&) { return {1.}; }
static genotype_comb_type reentrant_GLC(size_t, const genotype_comb_type&) { return {1.}; } static genotype_comb_type reentrant_GLC(size_t, const genotype_comb_type&) { return {1.}; }
static symmetry_table_type reentrant_sym(size_t, const symmetry_table_type& S) { return {permutation_type::identity(1), S.letters}; } // static symmetry_table_type reentrant_sym(size_t, const symmetry_table_type& S) { return {permutation_type::identity(1), S.letters}; }
template <typename VALUE_TYPE> template <typename VALUE_TYPE>
void init_iterators_rec(size_t node, const std::vector<bool>& recompute, std::vector<bool>& visited, void init_iterators_rec(size_t node, const std::vector<bool>& recompute, std::vector<bool>& visited,
......
...@@ -570,7 +570,7 @@ struct pedigree_tree_type { ...@@ -570,7 +570,7 @@ struct pedigree_tree_type {
: count(c), sub1(), sub2(), LR(true), RL(false) : count(c), sub1(), sub2(), LR(true), RL(false)
{} {}
bool is_null() const { return this == NULL; } bool is_null() const { return static_cast<const void*>(this) == NULL; }
bool operator == (const tree_descr_struc& other) const bool operator == (const tree_descr_struc& other) const
{ {
...@@ -587,10 +587,10 @@ struct pedigree_tree_type { ...@@ -587,10 +587,10 @@ struct pedigree_tree_type {
bool operator < (const tree_descr_struc& other) const bool operator < (const tree_descr_struc& other) const
{ {
if (!this) { if (!static_cast<const void*>(this)) {
return !!(&other); return !!static_cast<const void*>(&other);
} }
if (!(&other)) { if (!static_cast<const void*>(&other)) {
return false; return false;
} }
return count < other.count return count < other.count
...@@ -622,7 +622,7 @@ struct pedigree_tree_type { ...@@ -622,7 +622,7 @@ struct pedigree_tree_type {
bool next(tree_descr_bool_map_type& cursors, tree_descr_bool_map_type& visited) bool next(tree_descr_bool_map_type& cursors, tree_descr_bool_map_type& visited)
{ {
if (!this || visited[this]) { if (!static_cast<const void*>(this) || visited[this]) {
return true; return true;
} }
visited[this] = true; visited[this] = true;
...@@ -721,7 +721,7 @@ struct pedigree_tree_type { ...@@ -721,7 +721,7 @@ struct pedigree_tree_type {
std::ostream& operator << (std::ostream& os, const tree_descr_struc& td) std::ostream& operator << (std::ostream& os, const tree_descr_struc& td)
{ {
static int recurse = 0; static int recurse = 0;
if (!(&td)) { if (!static_cast<const void*>(&td)) {
return os; return os;
} }
os << td.count; os << td.count;
...@@ -944,9 +944,10 @@ compute_tree_mapping(const std::vector<int>& o1, const std::vector<int>& o2) ...@@ -944,9 +944,10 @@ compute_tree_mapping(const std::vector<int>& o1, const std::vector<int>& o2)
return ret; return ret;
} }
#if 0
struct expression_symmetries_type { struct expression_symmetries_type {
std::vector<std::vector<bool>> symmetry_cursors; // std::vector<std::vector<bool>> symmetry_cursors;
std::vector<std::vector<bool>> latent_symmetry_cursors; // std::vector<std::vector<bool>> latent_symmetry_cursors;
std::vector<std::vector<std::vector<int>>> input_cliques; std::vector<std::vector<std::vector<int>>> input_cliques;
std::vector<std::vector<std::vector<int>>> latent_input_cliques; std::vector<std::vector<std::vector<int>>> latent_input_cliques;
...@@ -954,7 +955,7 @@ struct expression_symmetries_type { ...@@ -954,7 +955,7 @@ struct expression_symmetries_type {
const pedigree_tree_type* tree; const pedigree_tree_type* tree;
expression_symmetries_type(const pedigree_tree_type& t) expression_symmetries_type(const pedigree_tree_type& t)
: symmetry_cursors(), latent_symmetry_cursors(), input_cliques(), latent_input_cliques(), tree(&t) : /*symmetry_cursors(), latent_symmetry_cursors(),*/ input_cliques(), latent_input_cliques(), tree(&t)
{ {
MSG_DEBUG_INDENT_EXPR("[EXP SYM TYPE] "); MSG_DEBUG_INDENT_EXPR("[EXP SYM TYPE] ");
auto descr = tree->root_descr(); auto descr = tree->root_descr();
...@@ -1786,7 +1787,7 @@ struct symmetry_propagator { ...@@ -1786,7 +1787,7 @@ struct symmetry_propagator {
return compute_propagated_(true, exp_sym.latent_input_cliques, lat_sym_propagators, get_input_descr, get_sym_group, get_lat_sym_group, labels, inf_mat, get_node_lumper); return compute_propagated_(true, exp_sym.latent_input_cliques, lat_sym_propagators, get_input_descr, get_sym_group, get_lat_sym_group, labels, inf_mat, get_node_lumper);
} }
}; };
#endif
#if 0 #if 0
......
...@@ -21,18 +21,7 @@ namespace output { ...@@ -21,18 +21,7 @@ namespace output {
}; };
template <typename _Stream, typename _Container> template <typename _Stream, typename _Container>
_Stream& output_value_container(_Stream& os, _Container&& container) { _Stream& output_value_container(_Stream& os, _Container&& container);
typedef traits<_Container> _Traits;
auto beg = container.begin(), end = container.end();
os << _Traits::prefix;
if (beg != end) {
os << (*beg);
for (++beg; beg != end; ++beg) {
os << _Traits::separator << (*beg);
}
}
return os << _Traits::suffix;
};
namespace detail { namespace detail {
template<typename T> template<typename T>
......
...@@ -205,6 +205,7 @@ protected: ...@@ -205,6 +205,7 @@ protected:
}; };
#endif #endif
#if 0
struct letter_permutation_type { struct letter_permutation_type {
std::map<label_type::letter_type, label_type::letter_type> table; std::map<label_type::letter_type, label_type::letter_type> table;
...@@ -407,8 +408,8 @@ struct letter_permutation_type { ...@@ -407,8 +408,8 @@ struct letter_permutation_type {
}; };
struct symmetry_table_type { struct symmetry_table_type {
//std::vector<std::pair<label_type::letter_type, label_type::letter_type>> switches; /* over letters */ std::vector<std::pair<label_type::letter_type, label_type::letter_type>> switches; /* over letters */
//std::map<int, int> table; /* over states */ std::map<int, int> table; /* over states */
permutation_type table; permutation_type table;
letter_permutation_type letters; letter_permutation_type letters;