probabilities.cc 18.5 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
#include "computations/probabilities.h"
Damien Leroux's avatar
Damien Leroux committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#include <bitset>

static inline
size_t next_power_of_two(size_t n)
{
    --n;
    n|=n>>1;
    n|=n>>2;
    n|=n>>4;
    n|=n>>8;
    n|=n>>16;
    /*n|=n>>32;*/
    ++n;
    return n;
}

#if 0
static inline
size_t transform_index(size_t index, size_t shift_left_mask, size_t shift_right_mask, size_t parent_bits, size_t no_shift_mask)
{
    return ((index & shift_left_mask) << order)
         | ((index & shift_right_mask) >> parent_bits)
         | (index & no_shift_mask);
}
#endif

#define DUMPHEX(_x_) MSG_DEBUG(std::setw(30) << std::left << #_x_ << std::bitset<sizeof(_x_) * 8>(_x_))
#define DUMP(_x_) MSG_DEBUG(std::setw(30) << std::left << #_x_ << (_x_))
struct transform_index {
    /* the number of parents rounded to the next power of two */
    size_t n_parents_pow2;
    /* the number of bits required to encode the parent indices */
    size_t parent_bits;
    /* the number of columns in the resulting matrix (2^(n_loci * n_parents)) */
    size_t n_col;
    /* the amount of bits to shift the new locus left */
    size_t shift_left;
    /* mask for the bits to shift right */
    size_t shift_right_mask;
    /* mask for the bits to shift left */
    size_t shift_left_mask;
    /* mask for the bits to NOT shift */
    size_t no_shift_mask;

    transform_index(size_t n_parents, size_t n_loci, size_t n_order)
        : n_parents_pow2(next_power_of_two(n_parents))
        , parent_bits(__builtin_ctz(n_parents_pow2))
        , n_col(1 << (parent_bits * n_loci))
        , shift_left(parent_bits * n_order)
        , shift_right_mask((1 << (shift_left + parent_bits)) - 1)
        , shift_left_mask(n_parents_pow2 - 1)
        , no_shift_mask(~shift_right_mask)
    {
        /*DUMP(n_parents_pow2);*/
        /*DUMP(n_col);*/
        /*DUMP(parent_bits);*/
        /*DUMPHEX(shift_left);*/
        /*DUMPHEX(shift_left_mask);*/
        /*DUMPHEX(shift_right_mask);*/
        /*DUMPHEX(no_shift_mask);*/
    }

    size_t operator () (size_t index)
    {
        size_t loc_to_shift_left = index & shift_left_mask;
        size_t loc_shifted_left = loc_to_shift_left << shift_left;
        size_t loc_to_shift_right = index & shift_right_mask;
        size_t loc_shifted_right = loc_to_shift_right >> parent_bits;
        size_t loc_unchanged = index & no_shift_mask;
        /*DUMPHEX(index);*/
        /*DUMPHEX(loc_to_shift_left);*/
        /*DUMPHEX(loc_shifted_left);*/
        /*DUMPHEX(loc_to_shift_right);*/
        /*DUMPHEX(loc_shifted_right);*/
        /*DUMPHEX(loc_unchanged);*/
        size_t ret = loc_shifted_left | loc_shifted_right | loc_unchanged;
        /*DUMPHEX(ret);*/
        /*MSG_DEBUG("transform(" << index << ") = " << ret);*/
        return ret;
    }
};


MatrixXd permutation_matrix(size_t n_parents, size_t n_loci, size_t order)
{
    transform_index ti(n_parents, n_loci, order);
    MatrixXd ret = MatrixXd::Zero(ti.n_col, ti.n_col);
    for (size_t index = 0; index < ti.n_col; ++index) {
        ret(ti(index), index) = 1;
    }
    if (n_parents < ti.n_parents_pow2) {
        MatrixXd reduction_matrix_0 = MatrixXd::Identity(n_parents, ti.n_parents_pow2);
        MatrixXd reduction_matrix = reduction_matrix_0;
        for (size_t i = 1; i < n_loci; ++i) {
            MatrixXd tmp = kroneckerProduct(reduction_matrix, reduction_matrix_0);
            reduction_matrix = tmp;
        }
        /*MSG_DEBUG("reduction_matrix_0" << std::endl << reduction_matrix_0);*/
        /*MSG_DEBUG("reduction_matrix" << std::endl << reduction_matrix);*/
        ret = reduction_matrix * ret * reduction_matrix.transpose();
    }
120
    /*MSG_DEBUG("Permutation matrix (n_parents=" << n_parents << ", n_loci=" << n_loci << ", order=" << order << ')' << std::endl << ret);*/
Damien Leroux's avatar
Damien Leroux committed
121
122
123
    return ret;
}

