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

Fixed a lot of bugs. Cache is working properly.

parent fb382a1c
<breeding-design>
<ancestor-profile name="diallelic-homozygous">
<genotypes>
<genotype mater="0" pater="0"/>
<genotype mater="2" pater="2"/>
</genotypes>
</ancestor-profile>
<generation name="A">
<ancestor haplotype="a"/>
<ancestor haplotype="a" profile="diallelic-homozygous"/>
</generation>
<generation name="B">
<ancestor haplotype="b"/>
<ancestor haplotype="b" profile="diallelic-homozygous"/>
</generation>
<generation name="C">
<ancestor haplotype="c"/>
<ancestor haplotype="c" profile="diallelic-homozygous"/>
</generation>
<generation name="D">
<ancestor haplotype="d"/>
<ancestor haplotype="d" profile="diallelic-homozygous"/>
</generation>
<generation name="E">
<ancestor haplotype="e"/>
<ancestor haplotype="e" profile="diallelic-homozygous"/>
</generation>
<generation name="F">
<ancestor haplotype="f"/>
<ancestor haplotype="f" profile="diallelic-homozygous"/>
</generation>
<generation name="G">
<ancestor haplotype="g"/>
<ancestor haplotype="g" profile="diallelic-homozygous"/>
</generation>
<generation name="H">
<ancestor haplotype="h"/>
<ancestor haplotype="h" profile="diallelic-homozygous"/>
</generation>
<generation name="AB">
<cross p1="A" p2="B"/>
<generation name="BA">
<cross p1="B" p2="A"/>
</generation>
<generation name="CD">
<cross p1="C" p2="D"/>
<generation name="DC">
<cross p1="D" p2="C"/>
</generation>
<generation name="EF">
<cross p1="E" p2="F"/>
<generation name="HE">
<cross p1="H" p2="E"/>
</generation>
<generation name="GH">
<cross p1="G" p2="H"/>
<generation name="GF">
<cross p1="G" p2="F"/>
</generation>
<generation name="ABCD">
<cross p1="AB" p2="CD"/>
<generation name="DCBA">
<cross p1="DC" p2="BA"/>
</generation>
<generation name="EFGH">
<cross p1="EF" p2="GH"/>
<generation name="HEGF">
<cross p1="HE" p2="GF"/>
</generation>
<generation name="IC8">
<cross p1="ABCD" p2="EFGH"/>
<cross p1="HEGF" p2="DCBA"/>
</generation>
<generation name="MAGIC8-RIL">
<ril p1="IC8"/>
......
......@@ -20,7 +20,12 @@
<generation name="CP">
<cross p1="AB" p2="CD"/>
</generation>
<generation name="CP-hd">
<dh p1="CP"/>
</generation>
<!--
<generation name="self-CP">
<self p1="CP"/>
</generation>
-->
</breeding-design>
<format-specification>
<marker-observations generation="MAGIC8-RIL">
<marker-observations generation="MAGIC8-RIL" format="mapmaker">
<observation symbol="0">
<genotype from-mother="0" from-father="0" />
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<genotype from-mother="2" from-father="2" />
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<genotype from-mother="0" from-father="0" />
<genotype from-mother="2" from-father="2" />
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
<marker-observations generation="A" format="mapmaker">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
<marker-observations generation="B" format="mapmaker">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
<marker-observations generation="C" format="mapmaker">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
<marker-observations generation="D" format="mapmaker">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
<marker-observations generation="E" format="mapmaker">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
<marker-observations generation="F" format="mapmaker">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
<marker-observations generation="G" format="mapmaker">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
<marker-observations generation="H" format="mapmaker">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
</observation>
</marker-observations>
</format-specification>
<format-specification>
<marker-observations generation="F2" format="mapmaker">
<observation symbol="A">
<genotype>aa</genotype>
<genotype from-mother="a" from-father="a"/>
</observation>
<observation symbol="B">
<genotype>bb</genotype>
<genotype from-mother="b" from-father="b"/>
</observation>
<observation symbol="H">
<genotype>ab</genotype>
<genotype>ba</genotype>
<genotype from-mother="a" from-father="b"/>
<genotype from-mother="b" from-father="a"/>
</observation>
<observation symbol="C">
<genotype>aa</genotype>
<genotype>ba</genotype>
<genotype>ab</genotype>
<genotype from-mother="a" from-father="b"/>
<genotype from-mother="b" from-father="a"/>
<genotype from-mother="a" from-father="a"/>
</observation>
<observation symbol="D">
<genotype>bb</genotype>
<genotype>ba</genotype>
<genotype>ab</genotype>
<genotype from-mother="a" from-father="b"/>
<genotype from-mother="b" from-father="a"/>
<genotype from-mother="b" from-father="b"/>
</observation>
<observation symbol="-">
<genotype>aa</genotype>
<genotype>bb</genotype>
<genotype>ab</genotype>
<genotype>ba</genotype>
<genotype from-mother="a" from-father="b"/>
<genotype from-mother="b" from-father="a"/>
<genotype from-mother="a" from-father="a"/>
<genotype from-mother="b" from-father="b"/>
</observation>
</marker-observations>
<copy generation="F3" from="F2"/>
......@@ -34,10 +34,10 @@
<copy generation="F7" from="F2"/>
<marker-observations generation="BC" format="mapmaker">
<observation symbol="A">
<genotype>aa</genotype>
<genotype from-mother="a" from-father="a"/>
</observation>
<observation symbol="H">
<genotype>ab</genotype>
<genotype from-mother="a" from-father="b"/>
</observation>
</marker-observations>
</format-specification>
......@@ -348,7 +348,7 @@ namespace std {
for (auto& d: rs.P) {
accum ^= hd(d);
}
std::cout << "hashing " << rs << " => " << accum << std::endl;
/*std::cout << "hashing " << rs << " => " << accum << std::endl;*/
return accum;
}
};
......@@ -405,22 +405,22 @@ struct algebraic_genotype { /* must fit in a QWORD */
: genotype({{0, 0}, {0, 0}}), type(Type::Polynom), poly(rs)
{}
static algebraic_genotype Genotype(char ff, char fs, char sf, char ss, const fast_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};
}
static algebraic_genotype Cross(char ff, char fs, char sf, char ss, const fast_polynom& rs)
static algebraic_genotype Cross(ancestor_allele ff, ancestor_allele fs, ancestor_allele sf, ancestor_allele ss, const fast_polynom& rs)
{
return {{{ff, fs}, {sf, ss}}, Type::Cross, rs};
}
static algebraic_genotype Gamete(char ff, char sf, const fast_polynom& rs)
static algebraic_genotype Gamete(ancestor_allele ff, ancestor_allele sf, const fast_polynom& rs)
{
return {{{ff, 0}, {sf, 0}}, Type::Gamete, rs};
}
static algebraic_genotype Haplotype(char ff, char sf, const fast_polynom& rs)
static algebraic_genotype Haplotype(ancestor_allele ff, ancestor_allele sf, const fast_polynom& rs)
{
return {{{ff, 0}, {sf, 0}}, Type::Haplotype, rs};
}
......@@ -469,6 +469,8 @@ struct algebraic_genotype { /* must fit in a QWORD */
throw error("can't + two non-null algebraic_genotype");
}
bool operator == (const algebraic_genotype ag) const { return type == ag.type && genotype == ag.genotype && poly == ag.poly; }
friend std::ostream& operator << (std::ostream& os, algebraic_genotype::Type t)
{
switch(t) {
......@@ -493,16 +495,16 @@ struct algebraic_genotype { /* must fit in a QWORD */
os << ag.genotype.first.first << ',' << ag.genotype.second.first;
break;
case algebraic_genotype::Type::Gamete:
os << (char)('0' + ag.genotype.first.first);
os << (char)('0' + ag.genotype.first.first.tag);
os << ',';
os << (char)('0' + ag.genotype.second.first);
os << (char)('0' + ag.genotype.second.first.tag);
break;
case algebraic_genotype::Type::Cross:
os << (char)('0' + ag.genotype.first.first);
os << (char)('0' + ag.genotype.first.second);
os << (char)('0' + ag.genotype.first.first.tag);
os << (char)('0' + ag.genotype.first.second.tag);
os << ',';
os << (char)('0' + ag.genotype.second.first);
os << (char)('0' + ag.genotype.second.second);
os << (char)('0' + ag.genotype.second.first.tag);
os << (char)('0' + ag.genotype.second.second.tag);
break;
case Type::RilPolynom:
case Type::Polynom:
......@@ -548,10 +550,10 @@ struct algebraic_genotype { /* must fit in a QWORD */
case Type::DoublingGamete:
switch (other.type) {
case Type::Genotype:
return Genotype(other.genotype.first % (bool)genotype.first.first,
other.genotype.first % (bool)genotype.first.first,
other.genotype.second % (bool)genotype.second.first,
other.genotype.second % (bool)genotype.second.first,
return Genotype(other.genotype.first % (bool)genotype.first.first.tag,
other.genotype.first % (bool)genotype.first.first.tag,
other.genotype.second % (bool)genotype.second.first.tag,
other.genotype.second % (bool)genotype.second.first.tag,
poly * other.poly);
break;
case Type::Gamete:
......@@ -573,21 +575,21 @@ struct algebraic_genotype { /* must fit in a QWORD */
throw error("Unsupported Genotype*Haplotype");
break; /* make SJ happy */
case Type::DoublingGamete:
return Genotype(genotype.first % (bool)other.genotype.first.first,
genotype.first % (bool)other.genotype.first.first,
genotype.second % (bool)other.genotype.second.first,
genotype.second % (bool)other.genotype.second.first,
return Genotype(genotype.first % (bool)other.genotype.first.first.tag,
genotype.first % (bool)other.genotype.first.first.tag,
genotype.second % (bool)other.genotype.second.first.tag,
genotype.second % (bool)other.genotype.second.first.tag,
poly * other.poly);
case Type::Gamete:
return Haplotype(genotype.first % (bool)other.genotype.first.first,
genotype.second % (bool)other.genotype.second.first,
return Haplotype(genotype.first % (bool)other.genotype.first.first.tag,
genotype.second % (bool)other.genotype.second.first.tag,
poly * other.poly);
/*break;*/
case Type::Cross:
return Genotype(genotype.first % (bool)other.genotype.first.first,
genotype.first % (bool)other.genotype.first.second,
genotype.second % (bool)other.genotype.second.first,
genotype.second % (bool)other.genotype.second.second,
return Genotype(genotype.first % (bool)other.genotype.first.first.tag,
genotype.first % (bool)other.genotype.first.second.tag,
genotype.second % (bool)other.genotype.second.first.tag,
genotype.second % (bool)other.genotype.second.second.tag,
poly * other.poly);
/*break;*/
case Type::Null:
......@@ -603,8 +605,8 @@ struct algebraic_genotype { /* must fit in a QWORD */
genotype.second.first, other.genotype.second.first,
poly * other.poly);
case Type::Genotype:
return Haplotype(other.genotype.first % (bool)genotype.first.first,
other.genotype.second % (bool)genotype.second.first,
return Haplotype(other.genotype.first % (bool)genotype.first.first.tag,
other.genotype.second % (bool)genotype.second.first.tag,
poly * other.poly);
case Type::DoublingGamete:
if (genotype == other.genotype) {
......@@ -621,6 +623,8 @@ struct algebraic_genotype { /* must fit in a QWORD */
case Type::Haplotype:
switch (other.type) {
case Type::Haplotype:
/*std::cout << "H1" << std::endl << (*this) << std::endl;*/
/*std::cout << "H2" << std::endl << other << std::endl;*/
return Genotype(genotype.first.first, other.genotype.first.first,
genotype.second.first, other.genotype.second.first,
poly * other.poly);
......@@ -636,10 +640,10 @@ struct algebraic_genotype { /* must fit in a QWORD */
case Type::Cross:
switch (other.type) {
case Type::Genotype:
return Genotype(other.genotype.first % (bool)genotype.first.first,
other.genotype.first % (bool)genotype.first.second,
other.genotype.second % (bool)genotype.second.first,
other.genotype.second % (bool)genotype.second.second,
return Genotype(other.genotype.first % (bool)genotype.first.first.tag,
other.genotype.first % (bool)genotype.first.second.tag,
other.genotype.second % (bool)genotype.second.first.tag,
other.genotype.second % (bool)genotype.second.second.tag,
poly * other.poly);
case Type::Haplotype: throw error("Unsupported Cross*Haplotype");
case Type::Gamete: throw error("Unsupported Cross*Gamete");
......@@ -833,6 +837,8 @@ struct eval_poly_optim {
Matrix<int, Dynamic, Dynamic> rebuild;
size_t pow_max;
/*GenoMatrix debug_G;*/
eval_poly_optim()
: indices(), polynoms(), rebuild(), pow_max(0)
{}
......@@ -858,6 +864,8 @@ struct eval_poly_optim {
}
pow_max = polynoms[0].degree() + 1;
/*std::cout << "There are " << uniq.size() << " dictinct polynoms in " << (G.innerSize() * G.outerSize()) << " cells. pow_max=" << pow_max << std::endl;*/
/*debug_G = G;*/
/*std::cout << debug_G << std::endl;*/
/*std::cout << rebuild << std::endl;*/
}
......@@ -870,6 +878,7 @@ struct eval_poly_optim {
cache(const eval_poly_optim* et)
: m_epo(et), m_cache(), r_pow(), s_pow(), values()
{
/*std::cout << "new cache for" << std::endl << m_epo->debug_G << std::endl;*/
values.resize(m_epo->polynoms.size());
r_pow.resize(m_epo->pow_max);
s_pow.resize(m_epo->pow_max);
......@@ -890,6 +899,7 @@ struct eval_poly_optim {
cache& operator = (cache&& tmp)
{
m_epo = tmp.m_epo;
m_cache.swap(tmp.m_cache);
r_pow.swap(tmp.r_pow);
s_pow.swap(tmp.s_pow);
values.swap(tmp.values);
......@@ -933,6 +943,10 @@ struct eval_poly_optim {
}
}
/*std::cout << "Evaluating" << std::endl << m_epo->debug_G << std::endl;*/
/*std::cout << "with X=" << dist << std::endl;*/
/*std::cout << ret << std::endl;*/
return ret;
}
};
......
......@@ -160,12 +160,18 @@ operator & (cache_output& c, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _Max
sz = mat.outerSize();
/*MSG_DEBUG("outer size " << sz);*/
c & sz;
_Scalar* data = mat.data();
for (size_t i = 0; i < mat.size(); ++i) {
c & data[i];
}
#if 0
for (int j = 0; j < mat.outerSize(); ++j) {
for (int i = 0; i < mat.innerSize(); ++i) {
/*MSG_DEBUG('(' << i << ',' << j << ") " << mat(i, j));*/
c & mat(i, j);
}
}
#endif
return c;
}
......@@ -174,11 +180,16 @@ cache_input&
operator & (cache_input& c, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& mat)
{
int sz1, sz2;
_Scalar s;
c & sz1;
c & sz2;
/*MSG_DEBUG(__FILE__ << ':' << __LINE__ << '\t' << __func__);*/
mat.resize(sz1, sz2);
_Scalar* data = mat.data();
for (size_t i = 0; i < mat.size(); ++i) {
c & data[i];
}
#if 0
_Scalar s;
/*MSG_DEBUG(__FILE__ << ':' << __LINE__ << '\t' << __func__);*/
for (int j = 0; j < mat.outerSize(); ++j) {
for (int i = 0; i < mat.innerSize(); ++i) {
......@@ -187,6 +198,7 @@ operator & (cache_input& c, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxR
mat(i, j) = s;
}
}
#endif
return c;
}
......@@ -219,5 +231,20 @@ cache_file<DIRECTION>& operator & (cache_file<DIRECTION>& cf, labelled_matrix<M,
template <typename DIRECTION>
cache_file<DIRECTION>& operator & (cache_file<DIRECTION>& cf, ancestor_allele& a)
{
return cf & a.ancestor & a.tag;
}
template <typename DIRECTION>
cache_file<DIRECTION>& operator & (cache_file<DIRECTION>& cf, allele_pair& a)
{
return cf & a.first & a.second;
}
#endif
......@@ -9,10 +9,11 @@
#include <tuple>
#include <future>
#include <type_traits>
#include <openssl/evp.h>
#include <vector>
#include <map>
#include "labelled_matrix.h"
namespace detail {
template<typename T>
struct has_const_iterator
......@@ -56,6 +57,10 @@ namespace detail {
&& has_begin_end<T>::end_value> {};
}
#ifdef USE_OPENSSL_MD5
#include <openssl/evp.h>
struct md5_digest {
EVP_MD_CTX md_ctx;
bool blocked;
......@@ -188,6 +193,109 @@ struct md5_digest {
}
};
#else
struct md5_digest {
typedef size_t context_type;
size_t context;
bool blocked;
md5_digest()
{
context = 0xdeadb33f;
blocked = false;
}
operator std::string ()
{
std::stringstream s;
s << std::hex << context;
return s.str();
static unsigned char hex[16] = {'0', '1', '2', '3', '4', '5', '6', '7',
'8', '9', 'a', 'b', 'c', 'd', 'e', 'f'};
std::string buf;
buf.resize(sizeof(context_type) * 2);
auto tmp = context;
for (size_t i = 0; i < (2 * sizeof(context_type)); ++i) {
buf[i] = hex[tmp & 0xF];
tmp >>= 4;
}
return buf;
}
#define MAGIC_CONSTANT 0x0d853211
md5_digest& blend(size_t other_context)
{
#define HALFBITS (sizeof(context_type) * 4)
auto tmp = (context * MAGIC_CONSTANT);
context = (tmp >> HALFBITS) | (tmp << HALFBITS);
/*context |= tmp;*/
context ^= other_context;
return *this;
}
template <typename Whatever>
md5_digest& update(Whatever x)
{
return blend(std::hash<Whatever>()(x));
}
template <typename _Iterator>
md5_digest& update(_Iterator begin, _Iterator end)
{