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

Implementation of multi-pop in progress.

parent 2d03b15a
......@@ -15,43 +15,65 @@
#include "computations/probabilities.h"
#endif
template <typename R, typename C>
model extend(const model& M, const labelled_matrix<MatrixXd, R, C>& bloc)
{
return {const_cast<model*>(&M)->extend(bloc.data)};
}
template <typename L>
struct bloc_builder {
struct element {
int first_row;
std::vector<int> columns;
};
std::vector<element> elements;
std::vector<L> labels;
int n_rows;
int n_columns;
};
struct pop_into_block {
int first_row;
std::vector<int> columns;
void build(MatrixXd& block, const parental_origin_per_locus_type& popl)
{
auto slice = block.middleRows(first_row, popl.innerSize());
for (size_t i = 0; i < popl.outerSize(); ++i) {
slice.col(columns[i]) = popl.data.col(i);
}
}
};
std::vector<pop_into_block>
init_bloc_pop_builder_disconnected(int n_qtl);
std::vector<pop_into_block>
init_bloc_pop_builder_connected(int n_qtl);
/*template <typename R, typename C>*/
/*model extend(const model& M, const labelled_matrix<MatrixXd, R, C>& bloc)*/
/*{*/
/*return {const_cast<model*>(&M)->extend(bloc.data)};*/
/*}*/
static inline
model extend(const model& M, const MatrixXd& bloc)
model extend(const model& M, const value<model_block_type>& bloc)
{
return {const_cast<model*>(&M)->extend(bloc)};
}
template <typename R, typename C>
model extend(const model& M, const value<labelled_matrix<MatrixXd, R, C>>& bloc)
{
return {const_cast<model*>(&M)->extend(bloc->data)};
}
static inline
model extend(const model& M, const value<MatrixXd>& bloc)
{
return {const_cast<model*>(&M)->extend(*bloc)};
}
MatrixXd
traits_to_matrix();
std::vector<double>
trait_by_name(population_value pop, const std::string& name);
MatrixXd
value<MatrixXd>
trait_permutations(const std::vector<double>& trait, int n);
VectorXd
f_tester(const model& M, const model& M0, const parental_origin_per_locus_type& popl);
MatrixXd
model_block_type
cross_indicator();
model
......@@ -71,7 +93,7 @@ int max_col(const MatrixXd& mat)
}
model
init_model(const MatrixXd& Y, const MatrixXd& indicator);
init_model(const value<MatrixXd>& Y, const value<model_block_type>& indicator);
VectorXd
row_max(const MatrixXd& mat);
......
......@@ -6,6 +6,7 @@
/*#include "algebraic_genotype.h"*/
#include "data/genoprob_computer.h"
#include "data/qtl_chromosome.h"
#include "lazy.h"
#include "lumping.h"
#include <Eigen/SparseCore>
#include <unsupported/Eigen/MatrixFunctions>
......@@ -431,9 +432,11 @@ struct generation_rs {
eval_poly_optim epo;
MatrixXd redux_matrix;
spin_lock redux_lock;
const MatrixXd& get_redux()
{
scoped_lock redux_guard(redux_lock);
if (redux_matrix.size() == 0) {
redux_matrix = make_redux(matrix_labels_to_cliques(main_process().column_labels));
}
......@@ -441,9 +444,11 @@ struct generation_rs {
}
std::vector<allele_pair> unique_labels;
spin_lock ul_lock;
const std::vector<allele_pair>& get_unique_labels()
{
scoped_lock ul_guard(ul_lock);
if (unique_labels.size() == 0) {
const MatrixXd& redux = get_redux();
const auto& breakmat = main_process();
......@@ -464,9 +469,14 @@ struct generation_rs {
std::vector<MatrixXb> state_to_parental_origin_haplo2;
std::vector<labelled_matrix<MatrixXd, std::vector<char>, allele_pair>> state_to_parental_origin;
/*spin_lock stpo_lock;*/
std::recursive_mutex stpo_lock;
labelled_matrix<MatrixXd, std::vector<char>, allele_pair>&
compute_state_to_parental_origin(size_t n_selected_qtl)
{
/*scoped_lock stpo_guard(stpo_lock);*/
std::lock_guard<std::recursive_mutex> stpo_guard(stpo_lock);
if (state_to_parental_origin.size() > n_selected_qtl) {
return state_to_parental_origin[n_selected_qtl];
}
......@@ -527,8 +537,8 @@ struct generation_rs {
ap_labels.begin(), ap_labels.end(),
h1n.cast<double>() + h2n.cast<double>());
}
std::cout << "Generation " << name << " with " << n_selected_qtl << " pop matrix:" << std::endl;
std::cout << state_to_parental_origin[n_selected_qtl] << std::endl;
/*std::cout << "Generation " << name << " with " << n_selected_qtl << " pop matrix:" << std::endl;*/
/*std::cout << state_to_parental_origin[n_selected_qtl] << std::endl;*/
return state_to_parental_origin[n_selected_qtl];
}
......@@ -910,9 +920,12 @@ struct generation_rs {
|| (a.first == b.second && a.second == b.first);
}
spin_lock ulv_lock;
/* NOT CONST! only for the sake of lazy evaluation */
const std::vector<VectorXd>& unphased_LV() const
{
scoped_lock ulv_guard(*const_cast<spin_lock*>(&ulv_lock));
if (!unique_unphased_LV.size()) {
std::vector<VectorXd>& uulv = *const_cast<std::vector<VectorXd>*>(&unique_unphased_LV);
#if 0
......
......@@ -10,6 +10,8 @@
#include "proba_config.h"
#include "generation_rs.h"
#include "../3rd-party/ThreadPool/ThreadPool.h"
struct bad_settings_exception : public std::exception {
std::string m;
bad_settings_exception(const std::string& msgs)
......@@ -166,6 +168,9 @@ struct settings_t {
detection_method_t detection_method;
double detection_window;
ThreadPool* pool;
static std::thread::id main_thread;
settings_t()
: notes()
, design(0)
......@@ -188,6 +193,7 @@ struct settings_t {
, cofactor_threshold(0)
, detection_method(detection_method_t::Undef)
, detection_window(10.)
, pool(0)
{}
std::vector<std::pair<const chromosome*, double>>
......@@ -233,7 +239,25 @@ struct settings_t {
for (auto& kv: populations) {
kv.second.name = kv.first;
}
if (parallel > 1) {
pool = new ThreadPool(parallel);
} else {
pool = new ThreadPool(1);
}
}
template <class F, class... Args>
auto enqueue(F&& f, Args&&... args)
-> std::future<typename std::result_of<F(Args...)>::type>
{
if (std::this_thread::get_id() != main_thread) {
return std::async(std::launch::deferred,
std::forward<F>(f),
std::forward<Args>(args)...);
}
return pool->enqueue(std::forward<F>(f),
std::forward<Args>(args)...);
}
};
std::ostream& operator << (std::ostream&, const settings_t&);
......
TARGET=spell-qtl
GCC_VERSION=-4.9
CXX=g++$(GCC_VERSION)
#CXX=clang++
COV=gcov$(GCC_VERSION)
INC=-I../include -I../include/input -I/usr/local/include/eigen3
LD=$(CXX)
LIBS=-lexpat -lssl
LIBS=-rdynamic -lexpat -lssl -ldl
VERSION_DEFINES=$(shell git tag | sed 's/\([0-9]*\)[.]\([0-9]*\)[.]\([0-9]*\)/-DVERSION_MAJOR=\\\"\1\\\" -DVERSION_MINOR=\\\"\2\\\" -DVERSION_PATCH=\\\"\3\\\"/')
CXXARGS=--std=gnu++0x -Wextra -Wall -pthread -Wno-unused-local-typedefs -Wextra -DEIGEN_NO_DEPRECATED_WARNING $(VERSION_DEFINES)
......@@ -18,9 +20,10 @@ COMPUTATIONS_SRC=basic_data.cc probabilities.cc model.cc
SRC=$(addprefix input/,$(INPUT_SRC)) $(addprefix computations/,$(COMPUTATIONS_SRC)) $(MAIN_SRC)
OBJ=$(subst .cc,.o,$(SRC))
DEP=$(subst .cc,.d,$(SRC))
COV_OBJ=$(subst .cc,.cov.o,$(SRC))
#DEBUG_OPTS=-Winvalid-pch -ggdb #-H
#DEBUG_OPTS=-Winvalid-pch -ggdb -O #-H
OPT_OPTS=-O3 -DEIGEN_NO_DEBUG -DNDEBUG
#COV_OPTS=$(DEBUG_OPTS) --coverage
#PROF_OPTS=-pg
......@@ -29,7 +32,8 @@ INFO_FILE=test-coverage.info
C=$(CXX) $(CXXARGS) $(INC) $(OPT_OPTS) $(DEBUG_OPTS) $(PROF_OPTS)
all: spell-qtl
all: $(TARGET)
include .depend
......@@ -52,19 +56,25 @@ test-coverage: test_fast_polynom
genhtml --highlight --legend --output-directory TestCodeCoverage $(INFO_FILE)
.depend: $(SRC) Makefile
$C --depend $(SRC) > $@
#.depend: $(SRC) Makefile
# $C --depend $(SRC) > $@
.depend: $(DEP)
cat $(DEP) > $@
probabilities.o test_selfing.o markov_population.o test_outbred.o outbred.o test_settings.o $(OBJ):%.o: %.cc
$C -c $< -o $@
$(DEP):%.d: %.cc
$C -MM $< -MT $(subst .d,.o,$@) -o $@
polynom.o: polynom.cc
$C -DTEST_POLY -c $(DEBUG_OPTS) $< -o $@ -DNDEBUG -g
clean:
rm -f $(OBJ) $(PCH)
#rm -f $(OBJ) $(PCH)
rm -f $(DEP) .depend
rm -f $(OBJ)
rm -f $(TARGET)
.FORCE:
......@@ -80,7 +90,7 @@ test.poly: polynom.o input/read_mark.o input/read_map.o input/design.o input/rea
test.cp: outbred.cc input/read_mark.o
$C -DTEST_OUTBRED $< input/read_mark.o -o $@ $(DEBUG_OPTS) $(LIBS)
spell-qtl: $(OBJ)
$(TARGET): $(OBJ)
$C $(OBJ) $(LIBS) -o $@
test_selfing: test_selfing.o input/read_mark.o input/read_map.o input/design.o input/read_trait.o input/read_settings.o
......
#include "computations/model.h"
template <typename L>
bloc_builder<L>
init_connected_bloc_builder(
std::function<int (const population&)> get_pop_size,
std::function<const std::vector<L>& (const population&)> get_labels)
{
int row_id = 0;
std::map<L, int> index;
for (auto& kv: active_settings->populations) {
for (const L& l: get_labels(kv.second)) {
int sz = index.size();
index.insert({l, sz});
}
}
bloc_builder<L> ret;
ret.labels.reserve(index.size());
ret.elements.reserve(active_settings->populations.size());
for (auto& kv: index) { ret.labels.push_back(kv.first); }
for (auto& kv: active_settings->populations) {
ret.elements.emplace_back();
auto& e = ret.elements.back();
const auto& lv = get_labels(kv.second);
e.first_row = row_id;
e.columns.reserve(lv.size());
for (const L& l: lv) {
e.columns.push_back(index[l]);
}
row_id += get_pop_size(kv.second);
}
return ret;
}
template <typename L>
bloc_builder<L>
init_disconnected_bloc_builder(
std::function<int (const population&)> get_pop_size,
std::function<const std::vector<L>& (const population&)> get_labels)
{
int row_id = 0;
int col_id = 0;
bloc_builder<L> ret;
size_t n_labels = 0;
for (auto& kv: active_settings->populations) {
n_labels += get_labels(kv.second).size();
}
ret.labels.reserve(n_labels);
ret.elements.reserve(active_settings->populations.size());
for (auto& kv: active_settings->populations) {
ret.elements.emplace_back();
auto& e = ret.elements.back();
const auto& lv = get_labels(kv.second);
e.first_row = row_id;
e.columns.reserve(lv.size());
for (const L& l: lv) {
e.columns.push_back(col_id++);
ret.labels.push_back(l);
}
row_id += get_pop_size(kv.second);
}
return ret;
}
#if 0
template <typename LABEL_TYPE>
struct unique_ordered_labels {
std::map<LABEL_TYPE, int> index;
std::vector<LABEL_TYPE> labels;
/*unique_ordered_labels(*/
};
std::vector<pop_into_block>
init_bloc_pop_builder_connected(int n_qtl)
{
/*int col_id = 0;*/
int row_id = 0;
std::set<std::vector<char>> uniq;
std::map<std::vector<char>, int> index;
for (auto& kv: active_settings->populations) {
const std::vector<std::vector<char>>& labels
= generation_by_name(kv.second.qtl_generation_name)->compute_state_to_parental_origin(n_qtl).row_labels;
for (auto& l: labels) {
if (uniq.find(l) == uniq.end()) {
index[l] = uniq.size();
uniq.insert(l);
}
}
}
std::vector<pop_into_block> ret;
ret.reserve(active_settings->populations.size());
for (auto& kv: active_settings->populations) {
const std::vector<std::vector<char>>& labels
= generation_by_name(kv.second.qtl_generation_name)->compute_state_to_parental_origin(n_qtl).row_labels;
ret.push_back({row_id, {}});
ret.back().columns.resize(labels.size());
for (auto& l: labels) {
ret.back().columns.push_back(index[l]);
}
row_id += kv.second.observed_traits.front().values.size();
}
return ret;
}
#endif
MatrixXd
traits_to_matrix()
{
......@@ -22,6 +137,57 @@ traits_to_matrix()
}
std::vector<pop_into_block>
init_bloc_pop_builder_disconnected(int n_qtl)
{
int col_id = 0;
int row_id = 0;
std::function<int()> incr = [&] () { return col_id++; };
std::vector<pop_into_block> ret;
ret.reserve(active_settings->populations.size());
for (auto& kv: active_settings->populations) {
const std::vector<std::vector<char>>& labels
= generation_by_name(kv.second.qtl_generation_name)->compute_state_to_parental_origin(n_qtl).row_labels;
ret.push_back({row_id, {}});
ret.back().columns.resize(labels.size());
std::generate_n(ret.back().columns.begin(), labels.size(), incr);
row_id += kv.second.observed_traits.front().values.size();
}
return ret;
}
std::vector<pop_into_block>
init_bloc_pop_builder_connected(int n_qtl)
{
int row_id = 0;
std::set<std::vector<char>> uniq;
std::map<std::vector<char>, int> index;
for (auto& kv: active_settings->populations) {
const std::vector<std::vector<char>>& labels
= generation_by_name(kv.second.qtl_generation_name)->compute_state_to_parental_origin(n_qtl).row_labels;
for (auto& l: labels) {
if (uniq.find(l) == uniq.end()) {
index[l] = uniq.size();
uniq.insert(l);
}
}
}
std::vector<pop_into_block> ret;
ret.reserve(active_settings->populations.size());
for (auto& kv: active_settings->populations) {
const std::vector<std::vector<char>>& labels
= generation_by_name(kv.second.qtl_generation_name)->compute_state_to_parental_origin(n_qtl).row_labels;
ret.push_back({row_id, {}});
ret.back().columns.resize(labels.size());
for (auto& l: labels) {
ret.back().columns.push_back(index[l]);
}
row_id += kv.second.observed_traits.front().values.size();
}
return ret;
}
std::vector<double>
trait_by_name(population_value pop, const std::string& name)
{
......@@ -34,12 +200,13 @@ trait_by_name(population_value pop, const std::string& name)
}
MatrixXd
value<MatrixXd>
trait_permutations(const std::vector<double>& trait, int n)
{
MatrixXd ret(trait.size(), n);
value<MatrixXd> ret = MatrixXd();
ret->resize(trait.size(), n);
for (int j = 0; j < n; ++j) {
auto col = ret.col(j);
auto col = ret->col(j);
for (size_t i = 0; i < trait.size(); ++i) {
col(i) = trait[i];
}
......@@ -53,10 +220,11 @@ VectorXd
f_tester(const model& M, const model& M0, const parental_origin_per_locus_type& popl)
{
model M1 = extend(M0, popl);
M1.compute();
return f_test(*const_cast<model*>(&M), *const_cast<model*>(&M1));
}
MatrixXd
model_block_type
cross_indicator()
{
int n_cols = active_settings->populations.size();
......@@ -64,12 +232,13 @@ cross_indicator()
for (auto& pop: active_settings->populations) {
n_rows += pop.second.observed_traits.front().values.size();
}
MatrixXd val = MatrixXd::Zero(n_rows, n_cols);
model_block_type val;
val.data = MatrixXd::Zero(n_rows, n_cols);
n_rows = 0;
int p = 0;
for (auto& pop: active_settings->populations) {
int sz = pop.second.observed_traits.front().values.size();
val.block(n_rows, p, sz, 1) = VectorXd::Ones(sz);
val.data.block(n_rows, p, sz, 1) = VectorXd::Ones(sz);
n_rows += sz;
++p;
}
......@@ -80,7 +249,8 @@ cross_indicator()
model
basic_model()
{
model M(traits_to_matrix());
value<MatrixXd> traits = make_value<mem>(traits_to_matrix);
model M(*traits);
M.add_bloc(cross_indicator());
const char* decomp_method = getenv("DECOMP_METHOD");
std::string s(decomp_method ? decomp_method : "");
......@@ -95,6 +265,7 @@ basic_model()
/*std::cout << "coefs" << std::endl << M.coefficients() << std::endl;*/
/*auto r = M.residuals().array();*/
/*std::cout << "residuals" << std::endl << (r * r).colwise().sum() << std::endl;*/
M.compute();
return M;
}
......@@ -116,7 +287,7 @@ f_test_along_chromosome(const value<model>& Mcurrent, const value<model>& M0,
model
init_model(const MatrixXd& Y, const MatrixXd& indicator)
init_model(const value<MatrixXd>& Y, const value<model_block_type>& indicator)
{
model M(Y);
const char* decomp_method = getenv("DECOMP_METHOD");
......@@ -127,12 +298,13 @@ init_model(const MatrixXd& Y, const MatrixXd& indicator)
M.use_QR();
}
M.add_bloc(indicator);
std::cout << "Xt" << std::endl << M.X().transpose() << std::endl;
std::cout << "Yt" << std::endl << M.Y().transpose() << std::endl;
std::cout << "XtX^-1" << std::endl << M.XtX_pseudo_inverse() << std::endl;
std::cout << "coefs" << std::endl << M.coefficients() << std::endl;
auto r = M.residuals().array();
std::cout << "residuals" << std::endl << (r * r).colwise().sum() << std::endl;
/*std::cout << "Xt" << std::endl << M.X().transpose() << std::endl;*/
/*std::cout << "Yt" << std::endl << M.Y().transpose() << std::endl;*/
/*std::cout << "XtX^-1" << std::endl << M.XtX_pseudo_inverse() << std::endl;*/
/*std::cout << "coefs" << std::endl << M.coefficients() << std::endl;*/
/*auto r = M.residuals().array();*/
/*std::cout << "residuals" << std::endl << (r * r).colwise().sum() << std::endl;*/
M.compute();
return M;
}
......
......@@ -32,6 +32,18 @@ parental_origin(const locus_probabilities_type& lp, generation_value gen, const
return s * lp;
}
parental_origin_type
parental_origin2(qtl_chromosome_value qtl_chr,
population_value pop,
generation_value gen,
const state_to_parental_origin_matrix_type& stfopom,
const multi_generation_observations& mgo,
const selected_qtls_on_chromosome_type& sqoc)
{
/*auto& s = const_cast<generation_rs*>(gen)->compute_state_to_parental_origin(sqoc.size());*/
return stfopom * locus_probabilities(qtl_chr, pop, gen, mgo, sqoc);
}
collection<parental_origin_per_locus_type>
parental_origin_per_locus(const collection<parental_origin_type>& apo)
......@@ -91,11 +103,14 @@ compute_parental_origins(const value<population_value>& pop,
= make_collection(population_marker_obs,
pop, chr, ap);
collection<locus_probabilities_type> alp
= make_collection(locus_probabilities, qtl_chr, pop, gen, apmo, sqoc);
/*collection<locus_probabilities_type> alp*/
/*= make_collection(locus_probabilities, qtl_chr, pop, gen, apmo, sqoc);*/
value<state_to_parental_origin_matrix_type>
stfopom = make_value(state_to_parental_origin_matrix, gen, qtl_chr, sqoc);
collection<parental_origin_type> apo
= make_collection(parental_origin, alp, gen, sqoc);
= make_collection(parental_origin2, qtl_chr, pop, gen, stfopom, apmo, sqoc);
return parental_origin_per_locus(apo);
}
......
......@@ -14,6 +14,10 @@ int main(int argc, const char** argv)
/*std::string digest = test;*/
/*std::cout << NORMAL << "digest " << digest << std::endl;*/
/*std::cout << get_func_name(f_tester) << std::endl;*/
/*return 0;*/
active_settings = settings_t::from_args(argc, argv);
if (!active_settings) {
exit(0);
......@@ -25,8 +29,10 @@ int main(int argc, const char** argv)
active_settings->finalize();
std::cout << "Cross Indicator" << std::endl << cross_indicator() << std::endl;
std::vector<double> t = trait_by_name(&active_settings->populations["toto"], "trait1");
std::cout << trait_permutations(t, 50) << std::endl;
/*std::cout << trait_permutations(t, 50) << std::endl;*/
active_settings->working_set.emplace_back(&active_settings->map[0]);
value<chromosome_value> chr = &active_settings->map[0];
......@@ -39,16 +45,17 @@ int main(int argc, const char** argv)
/*computations::f_test_along_chromosome ftac(M0, M0, pop, chr, gen, computations::selected_qtls_on_chromosome(qtl_chr));*/