Damien Leroux's avatar
Damien Leroux committed
124
125


126
#if 0
127
128
locus_probabilities_type
locus_prob(qtl_chromosome_value qtl_chr, double noise, generation_value gen,
129
           const std::vector<double>& loci, const MatrixXd& LV)
130
{
Damien Leroux's avatar
Damien Leroux committed
131
    segment_computer_t gc(*gen, qtl_chr, noise, loci);
132
133
134
135
136
	/*MSG_DEBUG("got gc " << gc);*/
    return gc.compute(LV);
}


Damien Leroux's avatar
Damien Leroux committed
137
locus_probabilities_type
Damien Leroux's avatar
Damien Leroux committed
138
locus_probabilities(qtl_chromosome_value qtl_chr, population_value pop,
Damien Leroux's avatar
Damien Leroux committed
139
                    generation_value gen,
140
                    const MatrixXd& mgo,
141
142
                    const selected_qtls_on_chromosome_type&,
                    const std::vector<double>& loci)
Damien Leroux's avatar
Damien Leroux committed
143
{
144
145
146
147
148
149
150
151
152
153
#if 0
    value<MatrixXb> LV = mgo.LV;
    return *make_value</*Disk|*/Sync>(locus_prob,
                                  value<qtl_chromosome_value>{qtl_chr},
                                  value<double>{pop->noise},
                                  value<generation_value>{gen},
                                  value<std::vector<double>>{loci},
                                  LV);
    /*MSG_DEBUG("LV" << std::endl << mgo.LV);*/
    /*locus_probabilities_type ret = gc.compute(mgo.LV);*/
154
    /*MSG_DEBUG("lp" << std::endl << ret);*/
155
156
157
158
159
160
    /*return ret;*/
#else
    generation_rs::segment_computer_t gc = gen->segment_computer(qtl_chr, pop->noise, loci);
	/*MSG_DEBUG("got gc " << gc);*/
    return gc.compute(mgo.LV);
#endif
Damien Leroux's avatar
Damien Leroux committed
161
}
162
#endif
Damien Leroux's avatar
Damien Leroux committed
163
164


165
optim_segment_computer_t
166
make_segment_computer(const chromosome* chr, population_value pop, const locus_key& lk, int ind, size_t)
167
168
169
170
171
172
{
    /*qtl_chromosome qc(ck->chr, lk);*/
    const geno_matrix* gen = pop->gen.get();
    return {*gen, chr, lk, 0/* FIXME: where is the noise parameter? */, pop->get_LV(chr->name, ind)};
}

173
/* FIXME! */
174
175
locus_probabilities_type
locus_probabilities(const context_key& ck, const locus_key& lk,
176
177
                    /*const MatrixXd& mgo,*/
                    int ind,
Damien Leroux's avatar
Damien Leroux committed
178
                    const std::vector<double>& loci)
179
{
180
181
    /*MSG_DEBUG("locus_probabilities(" << ck << ", " << lk << ", " << ind << ", " << loci << ")");*/
    /*MSG_QUEUE_FLUSH();*/
182
    /*const MatrixXd& LV = ck->pop->get_LV(ck->chr->name, ind);*/
183
    auto ret = make_value</*Mem|*/Sync>(make_segment_computer, as_value(ck->chr), as_value(ck->pop), as_value(lk), as_value(ind), as_value(ck->pop->indices[ind]))->compute(loci);
184
    /*MSG_DEBUG("locus_probabilities(" << ck << ", " << lk << ", " << ind << ", " << loci << ")" << std::endl << ret);*/
185
    return ret;
186
187
188
189
190
191
}




#if 0
Damien Leroux's avatar
Damien Leroux committed
192
193
state_to_parental_origin_matrix_type
state_to_parental_origin_matrix(generation_value gen,
194
                                int n_qtl)
Damien Leroux's avatar
Damien Leroux committed
195
{
196
    /*MSG_DEBUG("state_to_parental_origin_matrix " << ((void*)gen));*/
197
    return const_cast<generation_rs*>(gen)->get_state_to_parental_origin(n_qtl);
Damien Leroux's avatar
Damien Leroux committed
198
199
200
201
202
}

