jobs.cc 24.1 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/>.
 */

18
19
#include "dispatch.h"
#include "output.h"
20
#include "pedigree.h"
21
/*#include "factor_var4.h"*/
22
#include "graphnode2.h"
23
#include "tensor.h"
24

25
bool job_check_fourcc(ifile& ifs, const char* fourcc)
26
27
28
29
30
31
32
33
34
{
    if (check_fourcc(ifs, fourcc)) {
        MSG_ERROR("NOT A JOB DATA FILE OR HAS BEEN CORRUPTED", "Re-run this tool.");
        return true;
    }
    return false;
}


35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
inline
std::ostream&
operator << (std::ostream& os, const std::pair<char, char>& p) { return os << ((int) p.first) << ((int) p.second); }


inline
std::ostream&
operator << (std::ostream& os, const std::map<bn_label_type, double>& obs)
{
    os << '{';
    auto i = obs.begin(), j = obs.end();
    if (i != j) {
        os << i->first << '=' << i->second;
        for (++i; i != j; ++i) {
            os << ' ' << i->first << '=' << i->second;
        }
    }
    return os << '}';
}


56
57
void
dispatch_geno_probs(
58
59
60
61
62
        const std::vector<gencomb_type>& lc,  /* table parent states -> spawnling states */
        const std::vector<label_type>& labels, /* spawnling labels */
        const std::map<label_type, double>& geno_probs, /* spawnling genotype probabilities */
        size_t ind, /* spawnling number */
        std::map<size_t, std::vector<double>>& state_prob) /* all computed locus vectors to fetch the parents' states */
63
{
64
/*
65
    scoped_indent _(MESSAGE("[dispatchGP] "));
66
    MSG_DEBUG("ind " << ind);
67
68
69
70
    MSG_DEBUG("lc " << lc);
    MSG_DEBUG("labels " << labels);
    MSG_DEBUG("geno_probs " << geno_probs);
    MSG_DEBUG("state_prob " << state_prob);
71
*/
72
//    MSG_QUEUE_FLUSH();
73
74
75
76
77
78
79
80
81
82
83
84
85
    std::vector<double> ret(lc.size());
    std::map<label_type, double> norms;
    for (size_t i = 0; i < ret.size(); ++i) {
        ret[i] = 0;
        for (const auto& e: lc[i]) {
            double coef = e.coef;
            for (const auto& k: e.keys) {
                coef *= (state_prob[k.parent][k.state] > 0 ? 1. : 0.);
            }
            ret[i] += coef;
        }
        norms[labels[i]] += ret[i];
    }
86
87
88
89
90
91
92
93
    for (auto& kv: norms) if (kv.second) {
        auto it = geno_probs.find(kv.first);
        if (it == geno_probs.end()) {
            kv.second = 0;
        } else {
            kv.second = it->second / kv.second;
        }
    }
94
95
96
97
98
99
    for (size_t i = 0; i < ret.size(); ++i) { ret[i] *= norms[labels[i]]; }
    state_prob[ind].swap(ret);
}



100
101
102
103
104
105
106
107
108
109
110
111
112
113
std::vector<bn_label_type>
get_domain(std::map<var_vec, genotype_comb_type>& domains, variable_index_type var)
{
    std::vector<bn_label_type> ret;
    const auto& dom = domains[{var}];
    ret.reserve(dom.size());
    for (const auto& e: dom) {
        ret.push_back(e.keys.keys.front().state);
    }
    return ret;
}



114
115
116
117
118
119
120
121
122
genotype_comb_type
build_evidence(variable_index_type id, const std::map<bn_label_type, double>& obs)
{
    genotype_comb_type tmp;
    for (const auto& kv: obs) {
        tmp.m_combination.emplace_back(genotype_comb_type::element_type{{{id, kv.first}}, kv.second});
    }
    return tmp;
}
123

