model.cc 45.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/* Spell-QTL  Software suite for the QTL analysis of modern datasets.
 * Copyright (C) 2016,2017  Damien Leroux <damien.leroux@inra.fr>, Sylvain Jasson <sylvain.jasson@inra.fr>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

Damien Leroux's avatar
Damien Leroux committed
18
19
#include "computations/model.h"

20
#if 0
Damien Leroux's avatar
Damien Leroux committed
21
22
23
MatrixXd
traits_to_matrix()
{
24
25
	std::function<std::vector<std::string>(const qtl_pop_type&)>
		get_labels = [](const qtl_pop_type& pop)
26
27
28
29
30
31
32
33
		{
			std::vector<std::string> labels;
			labels.reserve(pop.observed_traits.size());
			for (auto& ot: pop.observed_traits) {
				labels.push_back(ot.name);
			}
			return labels;
		};
34
35
	std::function<int(const qtl_pop_type&)>
		get_pop_size = [](const qtl_pop_type& pop)
36
37
38
		{
			return pop.observed_traits.front().values.size();
		};
39
	block_builder<std::string> bb;
Damien Leroux's avatar
Damien Leroux committed
40
	init_connected_block_builder(bb, get_pop_size, get_labels, all_populations());
41
42
43
44
45
46
	MatrixXd ret(bb.n_rows, bb.n_columns);
	auto builder = bb.begin(ret);
	for (auto& pop: active_settings->populations) {
		auto& observed_traits = pop.second.observed_traits;
		MatrixXd tmp(observed_traits.front().values.size(),
					 observed_traits.size());
Damien Leroux's avatar
Damien Leroux committed
47
48
        for (size_t i = 0; i < observed_traits.front().values.size(); ++i) {
            for (size_t j = 0; j < observed_traits.size(); ++j) {
49
                tmp(i, j) = observed_traits[j].values[i];
Damien Leroux's avatar
Damien Leroux committed
50
51
            }
        }
52
53
54
55
		builder += tmp;
	}
	/*MSG_DEBUG("traits" << std::endl << ret);*/
	return ret;
Damien Leroux's avatar
Damien Leroux committed
56
}
57
#endif
Damien Leroux's avatar
Damien Leroux committed
58
59


60

61
const MatrixXd&
62
63
trait_by_name(population_value pop, const std::string& name)
{
64
	static MatrixXd none;
65
66
67
68
69
    for (auto& ot: pop->observed_traits) {
        if (ot.name == name) {
            return ot.values;
        }
    }
70
    return none;
71
72
73
}


74
75
76
77
78
79
80
81
82
83
template <typename M>
void
random_permutation(M&& block, const MatrixXd& values) {
    PermutationMatrix<Dynamic,Dynamic> perm(block.rows());
    perm.setIdentity();
    std::random_shuffle(perm.indices().data(), perm.indices().data()+perm.indices().size());
    block = perm * values; // permute rows
}


84
value<MatrixXd>
85
trait_permutations(const std::string& trait_name, int n)
86
{
87
	size_t n_rows = 0;
88
	std::vector<MatrixXd> traits;
89
90
91
92
    const auto& populations = active_settings->linked_pops[trait_name];

    /* FIXME: are linked pops inter-permutable? */
    for (const auto& pop_vec: populations) {
93
//        traits.emplace_back();
94
95
        for (const auto& pop_ptr: pop_vec) {
            n_rows += pop_ptr->size();
96
97
98
//            const auto& trait = pop_ptr->observed_traits.front().values;
//            traits.back().insert(traits.back().end(), trait.begin(), trait.end());
            traits.emplace_back(pop_ptr->observed_traits.front().values);
99
100
101
        }
    }

102
    value<MatrixXd> ret = MatrixXd();
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
    int block_cols = traits[0].cols();
	ret->resize(n_rows, block_cols * n);
    for (int j = 0; j < ret->cols(); j += block_cols) {
        n_rows = 0;
    	for (size_t t = 0; t < traits.size(); ++t) {
            random_permutation(ret->block(n_rows, j, traits[t].rows(), block_cols), traits[t]);
            n_rows += traits[t].rows();
        }
//		size_t ts = traits[t].size();
//		auto slice = ret->middleRows(n_rows, ts);
//		for (int j = 0; j < n; ++j) {
//			auto col = slice.col(j);
//			for (size_t i = 0; i < traits[t].size(); ++i) {
//				col(i) = traits[t][i];
//			}
//			std::random_shuffle(col.data(), col.data() + ts);
//		}
//		n_rows += ts;
121
	}
122
123
124
125
    return ret;
}


126
127
128
129
value<MatrixXd>
trait_matrix(const std::string& trait_name, const collection<population_value>& all_pops)
{
	size_t n_rows = 0;
130
	std::vector<const MatrixXd*> traits;
131
132
133
	traits.reserve(all_pops.size());
	for (auto& pop: all_pops) {
		traits.push_back(&trait_by_name(*pop, trait_name));
134
135
        MSG_DEBUG("have single_trait from sub-pop with " << traits.back()->rows() << " values");
		n_rows += traits.back()->rows();
136
	}
137
    MSG_DEBUG("have a total of " << n_rows << " rows");
138
    value<MatrixXd> ret = MatrixXd();
139
	ret->resize(n_rows, traits[0]->cols());
140
141
142
	n_rows = 0;
	int t = 0;
	for (auto& pop: all_pops) {
143
144
		size_t ts = traits[t]->rows();
		ret->middleRows(n_rows, ts) = *traits[t];
145
        n_rows += ts;
146
        ++t;
147
148
149
150
151
152
        (void)pop;
	}
    return ret;
}


