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

Pre-0.1.

parent 00113d3e
......@@ -27,7 +27,7 @@ public:
void set_title(const std::string& s)
{
std::unique_lock<msg_handler_t::lock_type> lock(msg_handler_t::mutex);
std::lock_guard<std::mutex> lock(title_mutex);
title = s;
}
......@@ -47,6 +47,7 @@ private:
std::atomic_ulong done;
std::string title;
std::mutex title_mutex;
void display_progress();
};
......@@ -71,7 +72,12 @@ inline void ThreadPool::display_progress()
msg << SAVE_CURSOR;
msg /*<< HIDE_CURSOR*/ << GOTO_0_0;
msg << msg_handler_t::i();
msg << "[SPELL-QTL] Current: " << msg_handler_t::w() << title << msg_handler_t::i() << ERASE_TO_EOL << std::endl;
msg << "[SPELL-QTL] Current: " << msg_handler_t::w();
{
std::lock_guard<std::mutex> lock(title_mutex);
msg << title;
}
msg << msg_handler_t::i() << ERASE_TO_EOL << std::endl;
msg << "[SPELL-QTL] " << msg_handler_t::n() << workers.size() << " threads, " << total_queued << " tasks queued in previous batches" << msg_handler_t::i() << ERASE_TO_EOL << std::endl;
if (queued) {
msg << "[SPELL-QTL] Task info: ";
......
......@@ -439,9 +439,9 @@ register_task_in_progress(std::shared_ptr<async_computation<Ret(Args...)>> v,
Ret (*f) (Args...),
const value<typename clean_type<Args>::type>&... args)
{
__get_in_progress_mutex<Ret, Args...>().lock();
/*__get_in_progress_mutex<Ret, Args...>().lock();*/
__get_in_progress_registry<Ret, Args...>().get(f, args...) = v;
__get_in_progress_mutex<Ret, Args...>().unlock();
/*__get_in_progress_mutex<Ret, Args...>().unlock();*/
return v;
}
......@@ -480,10 +480,15 @@ template <typename Ret, typename... Args>
value_type& __get_noconst()
{
std::unique_lock<std::mutex> read_guard(mutex);
/*std::lock_guard<decltype(mutex)> read_guard(mutex);*/
/*if (m_future.valid()) {*/
/*return m_storage = m_future.get();*/
/*}*/
while (!mutex.try_lock());
if (m_future.valid()) {
return m_storage = m_future.get();
m_storage = m_future.get();
}
mutex.unlock();
return m_storage;
}
......@@ -496,7 +501,7 @@ template <typename Ret, typename... Args>
std::tuple<value<typename clean_type<Args>::type>...> dependencies;
value_type m_storage;
std::future<Ret> m_future;
std::mutex mutex;
std::recursive_mutex mutex;
};
template <typename Ret, typename... Args>
......@@ -505,6 +510,10 @@ std::shared_ptr<async_computation<Ret(Args...)>>
Ret (*func) (Args...),
const value<typename clean_type<Args>::type>&... args)
{
struct _mac_guard {
_mac_guard() { __get_in_progress_mutex<Ret, Args...>().lock(); }
~_mac_guard() { __get_in_progress_mutex<Ret, Args...>().unlock(); }
} _;
auto& r = __get_in_progress_registry<Ret, Args...>();
auto exists = r.find(func, args...);
if (exists) {
......@@ -524,6 +533,10 @@ std::shared_ptr<async_computation<Ret(Args...)>>
std::function<Ret(Args...)>& proxy,
const value<typename clean_type<Args>::type>&... args)
{
struct _mac_guard {
_mac_guard() { __get_in_progress_mutex<Ret, Args...>().lock(); }
~_mac_guard() { __get_in_progress_mutex<Ret, Args...>().unlock(); }
} _;
auto& r = __get_in_progress_registry<Ret, Args...>();
auto exists = r.find(func, args...);
if (exists) {
......@@ -736,29 +749,29 @@ template <typename T>
: m_min(min), m_max(max), m_step(step)
{}
struct iterator {
value<T> m_data;
T m_data;
T m_step;
T m_max;
iterator(T value, T max, T step)
: m_data(new immediate_value<T>(value)), m_step(step), m_max(max)
: m_data(value), m_step(step), m_max(max)
{}
iterator& operator ++ ()
{
*m_data += m_step;
if (*m_data > m_max) {
*m_data = m_max;
m_data += m_step;
if (m_data > m_max) {
m_data = m_max;
}
return *this;
}
bool operator == (const iterator& i) const
{
return *m_data == *i.m_data;
return m_data == i.m_data;
}
bool operator != (const iterator& i) const
{
return *m_data != *i.m_data;
return m_data != i.m_data;
}
const value<T>& operator * () const { return m_data; }
value<T> operator * () const { return {m_data}; }
};
iterator begin() const { return {m_min, m_max, m_step}; }
iterator end() const { return {m_max, m_max, 0}; }
......@@ -866,8 +879,8 @@ struct with_mem_cache_traits {
return_type
create(CachingPolicy _Sync, Ret (&f) (Args...), const clean_value_type<Args>&... x)
{
/*static std::mutex _;*/
/*std::unique_lock<std::mutex> lock_guard(_);*/
static std::mutex _;
std::lock_guard<std::mutex> lock_guard(_);
/*value<Ret>* ret_in_progress = __get_in_progress_registry<Ret, Args...>().find(&f, x...);*/
......@@ -893,6 +906,8 @@ struct without_mem_cache_traits {
return_type
create(CachingPolicy _Sync, Ret (&f) (Args...), const clean_value_type<Args>&... x)
{
static std::mutex _;
std::lock_guard<std::mutex> lock_guard(_);
/*value<Ret>* ret_in_progress = __get_in_progress_registry<Ret, Args...>().find(&f, x...);*/
/*if (ret_in_progress) {*/
......
......@@ -17,6 +17,23 @@ typedef labelled_matrix<MatrixXd, allele_pair, int> geno_prob_type;
typedef labelled_matrix<MatrixXd, std::vector<char>, allele_pair> state_to_parental_origin_matrix_type;
typedef labelled_matrix<MatrixXd, std::vector<char>, double> parental_origin_type;
typedef labelled_matrix<MatrixXd, int, std::vector<char>> parental_origin_per_locus_type;
typedef Eigen::Matrix<bool, 1, Eigen::Dynamic> bool_row_type;
struct compare_row {
bool operator () (const bool_row_type& r1, const bool_row_type& r2) const
{
for (int i = 0; i < r1.innerSize(); ++i) {
if (r1(i) > r2(i)) {
return true;
} else if (r1(i) < r2(i)) {
return false;
}
}
return false;
}
};
typedef std::set<bool_row_type, compare_row> bool_row_set;
std::vector<std::pair<const chromosome*, double>>
selected_qtls();
......@@ -61,6 +78,10 @@ marker_observation_specifications(generation_value);
labelled_matrix<MatrixXd, std::vector<char>, allele_pair>
compute_state_to_parental_origin(const context_key& ck, const locus_key& lk);
MatrixXb
get_contrast_groups(generation_value gen, const locus_key& lk);
/*observation_vectors_type*/
/*observation_vectors(const generation_rs*);*/
......
This diff is collapsed.
......@@ -38,13 +38,14 @@ struct block_builder {
return *this;
}
template <typename MATRIX>
struct iterator {
element_iterator_type ei;
MatrixXd* output;
MATRIX* output;
/* assemble */
iterator&
operator += (const MatrixXd& block)
operator += (const MATRIX& block)
{
/*MSG_DEBUG("output " << output);*/
/*MSG_DEBUG(output->innerSize() << ',' << output->outerSize());*/
......@@ -58,10 +59,16 @@ struct block_builder {
++ei;
return *this;
}
};
template <typename MATRIX>
struct disasm_iterator {
element_iterator_type ei;
MATRIX* output;
/* disassemble */
iterator&
operator -= (MatrixXd& block)
disasm_iterator&
operator -= (MATRIX& block)
{
/*MSG_DEBUG("output " << output);*/
/*MSG_DEBUG(output->innerSize() << ',' << output->outerSize());*/
......@@ -79,10 +86,17 @@ struct block_builder {
}
};
iterator begin(MatrixXd& output) const
template <typename MATRIX>
iterator<MATRIX> begin(MATRIX& output) const
{
/*output.resize(n_rows, n_columns);*/
output = MatrixXd::Zero(n_rows, n_columns);
output = MATRIX::Zero(n_rows, n_columns);
return {elements.begin(), &output};
}
template <typename MATRIX>
disasm_iterator<MATRIX> disasm_begin(MATRIX& output) const
{
return {elements.begin(), &output};
}
};
......@@ -146,12 +160,13 @@ void
init_connected_block_builder(
block_builder<L>& ret,
std::function<int (const population&)> get_pop_size,
std::function<std::vector<L> (const population&)> get_labels)
std::function<std::vector<L> (const population&)> get_labels,
const collection<population_value>& all_pops)
{
std::map<L, int> index;
for (auto& kv: active_settings->populations) {
for (const L& l: get_labels(kv.second)) {
for (auto& pop: all_pops) {
for (const L& l: get_labels(**pop)) {
int sz = index.size();
index.insert({l, sz});
}
......@@ -173,12 +188,12 @@ init_connected_block_builder(
ret.elements.reserve(active_settings->populations.size());
for (auto& kv: index) { ret.labels.push_back(kv.first); }
for (auto& kv: active_settings->populations) {
for (auto& kv: all_pops) {
ret.elements.emplace_back();
auto& e = ret.elements.back();
const auto& lv = get_labels(kv.second);
const auto& lv = get_labels(**kv);
e.first_row = ret.n_rows;
e.n_rows = get_pop_size(kv.second);
e.n_rows = get_pop_size(**kv);
e.columns.reserve(lv.size());
for (const L& l: lv) {
e.columns.push_back(index[l]);
......@@ -192,24 +207,25 @@ void
init_disconnected_block_builder(
block_builder<L>& ret,
std::function<int (const population&)> get_pop_size,
std::function<std::vector<L>(const population&)> get_labels)
std::function<std::vector<L>(const population&)> get_labels,
const collection<population_value>& all_pops)
{
ret.n_rows = 0;
int col_id = 0;
ret.n_columns = 0;
for (auto& kv: active_settings->populations) {
ret.n_columns += get_labels(kv.second).size();
for (auto& kv: all_pops) {
ret.n_columns += get_labels(**kv).size();
}
ret.labels.reserve(ret.n_columns);
ret.elements.reserve(active_settings->populations.size());
for (auto& kv: active_settings->populations) {
for (auto& kv: all_pops) {
ret.elements.emplace_back();
auto& e = ret.elements.back();
const auto& lv = get_labels(kv.second);
const auto& lv = get_labels(**kv);
e.first_row = ret.n_rows;
e.n_rows = get_pop_size(kv.second);
e.n_rows = get_pop_size(**kv);
e.columns.reserve(lv.size());
for (const L& l: lv) {
e.columns.push_back(col_id++);
......@@ -224,10 +240,14 @@ collection<model_block_type>
compute_parental_origins_multipop(const collection<population_value>& all_pops,
const value<chromosome_value>& chr,
const value<locus_key>& lk,
const value<std::vector<double>>& loci);
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci);
collection<std::vector<char>>
contrast_groups_connected();
/*collection<std::vector<char>>*/
/*contrast_groups_connected();*/
MatrixXd
contrast_groups(const collection<population_value>& all_pops, const locus_key& lk);
VectorXd
f_test_with_permutations(const std::string& trait, int n_permut, chromosome_value chr, const locus_key& lk);
......@@ -248,13 +268,15 @@ collection<model_block_type>
assemble_parental_origins_multipop(
const value<chromosome_value>& chr,
const value<locus_key>& lk,
const std::vector<collection<parental_origin_per_locus_type>>& all_popl);
const std::vector<collection<parental_origin_per_locus_type>>& all_popl,
const collection<population_value>& all_pops);
collection<parental_origin_per_locus_type>
disassemble_parental_origins_multipop(
const chromosome_value& chr,
const locus_key& lk,
model_block_type& block);
model_block_type& block,
const collection<population_value>& all_pops);
struct computation_along_chromosome {
MatrixXd coefficients;
......@@ -274,12 +296,13 @@ struct computation_along_chromosome {
}
};
enum ComputationType { FTest=1, FTestLOD=2, R2=4, Chi2=8, Mahalanobis=16 };
enum ComputationType { NoTest=0, FTest=1, FTestLOD=2, R2=4, Chi2=8, Mahalanobis=16 };
enum ComputationResults { Test=0, RSS=1, Residuals=2, Coefficients=4, Rank=8 };
constexpr ComputationType operator | (ComputationType a, ComputationType b) { return (ComputationType) (((int)a) | ((int)b)); }
constexpr ComputationResults operator | (ComputationResults a, ComputationResults b) { return (ComputationResults) (((int)a) | ((int)b)); }
namespace std {
template <>
struct hash<computation_along_chromosome> {
......@@ -314,10 +337,13 @@ void
compute_along_chromosome(computation_along_chromosome& ret,
ComputationType ct, ComputationResults cr,
const value<model>& Mcurrent, const value<model>& M0,
const model_block_key& base_key,
chromosome_value chr,
const std::vector<double>& loci,
const collection<parental_origin_per_locus_type>& popl)
{
if (cr & RSS) {
ret.rss.resize(Mcurrent->Y().innerSize(), Mcurrent->Y().outerSize() * popl.size());
ret.rss.resize(Mcurrent->Y().outerSize(), popl.size());
}
if (cr & Residuals) {
ret.residuals.resize(Mcurrent->Y().innerSize(), Mcurrent->Y().outerSize() * popl.size());
......@@ -352,15 +378,18 @@ compute_along_chromosome(computation_along_chromosome& ret,
collection<int> joiner;
joiner.reserve(popl.size());
/* FIXME: create correct keys? */
value<model_block_key> vmbk = model_block_key();
/* FIXME: create correct keys? FIXED */
for (auto i: range<int>(0, popl.size(), 1)) {
DUMP_FILE_LINE();
value<model_block_key> vmbk = model_block_key(base_key);
*vmbk += std::make_pair(chr, loci[*i]);
/*MSG_DEBUG("new block key " << vmbk);*/
joiner.push_back(make_value<_Policy>(compute_one,
i, pcac, /* flower */
vct, vcr,
Mcurrent, M0,
vmbk,
/*as_value(base_key + std::make_pair(chr, loci[i])),*/
popl[i]));
}
......
......@@ -14,9 +14,12 @@
/*const multi_generation_observations&,*/
/*const selected_qtls_on_chromosome_type&);*/
MatrixXd permutation_matrix(size_t n_parents, size_t n_loci, size_t order);
locus_probabilities_type
locus_probabilities(const context_key&, const locus_key&,
const multi_generation_observations&);
const multi_generation_observations&,
const std::vector<double>&);
state_to_parental_origin_matrix_type
state_to_parental_origin_matrix(generation_value, int);
......@@ -25,7 +28,7 @@ parental_origin_type
parental_origin(const locus_probabilities_type&, generation_value, const selected_qtls_on_chromosome_type&);
parental_origin_per_locus_type
joint_parental_origin_at_locus(const context_key& ck, const locus_key& lk, double locus);
joint_parental_origin_at_locus(const context_key& ck, const locus_key& lk);
collection<parental_origin_per_locus_type>
parental_origin_per_locus(const collection<parental_origin_type>&);
......
......@@ -47,16 +47,17 @@ struct chromosome {
void compute_haplo_sizes()
{
double prev = -1.;
double prev = 0;
int hs = 0;
haplo_sizes.reserve(raw.marker_locus.size());
for (double l: raw.marker_locus) {
++hs;
if (l != prev) {
haplo_sizes.push_back(hs);
hs = 0;
hs = 1;
prev = l;
} else {
++hs;
}
prev = l;
}
if (hs) {
haplo_sizes.push_back(hs);
......
......@@ -50,7 +50,19 @@ struct locus_key_struc {
return !this || locus == no_locus;
}
int depth() const
size_t locus_order(double loc) const
{
if (is_empty()) {
/* FIXME: should whine. */
return 0;
}
if (locus <= loc) {
return 0;
}
return 1 + parent->locus_order(loc);
}
size_t depth() const
{
if (!this || locus == no_locus) {
return 0;
......@@ -67,6 +79,21 @@ struct locus_key_struc {
|| (parent && parent->has(loc));
}
void _fill_vec(std::vector<double>& vec) const
{
if (parent) {
parent->_fill_vec(vec);
}
vec.push_back(locus);
}
std::vector<double> to_vec() const
{
std::vector<double> ret;
_fill_vec(ret);
return ret;
}
/*locus_key_struc* operator + (const locus_key& k) = delete;*/
/*locus_key_struc* operator - (const locus_key& k) = delete;*/
......@@ -92,6 +119,8 @@ struct locus_key_struc {
return k;
} else if (l == k->locus) {
return k->parent;
} else if (l > k->locus) {
return k;
} else {
locus_key k2 = k->parent - l;
return locus_key(new locus_key_struc(k2, k->locus));
......@@ -121,10 +150,12 @@ struct locus_key_struc {
{
bool bk1 = !!k1;
bool bk2 = !!k2;
return (bk1 && bk2
&& k1->locus == k2->locus
&& (k1->parent == k2->parent))
|| !(bk1 || bk2);
if (bk1 && bk2) {
return k1->locus == k2->locus
&& (k1->parent == k2->parent);
} else {
return !(bk1 || bk2);
}
}
friend
......@@ -139,6 +170,37 @@ struct locus_key_struc {
return k2 ? k1 + k2->parent + k2->locus : k1;
}
friend
bool operator < (const locus_key& k1, const locus_key& k2)
{
if (k1->is_empty()) {
return !k2->is_empty();
} else if (k2->is_empty()) {
return false;
}
std::vector<double> v1 = k1->to_vec();
std::vector<double> v2 = k2->to_vec();
auto i1 = v1.begin();
auto j1 = v1.end();
auto i2 = v2.begin();
auto j2 = v2.end();
while (i1 != j1 && i2 != j2 && *i1 == *i2) { ++i1; ++i2; }
if (i1 == j1) {
return i2 != j2;
} else if (i2 == j2) {
return false;
} else {
return *i1 < *i2;
}
}
friend
bool operator > (const locus_key& k1, const locus_key& k2)
{
return k2 < k1;
}
friend
locus_key operator & (const locus_key& locus_key1, const locus_key& locus_key2)
{
......@@ -184,9 +246,15 @@ struct locus_key_struc {
}
if (parent) {
return kroneckerProduct(
#ifdef LK_REDUCE_R2L
ret
,
#endif
parent->reduce(n_geno, loc_remove)
#ifndef LK_REDUCE_R2L
,
ret
#endif
);
} else {
return ret;
......@@ -211,7 +279,7 @@ struct locus_key_struc {
namespace std {
template <>
struct hash<locus_key> {
size_t operator () (const locus_key& k)
size_t operator () (const locus_key& k) const
{
hash<double> h;
return k
......@@ -221,6 +289,29 @@ namespace std {
: 0;
}
};
inline locus_key __empty_locus_key() { static locus_key _(new locus_key_struc()); return _; }
inline
locus_key_struc::iterator
begin(const locus_key& lk)
{
if (!lk) {
return __empty_locus_key()->begin();
}
return lk->begin();
}
inline
locus_key_struc::iterator
end(const locus_key& lk)
{
if (!lk) {
return __empty_locus_key()->end();
}
return lk->end();
}
}
/* two responsibilities :
......
......@@ -36,7 +36,22 @@ typedef std::shared_ptr<message_struc> message_handle;
/*#define MAKE_MESSAGE(_dest_, _expr_) do { std::stringstream __s; __s << _expr_; _dest_ = __s.str(); } while (0)*/
#define MESSAGE(_expr_) dynamic_cast<std::stringstream*>(&(std::stringstream() << _expr_))->str()
/*#define CREATE_MESSAGE(_var_, _channel_, _expr_) message_handle _var_(new message_struc {_channel_, MESSAGE(_expr_)});*/
#ifndef SPELL_UNSAFE_OUTPUT
#define CREATE_MESSAGE(_channel_, _what_) msg_handler_t::enqueue(message_handle{new message_struc {_channel_, _what_}});
#else
#define CREATE_MESSAGE(_channel_, _what_) switch(_channel_) { \
case msg_channel::Out: \
std::cout << _what_; \
break; \
case msg_channel::Err: \
std::cerr << _what_; \
break; \
case msg_channel::Log: \
std::clog << _what_; \