parental_origin_type
parental_origin(const locus_probabilities_type& lp, generation_value gen, const selected_qtls_on_chromosome_type& sqoc)
{
203
    auto& s = const_cast<generation_rs*>(gen)->get_state_to_parental_origin(sqoc.size());
Damien Leroux's avatar
Damien Leroux committed
204
205
    return s * lp;
}
206
#endif
Damien Leroux's avatar
Damien Leroux committed
207

208
209
210
211
212

template <typename T>
void dump_hash(const std::string& n, const T& x)
{
    std::hash<T> h;
213
    MSG_DEBUG("hash(" << n << ") = " << h(x));
214
215
}

216
#if 0
217
218
219
220
221
parental_origin_type
parental_origin2(qtl_chromosome_value qtl_chr,
			 	 population_value pop,
			 	 generation_value gen,
				 const state_to_parental_origin_matrix_type& stfopom,
222
		 		 const MatrixXd& mgo,
223
224
		 		 const selected_qtls_on_chromosome_type& sqoc,
                 const std::vector<double>& loci)
225
226
{
    /*auto& s = const_cast<generation_rs*>(gen)->compute_state_to_parental_origin(sqoc.size());*/
227
    /*MSG_DEBUG("stfopom: " << stfopom.innerSize() << "," << stfopom.outerSize());*/
228
    auto lp = l ocus_probabilities(qtl_chr, pop, gen, mgo, sqoc, loci);
229
230
231
232
233
234
235
236
237
238
239
240
    /*MSG_DEBUG("lp: " << lp.innerSize() << "," << lp.outerSize());*/
    if (0) {
        md5_digest md5;
        md5 << mgo;
        std::string m = md5;
        MSG_DEBUG("mgo: MD5 = " << m << std::endl << mgo.LV);
        dump_hash("stfopom", stfopom);
        dump_hash("mgo", mgo);
        dump_hash("sqoc", sqoc);
        dump_hash("loci", loci);
    }
    return stfopom * lp;
241
}
242
243
244
245

parental_origin_type
parental_origin2(const context_key& ck,
                 const locus_key& lk,
246
                 const MatrixXd& mgo)
247
{
248
    auto lp = make_value<Disk>(l ocus_probabilities, value<context_key>{ck}, value<locus_key>{lk}, value<MatrixXd>{mgo}, as_value(ck->loci));
249
250
251
    auto stfopom = make_value<Mem>(compute_state_to_parental_origin, value<context_key>{ck}, value<locus_key>{lk});
    return (*stfopom) * (*lp);
}
252
#endif
253
254
255
256