153
#if 0
154
155
value<MatrixXd>
residuals_permutations(const model& M, int n)
Damien Leroux's avatar
Damien Leroux committed
156
{
157
158
159
160
161
162
163
164
165
166
	size_t n_rows = 0;
	std::vector<size_t> sizes;
	sizes.reserve(active_settings->populations.size());
	for (auto& pop: active_settings->populations) {
		sizes.push_back(pop.second.observed_traits.front().values.size());
		n_rows += sizes.back();
	}
    value<MatrixXd> ret = MatrixXd();
	ret->resize(n_rows, n);
	n_rows = 0;
167
	/*int t = 0;*/
168
    /* FIXME: only works for ONE single_trait at a time */
169
170
171
    /*std::cerr << "ret$dim(" << ret->innerSize() << ',' << ret->outerSize() << ')' << std::endl;*/
    /*std::cerr << "residuals$dim(" << M.residuals().innerSize() << ',' << M.residuals().outerSize() << ')' << std::endl;*/
    ret->colwise() = M.residuals().col(0);
172
173
	/*for (auto& pop: active_settings->populations) {*/
	for (size_t t = 0; t < active_settings->populations.size(); ++t) {
174
175
176
177
178
179
180
		size_t ts = sizes[t];
		auto slice = ret->middleRows(n_rows, ts);
		for (int j = 0; j < n; ++j) {
			auto col = slice.col(j);
			std::random_shuffle(col.data(), col.data() + ts);
		}
		n_rows += ts;
181
		/*++t;*/
182
183
	}
    return ret;
Damien Leroux's avatar
Damien Leroux committed
184
}
185
#endif
Damien Leroux's avatar
Damien Leroux committed
186

187

188
model_block_type
189
190
/*cross_indicator(const collection<population_value>& pops)*/
cross_indicator(const std::string& trait)
Damien Leroux's avatar
Damien Leroux committed
191
{
192
    const auto& pops = active_settings->linked_pops[trait];
193
    int n_cols = pops.size();
Damien Leroux's avatar
Damien Leroux committed
194
    int n_rows = 0;
195
196
197
198
    for (const auto& pop_vec: pops) {
        for (const auto& pop: pop_vec) {
            n_rows += pop->size();
        }
Damien Leroux's avatar
Damien Leroux committed
199
    }
200
	model_block_type val;
201
    /*int p = 0;*/
Damien Leroux's avatar
Damien Leroux committed
202
203
204
205
206
    std::vector<int> columns;
    const auto& lp = active_settings->linked_pops;
    std::vector<int> pop_col(lp.size(), -1);
    /*for (int i = 0; i < pop_col.size(); ++i) { pop_col[i] = i; }*/
    std::vector<int> block_size_by_col(lp.size(), 0);
207
    size_t col = 0, row = 0;
208
    val.data = MatrixXd::Zero(n_rows, n_cols);
209
210
211
212
213
214
    for (const auto& pop_vec: pops) {
        for (const auto& pop: pop_vec) {
            val.data.block(row, col, pop->size(), 1) = VectorXd::Ones(pop->size());
            row += pop->size();
        }
        ++col;
Damien Leroux's avatar
Damien Leroux committed
215
    }
216

Damien Leroux's avatar
Damien Leroux committed
217
    val.column_labels.clear();
218
219
220
221
222
223
224
    val.column_labels.reserve(val.data.outerSize());
    for (const auto& pop_vec: pops) {
        for (const auto& pop: pop_vec) {
            val.column_labels.emplace_back(pop->qtl_generation_name.begin(), pop->qtl_generation_name.end());
        }
    }

Damien Leroux's avatar
Damien Leroux committed
225
226
227
228
229
230
231
232
    /*for (const auto& pop: pops) {*/
        /*int sz = (*pop)->observed_traits.front().values.size();*/
        /*val.data.block(n_rows, p, sz, 1) = VectorXd::Ones(sz);*/
        /*val.column_labels.emplace_back();*/
        /*val.column_labels.back().push_back(' ');*/
        /*n_rows += sz;*/
        /*++p;*/
    /*}*/
233
	/*MSG_DEBUG("cross indicator " << MATRIX_SIZE(val) << std::endl << val);*/
Damien Leroux's avatar
Damien Leroux committed
234
    /*MSG_QUEUE_FLUSH();*/
Damien Leroux's avatar
Damien Leroux committed
235
236
237
238
    return val;
}


239
240
241
242
243
244
value<model_block_type>
make_covariables_block(const std::string& trait, const std::vector<std::string>& covar_names) {
    if (covar_names.size() == 0) {
        return model_block_type{};
    }
    collection<model_block_type> blocks;
245
//    MSG_DEBUG("make_covariables_block");
246
247
248
    for (const auto& vqp: active_settings->linked_pops[trait]) {
        for (const auto& qp: vqp) {
            blocks.emplace_back(model_block_type());
249
250
//            MSG_DEBUG("COVAR QP " << qp);
//            MSG_QUEUE_FLUSH();
251
252
253
254
255
            blocks.back()->data = qp->covariables.values;
            for (const auto& l: qp->covariables.dim_names) {
                blocks.back()->column_labels.emplace_back(l.begin(), l.end());
            }
            blocks.back()->row_labels.assign(qp->covariables.good_indices.begin(), qp->covariables.good_indices.end());
256
//            MSG_DEBUG("COVAR BLOCK " << MATRIX_SIZE(blocks.back()->data) << " CLAB " << blocks.back()->column_labels << " RLAB " << blocks.back()->row_labels);
257
258
        }
    }
259
//    MSG_DEBUG("make_covariables_block assemble");
260
261
262
263
    return assemble_block_multipop(blocks);
}