124
125
126
std::map<std::string, std::pair<count_jobs_type, job_type>>
job_registry = {
    {"dummy-one", {
127
128
        [](bn_settings_t*) { return 1; },
        [](bn_settings_t*, size_t i) { MSG_ERROR("running #" << i, ""); return true; }
129
130
    }},
    {"dummy", {
131
132
        [](bn_settings_t* settings) { return settings->count_markers(); },
        [](bn_settings_t*, size_t i) { MSG_ERROR("running #" << i, ""); return true; }
133
    }},
134
    {"compute-alleles-per-marker", {
135
136
        [](bn_settings_t* settings) { return settings->count_markers(); },
        [](bn_settings_t* settings, size_t mark)
137
        {
Damien Leroux's avatar
Damien Leroux committed
138
139
140
141
142
143
            std::string filename = settings->job_filename("compute-alleles-per-marker", mark);
            
            if (check_file(filename, false, true, false)) {
                return true;
            }

144
            /*const std::string& marker_name = settings->marker_names[mark];*/
145
146
147
148
149
150
151
152
            auto get_obs
                = [&] (const std::string& gen, size_t i)
                {
                    const auto& pop_obs = settings->observed_mark.find(gen)->second;
                    const marker_data& obs = pop_obs.observations;
                    auto ret = settings->marker_observation_specs
                        .find(pop_obs.format_name)->second
                        .score(obs.data.find(settings->marker_names[mark])->second[i]);
153
                    /*MSG_DEBUG("FETCHING OBS[" << marker_name << ':' << gen << ':' << i << "] => " << ret);*/
154
155
156
                    return ret;
                };

Damien Leroux's avatar
Damien Leroux committed
157
            ofile dump(filename);
158
159
160
161
162
163
164
165
166
167
168
169
170
171

            std::set<char> used_alleles;
            for (const auto& kv: settings->observed_mark) {
                if (settings->marker_observation_specs.find(kv.second.format_name)->second.domain != ODAllele) {
                    continue;
                }
                const std::vector<int>&
                    gen_ind = settings->pedigree.individuals_by_generation_name.find(kv.first)->second;
                for (size_t i = 0; i < gen_ind.size(); ++i) {
                    for (const auto& score: get_obs(kv.first, i)) {
                        used_alleles.insert(score.first);
                        used_alleles.insert(score.second);
                    }
                }
172
                /*MSG_DEBUG("Used alleles in generation " << kv.first << ": " << used_alleles.size());*/
173
            }
174
175
176
            if (used_alleles.size() == 0) {
                used_alleles.insert(1);
            }
177
178
179
180
181
            rw_base()(dump, used_alleles);
            return true;
        }
    }},
    {"compute-factor-graphs", {
182
183
        [](bn_settings_t* settings) { return settings->unique_n_alleles.size(); },
        [](bn_settings_t* settings, size_t n)
184
185
        {
            std::vector<size_t> unique_n_alleles;
186
            {
187
                ifile unafs(settings->job_filename("unique_n_alleles", 0));
188
189
                rw_base()(unafs, unique_n_alleles);
            }
Damien Leroux's avatar
Damien Leroux committed
190
191
192
193
            std::string filename = settings->job_filename("factor-graph", unique_n_alleles[n]);
            if (check_file(filename, false, true, false)) {
                return true;
            }
194
195
            /*MSG_DEBUG("Pedigree" << std::endl << settings->pedigree);*/
            /*MSG_DEBUG("COMPUTING FACTOR GRAPH FOR " << unique_n_alleles[n] << " ALLELES.");*/
196
            /*factor_graph fg(settings->pedigree, unique_n_alleles[n], settings->noise);*/
197
            auto fg = factor_graph_type::from_pedigree(settings->pedigree.items, unique_n_alleles[n], settings->output_mode & bn_settings_t::OutputGametes, settings->input_generations, settings->output_generations);
198
            /*MSG_DEBUG("COMPUTED FACTOR GRAPH" << std::endl << fg);*/
199
200
            /*fg.finalize();*/
            /*fg.save(settings->job_filename("factor-graph", unique_n_alleles[n]));*/
201
            /*fg->dump();*/
202
            //fg->to_image(MESSAGE("factor-graph-" << unique_n_alleles[n]), "png");
203
            auto I = instance(fg);
Damien Leroux's avatar
Damien Leroux committed
204
            ofile of(filename);
205
206
            I->file_io(of);
            rw_comb<int, bn_label_type>()(of, fg->domains);
207
208
209
            return true;
        }
    }},
210
    {"compute-LV", {
211
212
        [](bn_settings_t* settings) { return settings->count_markers(); },
        [](bn_settings_t* settings, size_t mark)
213
        {
Damien Leroux's avatar
Damien Leroux committed
214
215
216
217
218
            std::string filename = settings->job_filename("compute-LV", mark);
            if (check_file(filename, false, true, false)) {
                return true;
            }

219
            const std::string& marker_name = settings->marker_names[mark];
220
            std::set<char> used_alleles;
221
222
            std::map<char, char> allele_obs_to_idx;
            {
223
                ifile apm(settings->job_filename("compute-alleles-per-marker", mark));
224
225
226
227
                rw_base()(apm, used_alleles);
                char idx = 0;
                for (char c: used_alleles) {
                    allele_obs_to_idx[c] = idx++;
228
                    /*MSG_DEBUG("ALLELE CONVERT " << c << " => " << ((int) allele_obs_to_idx[c]));*/
229
230
                }
            }
231
232
233
            /*graph_type fg;*/
            /*fg.load(settings->job_filename("factor-graph", used_alleles.size()));*/
            ifile ifs(settings->job_filename("factor-graph", used_alleles.size()));
234
            std::unique_ptr<instance_type> I(new instance_type());
235
236
237
            std::map<var_vec, genotype_comb_type> domains;
            I->file_io(ifs);
            rw_comb<int, bn_label_type>()(ifs, domains);
238
239
240
241
242
243
244
245
246

            std::map<char, int> allele_obs2bn;
            for (char c: used_alleles) { allele_obs2bn[c] = allele_obs2bn.size(); }

            auto make_obs
                = [&] (const population_marker_observation& pop_obs, size_t i, size_t n)
                {
                    std::map<bn_label_type, double> ind_obs;
                    const auto& obs_spec = settings->marker_observation_specs.find(pop_obs.format_name)->second;
247
                    const auto& obsdat = pop_obs.observations.data.find(settings->marker_names[mark])->second[i];
248
                    std::vector<std::pair<char, char>>
249
                        obs_vec = obs_spec.score(obsdat);
250
251
                    genotype_comb_type::key_type key;
                    key.parent = n;
252
                    for (const auto& label: get_domain(domains, n)) {
253
254
                        ind_obs[label] = 0;
                    }
255
256
257
258
259
//                     MSG_DEBUG("i " << i);
//                     MSG_DEBUG("obs " << obsdat);
//                     MSG_DEBUG("obs_vec " << obs_vec);
//                     MSG_DEBUG("ind_dom = " << get_domain(domains, n));
//                     MSG_DEBUG("ind_obs = " << ind_obs);
260
261
262
263
264
265
266
                    /*MSG_DEBUG("obs spec");*/
                    /*for (const auto& s: obs_spec.scores) {*/
                        /*MSG_DEBUG("* " << s.first);*/
                        /*for (const auto& p: s.second) {*/
                            /*MSG_DEBUG("  " << p.first << p.second);*/
                        /*}*/
                    /*}*/
267
//                     MSG_DEBUG("obs for #" << n << " '" << pop_obs.observations.data.find(settings->marker_names[mark])->second[i] << "': " << obs_vec);
268
269
                    if (obs_spec.domain == ODAllele) {
                        for (const auto& obs: obs_vec) {
270
                            for (const auto& label: get_domain(domains, n)) {
271
                                if (label.first_allele == allele_obs_to_idx[obs.first] && label.second_allele == allele_obs_to_idx[obs.second]) {
272
273
274
275
276
277
                                    ind_obs[label] = 1;
                                }
                            }
                        }
                    } else {
                        for (const auto& obs: obs_vec) {
278
                            for (const auto& label: get_domain(domains, n)) {
279
280
281
282
283
284
                                if (label.first == obs.first && label.second == obs.second) {
                                    ind_obs[label] = 1;
                                }
                            }
                        }
                    }
285
//                     MSG_DEBUG("ind_obs = " << ind_obs);
286
287
288
                    return ind_obs;
                };

289
290
            /*auto instance = fg.instance();*/
            /*instance.clear_evidence();*/
291
            /*std::vector<std::shared_ptr<message_type>> evidence;*/
292

293
            for (const auto& pop_obs: settings->observed_mark) {
294
                /*MSG_DEBUG("Obs in generation " << pop_obs.first);*/
295
296
                size_t n = 0;
                for (size_t i: settings->pedigree.individuals_by_generation_name.find(pop_obs.first)->second) {
297
298
                    int id = settings->pedigree.ind2id(i);
                    auto obs = make_obs(pop_obs.second, n, id);
299
300
                    /*evidence.emplace_back(std::make_shared<message_type>(message_type{build_evidence(id, obs)}));*/
                    I->add_evidence(id, build_evidence(id, obs));
301
                    MSG_DEBUG("observations for individual #" << n << " (id #" << id << ") = " << obs);
302
                    /*MSG_DEBUG("id=" << id << " id=" << id << " evidence " << evidence);*/
303
                    double accum = 0;
304
305
                    /*for (const auto& kv: obs) {*/
                        /*evidence.force_set({{node, kv.first}}, kv.second);*/
306
                        /*instance.evidence(node, kv.first, kv.second);*/
307
308
309
                        /*accum += kv.second;*/
                    /*}*/
                    if (0 && /* FIXME */ accum == 0) {
310
                        msg_handler_t::cout() << std::endl;
311
                        MSG_ERROR("Inconsistent observations for individual #" << n << " of generation " << pop_obs.first << " on marker " << marker_name, SPELL_STRING("Remove inconsistencies in your genotype data file for generation " << pop_obs.first));
312
313
314
315
316
317
                        return false;
                    }
                    ++n;
                }
            }

318
319
320
            /*instance.evidence(evidence);*/
            /*MSG_DEBUG("EVIDENCE");*/
            /*MSG_DEBUG("" << instance.evidence().dump());*/
321

322
323
324
325
326
327
            /*std::vector<graph_type::query_operation_type> marginals;*/
            std::map<edge_type, std::shared_ptr<message_type>> messages;
            /*for (const auto& i: settings->pedigree.items) {*/
                /*auto q = fg.build_query_operation({i.id});*/
                /*marginals.insert(marginals.end(), q.begin(), q.end());*/
            /*}*/
328
            /*std::map<int, genotype_comb_type> marginals;*/
329
330
331
            /*fg.compute_messages(evidence.begin(), evidence.end(), messages);*/
            /*fg.compute_full_factor_state(messages, marginals);*/
            message_type output = I->compute(0, domains);
332
            MSG_DEBUG("OUTPUT:");
333
            std::map<variable_index_type, std::map<label_type, double>> marginals;
334
335
            for (const auto& t: output) {
                MSG_DEBUG(" * " << t << std::endl);
336
337
338
339
340
341
342
                for (const auto& e: t) {
                    for (const auto& k: e.keys) {
                        label_type l = {k.state.first, k.state.second};
                        marginals[k.parent][l] += e.coef;
                    }
                }
            }
343

344
            if (settings->output_mode & bn_settings_t::OutputOnePoint) {
345
                std::string op_dirname = SPELL_STRING(settings->work_directory << '/' << settings->pop_name << ".1-point");
346
                ensure_directories_exist(op_dirname);
347
                std::string op_filename = SPELL_STRING(op_dirname << '/' << settings->pop_name << ".pedigree-and-probabilities." << marker_name << ".csv");
348
                std::ofstream ofs(op_filename);
349
350
351
352
353
                ofs << "Gen;Id;P1;P2;Prob" << std::endl;
                for (const auto& pi: settings->pedigree.items) {
                    ofs << pi.gen_name << ';' << pi.id << ';' << pi.p1 << ';' << pi.p2 << ";{";
                    bool first = true;
                    for (const auto& kv: marginals[pi.id]) {
354
                        ofs << (first ? "\"" : ", \"") << kv.first << "\": " << kv.second;
355
356
357
358
                        first = false;
                    }
                    ofs << '}' << std::endl;
                }
359
                settings->closing_messages.push_back(SPELL_STRING(GREEN << "1-point Parental Origin Probabilities for marker `" << marker_name << "' written in file " << WHITE << op_filename << NORMAL << std::endl));
360
361
            }

362
363
364
365
366
367
368
369
370
            if (settings->output_mode & bn_settings_t::OutputPopData) {
                MSG_DEBUG("PROJECTED MARGINALS");
                for (const auto &km: marginals) {
                    std::stringstream ss;
                    ss << " * " << km.first;
                    for (const auto &kv: km.second) {
                        ss << ' ' << kv.first << '=' << kv.second;
                    }
                    MSG_DEBUG(ss.str());
371
                }
372
                //            MSG_QUEUE_FLUSH();
373

374
                std::map<size_t, std::vector<double>> state_prob, output_prob;
Damien Leroux's avatar
WIP.    
Damien Leroux committed
375
                std::map<int, VectorXd> gamete_prob;
376

377
378
379
380
381
382
383
384
385
386
                for (size_t ind = 1; ind < settings->pedigree.tree.m_ind_number_to_node_number.size(); ++ind) {
                    size_t node = settings->pedigree.tree.m_ind_number_to_node_number[ind];
                    auto gen = settings->pedigree.get_gen(ind);
                    const auto &labels = gen->labels;
                    const auto &LC = settings->pedigree.LC[node];
                    int variable = settings->pedigree.ind2id(ind);
                    dispatch_geno_probs(LC, labels, marginals[variable], node, state_prob);
                    if (std::find(settings->output_generations.begin(), settings->output_generations.end(),
                                  gen->name) != settings->output_generations.end()) {
                        output_prob[node] = state_prob[node];
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
                        if ((settings->output_mode & bn_settings_t::OutputGametes) && !settings->pedigree.items[ind - 1].is_ancestor()) {
                            auto& mat = marginals[-variable];
                            MSG_DEBUG("marginals for gamete " << mat);
                            std::map<label_type, double> prob;
                            for (const auto& kv: mat) {
                                prob[{kv.first.first(), kv.first.second()}] = kv.second;
                            }
                            if (settings->pedigree.items[ind - 1].is_dh()) {
                                gamete_prob[ind].resize(2);
                                gamete_prob[ind] << prob[{GAMETE_L, 0}], prob[{GAMETE_R, 0}];
                            } else {
                                gamete_prob[ind].resize(4);
                                gamete_prob[ind] << prob[{GAMETE_L, GAMETE_L}], prob[{GAMETE_L, GAMETE_R}], prob[{GAMETE_R, GAMETE_L}], prob[{GAMETE_R, GAMETE_R}];
                            }
                            MSG_DEBUG("GAMETE #" << ind << " prob " << gamete_prob[ind - 1].transpose());
damien's avatar
damien committed
402
                        }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
403
                    }
404
405
                }

406
                MSG_DEBUG("LOCUS VECTORS");
407

408
409
410
                for (const auto &kv: state_prob) {
                    MSG_DEBUG(" * " << kv.first << ' ' << kv.second);
                }
411

Damien Leroux's avatar
Damien Leroux committed
412
                ofile ofs(filename);
413
                rw_base()(ofs, output_prob);
Damien Leroux's avatar
WIP.    
Damien Leroux committed
414
415
                
                if (settings->output_mode & bn_settings_t::OutputGametes) {
damien's avatar
damien committed
416
                    ofile og(settings->job_filename("gametes-LV", mark));
Damien Leroux's avatar
WIP.    
Damien Leroux committed
417
418
                    rw_base()(og, gamete_prob);
                }
419
            }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
420
            
421
422
            return true;

423
424
        }
    }},
425
    {"collect-LV", {
426
427
        [](bn_settings_t*) { return 1; },
        [](bn_settings_t* settings, size_t)
428
        {
429
430
431
432
            settings->load_full_pedigree_data();
            const pedigree_type& ped = settings->pedigree;
            std::map<std::string, std::map<size_t, std::vector<double>>> mark_state_prob;
            size_t n_mark = settings->count_markers();
Damien Leroux's avatar
WIP.    
Damien Leroux committed
433
434
            bool have_gametes = settings->output_mode & bn_settings_t::OutputGametes;
            gamete_LV_database gamete_LV;
Damien Leroux's avatar
WIP.    
Damien Leroux committed
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
//             for (size_t m = 0; m < n_mark; ++m) {
//                 ifile ifs(settings->job_filename("compute-LV", m));
//                 rw_base() (ifs, mark_state_prob[settings->marker_names[m]]);
//                 if (have_gametes) {
//                     ifile ig(settings->job_filename("gametes-LV", m));
//                     std::map<int, MatrixXd> gp;
//                     rw_base() (ig, gp);
//                     for (const auto& kv: gp) {
//                         if (kv.second.size() == 0) {
//                             MSG_DEBUG("HAVE EMPTY GAMETE LV " << kv.first);
//                         } else {
//                             if (kv.second.sum() == 0) {
//                                 MSG_DEBUG("HAVE ZERO GAMETE LV " << kv.first);
//                             } else {
//                                 gamete_LV.add_gametes(settings->marker_names[m], kv.first, kv.second, ped.items[kv.first - 1].is_dh());
//                             }
//                         }
//                     }
//                 }
//             }

            pop_data_type pop_data;
            pop_data.ancestor_names = settings->pedigree.ancestor_names;

            MSG_DEBUG("with_LC? " << std::boolalpha << settings->pedigree.with_LC);
            
461
            for (size_t m = 0; m < n_mark; ++m) {
Damien Leroux's avatar
WIP.    
Damien Leroux committed
462
                std::map<size_t, std::vector<double>> mark_states;
463
                ifile ifs(settings->job_filename("compute-LV", m));
Damien Leroux's avatar
WIP.    
Damien Leroux committed
464
465
466
467
468
469
470
471
472
473
474
                msg_handler_t::cout() << "\rReading marker data #" << m << '/' << n_mark << "...";
                rw_base() (ifs, mark_states);
                const std::string& mark = settings->marker_names[m];
                if (settings->pedigree.with_LC) {
                    msg_handler_t::cout() << " LC...";
                    for (auto& ind_sp: mark_states) {
                        const std::string& gen_name = ped.generations[ped.node_generations[ind_sp.first]]->name;
                        Eigen::Map<Eigen::VectorXd> sp(ind_sp.second.data(), ind_sp.second.size());
                        pop_data.LV.data_by_marker[gen_name][mark].emplace_back(sp);
                    }
                }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
475
                if (have_gametes) {
Damien Leroux's avatar
WIP.    
Damien Leroux committed
476
                    msg_handler_t::cout() << " gametes...";
damien's avatar
damien committed
477
                    ifile ig(settings->job_filename("gametes-LV", m));
Damien Leroux's avatar
WIP.    
Damien Leroux committed
478
479
480
                    std::map<int, MatrixXd> gp;
                    rw_base() (ig, gp);
                    for (const auto& kv: gp) {
damien's avatar
damien committed
481
482
                        if (kv.second.size() == 0) {
                            MSG_DEBUG("HAVE EMPTY GAMETE LV " << kv.first);
Damien Leroux's avatar
WIP.    
Damien Leroux committed
483
484
                        } else if (kv.second.sum() == 0) {
                            MSG_DEBUG("HAVE ZERO GAMETE LV " << kv.first);
damien's avatar
damien committed
485
                        } else {
Damien Leroux's avatar
WIP.    
Damien Leroux committed
486
                            gamete_LV.add_gametes(settings->marker_names[m], kv.first, kv.second, ped.items[kv.first - 1].is_dh());
damien's avatar
damien committed
487
                        }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
488
489
                    }
                }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
490
                msg_handler_t::cout() << "\x1b[K";
491
            }
Damien Leroux's avatar
WIP.    
Damien Leroux committed
492
493

            if (settings->output_mode & bn_settings_t::OutputGametes) {
494
                std::string filename = SPELL_STRING(settings->work_directory << "/" << settings->pop_name << ".cache/gamete.data");
Damien Leroux's avatar
WIP.    
Damien Leroux committed
495
496
                ofile og(filename);
                rw_base() (og, gamete_LV);
damien's avatar
damien committed
497
498
                MSG_DEBUG("Gametes: " << gamete_LV.data);
//                 abort();
Damien Leroux's avatar
WIP.    
Damien Leroux committed
499
500
            }

501
502
503
504
505
506

            pop_data.name = settings->pop_name;
            for (const auto& om: settings->observed_mark) {
                pop_data.marker_observation_filenames[om.first] = om.second.filename;
            }

Damien Leroux's avatar
Damien Leroux committed
507
//             settings->load_full_pedigree_data();
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523

            /*pop_data.genetic_map_filename = settings->map_filename;*/
            pop_data.pedigree_filename = settings->pedigree_filename;
            /*pop_data.qtl_generation_name = settings->qtl_generation_name;*/
            pop_data.generations = ped.generations;
            pop_data.families.clear();
            for (size_t ind_node: ped.tree.m_ind_number_to_node_number) {
                if (ind_node == (size_t) -1) {
                    continue;
                }
                auto gen = pop_data.generations[ped.node_generations[ind_node]];
                pop_data.families[gen->name].push_back(ind_node);
                pop_data.gen_by_id[ind_node] = gen;
            }
            MSG_DEBUG("POP-DATA:");
            MSG_DEBUG("" << pop_data);
524
            std::string filename = pop_data.save(settings->work_directory + "/" + settings->pop_name + ".cache");
525
            MSG_DEBUG("Saved output in " << filename);
Damien Leroux's avatar
Damien Leroux committed
526
            {
527
                ofile output(settings->job_filename("output", 0));
Damien Leroux's avatar
Damien Leroux committed
528
529
                rw_base() (output, filename);
            }
530
531
532
            return true;
        }
    }},
533
534
};

535
#include "output_impl.h"