/* Compute the joint probability of selected loci UNION new locus */
geno_prob_type
Damien Leroux's avatar
Damien Leroux committed
257
joint_geno_prob_at_locus(const context_key& ck, const locus_key& lk)
258
{
Damien Leroux's avatar
Damien Leroux committed
259
260
261
//    DUMP_FILE_LINE();
//    MSG_DEBUG("* Computing geno prob in " << ck << " given " << lk);
//    MSG_DEBUG("=====================================================");
Damien Leroux's avatar
Damien Leroux committed
262
    auto it = ck->locus_indices.find(lk->locus);
263
264
    if (it == ck->locus_indices.end()) {
        /* FIXME: whine */
Damien Leroux's avatar
Damien Leroux committed
265
266
//        MSG_ERROR("THIS LOCUS DOESN'T EXIST IN THIS CONTEXT! ck.loci=[" << ck->loci << "] locus=" << lk->locus, "");
//        MSG_QUEUE_FLUSH();
Damien Leroux's avatar
Damien Leroux committed
267
        throw 0;
268
269
        return {};
    }
Damien Leroux's avatar
Damien Leroux committed
270
//    MSG_DEBUG_INDENT;
271
272
    size_t loc_idx = it->second;
    value<context_key> vck = ck;
Damien Leroux's avatar
Damien Leroux committed
273
    value<locus_key> vlk = lk->parent;
274
275
276
    /*collection<MatrixXd>*/
        /*vmgo = make_collection<Disk>(population_marker_obs,*/
                                     /*vck, range<int>(0, ck->pop->size(), 1));*/
277
    collection<locus_probabilities_type>
278
279
280
        /*alp = make_collection<Disk>(locus_probabilities, vck, vlk, vmgo, as_value((*vck)->loci));*/
        alp = make_collection<Disk>(locus_probabilities,
                                    vck, vlk, range<int>(0, ck->pop->size(), 1), as_value((*vck)->loci));
281
    for (auto& x: alp) {
Damien Leroux's avatar
Damien Leroux committed
282
283
        /*MSG_DEBUG(x);*/
        (void)*x;
284
    }
285
286
287
288
289

    geno_prob_type ret;

    ret.row_labels = alp[0]->row_labels;
    ret.data.resize(alp[0]->innerSize(), alp.size());
Damien Leroux's avatar
Damien Leroux committed
290
//    DUMP_FILE_LINE();
291
    for (size_t ind = 0; ind < ck->pop->size(); ++ind) {
Damien Leroux's avatar
Damien Leroux committed
292
293
294
//        MSG_DEBUG("n.cols = " << ret.data.outerSize() << " vs " << alp[ind]->data.outerSize());
//        MSG_DEBUG("col.size = " << ret.data.innerSize() << " vs " << alp[ind]->data.innerSize());
//        MSG_QUEUE_FLUSH();
295
296
297
        ret.data.col(ind) = alp[ind]->data.col(loc_idx);
    }

Damien Leroux's avatar
Damien Leroux committed
298
299
300
301
//    DUMP_FILE_LINE();
//    MSG_DEBUG("computed " << ck << ' ' << lk);
//    MSG_DEBUG(MATRIX_SIZE(ret));
//    MSG_DEBUG("-----------------------------------------------------");
302

Damien Leroux's avatar
Damien Leroux committed
303
304
305
306
    /*if (lk && lk->locus != locus_key_struc::no_locus) {*/
    if (!(*vlk)->is_empty()) {
        /*std::vector<double> test_loci = {lk->locus};*/
        /*context_key ck2(new context_key_struc(ck->pop, ck->chr, test_loci));*/
307
308
309
        auto pop_pred = make_value<Disk|Sync>(joint_geno_prob_at_locus,
                                              as_value(ck), vlk);
        /*(void)*pop_pred;*/
Damien Leroux's avatar
Damien Leroux committed
310
311
312
313
//        DUMP_FILE_LINE();
//        MSG_DEBUG("previous geno_prob " << ck << ' ' << lk->parent << '+' << lk->locus);
//        MSG_DEBUG(MATRIX_SIZE((*pop_pred)));
//        MSG_DEBUG("-----------------------------------------------------");
314
315
        /*MatrixXd ones = MatrixXd::Ones(ck->gen->get_unique_labels().size(), 1);*/
        MatrixXd ones = MatrixXd::Ones(ret.data.rows() / pop_pred->data.rows(), 1);
Damien Leroux's avatar
Damien Leroux committed
316
        MatrixXd hada_mult = kroneckerProduct(pop_pred->data, ones);
Damien Leroux's avatar
Damien Leroux committed
317
318
//        MSG_DEBUG("Previous geno_prob after kronecker to prepare for the Hadamard product:" << std::endl << hada_mult);
//        MSG_DEBUG("-----------------------------------------------------");
Damien Leroux's avatar
Damien Leroux committed
319
320
        MatrixXd tmp = (ret.data.array() * hada_mult.array()).matrix();
        ret.data = tmp;
321
322
    }

Damien Leroux's avatar
Damien Leroux committed
323
324
325
326
//    DUMP_FILE_LINE();
//    MSG_DEBUG("joint probabilities " << ck << ' ' << lk);
//    MSG_DEBUG(MATRIX_SIZE(ret));
//    MSG_DEBUG("=====================================================");
327

Damien Leroux's avatar
Damien Leroux committed
328
329
330
//    DUMP_FILE_LINE();
//    MSG_DEBUG_DEDENT;
//    MSG_QUEUE_FLUSH();
331
332
333
334
    return ret;
}


335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
/* Compute the joint probability of selected loci UNION new locus */
parental_origin_per_locus_type
joint_dominance_at_locus(const context_key& ck, const locus_key& lk, bool)
{
    value<context_key> vck{ck};
    value<locus_key> vlk{lk};
    value<locus_key> vlkp{lk->parent};
    auto pop = make_value<Disk>(joint_geno_prob_at_locus, vck, vlk);
    auto stdom = make_value<Mem>(compute_state_to_dominance, vck, vlkp);

    /*(void) *pop;*/
    /*(void) *stfopom;*/
    /*MSG_DEBUG("joint_dominance_at_locus");*/
    /*MSG_DEBUG("******* " << MATRIX_SIZE(*pop));*/
    /*MSG_DEBUG(pop);*/
    /*MSG_DEBUG("******* " << MATRIX_SIZE(*stdom));*/
    /*MSG_DEBUG(stdom);*/
    /*MSG_QUEUE_FLUSH();*/

    return ((*stdom) * (*pop)).transpose();
}