264
#if 0
Damien Leroux's avatar
Damien Leroux committed
265
266
267
model
basic_model()
{
268
	value<MatrixXd> traits = make_value<Sync|Mem>(traits_to_matrix);
Damien Leroux's avatar
Damien Leroux committed
269
    model M(*traits, all_populations());
270
    M.add_block({}, cross_indicator(all_populations()));
Damien Leroux's avatar
Damien Leroux committed
271
272
273
274
275
276
277
278
    const char* decomp_method = getenv("DECOMP_METHOD");
    std::string s(decomp_method ? decomp_method : "");
    /*std::cout << "Xt" << std::endl << M.X().transpose() << std::endl;*/
    /*std::cout << "Yt" << std::endl << M.Y().transpose() << std::endl;*/
    /*std::cout << "XtX^-1" << std::endl << M.XtX_pseudo_inverse() << std::endl;*/
    /*std::cout << "coefs" << std::endl << M.coefficients() << std::endl;*/
    /*auto r = M.residuals().array();*/
    /*std::cout << "residuals" << std::endl << (r * r).colwise().sum() << std::endl;*/
279
	M.compute();
Damien Leroux's avatar
Damien Leroux committed
280
281
    return M;
}
282
#endif
Damien Leroux's avatar
Damien Leroux committed
283

284
285
286
287
288
289
290
291
#if 0
VectorXd
f_tester(const model& M, const model& M0, const parental_origin_per_locus_type& popl)
{
    model M1 = extend(M0, popl);
	M1.compute();
    return f_test(*const_cast<model*>(&M), *const_cast<model*>(&M1));
}
Damien Leroux's avatar
Damien Leroux committed
292
293
294
295
296

MatrixXd
f_test_along_chromosome(const value<model>& Mcurrent, const value<model>& M0,
                                 const collection<parental_origin_per_locus_type>& popl)
{
297
                DUMP_FILE_LINE();
Damien Leroux's avatar
Damien Leroux committed
298
299
300
    collection<VectorXd> aft
        = make_collection(f_tester, Mcurrent, M0, popl);

301
                DUMP_FILE_LINE();
Damien Leroux's avatar
Damien Leroux committed
302
303
304
    MatrixXd ret(aft[0]->innerSize(), aft.size());

    for (size_t i = 0; i < aft.size(); ++i) {
305
                DUMP_FILE_LINE();
Damien Leroux's avatar
Damien Leroux committed
306
307
308
309
        ret.col(i) = *aft[i];
    }
    return ret;
}
310
311
#endif

312
model extend_model(const model& m0, const model_block_key& k, const parental_origin_per_locus_type& block)
313
{
314
    model ret = m0.extend(k, block);
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
    ret.compute();
    return ret;
}

#if 0
VectorXd compute_ftest(const model& m0, const model& m1)
{
    return f_test(m0, m1);
}

VectorXd compute_r2(const model& m0, const model& m1) { return r2(m0, m1); }
VectorXd compute_chi2(const model& m0, const model& m1)
{
    VectorXd ret(1);
    ret(0) = chi2(m0, m1);
    return ret;
}

void
consolidate_computation_result(MatrixXd& ret,
                               const MatrixXd& (model::* func) () const,
                               const collection<model>& models)
{
    ret.resize(((*models[0]).*func)().innerSize(),
               ((*models[0]).*func)().outerSize() * models.size());
    size_t bl_size = ((*models[0]).*func)().outerSize();
    if (bl_size > 1) {
        size_t ofs = 0;
        for (size_t i = 0; i < models.size(); ++i) {
            ret.middleCols(ofs, bl_size) = ((*models[i]).*func)();
            ofs += bl_size;
        }
    } else {
        for (size_t i = 0; i < models.size(); ++i) {
            ret.col(i) = ((*models[i]).*func)().col(0);
        }
    }
}

void
consolidate_computation_result(MatrixXd& ret,
                               const VectorXd& (model::* func) () const,
                               const collection<model>& models)
{
    ret.resize(((*models[0]).*func)().innerSize(), models.size());
    for (size_t i = 0; i < models.size(); ++i) {
        ret.col(i) = ((*models[i]).*func)();
    }
}

void
consolidate_computation_result(VectorXd& ret,
                               int (model::* func) () const,
                               const collection<model>& models)
{
    ret.resize(models.size());
    for (size_t i = 0; i < models.size(); ++i) {
        ret(i) = ((*models[i]).*func)();
    }
}
#endif

void
one_computation_result(int i, MatrixXd& ret, const MatrixXd& model_member)
{
    size_t bl_size = model_member.outerSize();
    if (bl_size > 1) {
        ret.middleCols(bl_size * i, bl_size) = model_member;
    } else {
        ret.col(i) = model_member.col(0);
    }
}

void
one_computation_result(int i, MatrixXd& ret, const VectorXd& model_member)
{
    ret.col(i) = model_member.col(0);
}

void
one_computation_result(int i, VectorXd& ret, int model_member)
{
    ret(i) = model_member;
}


Damien Leroux's avatar
Damien Leroux committed
401
402
403
404
405
406
407
408
409
410
411
412
413
struct dump_model {
    std::string name;
    const model& m;

    friend
        std::ostream&
        operator << (std::ostream& os, const dump_model& dm)
        {
            return os << '[' << dm.name << " {" << dm.m.keys() << "} rank=" << dm.m.rank() << " rss=" << dm.m.rss().transpose() << ']';
        }
};


