Commit 9fa68f57 authored by Damien Leroux's avatar Damien Leroux
Browse files

WIP.

parent ba677a50
......@@ -401,7 +401,7 @@ struct algebraic_genotype { /* must fit in a QWORD */
type = Type::Null;
/*r_exp = s_exp = 0;*/
/*poly = {0, 0, {1}};*/
poly = fast_polynom::one;
poly = fast_polynom::zero;
}
algebraic_genotype(locus_pair lp, Type t, const fast_polynom& rs)
......@@ -413,6 +413,11 @@ struct algebraic_genotype { /* must fit in a QWORD */
: genotype({{0, 0}, {0, 0}}), type(Type::Polynom), poly(rs)
{}
static algebraic_genotype Polynom(ancestor_allele ff, ancestor_allele fs, ancestor_allele sf, ancestor_allele ss, const fast_polynom& rs)
{
return {{{ff, fs}, {sf, ss}}, Type::Polynom, rs};
}
static algebraic_genotype Genotype(ancestor_allele ff, ancestor_allele fs, ancestor_allele sf, ancestor_allele ss, const fast_polynom& rs)
{
return {{{ff, fs}, {sf, ss}}, Type::Genotype, rs};
......@@ -515,8 +520,8 @@ struct algebraic_genotype { /* must fit in a QWORD */
os << (char)('0' + ag.genotype.second.second.tag);
break;
case Type::RilPolynom:
case Type::Polynom:
break;
case Type::Polynom:
default:
os << ag.genotype;
};
......@@ -574,7 +579,15 @@ struct algebraic_genotype { /* must fit in a QWORD */
};
break;
case Type::Polynom:
return {other.genotype, other.type, poly * other.poly};
return Polynom(genotype.first.first, genotype.first.second,
other.genotype.second.first, other.genotype.second.second,
poly * other.poly);
/*switch (other.type) {*/
/*case Type::Polynom:*/
/*default:*/
/*return {other.genotype, other.type, poly * other.poly};*/
/*};*/
break;
case Type::Genotype:
switch (other.type) {
case Type::Genotype:
......@@ -602,7 +615,11 @@ struct algebraic_genotype { /* must fit in a QWORD */
/*break;*/
case Type::Null:
throw error("Unsupported Genotype*Null");
case Type::Polynom: break; /* can't happen */
case Type::Polynom: // break; /* can't happen */
/*return Polynom(genotype.first.first, genotype.first.second,*/
/*other.genotype.second.first, other.genotype.second.second,*/
/*poly * other.poly);*/
return {genotype, type, poly * other.poly};
case Type::RilPolynom: break; /* can't happen */
};
break;
......@@ -937,6 +954,10 @@ struct eval_poly_optim {
for (size_t p = 0; p < m_epo->polynoms.size(); ++p) {
const impl::f_polynom& poly = m_epo->polynoms[p];
if (poly.P.size() == 0) {
values[p] = 0;
continue;
}
double accum = r_pow[poly.r_exp] * s_pow[poly.s_exp];
double tmp = poly.P[0];
for (size_t i = 1; i < poly.P.size(); ++i) {
......
......@@ -53,6 +53,52 @@ ftest_along_chromosome(chromosome_value chr, const locus_key& lk,
*/
struct signal_display {
static const char* tick(double x)
{
static const char* ticks[9] = { " ", "\u2581", "\u2582", "\u2583", "\u2584", "\u2585", "\u2586", "\u2587", "\u2588" };
return ticks[x < 0. ? 0
: x >= 1. ? 8
: ((int) floor(x * 9))];
}
VectorXd values;
int imax_;
bool above_;
signal_display(const VectorXd& v, int imax, bool above)
: values(v.innerSize()), imax_(imax), above_(above)
{
int sig_cols = msg_handler_t::termcols() - 3;
values = v;
/*MSG_DEBUG("values.innerSize = " << values.innerSize());*/
while (values.innerSize() >= sig_cols) {
int i = values.innerSize() >> 1;
values = values.transpose() * kroneckerProduct(MatrixXd::Identity(i, i), MatrixXd::Constant(1, 2, .5));
/*MSG_DEBUG("values.innerSize = " << values.innerSize());*/
}
double vmin = values.minCoeff();
double vmax = values.maxCoeff();
values = ((values.array() - vmin) / (vmax - vmin)).matrix();
}
friend std::ostream& operator << (std::ostream& os, const signal_display& sd)
{
os << _WHITE << '[';
for (int i = 0; i < sd.values.innerSize(); ++i) {
if (i == sd.imax_) {
os << (sd.above_ ? _GREEN : _RED);
}
os << tick(sd.values(i));
if (i == sd.imax_) {
os << _WHITE;
}
}
return os << ']' << _NORMAL;
}
};
enum probability_mode { Joint, Single };
......@@ -505,6 +551,7 @@ struct model_manager {
default:
last_computation = NULL;
};
return find_max();
}
......@@ -562,6 +609,10 @@ struct model_manager {
model_block_key mbk = locus_base_key;
mbk += std::make_pair(chrom_under_study, (*loci)[i_max]);
/*MSG_DEBUG("locus_base_key " << locus_base_key << " mbk " << mbk);*/
signal_display sd(*last_computation, i_max, max > threshold);
MSG_DEBUG("[COMPUTATION] " << loci->front() << sd << loci->back() << " max=" << max << " at " << (*loci)[i_max]);
return {chrom_under_study,
testpos[i_max], max, i_max, max > threshold,
mbk,
......
......@@ -14,10 +14,12 @@
#include <deque>
extern "C" {
#include <unistd.h>
#include <sys/ioctl.h>
}
#define _WHITE "\x1b[37;1m"
#define _RED "\x1b[31;1m"
#define _GREEN "\x1b[32;1m"
#define _YELLOW "\x1b[33;1m"
#define _CYAN "\x1b[36;1m"
#define _NORMAL "\x1b[0;m"
......@@ -259,6 +261,13 @@ struct msg_handler_t {
static void enqueue(const message_handle& mh) { instance().queue.enqueue(mh); }
static void wait_for_flush() { instance().queue.wait_for_flush(); }
static int termcols()
{
struct winsize w;
ioctl(0, TIOCGWINSZ, &w);
return w.ws_col;
}
static lock_type mutex;
};
......
......@@ -38,7 +38,7 @@ inline void output_polynomial(std::ostream& os, const std::vector<coef_t>& coefs
int exponent = coefs.size() - 1;
bool is_first = true;
if (exponent == -1) {
os << "nil";
os << "0";
}
auto c = coefs.rbegin();
auto e = coefs.rend();
......
......@@ -123,19 +123,18 @@ inline
GenerationRS convert(const GenoMatrix& gm)
{
std::vector<allele_pair> labels;
/*std::map<allele_pair, int> label_counts;*/
std::map<allele_pair, int> label_counts;
for (int i = 0; i < gm.innerSize(); ++i) {
labels.push_back(gm(i, 0).genotype.first);
/*label_counts[labels.back()]++;*/
label_counts[labels.back()]++;
}
/*for (auto kv: label_counts) {*/
/*std::cout << kv.first << ':' << kv.second << ' ';*/
/*CREATE_MESSAGE(msg_channel::Log, MESSAGE(kv.first << ':' << kv.second << ' '));*/
/*}*/
/*std::cout << std::endl;*/
/*std::cout << labels << std::endl;*/
/*MSG_DEBUG(std::endl << labels << std::endl);*/
for (int i = 0; i < gm.outerSize(); ++i) {
if (!(labels[i] == gm(0, i).genotype.second)) {
std::cerr << "oops label mismatch @" << i << ": " << labels[i] << ' ' << gm(0, i).genotype.second << std::endl;
MSG_ERROR("oops label mismatch @" << i << ": " << labels[i] << ' ' << gm(0, i).genotype.second, "");
}
}
GenerationRS ret(labels);
......@@ -281,8 +280,12 @@ typedef std::vector<std::pair<char, std::vector<allele_pair>>> marker_observatio
static inline std::ostream& operator << (std::ostream& os, const marker_observation_spec& mos)
{
os << "MARKER OBS[" << std::endl;
for (auto& p: mos) {
os << " * " << p.first << ": " << p.second << std::endl;
if (mos.front().first) {
for (auto& p: mos) {
os << " * " << p.first << ": " << p.second << std::endl;
}
} else {
os << "RELATIVE unknown='" << mos.back().first << '\'';
}
return os << ']';
}
......@@ -700,6 +703,59 @@ struct generation_rs {
};
struct generation_design_broman : public generation_design_base {
const generation_rs* parent1;
const generation_rs* parent2;
const generation_rs* parent3;
const generation_rs* parent4;
const generation_rs* parent5;
const generation_rs* parent6;
const generation_rs* parent7;
const generation_rs* parent8;
generation_design_broman(const generation_rs* my_gen,
const generation_rs* p1, const generation_rs* p2,
const generation_rs* p3, const generation_rs* p4,
const generation_rs* p5, const generation_rs* p6,
const generation_rs* p7, const generation_rs* p8)
: generation_design_base(my_gen)
, parent1(p1), parent2(p2), parent3(p3), parent4(p4)
, parent5(p5), parent6(p6), parent7(p7), parent8(p8)
{}
MatrixXb compute_mask(const chromosome* chr, multi_generation_observations& mgo, int i) const
{
auto p1 = parent1->design->compute_possibles(chr, mgo, mgo.get_mother(g, i));
auto p2 = parent2->design->compute_possibles(chr, mgo, mgo.get_mother(g, i));
auto p3 = parent3->design->compute_possibles(chr, mgo, mgo.get_mother(g, i));
auto p4 = parent4->design->compute_possibles(chr, mgo, mgo.get_mother(g, i));
auto p5 = parent5->design->compute_possibles(chr, mgo, mgo.get_mother(g, i));
auto p6 = parent6->design->compute_possibles(chr, mgo, mgo.get_mother(g, i));
auto p7 = parent7->design->compute_possibles(chr, mgo, mgo.get_mother(g, i));
auto p8 = parent8->design->compute_possibles(chr, mgo, mgo.get_mother(g, i));
auto nrow = p1.innerSize() + p2.innerSize() + p3.innerSize() + p4.innerSize()
+ p5.innerSize() + p6.innerSize() + p7.innerSize() + p8.innerSize();
MatrixXb ret(nrow, p1.outerSize());
nrow = 0;
ret.topRows(nrow = p1.innerSize()) = p1;
ret.middleRows(nrow, p1.innerSize()) = p2;
ret.middleRows(nrow += p1.innerSize(), p2.innerSize()) = p2;
ret.middleRows(nrow += p2.innerSize(), p3.innerSize()) = p3;
ret.middleRows(nrow += p3.innerSize(), p4.innerSize()) = p4;
ret.middleRows(nrow += p4.innerSize(), p5.innerSize()) = p5;
ret.middleRows(nrow += p5.innerSize(), p6.innerSize()) = p6;
ret.middleRows(nrow += p6.innerSize(), p7.innerSize()) = p7;
ret.middleRows(nrow += p7.innerSize(), p8.innerSize()) = p8;
return ret;
}
generation_design_base* _clone(const generation_rs* my_gen) const { return new generation_design_cross(my_gen, parent1, parent2); }
void to_str(std::ostream& os) const
{
os << "(" << (*parent1->design) << " x " << (*parent2->design) << (*parent3->design) << " x " << (*parent4->design)
<< (*parent5->design) << " x " << (*parent2->design) << (*parent1->design) << " x " << (*parent2->design) << ')';
}
virtual std::pair<const generation_rs*, const generation_rs*> get_parents() const { return {parent1, parent2}; }
};
struct generation_design_double_haploid : public generation_design_base {
const generation_rs* parent;
generation_design_double_haploid(const generation_rs* my_gen, const generation_rs* p)
......@@ -1711,19 +1767,6 @@ void population::init_observations(
ret = MatrixXb::Ones(g->main_process().innerSize(), chr->raw.marker_name.size());
}
}
#if 0
const auto& ovmap = g->observation_vectors;
if (i == -1 || !ovmap.size()) {
ret = MatrixXb::Ones(g->main_process().innerSize(), chr->marker_name.size());
} else {
std::vector<char> obs = get_observations(chr, g, i);
if (obs.size()) {
ret = g->raw_observations(ovmap, obs);
} else {
ret = MatrixXb::Ones(g->main_process().innerSize(), chr->marker_name.size());
}
}
#endif
}
......@@ -2057,6 +2100,111 @@ public:
}
#endif
static
generation_rs* broman8(const std::string& name,
generation_rs* A, generation_rs* B, generation_rs* C, generation_rs* D,
generation_rs* E, generation_rs* F, generation_rs* G, generation_rs* H)
{
generation_rs* ret = fetch_impl(name);
if (ret->design != NULL) {
return ret;
}
ret->design = new generation_design_broman(ret, A, B, C, D, E, F, G, H);
double v = .25, w = .5;
std::vector<generation_rs*> parents = {A, B, C, D, E, F, G, H};
/*std::vector<size_t> sizes = {*/
/*A->main_process().innerSize(), B->main_process().innerSize(),*/
/*C->main_process().innerSize(), D->main_process().innerSize(),*/
/*E->main_process().innerSize(), F->main_process().innerSize(),*/
/*G->main_process().innerSize(), H->main_process().innerSize()*/
/*};*/
/*size_t size = std::accumulate(sizes.begin(), sizes.end(), 0);*/
fast_polynom quarter = {.25};
fast_polynom two = {-2.};
MSG_DEBUG("quarter=" << quarter);
MSG_DEBUG("two=" << two);
algebraic_genotype
ag_quarter
= {{{0, 0}, {0,0}},
algebraic_genotype::Type::Polynom,
quarter},
ag_one
= {{{0, 0}, {0,0}},
algebraic_genotype::Type::Polynom,
fast_polynom::one},
ag_two
= {{{0, 0}, {0,0}},
algebraic_genotype::Type::Polynom,
two},
ag_minus_one
= {{{0, 0}, {0,0}},
algebraic_genotype::Type::Polynom,
fast_polynom::zero - fast_polynom::one};
MSG_DEBUG("ag_quarter=" << ag_quarter);
MSG_DEBUG("ag_one=" << ag_one);
MSG_DEBUG("ag_two=" << ag_two);
MSG_DEBUG("ag_minus_one=" << ag_minus_one);
GenoMatrix I4 = GenoMatrix::Constant(4, 1, ag_one).asDiagonal();
GenoMatrix I8 = GenoMatrix::Constant(8, 1, ag_two).asDiagonal();
GenoMatrix M2seed(2, 2);
M2seed << ag_minus_one, ag_one, ag_one, ag_minus_one;
GenoMatrix M1 = GenoMatrix::Constant(8, 8, ag_quarter) + I8,
M2 = kroneckerProduct(I4, M2seed),
M_I = GenoMatrix::Constant(8, 1, ag_one).asDiagonal();
for (int i = 0; i < 8; ++i) {
for (int j = 0; j < 8; ++j) {
M1(i, j).type = algebraic_genotype::Type::Genotype;
M2(i, j).type = algebraic_genotype::Type::Genotype;
M_I(i, j).type = algebraic_genotype::Type::Genotype;
M1(i, j).genotype = {parents[i]->main_process().column_labels[0],
parents[j]->main_process().column_labels[0]};
M2(i, j).genotype = M1(i, j).genotype;
M_I(i, j).genotype = M1(i, j).genotype;
}
}
MSG_DEBUG(MATRIX_SIZE(M1) << std::endl << M1);
MSG_DEBUG(MATRIX_SIZE(M2) << std::endl << M2);
MSG_DEBUG(MATRIX_SIZE(M_I) << std::endl << M_I);
fast_polynom ril_dl = {1, -2};
fast_polynom ril = ril_dl;
/*GenoMatrix tmp(1, 1);*/
algebraic_genotype tmp = {{{0, 0}, {0,0}},
algebraic_genotype::Type::Polynom,
fast_polynom::zero};
for (int i = 1; i < BROMAN_ORDER; ++i) {
tmp.poly = (fast_polynom::one - ril) / (coef_t)2;
ret->P.emplace_back((process){convert(M1 * tmp + M_I), v});
ret->P.emplace_back((process){convert(M2 * tmp + M_I), w});
ril *= ril_dl;
v *= .5;
w *= .5;
}
/* the rest is approximately ril_dl[BROMAN_ORDER+1] with the same weight */
ril *= ril_dl;
tmp.poly = (fast_polynom::one - ril) / (coef_t)2;
ret->P.emplace_back((process){convert(M1 * tmp + M_I), v});
ret->P.emplace_back((process){convert(M2 * tmp + M_I), w});
ret->P[0].weight *= -1; /* fix first weight */
ret->lump();
ret->reduce_processes();
ret->precompute();
MSG_DEBUG(ret);
return ret;
}
generation_rs* reduce_alleles(const std::string& new_name) const
{
generation_rs* ret = fetch_impl(new_name);
......
......@@ -17,6 +17,7 @@
#include "fmin.h"
#define RIL_ORDER 10
#define BROMAN_ORDER 5
using namespace Eigen;
......
......@@ -8,11 +8,12 @@
namespace script {
enum Type { type_none, type_bool, type_model, type_model_manager, type_locus, type_chromosome, type_result, type_ct };
enum Type { type_none, type_filename, type_bool, type_model, type_model_manager, type_locus, type_chromosome, type_result, type_ct };
static inline Type str_to_type(const std::string& s)
{
static const std::map<std::string, Type> _ = {
{"filename", type_filename},
{"bool", type_bool},
{"model", type_model_manager},
{"locus", type_locus},
......@@ -27,6 +28,7 @@ static inline const char* type_to_str(Type t)
{
static const std::map<Type, const char*> _ = {
{type_none, "none"},
{type_filename, "filename"},
{type_bool, "bool"},
{type_model_manager, "model"},
{type_locus, "locus"},
......@@ -47,10 +49,10 @@ struct DuplicateName : public std::runtime_error {
{}
};
struct Environment;
struct UndefinedVariable : public std::runtime_error {
UndefinedVariable(const std::string& name)
: std::runtime_error(MESSAGE("Name \"" << name << "\" is not defined."))
{}
UndefinedVariable(const Environment& env, const std::string& name);
};
static inline ComputationType str_to_test(const std::string& s)
......@@ -101,6 +103,7 @@ public:
const chromosome* chrom;
bool b;
model M;
std::string s;
model_manager* mm;
model_manager::test_result res;
ComputationType ct;
......@@ -125,6 +128,7 @@ public:
operator model_manager::test_result& () { if (type != type_result) { throw WrongType(*this, type_result); } return data.res; }
operator bool () { if (type != type_bool) { throw WrongType(*this, type_bool); } return data.b; }
operator ComputationType () { if (type != type_ct) { throw WrongType(*this, type_ct); } return data.ct; }
operator const std::string& () { if (type != type_filename) { throw WrongType(*this, type_filename); } return data.s; }
Variable& operator = (const model& m) { reset(m); return *this; }
Variable& operator = (model_manager* m) { reset(type_model_manager); data.mm = m; return *this; }
......@@ -133,17 +137,20 @@ public:
Variable& operator = (double l) { reset(type_locus); data.loc = l; return *this; }
Variable& operator = (bool b) { reset(type_bool); data.b = b; return *this; }
Variable& operator = (ComputationType ct) { reset(type_ct); data.ct = ct; return *this; }
Variable& operator = (const std::string& s) { reset(s); return *this; }
void reset(Type new_type)
{
if (type == new_type) { return; }
switch (type) {
case type_filename: data.s.~basic_string(); break;
case type_model: data.M.~model(); break;
case type_result: data.res.~test_result(); break;
default:;
};
type = new_type;
switch (type) {
case type_filename: new(&data.s) std::string(); break;
case type_model: new(&data.M) model(); break;
case type_result: new (&data.res) model_manager::test_result(); break;
default:;
......@@ -162,12 +169,19 @@ public:
new (&data.res) model_manager::test_result(r);
}
void reset(const std::string& s)
{
reset(type_filename);
new (&data.s) std::string(s);
}
Variable& operator = (const Variable& v)
{
#define CASE_ASSIGN(__type, __field) case type_##__type: data.__field = v.data.__field; break;
reset(v.type);
switch (type) {
CASE_ASSIGN(locus, loc);
CASE_ASSIGN(filename, s);
CASE_ASSIGN(chromosome, chrom);
CASE_ASSIGN(bool, b);
CASE_ASSIGN(model, M);
......@@ -180,6 +194,7 @@ public:
return *this;
}
#if 0
Variable& operator = (const std::string& s)
{
std::stringstream ss(s);
......@@ -198,6 +213,24 @@ public:
};
return *this;
}
#endif
friend std::ostream& operator << (std::ostream& os, const Variable& v)
{
switch (v.type) {
case type_locus: os << "(locus) " << v.data.loc; break;
case type_filename: os << "(filename) " << v.data.s; break;
case type_chromosome: os << "(chromosome) " << v.data.chrom; break;
case type_bool: os << "(bool) " << (v.data.b ? "true" : "false"); break;
case type_model: os << "(model) " << v.data.M; break;
case type_model_manager: os << "(model_manager) " << v.data.mm->Mcurrent.keys(); break;
case type_result: os << "(result) " << v.data.res; break;
case type_ct: os << "(ComputationType) " << test_to_str(v.data.ct); break;
case type_none: os << "<none>"; break;
default:;
};
return os;
}
};
......@@ -227,15 +260,35 @@ struct Environment {
} else if (parent != NULL) {
return parent->resolve(name);
}
throw UndefinedVariable(name);
throw UndefinedVariable(*this, name);
}
bool has(const std::string& name)
{
auto it = vars.find(name);
return it != vars.end() || (parent && parent->has(name));
}
friend std::ostream& operator << (std::ostream& os, const Environment& env)
{
if ((&env) == NULL) {
return os << "NULL ENVIRONMENT";
}
os << "Variables defined:" << std::endl;
for (const auto& kv: env.vars) {
os << " - " << kv.first << " = " << kv.second << std::endl;
}
if (env.parent != NULL) {
os << *env.parent;
}
return os;
}
};
inline UndefinedVariable::UndefinedVariable(const Environment& env, const std::string& name)
: std::runtime_error(MESSAGE("Name \"" << name << "\" is not defined." << std::endl << env))
{}
/* EXPRESSIONS */
struct Expression {
......@@ -577,14 +630,14 @@ struct PrintScores : public Printable {
virtual std::vector<std::vector<double>> columns(Environment& env) const
{
model_manager& mm = env.resolve("");
return matrix_to_columns(*mm.last_computation);
return matrix_to_columns(mm.last_computation->transpose());
}
virtual std::vector<std::string> colnames(Environment& env) const
{
std::vector<std::string> ret;
model_manager& mm = env.resolve("");
ret.resize(mm.last_computation->outerSize(), test_to_str(env.resolve("Last Computation")));
ret.resize(mm.last_computation->innerSize(), test_to_str(env.resolve("Last computation")));
return ret;
}
};
......@@ -834,10 +887,13 @@ struct Print : public Output {
struct TableOutput : public Output {
std::vector<std::string> filename_