358
359
/* Compute the joint probability of selected loci UNION new locus */
parental_origin_per_locus_type
Damien Leroux's avatar
Damien Leroux committed
360
joint_parental_origin_at_locus(const context_key& ck, const locus_key& lk)
361
362
363
{
    value<context_key> vck{ck};
    value<locus_key> vlk{lk};
Damien Leroux's avatar
Damien Leroux committed
364
365
366
    value<locus_key> vlkp{lk->parent};
    auto pop = make_value<Disk>(joint_geno_prob_at_locus, vck, vlk);
    auto stfopom = make_value<Mem>(compute_state_to_parental_origin, vck, vlkp);
367

368
369
    /*(void) *pop;*/
    /*(void) *stfopom;*/
Damien Leroux's avatar
Damien Leroux committed
370
371
372
373
374
375
    /*MSG_DEBUG("joint_parental_origin_at_locus");*/
    /*MSG_DEBUG("******* " << MATRIX_SIZE(*pop));*/
    /*MSG_DEBUG(pop);*/
    /*MSG_DEBUG("******* " << MATRIX_SIZE(*stfopom));*/
    /*MSG_DEBUG(stfopom);*/
    /*MSG_QUEUE_FLUSH();*/
376
377
378
379
380

    return ((*stfopom) * (*pop)).transpose();
}


381

Damien Leroux's avatar
Damien Leroux committed
382
383
384
385

collection<parental_origin_per_locus_type>
parental_origin_per_locus(const collection<parental_origin_type>& apo)
{
386
387
388
389
390
391
392
393
    /* TODO :
     * - séparer la préparation (tout avant le premier for) (sync|mem)
     * - éliminer le for externe
     * - retourner un SEUL parental_origin_per_locus_type
     * - virer LD de la fin de la fonction, il faut que ça soit un calcul à part
     * - remplacer l'argument par context_key, locus_key, locus_index
     */

394
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
395
396
397
398
    auto& tmp = apo[0];
    size_t n_ind = apo.size();
    size_t n_par = tmp->innerSize();
    size_t n_loc = tmp->outerSize();
399
400
401
    /*MSG_DEBUG("n_ind " << n_ind);*/
    /*MSG_DEBUG("n_par " << n_par);*/
    /*MSG_DEBUG("n_loc " << n_loc);*/
402
    /*MSG_DEBUG("apo" << std::endl << apo);*/
Damien Leroux's avatar
Damien Leroux committed
403

404
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
405
406
407
408
409
    /* labels for individuals */
    std::vector<int> individuals;
    individuals.resize(n_ind);
    size_t index = 0;
    std::generate(individuals.begin(), individuals.end(), [&] () { return index++; });
410
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
411
412
    /* labels for parents */
    std::vector<std::vector<char>> parents = tmp->row_labels;
413
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
414
415
    /* loci */
    std::vector<double> loci = tmp->column_labels;
416
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
417
418
419
    /* return value */
    collection<parental_origin_per_locus_type> ret(loci.size());

420
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
421
422
423
    /* first round : init matrices */
    for (size_t l = 0; l < n_loc; ++l) {
        /*labelled_matrix<MatrixXd, int, std::vector<char>>& lm = ret[loci[l]];*/
424
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
425
426
427
428
429
        ret[l] = new unique_value<parental_origin_per_locus_type>();
        labelled_matrix<MatrixXd, int, std::vector<char>>& lm = *ret[l];
        lm.row_labels = individuals;
        lm.column_labels = parents;
        lm.data.resize(n_ind, n_par);
430
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
431
        for (int ind: individuals) {
432
            /*MSG_DEBUG("ind=" << ind << " l=" << l << " lm.data$size(" << lm.data.innerSize() << ',' << lm.data.outerSize() << ") apo[ind]->data$size(" << apo[ind]->data.innerSize() << ',' << apo[ind]->data.outerSize() << ')');*/
Damien Leroux's avatar
Damien Leroux committed
433
434
435
            lm.data.row(ind) = apo[ind]->data.col(l);
        }
    }
436
#ifdef NEED_TO_TRANSMIT_QTL_GENERATION_NAME_AND_CHROMOSOME_PTR_TO_HANDLE_LD_PROPERLY
437
438
439
440
441
442
443
444
445
446
    if (active_settings->ld_data.size()) {
        const LD_matrices& ld = active_settings->ld_data.begin()->second;
        /*assert(n_loc == ld.ld.size() && "WRONG LD SIZE");*/
        /*MSG_DEBUG("n_loc = " << n_loc << std::endl);*/
        /*MSG_DEBUG("ld.size = " << ld.ld.size() << std::endl);*/
        for (size_t l = 0; l < n_loc; ++l) {
            /*MSG_DEBUG("LD " << l << std::endl << ld.ld[l]);*/
            ret[l]->data = ret[l]->data * ld.ld[l].transpose();
        }
    }
447
#endif
448
    /*DUMP_FILE_LINE();*/
Damien Leroux's avatar
Damien Leroux committed
449
450
451
452
    return ret;
}