414
415
int compute_one(int i, computation_along_chromosome* ret,
                ComputationType ct, ComputationResults cr,
416
                size_t y_block_size,
417
                const model& Mcurrent, const model& M0,
418
                const model_block_key& k,
419
420
                const parental_origin_per_locus_type& popl)
{
421
    /*DUMP_FILE_LINE();*/
422
//    MSG_DEBUG("one popl" << std::endl << popl.transpose());
423
    model Mext = M0.extend(k, popl);
424
    Mext.compute();
Damien Leroux's avatar
Damien Leroux committed
425
426
//    MSG_DEBUG("Rank(Mext " << Mext.keys() << ") = " << Mext.rank());
//    MSG_DEBUG('[' << Mext.keys() << " rank=" << Mext.rank() << " rss=" << Mext.rss().transpose() << ']');
427

Damien Leroux's avatar
Damien Leroux committed
428
//    DUMP_FILE_LINE();
429
430
431
    if (cr & RSS) {
        one_computation_result(i, ret->rss, Mext.rss());
    }
Damien Leroux's avatar
Damien Leroux committed
432
//    DUMP_FILE_LINE();
433
434
435
    if (cr & Residuals) {
        one_computation_result(i, ret->residuals, Mext.residuals());
    }
Damien Leroux's avatar
Damien Leroux committed
436
//    DUMP_FILE_LINE();
437
438
439
440
    if (cr & Coefficients) {
        one_computation_result(i, ret->coefficients, Mext.coefficients());
    }

Damien Leroux's avatar
Damien Leroux committed
441
//    DUMP_FILE_LINE();
442
443
444
    if (cr & Rank) {
        one_computation_result(i, ret->rank, Mext.rank());
    }
Damien Leroux's avatar
Damien Leroux committed
445
//    DUMP_FILE_LINE();
446
447
448
449
    if (ct & (FTest|FTestLOD)) {
        f_test(Mcurrent, Mext, i,
               ct & FTest ? &ret->ftest_pvalue : NULL,
               ct & FTestLOD ? &ret->ftest_lod : NULL);
Damien Leroux's avatar
Damien Leroux committed
450

451
#if 0
452
        if ((ct & FTest) && Mcurrent.Y().outerSize() <= 5) {
453
            dump_model m1{"Mcurrent", Mcurrent}, m2{"Mext", Mext};
Damien Leroux's avatar
Damien Leroux committed
454
455
            MSG_DEBUG(m1 << " / " << m2 << " ftest=" << ret->ftest_pvalue.col(i).transpose());
        }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
456
#endif
457
    }
Damien Leroux's avatar
Damien Leroux committed
458
//    DUMP_FILE_LINE();
459
460
461
    if (ct & Mahalanobis) {
        mahalanobis(Mcurrent, Mext, i, &ret->mahalanobis);
    }
Damien Leroux's avatar
Damien Leroux committed
462
//    DUMP_FILE_LINE();
463
464
465
    if (ct & R2) {
        r2(Mcurrent, Mext, i, &ret->r2);
    }
Damien Leroux's avatar
Damien Leroux committed
466
//    DUMP_FILE_LINE();
467
    if (ct & Chi2) {
468
469
470
471
472
//        auto tmp = chi2(Mcurrent, Mext, y_block_size);
//        MSG_DEBUG("HAVE CHI2 COLUMN " << tmp.rows() << 'x' << tmp.cols());
//        MSG_DEBUG("HAVE recipient   " << ret->chi2.rows() << 'x' << ret->chi2.cols());
        ret->chi2.col(i) = chi2(Mcurrent, Mext, y_block_size);
//        ret->chi2.col(i) = chi2(Mcurrent, Mext, y_block_size).transpose();
473
    }
Damien Leroux's avatar
Damien Leroux committed
474
475
476
477
478
479
480
481
482
    if (ct & Chi2LOD) {
//        auto tmp = chi2(Mcurrent, Mext, y_block_size);
//        MSG_DEBUG("HAVE CHI2 COLUMN " << tmp.rows() << 'x' << tmp.cols());
//        MSG_DEBUG("HAVE recipient   " << ret->chi2.rows() << 'x' << ret->chi2.cols());
//        MSG_DEBUG("Chi2LOD");
        ret->chi2_lod.col(i) = chi2(Mcurrent, Mext, y_block_size, true);
//        ret->chi2.col(i) = chi2(Mcurrent, Mext, y_block_size).transpose();
    }
//    DUMP_FILE_LINE();
483
484
485
486

    return 0;
}

Damien Leroux's avatar
Damien Leroux committed
487
488

model
489
init_model(const value<MatrixXd>& Y, const value<model_block_type>& indicator, const value<model_block_type>& covars, double threshold)
Damien Leroux's avatar
Damien Leroux committed
490
{
491
    model M(Y, threshold, all_populations());
492
493
//    const char* decomp_method = getenv("DECOMP_METHOD");  /* FIXME get rid of this sh\bt. */
//    std::string s(decomp_method ? decomp_method : "");
494
    M.m_max_order = 1;
495
    M.add_block(model_block_key_struc::cross_indicator(), indicator);
Damien Leroux's avatar
Damien Leroux committed
496
//    MSG_DEBUG("covariables " << active_settings->with_covariables << " covars.cols=" << covars->data.cols());
497
    if (covars && covars->data.cols() > 0) {
498
499
        M.add_block(model_block_key_struc::covariables(), covars);
    }
500
	/*std::cout << "indicator" << std::endl << indicator->transpose() << std::endl;*/
501
502
503
504
505
506
507
    /*std::cout << "Xt" << std::endl << M.X().transpose() << std::endl;*/
    /*std::cout << "Yt" << std::endl << M.Y().transpose() << std::endl;*/
    /*std::cout << "XtX^-1" << std::endl << M.XtX_pseudo_inverse() << std::endl;*/
    /*std::cout << "coefs" << std::endl << M.coefficients() << std::endl;*/
    /*auto r = M.residuals().array();*/
    /*std::cout << "residuals" << std::endl << (r * r).colwise().sum() << std::endl;*/
	M.compute();
Damien Leroux's avatar
Damien Leroux committed
508
509
510
    return M;
}

