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

Added tests.

parent b6b6ba6e
GCC_VERSION=-4.9
CXX=g++$(GCC_VERSION)
CXXARGS=--std=gnu++0x -Wall -Wno-unused-local-typedefs -DEIGEN_NO_DEPRECATED_WARNING
CXXARGS=--std=gnu++0x -Wall -Wno-unused-local-typedefs -DEIGEN_NO_DEPRECATED_WARNING -pthread
COV=gcov$(GCC_VERSION)
INC=-I../include -I../include/input -I/usr/local/include/eigen3
LD=$(CXX)
LIBS=-lexpat
LIBS=-rdynamic -lexpat -ldl
VERSION_DEFINES=$(shell git tag | sed 's/\([0-9]*\)[.]\([0-9]*\)[.]\([0-9]*\)/-DVERSION_MAJOR=\\\"\1\\\" -DVERSION_MINOR=\\\"\2\\\" -DVERSION_PATCH=\\\"\3\\\"/')
C=$(CXX) $(CXXARGS) $(INC) -DTEST_POLY $(VERSION_DEFINES)
SRC=main.cc fast_polynom.cc markov_population.cc ../src/input/read_map.cc ../src/outbred.cc ../src/probapop_dtd.cc ../src/input/invoke_probapop.cc ../src/static_data.cc
SRC=main.cc cache.cc fast_polynom.cc markov_population.cc ril_properties.cc
EXT_SRC=../src/input/read_map.cc ../src/input/read_mark.cc ../src/input/read_trait.cc ../src/outbred.cc ../src/probapop_dtd.cc ../src/input/invoke_probapop.cc ../src/static_data.cc ../src/input/read_settings.cc ../src/input/xml/xml_design.cc ../src/input/xml/xml_format.cc ../src/input/xml/xml_settings.cc ../src/input/ld_matrices.cc ../src/input/pedigree.cc ../src/computations/basic_data.cc
#SRC=main.cc fast_polynom.cc ril_properties.cc ../src/input/read_map.cc ../src/outbred.cc ../src/probapop_dtd.cc ../src/input/invoke_probapop.cc ../src/static_data.cc
OBJ=$(subst .cc,.o,$(SRC))
EXT_OBJ=$(subst ../,,$(subst .cc,.o,$(EXT_SRC)))
DEBUG_OPTS=-ggdb
OPT_OPTS=-O3 -DEIGEN_NO_DEBUG -DNDEBUG
......@@ -25,19 +28,21 @@ INFO_FILE=test-coverage.info
#C=$(CXX) $(CXXARGS) $(INC) $(DEBUG_OPTS)
C=$(CXX) $(CXXARGS) $(INC) $(DEBUG_OPTS) $(PROF_OPTS) $(VERSION_DEFINES) $(COV_OPTS)
C=$(CXX) $(CXXARGS) $(INC) $(OPT_OPTS) $(VERSION_DEFINES)
C=$(CXX) $(CXXARGS) $(INC) $(OPT_OPTS) $(DEBUG_OPTS) $(VERSION_DEFINES)
all: test_suite
include .depend
PCH=catch.hpp.gch
pch: $(PCH)
$(PCH):%.h.gch: %.h Makefile
$(PCH):%.hpp.gch: %.hpp Makefile
$C $<
test_suite: $(OBJ)
$C $(OBJ) -o $@ $(LIBS)
test_suite: $(OBJ) $(EXT_OBJ)
$C $(OBJ) $(EXT_OBJ) -o $@ $(LIBS)
test: test_suite
./$<
......@@ -56,7 +61,10 @@ test-coverage: zerocounters test
.depend: $(SRC) Makefile
$C --depend $(SRC) > $@
$(OBJ):%.o: %.cc
$(OBJ):%.o: %.cc $(PCH)
$C -c $< -o $@
$(EXT_OBJ):%.o: ../%.cc $(PCH)
$C -c $< -o $@
clean:
......
#include "catch.hpp"
#include "error.h"
#include "settings.h"
#include "labelled_matrix.h"
#include "cache2.h"
#include "computations/base.h"
#include "computations/basic_data.h"
#include "computations/probabilities.h"
template <typename X>
size_t md_value(const X& x)
{
md5_digest md;
md << x;
std::string s = md;
/*std::cout << "hash(" << x << ") = " << s << std::endl;*/
return md.context;
}
static inline
std::ostream& operator << (std::ostream& os, const multi_generation_observations& mgo)
{
return os << "MGO[" << mgo.n_mark << ']' << std::endl << mgo.LV;
}
int mult(int a, int b) { MSG_DEBUG(a << " * " << b << " = " << (a * b)); return a * b; }
int mmult(int a, int b, const multi_generation_observations&) { MSG_DEBUG(a << " * " << b << " = " << (a * b)); return a * b; }
const char* bc_conf[] = {
"",
"-n", "test",
"-P", "2",
"-bds", "../data/design-BC.xml",
"-gm", "dataset/debug.map",
"-mos", "../data/format-ABHCD.xml",
"-p", "toto", "BC", "dataset/debug_BC.phen",
"-m", "BC", "dataset/debug_BC.gen",
"step", "1",
"-na",
"-wd", "./tmp",
NULL};
struct settings_fixture {
settings_fixture(const char** argv)
{
int argc = 0;
for (; argv[argc]; ++argc);
active_settings = settings_t::from_args(argc, argv);
/*REQUIRE(active_settings != NULL);*/
if (active_settings) {
active_settings->finalize();
} else {
/*throw bad_settings_exception("Couldn't initialize settings.");*/
MSG_ERROR("Couldn't initialize settings.", "");
}
}
~settings_fixture() { delete active_settings; active_settings = NULL; }
};
TEST_CASE("hash.1", "check hash of some single values")
{
settings_fixture _(bc_conf);
/*multi_generation_observations mgo1, mgo2;*/
/*mgo1.n_mark = 3;*/
/*mgo2.n_mark = 3;*/
/*mgo1.pop = &active_settings->populations[0];*/
/*mgo2.pop = &active_settings->populations[0];*/
CHECK(md_value(0) != md_value(1));
CHECK(md_value(std::string("toto")) != md_value(std::string("titi")));
MatrixXb m1(3, 3);
m1 << 1, 0, 0,
0, 1, 0,
1, 1, 1;
MatrixXb m2(3, 3);
m2 << 1, 0, 1,
1, 1, 0,
1, 1, 1;
CHECK(md_value(m1) != md_value(m2));
labelled_matrix<MatrixXb, char, char> lm1, lm2;
lm1.data = m1;
lm2.data = m2;
lm1.column_labels = {'a', 'b', 'c'};
lm2.column_labels = {'a', 'b', 'c'};
lm1.row_labels = {'a', 'b', 'c'};
lm2.row_labels = {'a', 'b', 'c'};
CHECK(md_value(lm1) != md_value(lm2));
/*mgo1.LV = m1; mgo2.LV = m2;*/
/*CHECK(md_value(mgo1) != md_value(mgo2));*/
lm2.data = m1;
CHECK(md_value(lm1) == md_value(lm2));
/*mgo1.LV = m1; mgo2.LV = m1;*/
/*CHECK(md_value(mgo1) != md_value(mgo2));*/
lm2.column_labels[1] = 'z';
CHECK(md_value(lm1) != md_value(lm2));
}
TEST_CASE("hash.2", "check hash of other single values")
{
settings_fixture _(bc_conf);
qtl_chromosome qc1(&active_settings->map[0]);
qtl_chromosome qc2(&active_settings->map[0]);
CHECK(md_value(qc1) == md_value(qc2));
qc2.add_qtl({0});
CHECK(md_value(qc1) != md_value(qc2));
qc1.add_qtl({0});
CHECK(md_value(qc1) == md_value(qc2));
qc1.remove_qtl({0});
qc1.add_qtl({2});
/*CHECK(md_value(qc1) != md_value(qc2));*/
CHECK(md_value(qc1) == md_value(qc2)); /* BY DESIGN! only the number of selected QTLs matter, not their positions! */
}
TEST_CASE("hash.3", "check hash of multiple arguments as in parental_origin2")
{
settings_fixture _(bc_conf);
qtl_chromosome qc(&active_settings->map[0]);
population* pop = &active_settings->populations.begin()->second;
generation_rs* gen = active_settings->design->generation["BC"];
labelled_matrix<Eigen::MatrixXd, std::vector<char>, allele_pair>& stf = gen->state_to_parental_origin[0];
multi_generation_observations mgo1(pop, qc.chr, gen, 0);
multi_generation_observations mgo2(pop, qc.chr, gen, 1);
std::vector<double> sqoc;
std::vector<double> loci = {0, 1};
CHECK(md_value(mgo1) != md_value(mgo2));
CHECK(new_redux::md5(mgo1) != new_redux::md5(mgo2));
CHECK(new_redux::md5(&qc, pop, gen, stf, mgo1, sqoc, loci) == new_redux::md5(qc, pop, gen, stf, mgo1, sqoc, loci));
CHECK(new_redux::md5(&qc, pop, gen, stf, mgo1, sqoc, loci) != new_redux::md5(&qc, pop, gen, stf, mgo2, sqoc, loci));
}
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,
const std::vector<double>& loci)
{
return {};
}
using cc = cached_computation<decltype(parental_origin2)>;
using ccv = cached_computed_value<decltype(parental_origin2)>;
TEST_CASE("hash.4", "check behaviour of cached_computation<parental_origin2>")
{
settings_fixture _(bc_conf);
qtl_chromosome qc_(&active_settings->map[0]);
value<qtl_chromosome_value> qc = &qc_;
value<population_value> pop {&active_settings->populations.begin()->second};
value<generation_value> gen {active_settings->design->generation["BC"]};
value<labelled_matrix<Eigen::MatrixXd, std::vector<char>, allele_pair>> stf = (*gen)->state_to_parental_origin[0];
value<multi_generation_observations> mgo1 = multi_generation_observations(*pop, (*qc)->chr, *gen, 0);
value<multi_generation_observations> mgo2 = multi_generation_observations(*pop, (*qc)->chr, *gen, 1);
value<selected_qtls_on_chromosome_type> sqoc = selected_qtls_on_chromosome_type{};
value<std::vector<double>> loci = std::vector<double>{0, 1};
CHECK(md_value(mgo1) != md_value(mgo2));
CHECK(new_redux::md5(mgo1) != new_redux::md5(mgo2));
CHECK(new_redux::md5(qc, pop, gen, stf, mgo1, sqoc, loci) != new_redux::md5(qc, pop, gen, stf, mgo2, sqoc, loci));
/*value<parental_origin_type> pot1 = make_value(parental_origin2, qc, pop, gen, stf, mgo1, sqoc, loci);*/
/*value<parental_origin_type> pot2 = make_value(parental_origin2, qc, pop, gen, stf, mgo2, sqoc, loci);*/
cc po1("parental_origin2", parental_origin2, qc, pop, gen, stf, mgo1, sqoc, loci);
cc po2("parental_origin2", parental_origin2, qc, pop, gen, stf, mgo2, sqoc, loci);
CHECK(po1.m_md5_hash.accum != po2.m_md5_hash.accum);
CHECK(po1.m_md5_hash.append != po2.m_md5_hash.append);
CHECK(po1.get_path() != po2.get_path());
MSG_DEBUG("po1 path: " << po1.get_path());
MSG_DEBUG("po2 path: " << po2.get_path());
ccv pov1(CachingPolicy::Disk, parental_origin2, qc, pop, gen, stf, mgo1, sqoc, loci);
ccv pov2(CachingPolicy::Disk, parental_origin2, qc, pop, gen, stf, mgo2, sqoc, loci);
CHECK(pov1.m_comp.m_md5_hash.accum != pov2.m_comp.m_md5_hash.accum);
CHECK(pov1.m_comp.m_md5_hash.append != pov2.m_comp.m_md5_hash.append);
CHECK(pov1.m_comp.get_path() != pov2.m_comp.get_path());
}
TEST_CASE("hash.5", "check hashing of computed values")
{
settings_fixture _(bc_conf);
value<population_value> pop {&active_settings->populations.begin()->second};
value<const pop_mgo_data*> pmd = new pop_mgo_data(*pop);
value<chromosome_value> chr = &active_settings->map[0];
value<multi_generation_observations>
mv1 = make_value<Sync|Disk>(population_marker_obs,
pop, chr, value<int>{0}, pmd);
value<multi_generation_observations>
mv2 = make_value<Sync|Disk>(population_marker_obs,
pop, chr, value<int>{1}, pmd);
auto tmp = mv1->LV;
tmp = mv2->LV;
/*std::cout << "mv1" << std::endl << (*mv1) << std::endl;*/
/*std::cout << "mv2" << std::endl << (*mv2) << std::endl;*/
CHECK(mv1 != mv2);
CHECK(mv1.hash() != mv2.hash());
CHECK(new_redux::md5(mv1) != new_redux::md5(mv2));
value<multi_generation_observations>
mv1bis = make_value<Disk>(population_marker_obs,
pop, chr, value<int>{0}, pmd);
CHECK(mv1.hash() == mv1bis.hash());
CHECK(new_redux::md5(mv1, mv2) != new_redux::md5(mv1, mv1));
}
/media/Stock/devel/MCQTL/v6/src/spell-qtl -n debug -P auto -wd debug_tmp_1384347995340803 -bds /media/Stock/devel/MCQTL/v6/data/design-BC.xml -gm debug.map -mos /media/Stock/devel/MCQTL/v6/data/format-ABHCD.xml -p toto BC debug_BC.phen -m BC debug_BC.gen step 1 estimate-at ch1:86.31 -na -wd ./tmp
2013-11-13 14:06:54.834235
POP A 1
POP F1 1
POP B 1
POP BC 100
QTL 86.31
trait trait1 qtl [('M_1_9', 86.31)] effect [{'aa': 1.7321, 'ba': 0.0}] epistasis [0.0] noise 2.0 0.0 1.0
*ch1 10 M_1_1 1.93 M_1_2 44.55 M_1_3 1.2 M_1_4 2.49 M_1_5 17.54 M_1_6 6.65 M_1_7 10.67 M_1_8 15.01 M_1_10 48.68 M_1_11
generation individual mother father
F1 0 0 0
BC 0 0 0
BC 1 0 0
BC 2 0 0
BC 3 0 0
BC 4 0 0
BC 5 0 0
BC 6 0 0
BC 7 0 0
BC 8 0 0
BC 9 0 0
BC 10 0 0
BC 11 0 0
BC 12 0 0
BC 13 0 0
BC 14 0 0
BC 15 0 0
BC 16 0 0
BC 17 0 0
BC 18 0 0
BC 19 0 0
BC 20 0 0
BC 21 0 0
BC 22 0 0
BC 23 0 0
BC 24 0 0
BC 25 0 0
BC 26 0 0
BC 27 0 0
BC 28 0 0
BC 29 0 0
BC 30 0 0
BC 31 0 0
BC 32 0 0
BC 33 0 0
BC 34 0 0
BC 35 0 0
BC 36 0 0
BC 37 0 0
BC 38 0 0
BC 39 0 0
BC 40 0 0
BC 41 0 0
BC 42 0 0
BC 43 0 0
BC 44 0 0
BC 45 0 0
BC 46 0 0
BC 47 0 0
BC 48 0 0
BC 49 0 0
BC 50 0 0
BC 51 0 0
BC 52 0 0
BC 53 0 0
BC 54 0 0
BC 55 0 0
BC 56 0 0
BC 57 0 0
BC 58 0 0
BC 59 0 0
BC 60 0 0
BC 61 0 0
BC 62 0 0
BC 63 0 0
BC 64 0 0
BC 65 0 0
BC 66 0 0
BC 67 0 0
BC 68 0 0
BC 69 0 0
BC 70 0 0
BC 71 0 0
BC 72 0 0
BC 73 0 0
BC 74 0 0
BC 75 0 0
BC 76 0 0
BC 77 0 0
BC 78 0 0
BC 79 0 0
BC 80 0 0
BC 81 0 0
BC 82 0 0
BC 83 0 0
BC 84 0 0
BC 85 0 0
BC 86 0 0
BC 87 0 0
BC 88 0 0
BC 89 0 0
BC 90 0 0
BC 91 0 0
BC 92 0 0
BC 93 0 0
BC 94 0 0
BC 95 0 0
BC 96 0 0
BC 97 0 0
BC 98 0 0
BC 99 0 0
data type BC
100 11 0 0
*M_1_1 AAHAHHAAAAAHHHHHHHAAHAHAHAAAHHAAHAHHHHHHAAAHAAAHAAAAHHHHAHAHAAHHAHHHAHAHHHHAAAAHHHHAAAAAAHHHAAHAHAAH
*M_1_2 AAHAHHAAAAAHHHHHHHAAHAHAHAAAHHAAHAHHHHHHAAAHAAAHAAAAHHHHAHAHAAHHAHHHAHAHHHHAAAAHHAHAAAAAAHHHAAHAHAAH
*M_1_3 HAAHHHAAAAHHAHHHAHHAAAAAHAAHHHHAHAHAHHHHAHAHHAHHAAHHHAHAHHAHAAHHAHAHHHAHHHHHAAAAHAHAAHAAAHHHAHAHAAAH
*M_1_4 HAAHHHAAAAHHAHHHAHHAAAAAHAAHHHHAHAHAHHHHAHAHHAHHAAHHHAHAHHAHAAHHAHAHHHAHHHHHAHAAHAHAAHAAAHHHAHAHAAAH
*M_1_5 HHAHHHAAAAHHAHHHAHHAAAAAHAAHHHHAHAHAHHHHAHAHHAHHAAHHHAHAHHAHAAHHAHAHHHAHHHHHAHAAHAHAAHAAAHHHAHAHAAAH
*M_1_6 HHHAHHAAAAHAAHHHAHHAAAAAHAAHHAHAHAAAAHAHAHAHHAHAAAHAHAAAHHAHAAHHAHAHHHAHHAHHAHAAHHHAAHAAAHAHAHAHAAHH
*M_1_7 HHHAHHAAAAHAHHHHAAHAAAAAHAAHHAHAHAAAAHAHAHAHHAHAAAHAHAAHHHAHAAHHAHAHHHAAHAHHAHAAHHHAAHAAAHAHAHAHHAHH
*M_1_8 HHHAHHAAAAHAHHHHAAAAAAAHHAAHAHHHAAAAAHHHAAAHHAHAAAHHHAAHHHAHAAHHAHAHHHAAHHHHAHAAHAHAHHAAAHAHHHAHHAAH
*M_1_9 HHHAHHAAAAAAHHHHAAAAAAAHHAAHAHHAAAAAAHHHAAAHHAHAAAHHHAAHHHAHAHHAAHAHHHAAHHHHAHAAHAHAHHAAAHHHHHAHHAAH
*M_1_10 HHHAAAAAAAAAHHHAAAAHAAAHHAAHAHHAHAAAHHAHAAAHHHAAAAHHHAHHAHAHAHHAAAAHHHHHHHHHAHAAHHHAHHAAAHAHAHHHHAAH
*M_1_11 HAHHHAHHAAAAAHHAAAAHAHHHHHHAAHAHHAHAAHAHHAAHHAAHAHHHHAHAHAAHAHHHAAAHHHAHHHAHAHAHHAHHHAAHHHAAAHHHHHAA
*trait1 1.18418792626 0.337700787437 2.28175624947 3.52322209027 -0.596558588708 -1.27862920958 0.938807892177 -1.24734101796 2.06076472265 2.09547233654 3.03461291362 0.778489982715 3.51622498085 3.1628990258 2.21428881317 1.12925136238 4.25315792042 1.50986086539 1.93185537356 2.44422225829 0.96569305948 0.550584219719 4.73423669691 1.96372287718 1.19501869 6.07932545183 5.60723843856 0.640709467987 1.85905547608 -2.51766234814 -1.24612905148 0.894894135603 3.18027875926 1.2576676374 3.3847955256 1.78456555031 0.756576254556 -0.899964726343 0.178917855055 1.01123237728 -0.470942872517 1.14406477304 -0.825023762995 -0.84515682456 0.197553608024 1.2600915796 1.54617811294 3.39266678402 5.18070564821 1.60608795457 0.793204409739 -1.34516397857 3.9927020874 6.53174458486 -1.32178471858 2.61298340274 -2.19115223868 2.02016169565 -0.639591460867 0.25743845654 2.24476865705 -0.540719847669 0.409164362921 3.79261998731 2.16615055057 -0.556558794589 3.90008017603 -2.95097889782 3.21518999608 0.24982658287 2.85539713678 -0.358847480116 -1.35620536322 2.52672122078 -2.18119457675 -1.52075080895 -3.9724640144 -0.314766291607 -0.0221562775252 -0.0599760297807 0.994936518994 1.81926934364 -5.51535527832 1.20469783476 -0.275270416186 1.27424725198 2.51467013001 2.64374873239 0.904250154527 -0.344742928516 -1.86072082471 -2.01093102907 -2.13110752049 -0.0152279233786 0.932986296956 1.60359650832 -3.49728195736 0.31245124448 1.40215725654 0.826546862988
......@@ -34,7 +34,7 @@ TEST_CASE("degree_and_valuation", "Check some wrappers")
CHECK(fzero.valuation() == 0);
CHECK(fr.valuation() == 1);
CHECK(fs.valuation() == 0);
std::cout << p << std::endl;
/*std::cout << p << std::endl;*/
}
......@@ -80,8 +80,8 @@ TEST_CASE("basic.add", "Check simple additions of axiomatic polynoms")
impl::polynom_tables::reset();
/*std::cout << "Add table" << std::endl << impl::polynom_tables::add_table() << std::endl;*/
fast_polynom _2s = {2, -2};
std::cout << ((const impl::d_polynom&) _2s) << std::endl;
std::cout << ((const impl::f_polynom&) _2s) << std::endl;
/*std::cout << ((const impl::d_polynom&) _2s) << std::endl;*/
/*std::cout << ((const impl::f_polynom&) _2s) << std::endl;*/
/*std::cout << "Add table" << std::endl << impl::polynom_tables::add_table() << std::endl;*/
fast_polynom _2r = {0, 2};
/*std::cout << "Add table" << std::endl << impl::polynom_tables::add_table() << std::endl;*/
......@@ -222,9 +222,9 @@ TEST_CASE("check.roots", "Check some polynom root-finding")
fast_polynom a = {0, 0, 0, 1};
fast_polynom b = {1, -1};
fast_polynom c = a * b;
std::cout << "c = " << c << std::endl;
/*std::cout << "c = " << c << std::endl;*/
impl::f_polynom fc = c;
std::cout << "fc = " << fc << std::endl;
/*std::cout << "fc = " << fc << std::endl;*/
CHECK(fc.r_exp == 3);
CHECK(fc.s_exp == 1);
CHECK(fc.P.size() == 1);
......
#include "catch.hpp"
#include "markov_population.h"
#include "input/read_map.h"
#include "input/read_mark.h"
#include "input.h"
#include "input/invoke_probapop.h"
#include "probapop_dtd.h"
#include <exception>
......@@ -12,7 +11,7 @@
#include "chrono.h"
#include "data/genoprob_computer.h"
#include "invoke_probapop.h"
/*#include "invoke_probapop.h"*/
Eigen::IOFormat genomatrix_fmt(/*precision*/StreamPrecision, /*flags*/DontAlignCols,
......@@ -123,12 +122,12 @@ std::pair<chromosome, marker_data> from_argv_outbred(int argc, const char* const
generation_rs* select_pop(const std::string& pop)
{
generation_rs* A = generation_rs::ancestor("A", 'a');
generation_rs* B = generation_rs::ancestor("B", 'b');
generation_rs* A = generation_rs::ancestor("A", 'a', {{0, 0}});
generation_rs* B = generation_rs::ancestor("B", 'b', {{0, 0}});
/*std::cout << "Pop type: " << pop << std::endl;*/
if (pop == "cp") {
generation_rs* C = generation_rs::ancestor("C", 'c');
generation_rs* D = generation_rs::ancestor("D", 'd');
generation_rs* C = generation_rs::ancestor("C", 'c', {{0, 0}});
generation_rs* D = generation_rs::ancestor("D", 'd', {{0, 0}});
return A->crossing("AB", B)->crossing("CP", C->crossing("CD", D));
}
generation_rs* F1 = A->crossing("F1", B);
......@@ -142,7 +141,7 @@ generation_rs* select_pop(const std::string& pop)
return F1->to_doubled_haploid("DH");
}
if (pop == "ril") {
return F1->to_doubled_haploid("DH")->to_ril("RIL");
return F1->to_ril("RIL");
}
if (pop[0] == 'f') {
std::stringstream s(pop);
......@@ -244,7 +243,7 @@ void check_delta(const std::string& pop_type, double locus, double error, double
int test_segment(int argc, const char* const* argv) {
std::cout << "Testing pop_type=" << argv[2] << ": " << argv2str(argc - 3, argv + 3) << std::endl;
/*std::cout << "Testing pop_type=" << argv[2] << ": " << argv2str(argc - 3, argv + 3) << std::endl;*/
std::vector<std::pair<char, std::vector<allele_pair>>> obs_spec = {
{'A', {{'a', 'a'}}},
{'B', {{'b', 'b'}}},
......@@ -294,7 +293,11 @@ int test_segment(int argc, const char* const* argv) {
/*auto output = pop.compute_over_segment(chr1.marker_locus, obs, obs_spec, step, 0).transpose();*/
qtl_chromosome qc(&chr1);
generation_rs::segment_computer_t sc = pop.segment_computer(&qc, step, 0);
auto output = sc.compute(obs, obs_spec).transpose();
/*auto output = sc.compute(obs, obs_spec).transpose();*/
pop.compute_observation_vectors(obs_spec);
auto LV = pop.raw_observations(pop.observation_vectors, obs);
auto output = sc.compute(LV).transpose();
/*std::cout << argv[2] << " LABELS " << output.column_labels << std::endl;*/
if (argv[2][0] == 'f' && argv[2][1] != '2') {
output.data *= (MatrixXd(4, 4) << 1, 0, 0, 0,
......@@ -319,7 +322,7 @@ int test_segment(int argc, const char* const* argv) {
int test_segment_outbred(int argc, const char* const* argv) {
std::cout << "Testing pop_type=" << argv[2] << ": " << argv2str(argc - 3, argv + 3) << std::endl;
/*std::cout << "Testing pop_type=" << argv[2] << ": " << argv2str(argc - 3, argv + 3) << std::endl;*/
allele_pair a0 = {'a', 'c'};
allele_pair a1 = {'a', 'd'};
allele_pair a2 = {'b', 'c'};
......@@ -363,8 +366,12 @@ int test_segment_outbred(int argc, const char* const* argv) {
/*auto output = pop.compute_over_segment(chr1.marker_locus, obs, obs_spec, step, 0).transpose();*/
qtl_chromosome qc(&chr1);
generation_rs::segment_computer_t sc = pop.segment_computer(&qc, step, 0);
auto output = sc.compute(obs, obs_spec).transpose();
std::cout << "cp LABELS " << output.column_labels << std::endl;
/*auto output = sc.compute(obs, obs_spec).transpose();*/
pop.compute_observation_vectors(obs_spec);
auto LV = pop.raw_observations(pop.observation_vectors, obs);
auto output = sc.compute(LV).transpose();
/*std::cout << "cp LABELS " << output.column_labels << std::endl;*/
/*std::cout << output.upperRows(10) << std::endl;*/
output.data *= (MatrixXd(4, 4) << 1, 0, 0, 0,
0, 0, 1, 0,
......@@ -483,18 +490,19 @@ TEST_CASE("markov_3V", "Check equality with old Probapop with 3V")
{'9', {{'c', 'c'}}}
};
generation_rs* A = generation_rs::ancestor("A", 'a');
generation_rs* B = generation_rs::ancestor("B", 'b');
generation_rs* C = generation_rs::ancestor("C", 'c');
generation_rs* F1 = A->crossing("AB_3V", B)->crossing("F1_3V", C);
generation_rs* F2 = F1->selfing("F2_3V");
/*generation_rs* A = generation_rs::ancestor("A", 'a', {{0, 0}});*/
/*generation_rs* B = generation_rs::ancestor("B", 'b', {{0, 0}});*/
/*generation_rs* C = generation_rs::ancestor("C", 'c', {{0, 0}});*/
/*generation_rs* F1 = A->crossing("AB_3V", B)->crossing("F1_3V", C);*/
/*generation_rs* F2 = F1->selfing("F2_3V");*/
/*std::cout << (*F1) << std::endl << (*F2) << std::endl;*/
std::cout << "A " << A << std::endl;
std::cout << "B " << B << std::endl;
std::cout << "C " << C << std::endl;
std::cout << "F1 " << F1 << std::endl;
std::cout << "F2 " << F2 << std::endl;
/*std::cout << "A " << A << std::endl;*/
/*std::cout << "B " << B << std::endl;*/
/*std::cout << "C " << C << std::endl;*/
/*std::cout << "F1 " << F1 << std::endl;*/
/*std::cout << "F2 " << F2 << std::endl;*/
#if 0
multi_generation_observations mgo(5);
chromosome chr1 = {
"chr1",
......@@ -540,6 +548,7 @@ TEST_CASE("markov_3V", "Check equality with old Probapop with 3V")
for (size_t i = 0; i < 65; ++i) {
check_delta("3V", i, delta(out.row(i), pp.row(i)), 1.e-4, 1, args);
}
#endif
}
......
#include "catch.hpp"
#include "error.h"
#include "input.h"
#include "generation_rs.h"
bool compare_processes(const generation_rs* g1, const generation_rs* g2)
{
const auto& p1 = g1->P;
const auto& p2 = g2->P;
/*MSG_DEBUG("g1: " << (*g1));*/
/*MSG_DEBUG("g2: " << (*g2));*/
REQUIRE(p1.size() == p2.size());
for (size_t i = 0; i < p1.size(); ++i) {
REQUIRE(p1[i].weight == p2[i].weight);
REQUIRE(p1[i].G == p2[i].G);
}
return true;
}
generation_rs*
rilify(const generation_rs* g, const char* n1, const char* n2)
{
return g->to_ril(n2);
(void)n1;
/*return g->to_doubled_haploid(n1)->to_ril(n2);*/
/*return g->to_ril(n2)->to_doubled_haploid(n1);*/
}
/* Line */
TEST_CASE("RILoA", "RIL(A) = A")
{
generation_rs::clear_dict();
generation_rs* A = generation_rs::ancestor("A", 'a', {{0, 0}});
generation_rs* RIL = rilify(A, "_", "RIL");
compare_processes(A, RIL);
}