453

454
455
456
457
458
459
460
461
462
/*
 * TODO
 * in cache2:
 * - when caching is Disk only, have a registry record all CURRENTLY COMPUTING tasks:
 *   task registers upon creation, unregister upon completion, so that two asynchronous
 *   Disk tasks don't overwrite each other.
 */

#if 0
Damien Leroux's avatar
Damien Leroux committed
463
464
collection<parental_origin_per_locus_type>
compute_parental_origins(const value<population_value>& pop,
465
466
                         const value<chromosome_value>& chr,
                         const value<std::vector<double>>& loci)
Damien Leroux's avatar
Damien Leroux committed
467
{
468
    DUMP_FILE_LINE();
Damien Leroux's avatar
Damien Leroux committed
469
470
471
    value<qtl_chromosome_value> qtl_chr
        = qtl_chromosome_by_name((*chr)->name);

472
    DUMP_FILE_LINE();
Damien Leroux's avatar
Damien Leroux committed
473
474
     value<selected_qtls_on_chromosome_type> sqoc = selected_qtls_on_chromosome(*qtl_chr);

475
    DUMP_FILE_LINE();
Damien Leroux's avatar
Damien Leroux committed
476
477
    value<generation_value> qtl_gen = qtl_generation(*pop);

478
479
480
    /*DUMP_FILE_LINE();*/
    /*collection<pedigree_type> ap*/
        /*= make_collection(pedigree, pop, qtl_gen, individual_range(*pop));*/
Damien Leroux's avatar
Damien Leroux committed
481

482
483
484
485
    DUMP_FILE_LINE();
    value<const pop_mgo_data*> pmd = new pop_mgo_data(*pop);

    DUMP_FILE_LINE();
486
    collection<MatrixXd> apmo
487
488
        = make_collection<Sync|Disk>(population_marker_obs,
                                    pop, chr, range<int>(0, (*pop)->size(), 1), pmd);
489
    /*MSG_DEBUG("apmo" << std::endl << apmo);*/
Damien Leroux's avatar
Damien Leroux committed
490

491
    /*collection<locus_probabilities_type> alp*/
492
        /*= make_collection(l ocus_probabilities, qtl_chr, pop, gen, apmo, sqoc);*/
493

494
    DUMP_FILE_LINE();
495
	value<state_to_parental_origin_matrix_type>
496
		stfopom = make_value<Disk>(state_to_parental_origin_matrix, qtl_gen, value<int>(sqoc->size()));
497

498
499
    /*MSG_DEBUG("nqtl=" << sqoc->size());*/
    /*MSG_DEBUG("[cpo]stfopom: " << stfopom->innerSize() << "," << stfopom->outerSize());*/
Damien Leroux's avatar
Damien Leroux committed
500

501
    DUMP_FILE_LINE();
Damien Leroux's avatar
Damien Leroux committed
502
    collection<parental_origin_type> apo
503
        = make_collection<Disk>(parental_origin2, qtl_chr, pop, qtl_gen, stfopom, apmo, sqoc, loci);
504
    /*MSG_DEBUG("apo" << std::endl << apo);*/
Damien Leroux's avatar
Damien Leroux committed
505

506
    DUMP_FILE_LINE();
Damien Leroux's avatar
Damien Leroux committed
507
508
    return parental_origin_per_locus(apo);
}
509
#endif
Damien Leroux's avatar
Damien Leroux committed
510

511

512
#include "output_impl.h"