Damien Leroux's avatar
Damien Leroux committed
511

512
513
514
515
516
517
VectorXd
row_max(const MatrixXd& mat)
{
    return mat.rowwise().maxCoeff();
}

Damien Leroux's avatar
Damien Leroux committed
518

519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
std::vector<double>
get_quantiles(const VectorXd& values, const std::vector<double>& quantiles)
{
    std::function<bool(double, double)> comp = [] (double a, double b) { return a > b; };
    /* quantiles should be in range [0;1) */
    size_t size = values.innerSize();
    size_t max_interesting_index
        = (size_t) (size * *std::max_element(quantiles.begin(),
                                             quantiles.end()));
    std::vector<double> sorted_values(values.data(),
                                      values.data() + size);
    std::partial_sort(sorted_values.begin(),
                      sorted_values.begin() + max_interesting_index + 1,
                      sorted_values.end(),
                      comp);
    std::vector<double> qvec;
    for (double q: quantiles) {
        qvec.push_back(sorted_values[(size_t) (q * size)]);
    }
538
    MSG_DEBUG("quantiles: " << qvec);
539
540
    return qvec;
}
Damien Leroux's avatar
Damien Leroux committed
541

542
/*
543
collection<block_builder<std::vector<char>>>
544
545
546
547
compute_block_builders(const collection<population_value>& all_pop,
					   const value<qtl_chromosome_value>& qtl_chr,
					   const value<MatrixXd>& LD_matrix)
{
548
	collection<block_builder<std::vector<char>>> ret;
549
550
551
	return ret;
}
*/
552

Damien Leroux's avatar
Damien Leroux committed
553

554
collection<model_block_type>
Damien Leroux's avatar
Damien Leroux committed
555
556
557
assemble_parental_origins_multipop(
        const value<chromosome_value>& chr,
        const value<locus_key>& lk,
Damien Leroux's avatar
Damien Leroux committed
558
        const std::vector<collection<parental_origin_per_locus_type>>& all_popl,
559
560
        const collection<population_value>& all_pop,
        bool block_is_POP)
561
{
Damien Leroux's avatar
Damien Leroux committed
562
    /* FIXME take into account the ACTUAL list of populations used! */
Damien Leroux's avatar
Damien Leroux committed
563
564
    /*MSG_DEBUG("ALL_POPL " << (*chr)->name << " " << lk);*/
    /*for (const auto& popl: all_popl) {*/
Damien Leroux's avatar
Damien Leroux committed
565
566
    /*MSG_DEBUG("======================================");*/
    /*MSG_DEBUG("" << popl);*/
Damien Leroux's avatar
Damien Leroux committed
567
    /*}*/
Damien Leroux's avatar
Damien Leroux committed
568
    /*MSG_QUEUE_FLUSH();*/
Damien Leroux's avatar
Damien Leroux committed
569
570
571
572
573
574
    static
    std::function<int(const qtl_pop_type&)>
            get_pop_size = [](const qtl_pop_type& pop)
    {
        return pop.size();
    };
575

576
    std::function<std::vector<std::vector<char>>(const qtl_pop_type&)>
Damien Leroux's avatar
Damien Leroux committed
577
578
            get_labels;

579
580
    if (block_is_POP) {
        get_labels
Damien Leroux's avatar
Damien Leroux committed
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
                = [=](const qtl_pop_type& pop)
        {
            context_key ck(new context_key_struc(&pop, *chr, std::vector<double>()));
            /*MSG_DEBUG(ck);*/
            return make_value<Sync|Mem>(compute_state_to_parental_origin,
                                        value<context_key>{ck},
                                        lk)->row_labels;
        };
    } else {
        get_labels
                = [=](const qtl_pop_type& pop)
        {
            context_key ck(new context_key_struc(&pop, *chr, std::vector<double>()));
            /*MSG_DEBUG(ck);*/
            return make_value<Sync|Mem>(compute_state_to_dominance,
                                        value<context_key>{ck},
                                        lk)->row_labels;
        };
    }

    /*DUMP_FILE_LINE();*/
    block_builder<std::vector<char>> bb;
    if (active_settings->connected) {
        /*MSG_DEBUG("CONNECTED!");*/
        init_connected_block_builder(bb, get_pop_size, get_labels, all_pop);
    } else {
        /*MSG_DEBUG("DISCONNECTED!");*/
        init_disconnected_block_builder(bb, get_pop_size, get_labels, all_pop);
    }
    /*DUMP_FILE_LINE();*/
    collection<model_block_type> all_blocks;
    all_blocks.reserve(all_popl.front().size());
    size_t n = all_popl.front().size();
    for (size_t i = 0; i < n; ++i) {
        /*DUMP_FILE_LINE();*/
        all_blocks.emplace_back(new immediate_value<model_block_type>(model_block_type()));
        /*DUMP_FILE_LINE();*/
        all_blocks.back()->column_labels = bb.labels;
        auto builder = bb.begin(all_blocks.back()->data);
        /*DUMP_FILE_LINE();*/
        /*MSG_DEBUG("building at #" << i);*/
        for (auto& popl: all_popl) {
            /*MSG_DEBUG("popl" << std::endl << popl);*/
            /*DUMP_FILE_LINE();*/
            /*MSG_DEBUG("  block " << std::endl << popl[i]->data);*/
            /*MSG_QUEUE_FLUSH();*/
            builder += popl[i]->data;
        }
        /*MSG_DEBUG('#' << i << std::endl << all_blocks.back());*/
    }

    /*MSG_DEBUG("ALL BLOCKS" << std::endl << all_blocks);*/

    return all_blocks;
}


