Commit 0e0d423d authored by Damien Leroux's avatar Damien Leroux
Browse files

Work in progress.

parent 7b81f94f
......@@ -55,7 +55,7 @@ inline std::ostream& operator << (std::ostream& os, const subset& s)
inline std::ostream& operator << (std::ostream& os, const std::set<subset>& S)
{
for (auto s: S) {
os << s << std::endl;
os << s << " #" << s.size() << std::endl;
}
return os;
}
......@@ -221,7 +221,7 @@ struct lumper {
/*}*/
/* MARK */
if (subsets.size() > 1) {
MSG_DEBUG("KEYS " << labels[subsets.begin()->second.front()] << ": " << subsets);
/*MSG_DEBUG("KEYS " << labels[subsets.begin()->second.front()] << ": " << subsets);*/
P.erase(C);
for (auto& s: subsets) {
auto io = P.insert(s.second);
......@@ -282,7 +282,7 @@ struct lumper {
/*}*/
/* MARK */
if (subsets.size() > 1) {
MSG_DEBUG("KEYS " << labels[subsets.begin()->second.front()] << ": " << subsets);
/*MSG_DEBUG("KEYS " << labels[subsets.begin()->second.front()] << ": " << subsets);*/
P.erase(C);
for (auto& s: subsets) {
auto io = P.insert(s.second);
......@@ -334,6 +334,62 @@ struct lumper {
/*accum.genotype = {reduced_labels[B.front()], reduced_labels[C.front()]};*/
scalar_type accum = compute_sum(*C.begin(), B);
#ifdef LUMP_BY_COLS
ret(j, i) = accum;
#else
ret(i, j) = accum;
#endif
++j;
}
++i;
}
tm.stop();
/*MSG_DEBUG("Created lumped matrix in " << tm.accum << " seconds");*/
return ret;
}
scalar_type compute_sum(size_t s, const subset& B, const VectorXd& abundances)
{
cs.start();
scalar_type accum = scalar_type();
/*fast_polynom accum = fast_polynom::zero;*/
/*std::cout << "accum = " << accum << std::endl;*/
/*for (size_t j = B.find_first(); j != subset::npos; j = B.find_next(j)) {*/
for (int j: B) {
/*std::cout << "accum + " << G.data(i, j);*/
#ifdef LUMP_BY_COLS
accum += G(j, s);
#else
accum += G(s, j);
#endif
/*std::cout << " = " << accum << std::endl;*/
}
cs.stop();
return accum;
}
MATRIX_TYPE to_matrix(const std::set<subset>& P, const VectorXd& abundances, std::vector<label_type>& reduced_labels)
{
tm.start();
reduced_labels.clear();
reduced_labels.reserve(P.size());
for (auto& S: P) {
reduced_labels.push_back(labels[*S.begin()]);
}
MATRIX_TYPE ret(reduced_labels.size(), reduced_labels.size());
size_t i = 0;
size_t j;
for (auto& C: P) {
j = 0;
for (auto& B: P) {
/*algebraic_genotype accum = compute_sum(*C.begin(), B);*/
/*accum.type = algebraic_genotype::Type::Genotype;*/
/*accum.genotype = {reduced_labels[B.front()], reduced_labels[C.front()]};*/
scalar_type accum = compute_sum(*C.begin(), B, abundances);
#ifdef LUMP_BY_COLS
ret(j, i) = accum;
#else
......
......@@ -25,14 +25,14 @@ COMPUTATIONS_SRC=basic_data.cc probabilities.cc model.cc frontends.cc
SRC=$(addprefix input/,$(INPUT_SRC)) $(addprefix computations/,$(COMPUTATIONS_SRC)) $(MAIN_SRC)
OBJ=$(subst .cc,.o,$(SRC))
DEP=$(subst .cc,.d,$(SRC)) test_script.d matrixexp.d
DEP=$(subst .cc,.d,$(SRC)) test_script.d matrixexp.d matrixexp3.d
COV_OBJ=$(subst .cc,.cov.o,$(SRC))
ALL_BUT_MAIN_OBJ=$(subst main.o ,,$(OBJ))
DEBUG_OPTS=-ggdb
#DEBUG_OPTS=-ggdb
#DEBUG_OPTS=-ggdb -DNDEBUG
#OPT_OPTS=-O3 -DEIGEN_NO_DEBUG -DNDEBUG
OPT_OPTS=-O3 -DEIGEN_NO_DEBUG -DNDEBUG
#OPT_OPTS=-O3 -DNDEBUG
#PROF_OPTS=-O2 -DEIGEN_NO_DEBUG -DNDEBUG -g -pg
......@@ -164,5 +164,5 @@ test_bayes_dai: ../sandbox/bayes.cc bayes.d $(ALL_BUT_MAIN_OBJ)
test_sib: pedigree_analysis.cc static_data.o
$C $< static_data.o -o $@
matrixexp: matrixexp.cc ../include/lumping2.h
$C $< -o $@
matrixexp: matrixexp3.cc ../include/lumping2.h
$C -DORDER=$(ORDER) $< -o $@
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment