Commit 39463cf5 authored by Damien Leroux's avatar Damien Leroux
Browse files

Spell-marker working. Integration with spell-qtl WIP.

parent a0c76b0a
......@@ -11,7 +11,7 @@
<generation name="BC1">
<cross p1="F1" p2="A"/>
</generation>
<generation name="BC-RIL">
<ril p1="BC1"/>
<generation name="BCS">
<self p1="BC1"/>
</generation>
</breeding-design>
<format-specification>
<marker-observations name="02">
<observation symbol="0">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="0" from-father="2"/>
<alleles from-mother="2" from-father="0"/>
</observation>
<observation symbol="2">
<alleles from-mother="2" from-father="2"/>
<alleles from-mother="0" from-father="2"/>
<alleles from-mother="2" from-father="0"/>
</observation>
<observation symbol="-">
<alleles from-mother="0" from-father="0"/>
<alleles from-mother="2" from-father="2"/>
<alleles from-mother="0" from-father="2"/>
<alleles from-mother="2" from-father="0"/>
</observation>
</marker-observations>
<format name="02">
<domain type="allele" alphabet-from="02"/>
<scores>
<score symbol="0">
<from mother="0" father="0"/>
</score>
<score symbol="2">
<from mother="2" father="2"/>
</score>
<score symbol="-">
<from mother="0" father="0"/>
<from mother="2" father="2"/>
</score>
</scores>
</format>
</format-specification>
......@@ -371,10 +371,10 @@ struct algebraic_genotype { /* must fit in a QWORD */
static const algebraic_genotype null;
locus_pair genotype;
enum class Type: int {Null, Genotype, Gamete, DoublingGamete, Haplotype, Cross, Polynom, RilPolynom} type;
enum class Type: short {Null, Genotype, Gamete, DoublingGamete, Haplotype, Cross, Polynom, RilPolynom} type;
/*rs_polynom poly;*/
fast_polynom poly;
int __align[2];
/*int __align[2];*/
#ifdef DEBUG_ALGEBRAIC_GENOTYPE
/*std::pair<const algebraic_genotype*, const algebraic_genotype*> history;*/
#endif
......@@ -522,16 +522,16 @@ struct algebraic_genotype { /* must fit in a QWORD */
break;
case algebraic_genotype::Type::DoublingGamete:
case algebraic_genotype::Type::Gamete:
os << (char)('0' + ag.genotype.first.first.tag);
os << (char)('0' + ag.genotype.first.first.ancestor);
os << ',';
os << (char)('0' + ag.genotype.second.first.tag);
os << (char)('0' + ag.genotype.second.first.ancestor);
break;
case algebraic_genotype::Type::Cross:
os << (char)('0' + ag.genotype.first.first.tag);
os << (char)('0' + ag.genotype.first.second.tag);
os << (char)('0' + ag.genotype.first.first.ancestor);
os << (char)('0' + ag.genotype.first.second.ancestor);
os << ',';
os << (char)('0' + ag.genotype.second.first.tag);
os << (char)('0' + ag.genotype.second.second.tag);
os << (char)('0' + ag.genotype.second.first.ancestor);
os << (char)('0' + ag.genotype.second.second.ancestor);
break;
case Type::RilPolynom:
break;
......@@ -589,10 +589,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.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,
return Genotype(other.genotype.first % (bool)genotype.first.first.ancestor,
other.genotype.first % (bool)genotype.first.first.ancestor,
other.genotype.second % (bool)genotype.second.first.ancestor,
other.genotype.second % (bool)genotype.second.first.ancestor,
/*other.poly.o(poly));*/
poly * other.poly);
break;
......@@ -646,22 +646,22 @@ 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.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,
return Genotype(genotype.first % (bool)other.genotype.first.first.ancestor,
genotype.first % (bool)other.genotype.first.first.ancestor,
genotype.second % (bool)other.genotype.second.first.ancestor,
genotype.second % (bool)other.genotype.second.first.ancestor,
/*poly.o(other.poly));*/
poly * other.poly);
case Type::Gamete:
return Haplotype(genotype.first % (bool)other.genotype.first.first.tag,
genotype.second % (bool)other.genotype.second.first.tag,
return Haplotype(genotype.first % (bool)other.genotype.first.first.ancestor,
genotype.second % (bool)other.genotype.second.first.ancestor,
poly * other.poly);
/*break;*/
case Type::Cross:
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,
return Genotype(genotype.first % (bool)other.genotype.first.first.ancestor,
genotype.first % (bool)other.genotype.first.second.ancestor,
genotype.second % (bool)other.genotype.second.first.ancestor,
genotype.second % (bool)other.genotype.second.second.ancestor,
poly * other.poly);
/*break;*/
case Type::Null:
......@@ -681,8 +681,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.tag,
other.genotype.second % (bool)genotype.second.first.tag,
return Haplotype(other.genotype.first % (bool)genotype.first.first.ancestor,
other.genotype.second % (bool)genotype.second.first.ancestor,
poly * other.poly);
case Type::DoublingGamete:
if (genotype == other.genotype) {
......@@ -716,19 +716,19 @@ 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.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,
return Genotype(other.genotype.first % (bool)genotype.first.first.ancestor,
other.genotype.first % (bool)genotype.first.second.ancestor,
other.genotype.second % (bool)genotype.second.first.ancestor,
other.genotype.second % (bool)genotype.second.second.ancestor,
poly * other.poly);
case Type::Haplotype: throw error("Unsupported Cross*Haplotype");
case Type::Gamete: throw error("Unsupported Cross*Gamete");
/*case Type::Cross: throw error("Unsupported Cross*Cross");*/
case Type::Cross:
return Cross(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,
return Cross(other.genotype.first % (bool)genotype.first.first.ancestor,
other.genotype.first % (bool)genotype.first.second.ancestor,
other.genotype.second % (bool)genotype.second.first.ancestor,
other.genotype.second % (bool)genotype.second.second.ancestor,
poly * other.poly);
case Type::Null: throw error("Unsupported Cross*Null");
case Type::DoublingGamete: throw error("Unsupported Cross*DoublingGamete");
......
......@@ -5,6 +5,7 @@
#include "generation_rs_fwd.h"
/*#include "bayes/graph.h"*/
#if 0
static inline
size_t idx_to_ofs(const std::vector<size_t>& dim, const std::vector<size_t>& idx)
{
......@@ -35,7 +36,13 @@ size_t size(const std::vector<size_t>& dim) {
}
return stride;
}
#endif
#include "pedigree.h"
#include "bayes/bn.h"
#include "bayes/cli.h"
#include "bayes/dispatch.h"
#endif
......@@ -494,11 +494,17 @@ struct pop_data_type {
prepend(os, "| ", (*kv.second));
}
os << "| Computed locus vectors:" << std::endl;
prepend(os, "| ", pop_data.LV);
/*for (const auto& kv: pop_data.LV) {*/
/*os << "| - " << kv.first << std::endl;*/
/*prepend(os, "| ", kv.second);*/
/*}*/
/*prepend(os, "| ", pop_data.LV);*/
for (const auto& kv: pop_data.LV.data) {
os << "| Chromosome " << kv.first << std::endl;
for (const auto& gen_lv: kv.second) {
os << "| Generation " << gen_lv.first << std::endl;
for (size_t i = 0; i < gen_lv.second.size(); ++i) {
os << "| #" << i << std::endl;
prepend(os, "| ", gen_lv.second[i]);
}
}
}
return os;
}
......
......@@ -234,7 +234,7 @@ 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;
return cf & a.ancestor;
}
......
......@@ -35,17 +35,6 @@ namespace std {
}
};
template <>
struct hash<multi_generation_observations>
/*: dummy_hash<multi_generation_observations> {};*/
{
size_t operator () (const multi_generation_observations& mgo)
{
hash<decltype(multi_generation_observations::LV)> h;
return h(mgo.LV);
}
};
/*template <typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>*/
/*struct hash<Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>>*/
/*: dummy_hash<Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>> {};*/
......
......@@ -147,10 +147,7 @@ static inline md5_digest& operator << (md5_digest& md5, const pop_mgo_data&)
return md5;
}
/*multi_generation_observations*/
/*population_marker_obs(population_value, chromosome_value, int, const pop_mgo_data*);*/
multi_generation_observations
MatrixXd
population_marker_obs(const context_key&, int);
collection<std::string>
......@@ -173,12 +170,5 @@ cache_output& operator & (cache_output& co, const population*& pop)
return co & pop->name;
}
template <typename CACHE_DIRECTION>
CACHE_DIRECTION& operator & (CACHE_DIRECTION& c, multi_generation_observations& mgo)
{
return c & mgo.n_mark & mgo.LV & mgo.pop;
}
#endif
......@@ -16,10 +16,10 @@
MatrixXd permutation_matrix(size_t n_parents, size_t n_loci, size_t order);
locus_probabilities_type
locus_probabilities(const context_key&, const locus_key&,
const multi_generation_observations&,
const std::vector<double>&);
/*locus_probabilities_type*/
/*locus_probabilities(const context_key&, const locus_key&,*/
/*const multi_generation_observations&,*/
/*const std::vector<double>&);*/
state_to_parental_origin_matrix_type
state_to_parental_origin_matrix(generation_value, int);
......
......@@ -4,8 +4,7 @@
static inline
md5_digest& operator << (md5_digest& md5, const ancestor_allele& aa)
{
md5.update(aa.ancestor);
return md5.update(aa.tag);
return md5.update(aa.ancestor);
}
......
......@@ -354,7 +354,7 @@ struct geno_prob_computer {
}
for (; (l + 1) < r && (*locus)[l + 1] <= loc; ++l);
for (; (*locus)[r - 1] > loc; --r);
for (; r > 0 && (*locus)[r - 1] > loc; --r);
return {l, r};
}
......@@ -426,6 +426,8 @@ struct geno_prob_computer {
/*MSG_DEBUG("LOCUS VECTORS" << std::endl << *LV);*/
labelled_matrix<MatrixXd, allele_pair, double>
ret(this_gen->get_unique_labels(), steps);
/*MSG_DEBUG(MATRIX_SIZE(ret));*/
/*MSG_DEBUG(this_gen->get_unique_labels());*/
size_t l = 0, last = steps.size() - 1, r = !!(LV->outerSize() - 1);
for (size_t i = 0; i < steps.size(); ++i) {
/* FIXME: optimize the CORRECT initialization of l & r */
......@@ -433,6 +435,7 @@ struct geno_prob_computer {
ret.data.col(i) = fast_compute_at(steps[i], l, r, order);
if ((ret.data.col(i).array() < 0).any()) {
MSG_DEBUG("Bad value at column #" << i << std::endl << ret.data.col(i).transpose());
MSG_DEBUG(this_gen->P[order].G);
#if 0
MSG_DEBUG("steps: " << steps);
MSG_DEBUG("LV:" << std::endl << (*LV));
......
......@@ -392,6 +392,8 @@ void message_queue::run()
#define MSG_DEBUG_INDENT CREATE_MESSAGE(msg_channel::Log, "\x1")
#define MSG_DEBUG_DEDENT CREATE_MESSAGE(msg_channel::Log, "\x2")
struct scoped_indent { scoped_indent() { MSG_DEBUG_INDENT; } ~scoped_indent() { MSG_DEBUG_DEDENT; } };
inline void msg_handler_t::state_t::check(bool fatal)
{
msg_handler_t::scoped_lock_type _(msg_handler_t::mutex);
......
......@@ -167,6 +167,7 @@ struct fast_polynom {
fast_polynom& operator = (std::initializer_list<coef_t>);
fast_polynom& operator = (const std::vector<coef_t>&);
bool operator < (const fast_polynom& f) const;
bool operator == (const fast_polynom& f) const;
bool operator == (int f) const;
bool inf(const fast_polynom& p) const;
......
......@@ -19,6 +19,11 @@ inline bool fast_polynom::operator == (const fast_polynom& f) const
return value == f.value;
}
inline bool fast_polynom::operator < (const fast_polynom& f) const
{DEBUG_SCOPE;
return value < f.value;
}
inline bool fast_polynom::operator == (int f) const
{DEBUG_SCOPE;
return value == f;
......@@ -55,11 +60,13 @@ inline fast_polynom fast_polynom::operator + (const fast_polynom& q) const
/*impl::d_polynom R = P.add(Q);*/
/*std::cout << "P+Q=[" << result.value << "] " << R << std::endl;*/
result = impl::polynom_tables::reg(P.add(Q));
/*MSG_DEBUG("COMP; " << (*this) << " + " << q << " = " << result);*/
} else {
/*const impl::d_polynom& P = *this;*/
/*const impl::d_polynom& Q = q;*/
/*std::cout << "(precomputed) P[" << value << "]=" << P << " Q[" << q.value << "]=" << Q << std::endl;*/
/*std::cout << "(precomputed) P+Q=[" << result.value << "] " << result << std::endl;*/
/*MSG_DEBUG("PRECOMP; " << (*this) << " + " << q << " = " << result);*/
}
return result;
}
......@@ -71,6 +78,9 @@ inline fast_polynom fast_polynom::operator - (const fast_polynom& q) const
const impl::d_polynom& P = *this;
const impl::d_polynom& Q = q;
result = impl::polynom_tables::reg(P.sub(Q));
/*MSG_DEBUG("COMP; " << (*this) << " - " << q << " = " << result);*/
} else {
/*MSG_DEBUG("PRECOMP; " << (*this) << " - " << q << " = " << result);*/
}
return result;
}
......@@ -82,6 +92,9 @@ inline fast_polynom fast_polynom::operator * (const fast_polynom& q) const
const impl::polynom& P = *this;
const impl::polynom& Q = q;
result = impl::polynom_tables::reg(P.mul(Q));
/*MSG_DEBUG("COMP; " << (*this) << " * " << q << " = " << result);*/
} else {
/*MSG_DEBUG("PRECOMP; " << (*this) << " * " << q << " = " << result);*/
}
return result;
}
......
......@@ -73,7 +73,16 @@ namespace std {
{
static hash<int> hi;
static hash<std::vector<coef_t>> hv;
return impl::ROTATE<7>(impl::ROTATE<7>(hi(f.r_exp)) ^ hi(f.s_exp)) ^ hv(f.P);
size_t ret = impl::ROTATE<7>(impl::ROTATE<7>(hi(f.r_exp)) ^ hi(f.s_exp)) ^ hv(f.P);
#if 0
{
std::stringstream debug;
debug << "f_poly::hash(" << f.r_exp << ", " << f.s_exp << ',';
for (auto c: f.P) { debug << ' ' << c; }
MSG_DEBUG(debug.str() << ")=" << ret);
}
#endif
return ret;
}
};
......@@ -91,7 +100,10 @@ namespace std {
size_t operator () (const impl::operands_t& o) const
{
static hash<int> hi;
return impl::ROTATE<7>(hi(o.first)) ^ hi(o.second);
/*size_t ret = impl::ROTATE<7>(hi(o.first)) ^ hi(o.second);*/
size_t ret = (hi(o.first) * 0xdeadb3ef) ^ hi(o.second);
/*MSG_DEBUG("operands::hash(" << o.first << ',' << o.second << ")=" << ret);*/
return ret;
}
};
}
......
......@@ -100,7 +100,18 @@ inline std::ostream& operator << (std::ostream& os, const fast_polynom& f)
inline std::ostream& operator << (std::ostream& os, const impl::polynom& p)
{
return os << "{f=" << p.to_f() << " d=" << p.to_d() << '}';
/*return os << "{f=" << p.to_f() << " d=" << p.to_d() << '}';*/
os << '{';
if (p.f_init) {
os << "f=" << p.f;
if (p.d_init) {
os << ' ';
}
}
if (p.d_init) {
os << "d=" << p.d;
}
return os << '}';
}
inline std::ostream& operator << (std::ostream& os, const impl::operation_table_t& ot)
......
......@@ -40,7 +40,8 @@ namespace impl {
/*std::cout << " vec.size " << second.size() << std::endl;*/
auto inserted = first.insert(p2i_t::value_type(p, (fast_polynom) fast_polynom::not_initialized));
if (inserted.second) {DEBUG_SCOPE;
inserted.first->second = (fast_polynom) second.size();
/*inserted.first->second = (fast_polynom) second.size();*/
inserted.first->second.value = second.size();
second.push_back(&inserted.first->first);
/*std::cout << " -> regg'ed " << p << " / " << inserted.first->first*/
/*<< " as " << ((int)inserted.first->second)*/
......@@ -48,6 +49,7 @@ namespace impl {
/*<< " vec.size " << second.size()*/
/*<< std::endl;*/
}
/*MSG_DEBUG("REG " << p << " VALUE " << inserted.first->second.value);*/
return inserted.first->second;
}
};
......
......@@ -240,17 +240,30 @@ GenerationRS double_haploid(const GenerationRS& M)
return convert(kroneckerProduct(M.data, doubling_gametes()));
}
inline
std::map<char, VectorXd> marker_observation_spec::compile() const
{
std::map<char, VectorXd> ret;
for (const auto& score: scores) {
VectorXd obs = VectorXd::Zero(domain_size * domain_size);
for (const auto& s: score.second) {
obs(s.first + s.second * domain_size) = 1.;
}
ret[score.first] = obs;
}
return ret;
}
static inline std::ostream& operator << (std::ostream& os, const marker_observation_spec& mos)
{
os << "MARKER OBS[" << std::endl;
if (mos.front().first) {
for (auto& p: mos) {
if (mos.scores.front().first) {
for (auto& p: mos.scores) {
os << " * " << p.first << ": " << p.second << std::endl;
}
} else {
os << "RELATIVE unknown='" << mos.back().first << '\'';
os << "RELATIVE unknown='" << mos.scores.back().first << '\'';
}
return os << ']';
}
......@@ -339,7 +352,7 @@ generation_rs::generation_design_ancestor::generation_design_ancestor(const gene
: generation_design_base(my_gen, true)
{}
inline
VectorLC generation_rs::generation_design_ancestor::raw_lincomb() const { return make_parent_comb(g, 1); }
VectorLC generation_rs::generation_design_ancestor::raw_lincomb() const { return make_parent_comb((size_t) -1, 1); }
inline
generation_rs::generation_design_base* generation_rs::generation_design_ancestor::_clone(const generation_rs* my_gen) const { return new generation_design_ancestor(my_gen); }
inline
......@@ -356,7 +369,7 @@ generation_rs::generation_design_adhoc::generation_design_adhoc(const generation
: generation_design_base(my_gen, true)
{}
inline
VectorLC generation_rs::generation_design_adhoc::raw_lincomb() const { return make_parent_comb(g, g->main_process().innerSize()); }
VectorLC generation_rs::generation_design_adhoc::raw_lincomb() const { return make_parent_comb(/*FIXME*/(size_t)g, g->main_process().innerSize()); }
inline
generation_rs::generation_design_base* generation_rs::generation_design_adhoc::_clone(const generation_rs* my_gen) const { return new generation_design_adhoc(my_gen); }
inline
......@@ -375,7 +388,7 @@ generation_rs::generation_design_self::generation_design_self(const generation_r
inline
VectorLC generation_rs::generation_design_self::raw_lincomb() const
{
auto ct = make_parent_comb(parent, parent->main_process().data.cols());
auto ct = make_parent_comb(/*FIXME*/(size_t)parent, parent->main_process().data.cols());
/*return g->lumper.cast<gencomb_type>() * kroneckerProduct(ct, Matrix<gencomb_type, 4, 1>::Constant(.25));*/
return kroneckerProduct(ct, Matrix<gencomb_type, 4, 1>::Constant(.25));
}
......@@ -421,8 +434,8 @@ generation_rs::generation_design_cross::generation_design_cross(const generation
inline
VectorLC generation_rs::generation_design_cross::raw_lincomb() const
{
auto ct1 = make_parent_comb(parent1, parent1->main_process().data.cols());
auto ct2 = make_parent_comb(parent2, parent2->main_process().data.cols());
auto ct1 = make_parent_comb(/*FIXME*/(size_t)parent1, parent1->main_process().data.cols());
auto ct2 = make_parent_comb(/*FIXME*/(size_t)parent2, parent2->main_process().data.cols());
VectorLC p1, p2, ret;
p1 = kroneckerProduct(ct1, Matrix<gencomb_type, 2, 1>::Constant(.5));
p2 = kroneckerProduct(ct2, Matrix<gencomb_type, 2, 1>::Constant(.5));
......@@ -494,7 +507,7 @@ generation_rs::generation_design_double_haploid::generation_design_double_haploi
inline
VectorLC generation_rs::generation_design_double_haploid::raw_lincomb() const
{
auto ct = make_parent_comb(parent, parent->main_process().data.cols());
auto ct = make_parent_comb(/*FIXME*/(size_t)parent, parent->main_process().data.cols());
/*return g->lumper.cast<gencomb_type>() * kroneckerProduct(ct, Matrix<gencomb_type, 2, 1>::Constant(.5));*/
return kroneckerProduct(ct, Matrix<gencomb_type, 2, 1>::Constant(.5));
}
......@@ -594,14 +607,13 @@ void generation_rs::precompute_unique_labels()
const auto& breakmat = main_process();
size_t i = 0;
size_t j = 0;
unique_labels.clear();
while (i < (size_t)redux.innerSize()) {
while (j < breakmat.column_labels.size() && redux(i, j) != 1) ++j;
unique_labels.push_back(breakmat.column_labels[j]);
while (j < breakmat.column_labels.size() && redux(i, j) == 1) ++j;
++i;
}
/*MSG_DEBUG("unique_labels(" << name << ") = " << unique_labels);*/
/*MSG_DEBUG("labels(" << name << ") = " << main_process().column_labels);*/
}
inline
......@@ -887,17 +899,17 @@ generation_rs::segment_computer_t generation_rs::segment_computer(
inline std::ostream& operator << (std::ostream& os, const generation_rs& g)
{
os << "== Generation: " << g.name << " ======================" << std::endl;
os << "State combination: " << g.this_lincomb << std::endl;
/*os << "State combination: " << g.this_lincomb << std::endl;*/
os << g.genotype_probabilities() << std::endl;
os << "Locus Transition" << std::endl;
os << " Process size: " << g.main_process().innerSize() << std::endl;
#define GENERATION_PRINT_ALL_PROCESSES
/*#define GENERATION_PRINT_ALL_PROCESSES*/
#ifdef GENERATION_PRINT_ALL_PROCESSES
for (auto& p: g.P) {
os << "* weight " << p.weight << std::endl;
os << p.G << std::endl;
os << "colsums " << p.G.data.colwise().sum() << std::endl;
os << "rowsums " << p.G.data.rowwise().sum().transpose() << std::endl;
/*os << "colsums " << p.G.data.colwise().sum() << std::endl;*/
/*os << "rowsums " << p.G.data.rowwise().sum().transpose() << std::endl;*/
}
#else
os << g.main_process() << std::endl;
......@@ -905,14 +917,15 @@ inline std::ostream& operator << (std::ostream& os, const generation_rs& g)
os << "(" << (g.P.size() - 1) << " other parallel processes not shown)" << std::endl;
}
#endif
os << "Weights:";
double sum = 0;
for (auto& p: g.P) {
os << ' ' << p.weight;
sum += p.weight;
}
os << std::endl;
os << "Total: " << sum << " in " << g.P.size() << " processes" << std::endl;
/*os << "Weights:";*/
/*double sum = 0;*/
/*for (auto& p: g.P) {*/
/*os << ' ' << p.weight;*/
/*sum += p.weight;*/
/*}*/
/*os << std::endl;*/
/*os << "Total: " << sum << " in " << g.P.size() << " processes" << std::endl;*/