collection<model_block_type>
assemble_parental_origins_multipop(
        const value<chromosome_value>& chr,
        const std::vector<locus_key>& lk,
        const std::vector<collection<parental_origin_per_locus_type>>& all_popl,
        const collection<population_value>& all_pop,
        bool block_is_POP)
{
    /* FIXME take into account the ACTUAL list of populations used! */
    /*MSG_DEBUG("ALL_POPL " << (*chr)->name << " " << lk);*/
    /*for (const auto& popl: all_popl) {*/
    /*MSG_DEBUG("======================================");*/
    /*MSG_DEBUG("" << popl);*/
    /*}*/
    /*MSG_QUEUE_FLUSH();*/
    static
    std::function<int(const qtl_pop_type&)>
            get_pop_size = [](const qtl_pop_type& pop)
    {
        return pop.size();
    };

    std::function<std::vector<std::vector<char>>(const qtl_pop_type&)>
            get_labels;

    /*DUMP_FILE_LINE();*/
    collection<model_block_type> all_blocks;
    all_blocks.reserve(all_popl.front().size());
    size_t n = all_popl.front().size();
    for (size_t i = 0; i < n; ++i) {
        if (block_is_POP) {
            get_labels
                    = [&](const qtl_pop_type& pop)
671
672
673
674
675
            {
                context_key ck(new context_key_struc(&pop, *chr, std::vector<double>()));
                /*MSG_DEBUG(ck);*/
                return make_value<Sync|Mem>(compute_state_to_parental_origin,
                                            value<context_key>{ck},
Damien Leroux's avatar
Damien Leroux committed
676
                                            as_value(lk[i]))->row_labels;
677
            };
Damien Leroux's avatar
Damien Leroux committed
678
679
680
        } else {
            get_labels
                    = [=](const qtl_pop_type& pop)
681
682
683
684
685
            {
                context_key ck(new context_key_struc(&pop, *chr, std::vector<double>()));
                /*MSG_DEBUG(ck);*/
                return make_value<Sync|Mem>(compute_state_to_dominance,
                                            value<context_key>{ck},
Damien Leroux's avatar
Damien Leroux committed
686
                                            as_value(lk[i]))->row_labels;
687
            };
Damien Leroux's avatar
Damien Leroux committed
688
        }
689

Damien Leroux's avatar
Damien Leroux committed
690
        /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
691
692
693
694
695
696
697
698
699
700
        block_builder<std::vector<char>> bb;
        if (active_settings->connected) {
            /*MSG_DEBUG("CONNECTED!");*/
            init_connected_block_builder(bb, get_pop_size, get_labels, all_pop);
        } else {
            /*MSG_DEBUG("DISCONNECTED!");*/
            init_disconnected_block_builder(bb, get_pop_size, get_labels, all_pop);
        }
        /*DUMP_FILE_LINE();*/
        all_blocks.emplace_back(new immediate_value<model_block_type>(model_block_type()));
Damien Leroux's avatar
Damien Leroux committed
701
        /*DUMP_FILE_LINE();*/
702
        all_blocks.back()->column_labels = bb.labels;
Damien Leroux's avatar
Damien Leroux committed
703
        auto builder = bb.begin(all_blocks.back()->data);
Damien Leroux's avatar
Damien Leroux committed
704
        /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
705
706
        /*MSG_DEBUG("building at #" << i);*/
        for (auto& popl: all_popl) {
707
            /*MSG_DEBUG("popl" << std::endl << popl);*/
708
            /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
709
            /*MSG_DEBUG("  block " << std::endl << popl[i]->data);*/
710
            /*MSG_QUEUE_FLUSH();*/
Damien Leroux's avatar
Damien Leroux committed
711
712
713
714
            builder += popl[i]->data;
        }
        /*MSG_DEBUG('#' << i << std::endl << all_blocks.back());*/
    }
715

Damien Leroux's avatar
Damien Leroux committed
716
    /*MSG_DEBUG("ALL BLOCKS" << std::endl << all_blocks);*/
717

Damien Leroux's avatar
Damien Leroux committed
718
    return all_blocks;
719
720
}

721

Damien Leroux's avatar
Damien Leroux committed
722
723
724
725
726

MatrixXd
contrast_groups(const collection<population_value>& all_pops, const locus_key& lk)
{
	static
727
728
        std::function<int(const qtl_pop_type&)>
		get_pop_size = [](const qtl_pop_type&)
Damien Leroux's avatar
Damien Leroux committed
729
730
731
732
		{
			return 2;
		};

733
734
    std::function<std::vector<std::vector<char>>(const qtl_pop_type&)>
		get_labels = [=](const qtl_pop_type& pop)
Damien Leroux's avatar
Damien Leroux committed
735
736
737
738
739
740
		{
            static chromosome dummy;
            context_key ck(new context_key_struc(&pop, NULL, std::vector<double>()));
            /*MSG_DEBUG(ck);*/
            return make_value<Mem>(compute_state_to_parental_origin,
                                   value<context_key>{ck},
741
                                   as_value(lk->parent))->row_labels;
Damien Leroux's avatar
Damien Leroux committed
742
743
		};

744
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
745
746
747
748
749
750
751
752
753
754
755
756
	block_builder<std::vector<char>> bb;
	if (active_settings->connected) {
		/*MSG_DEBUG("CONNECTED!");*/
		init_connected_block_builder(bb, get_pop_size, get_labels, all_pops);
	} else {
		/*MSG_DEBUG("DISCONNECTED!");*/
		init_disconnected_block_builder(bb, get_pop_size, get_labels, all_pops);
	}

    MatrixXb big;
    auto builder = bb.begin(big);
    for (const auto& p: all_pops) {
Damien Leroux's avatar
Damien Leroux committed
757
        builder += get_contrast_groups((*p)->gen.get(), lk);
Damien Leroux's avatar
Damien Leroux committed
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
    }
    /*MSG_DEBUG("big" << std::endl << big);*/

    MatrixXb tmp = big.transpose() * big;
    MatrixXb test = MatrixXb::Zero(tmp.innerSize(), tmp.outerSize());
    while (tmp != test) {
        test = tmp;
        tmp = test * test;
    }
    /*MSG_DEBUG("test" << std::endl << test);*/

    bool_row_set uniq;
    for (int i = 0; i < test.innerSize(); ++i) {
        uniq.insert(test.row(i));
    }

    test.resize(uniq.size(), test.outerSize());
    int r = 0;
    for (const auto& row: uniq) {
        test.row(r++) = row;
    }

    return test.cast<double>();
}



785
786
787
788
collection<parental_origin_per_locus_type>
disassemble_parental_origins_multipop(
        const chromosome_value& chr,
        const locus_key& lk,
Damien Leroux's avatar
Damien Leroux committed
789
790
        model_block_type& block,
        const collection<population_value>& all_pop)
791
792
{
	static
793
794
        std::function<int(const qtl_pop_type&)>
		get_pop_size = [](const qtl_pop_type& pop)
795
796
797
798
		{
			return pop.observed_traits.front().values.size();
		};

799
800
    std::function<std::vector<std::vector<char>>(const qtl_pop_type&)>
		get_labels = [=](const qtl_pop_type& pop)
801
802
803
804
805
		{
            context_key ck(new context_key_struc(&pop, chr, std::vector<double>()));
            /*MSG_DEBUG(ck);*/
            return make_value<Mem>(compute_state_to_parental_origin,
                                   value<context_key>{ck},
806
                                   as_value(lk))->row_labels;
807
808
		};

809
    /*DUMP_FILE_LINE();*/
810
811
812
	block_builder<std::vector<char>> bb;
	if (active_settings->connected) {
		/*MSG_DEBUG("CONNECTED!");*/
Damien Leroux's avatar
Damien Leroux committed
813
		init_connected_block_builder(bb, get_pop_size, get_labels, all_pop);
814
815
	} else {
		/*MSG_DEBUG("DISCONNECTED!");*/
Damien Leroux's avatar
Damien Leroux committed
816
		init_disconnected_block_builder(bb, get_pop_size, get_labels, all_pop);
817
818
819
	}
    collection<parental_origin_per_locus_type> ret;
    ret.reserve(bb.elements.size());
Damien Leroux's avatar
Damien Leroux committed
820
    auto builder = bb.disasm_begin(block.data);
821
822
823
824
825
826
827
828
    for (size_t i = 0; i < bb.elements.size(); ++i) {
        ret.emplace_back(new immediate_value<model_block_type>(parental_origin_per_locus_type()));
        builder -= ret.back()->data;
    }
    return ret;
}


Damien Leroux's avatar
Damien Leroux committed
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
std::vector<double>
compute_effective_loci(const locus_key& lk, const std::vector<double>& loci)
{
    if (lk && lk->locus != locus_key_struc::no_locus) {
        std::vector<double> test_loci;
        test_loci.reserve(loci.size());
        for (double k: loci) {
            if (!lk->has(k)) { test_loci.push_back(k); }
        }
        return test_loci;
    } else {
        return loci;
    }
}


845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
template <typename Ret, typename... Args>
collection<model_block_type>
_compute_parental_origins_multipop(Ret (&Joint_POP_at_locus)(Args...),
                                   const collection<population_value>& all_pop,
                                   const value<chromosome_value>& chr,
                                   const value<locus_key>& lk,
                                   const value<std::vector<double>>& loci,
                                   const std::vector<double>& test_loci)
{
    /*DUMP_FILE_LINE();*/
	std::vector<collection<parental_origin_per_locus_type>> all_popl;
    /*DUMP_FILE_LINE();*/
	all_popl.reserve(all_pop.size());
    /*DUMP_FILE_LINE();*/

    /*MSG_DEBUG("loci are " << loci);*/
    /*MSG_DEBUG("given " << lk << ", test loci are " << test_loci);*/

    collection<locus_key> lkeys;
    lkeys.reserve(test_loci.size());
    for (double l: test_loci) {
        lkeys.emplace_back(*lk + l);
    }

	for (auto& p: all_pop) {
        /*DUMP_FILE_LINE();*/
        context_key ck(new context_key_struc(*p, *chr, *loci));
        all_popl.emplace_back(make_collection<Disk>(Joint_POP_at_locus,
                                                    value<context_key>{ck}, lkeys));
        /*std::cout << all_popl << std::endl;*/
		/*all_popl.emplace_back(compute_parental_origins(p, (*qtl_chr)->chr, loci));*/

	}
    for (auto& pl: all_popl) {
        for (auto& l: pl) {
            (void) *l;
        }
    }
    return assemble_parental_origins_multipop(chr, lk, all_popl, all_pop, &Joint_POP_at_locus == &joint_parental_origin_at_locus);
}



collection<model_block_type>
compute_dominance_multipop(const collection<population_value>& all_pop,
Damien Leroux's avatar
Damien Leroux committed
890
891
892
893
894
                           const value<chromosome_value>& chr,
                           const value<locus_key>& lk,
        /*const value<qtl_chromosome_value>& qtl_chr,*/
                           const value<std::vector<double>>& loci,
                           const std::vector<double>& test_loci)
895
896
{
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
897
    std::vector<collection<parental_origin_per_locus_type>> all_popl;
898
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
899
    all_popl.reserve(all_pop.size());
900
901
902
903
904
905
906
907
908
909
910
    /*DUMP_FILE_LINE();*/

    /*MSG_DEBUG("loci are " << loci);*/
    /*MSG_DEBUG("given " << lk << ", test loci are " << test_loci);*/

    collection<locus_key> lkeys;
    lkeys.reserve(test_loci.size());
    for (double l: test_loci) {
        lkeys.emplace_back(*lk + l);
    }

Damien Leroux's avatar
Damien Leroux committed
911
    for (auto& p: all_pop) {
912
913
914
915
916
        /*DUMP_FILE_LINE();*/
        context_key ck(new context_key_struc(*p, *chr, *loci));
        all_popl.emplace_back(make_collection<Disk>(joint_dominance_at_locus,
                                                    value<context_key>{ck}, lkeys, as_value(false)));
        /*std::cout << all_popl << std::endl;*/
Damien Leroux's avatar
Damien Leroux committed
917
        /*all_popl.emplace_back(compute_parental_origins(p, (*qtl_chr)->chr, loci));*/
918

Damien Leroux's avatar
Damien Leroux committed
919
    }
920
921
922
923
924
925
926
927
928
929
    for (auto& pl: all_popl) {
        for (auto& l: pl) {
            (void) *l;
        }
    }
    return assemble_parental_origins_multipop(chr, lk, all_popl, all_pop, false);

}


Damien Leroux's avatar
Damien Leroux committed
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
collection<model_block_type>
compute_dominance_multipop(const collection<population_value>& all_pop,
                           const value<chromosome_value>& chr,
                           const std::vector<locus_key>& lk,
        /*const value<qtl_chromosome_value>& qtl_chr,*/
                           const value<std::vector<double>>& loci,
                           const std::vector<double>& test_loci)
{
    /*DUMP_FILE_LINE();*/
    std::vector<collection<parental_origin_per_locus_type>> all_popl;
    /*DUMP_FILE_LINE();*/
    all_popl.reserve(all_pop.size());
    /*DUMP_FILE_LINE();*/

    /*MSG_DEBUG("loci are " << loci);*/
    /*MSG_DEBUG("given " << lk << ", test loci are " << test_loci);*/

    collection<locus_key> lkeys;
    lkeys.reserve(test_loci.size());
    auto lki = lk.begin();
    for (double l: test_loci) {
        lkeys.emplace_back(*lki++ + l);
    }

    for (auto& p: all_pop) {
        /*DUMP_FILE_LINE();*/
        context_key ck(new context_key_struc(*p, *chr, *loci));
        all_popl.emplace_back(make_collection<Disk>(joint_dominance_at_locus,
                                                    value<context_key>{ck}, lkeys, as_value(false)));
        /*std::cout << all_popl << std::endl;*/
        /*all_popl.emplace_back(compute_parental_origins(p, (*qtl_chr)->chr, loci));*/

    }
    for (auto& pl: all_popl) {
        for (auto& l: pl) {
            (void) *l;
        }
    }
    return assemble_parental_origins_multipop(chr, lk, all_popl, all_pop, false);
}


Damien Leroux's avatar
Damien Leroux committed
972
973
974
975
collection<model_block_type>
compute_parental_origins_multipop(const collection<population_value>& all_pop,
                                  const value<chromosome_value>& chr,
                                  const value<locus_key>& lk,
Damien Leroux's avatar
Damien Leroux committed
976
        /*const value<qtl_chromosome_value>& qtl_chr,*/
Damien Leroux's avatar
Damien Leroux committed
977
978
                                  const value<std::vector<double>>& loci,
                                  const std::vector<double>& test_loci)
Damien Leroux's avatar
Damien Leroux committed
979
{
980
#if 1
Damien Leroux's avatar
Damien Leroux committed
981
982
983
    /*locus_key lk(new locus_key_struc((*qtl_chr)->qtl));*/
    /*int n_qtl = (*qtl_chr)->qtl.size();*/

984
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
985
    std::vector<collection<parental_origin_per_locus_type>> all_popl;
986
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
987
    all_popl.reserve(all_pop.size());
988
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
989

Damien Leroux's avatar
Damien Leroux committed
990
#if 0
Damien Leroux's avatar
Damien Leroux committed
991
    std::vector<double> test_loci;
992
993
    /*MSG_DEBUG("have lk=" << lk);*/
    if (lk && *lk && (*lk)->locus != locus_key_struc::no_locus) {
Damien Leroux's avatar
Damien Leroux committed
994
995
996
997
998
999
1000
        test_loci.reserve(loci->size());
        for (double k: *loci) {
            if (!(*lk)->has(k)) { test_loci.push_back(k); }
        }
    } else {
        test_loci = *loci;
    }
For faster browsing, not all history is shown. View entire blame