settings.h 22 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
#ifndef _SPEL_SETTINGS_H_
#define _SPEL_SETTINGS_H_
20

Damien Leroux's avatar
Damien Leroux committed
21
22
23
24
#include "input/read_mark.h"
#include "input/read_map.h"
#include "input/read_trait.h"
#include "input/design.h"
25
#include <exception>
26
27
/*#include <Eigen/Core>*/
#include "eigen.h"
28
#include "proba_config.h"
29
30
31
32
#include "data/genoprob_computer2.h"
#include "data/qtl_chromosome.h"
#include "input/read_mark.h"
#include "input/read_trait.h"
33

34
35
36
#include "bayes/output.h"

#include "input/observations.h"
Damien Leroux's avatar
Damien Leroux committed
37
#include <thread>
38

39
//#include "../3rd-party/ThreadPool/ThreadPool.h"
40

41
42
#include "basic_file_checks.h"

43
44
45
#include "linear_combination.h"


46
47
#include "task_pool.h"

Damien Leroux's avatar
Damien Leroux committed
48
49
50
51
52
53
54
enum ComputationType { NoTest=0, FTest=1, FTestLOD=2, R2=4, Chi2=8, Chi2LOD=16, Mahalanobis=32 };
enum ComputationResults { Test=0, RSS=1, Residuals=2, Coefficients=4, Rank=8 };

constexpr ComputationType operator | (ComputationType a, ComputationType b) { return  (ComputationType) (((int)a) | ((int)b)); }
constexpr ComputationResults operator | (ComputationResults a, ComputationResults b) { return  (ComputationResults) (((int)a) | ((int)b)); }


55

56
57
58
59
60
61
struct pleiotropy_descr {
    std::string name;
    double tolerance;
    std::vector<std::string> traits;
};

62
63
64
65
66

struct marker_observation_spec {
    ObservationDomain domain;
    size_t domain_size;
    std::vector<std::pair<char, std::vector<std::pair<char, char>>>> scores;
67
68
69
70
71
72
73
    /*std::map<char, VectorXd> compile() const;*/
    const std::vector<std::pair<char, char>>& score(char obs) const
    {
        static std::vector<std::pair<char, char>> empty;
        auto it = std::find_if(scores.begin(), scores.end(), [=] (const std::pair<char, std::vector<std::pair<char, char>>>& e) { return e.first == obs; });
        return it != scores.end() ? it->second : empty;
    }
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
};


struct population_marker_observation {
    std::string filename;
    std::string format_name;
    size_t process_size;
    marker_data observations;

    std::map<char, VectorXd> observation_vectors;

    MatrixXd raw_observations(const std::vector<char>& observations) const
    {
        MatrixXd ret = MatrixXd::Ones(process_size, observations.size());
        std::string debug_obs(observations.begin(), observations.end());
        /*MSG_DEBUG("raw_observations[start] " << name << " obs=" << debug_obs);*/
        /*MSG_DEBUG(main_process().column_labels);*/
        /*for (auto& kv: obs_vectors) {*/
        /*MSG_DEBUG(" * " << kv.first << "   " << kv.second.transpose());*/
        /*}*/
        for (size_t i = 0; i < observations.size(); ++i) {
            std::map<char, VectorXd>::const_iterator ovi = observation_vectors.find(observations[i]);
            if (ovi != observation_vectors.end()) {
                ret.col(i) = ovi->second;
            } else {
                /* error! */
                MSG_ERROR("OBSERVATION NOT FOUND! obs #" << i << " chr(" << ((int)observations[i]) << ')', "FIX IT");
            }
        }
        /*MSG_DEBUG("raw_observations " << name << " obs=" << debug_obs << std::endl << ret);*/
        /*MSG_DEBUG("raw_observations obs=" << debug_obs << std::endl << ret);*/
        return ret;
    }

    operator marker_data& () { return observations; }
    operator std::string () { return filename; }
};


struct LD_matrices {
    std::vector<double> loci;
    std::vector<MatrixXd> ld;
};


119
120
121
struct bad_settings_exception : public std::exception {
    std::string m;
    bad_settings_exception(const std::string& msgs)
122
123
124
125
126
127
128
129
130
        : m(MESSAGE("One or more errors were found." << std::endl << msgs))
    {}
    const char* what() const throw() { return m.c_str(); }
};


struct trait_metadata_type {
    std::vector<std::string> dim_names;
    MatrixXd P;
131
132
133
134
};



135
struct settings_t {
136
    std::string name;
137
    std::string notes;
138

139
    std::string map_filename;
140
    std::vector<chromosome> map;
141
    /*format_specification_t* marker_observation_specs;*/
142
143
144
    /*std::multimap<std::string, qtl_pop_type> populations;*/
    std::vector<std::shared_ptr<qtl_pop_type>> populations;
    std::map<std::string, std::vector<std::shared_ptr<qtl_pop_type>>> pops_by_trait;
145

146
147
    std::vector<qtl_chromosome> working_set;

148
    double step;
149
150
    int parallel;
    std::string work_directory;
151
    bool connected;
152
153
154
155
156
157
158
    bool epistasis;
    bool pleiotropy;
    double pleiotropy_tolerance;
    std::vector<std::string> chromosome_selection;
    std::string cofactor_marker_selection_filename;
    double cofactor_marker_selection_distance;

159
160
    std::string epistasis_qtl_selection_filename;

161
162
163
164
165
166
167
168
169
170
    std::map<std::string, double> qtl_thresholds;
    int n_permutations;
    double qtl_threshold_quantile;
    std::map<std::string, double> cofactor_thresholds;

    std::string skeleton_mode;
    std::vector<std::string> skeleton_markers;
    double skeleton_interval;

    std::string cofactor_algorithm;
Damien Leroux's avatar
Damien Leroux committed
171
    double cofactor_exclusion_window_size;
172
173
174

    std::string qtl_algorithm;

175
176
    std::string initial_selection;

177
#if 0
178
179
180
181
182
183
184
185
186
187
188
189
    enum class detection_method_t : char {
        Undef,
        iQTLm,
        CIM
    };

    bool P_required;
    bool T_required;
    bool E_required;
    bool D_required;
    detection_method_t detection_method;
    double detection_window;
190
#endif
191

192
//	ThreadPool* pool;
193
194
	static std::thread::id main_thread;

195
196
197
	double tolerance;
    std::map<const chromosome*, std::vector<double>> estimation_loci;

198
199
200
201
    /* FIXME LD data is global to the pedigree (to the ancestors), not local to a population */
    std::map<std::string, std::vector<std::string>> ld_parents;
    std::map<std::string, std::map<const chromosome*, LD_matrices>> ld_data;
    std::map<std::string, std::map<const chromosome*, std::map<double, int>>> parent_count_per_pop_per_chr;
202

203
204
    LV_database LV;

205
    std::map<std::string, std::vector<std::vector<std::shared_ptr<qtl_pop_type>>>> linked_pops;
Damien Leroux's avatar
Damien Leroux committed
206
207
208

    std::unordered_map<std::thread::id, std::vector<std::string>> thread_stacks;

209
210
211
    size_t max_order;
    bool cross_indicator_can_interact;

212
213
214
215
    std::vector<std::string> npoint_gen;

    bool output_npoint;

216
217
    std::unique_ptr<pop_data_type> pop_data;
    std::vector<std::pair<std::string, std::string>> datasets;
218
219
220
221
222
223
224
225
226
    std::vector<std::pair<std::string, std::string>> covariables;
    std::vector<pleiotropy_descr> pleiotropic_traits_descr;

    std::vector<std::string> with_covariables;
    std::vector<std::string> with_traits;
    std::vector<std::string> with_lg;

    std::vector<trait> covar_per_pop;

Damien Leroux's avatar
Damien Leroux committed
227
    double lod_support_inner, lod_support_outer;
228
229

    std::map<std::string, trait_metadata_type> trait_metadata;
230

Damien Leroux's avatar
Damien Leroux committed
231
232
    double Rjoin;
    double Rskip;
Damien Leroux's avatar
Damien Leroux committed
233
234
235
    
    std::string report_format_string;
    ComputationResults additional_table_output;
Damien Leroux's avatar
Damien Leroux committed
236

237
    settings_t()
238
239
        : notes()
        , map_filename(), map()
240
        /*, populations()*/
241
        , pops_by_trait()
242
        , working_set()
243
244
        , step(1)
        , parallel(0)
245
        , work_directory(".")
246
247
248
249
250
251
252
        , connected(false)
        , epistasis(false)
        , pleiotropy(false)
        , pleiotropy_tolerance(.001)
        , chromosome_selection()
        , cofactor_marker_selection_filename()
        , cofactor_marker_selection_distance(10)
253
        , epistasis_qtl_selection_filename()
254
        , qtl_thresholds()
255
        , n_permutations(10000)
256
        , qtl_threshold_quantile(0.05)
257
258
259
260
261
        , cofactor_thresholds()
        , skeleton_mode("auto")
        , skeleton_markers()
        , skeleton_interval(20.)
        , cofactor_algorithm("forward")
Damien Leroux's avatar
Damien Leroux committed
262
        , cofactor_exclusion_window_size(30.)
263
264
265
        , qtl_algorithm("iQTLm")
        /*, detection_method(detection_method_t::Undef)*/
        /*, detection_window(10.)*/
266
//		, pool(0)
267
268
		, tolerance(1.e-10)
        , estimation_loci()
269
270
        , ld_parents()
        , ld_data()
271
        , parent_count_per_pop_per_chr()
272
        , LV()
Damien Leroux's avatar
Damien Leroux committed
273
274
        , linked_pops()
        , thread_stacks()
275
276
        , max_order(1)
        , cross_indicator_can_interact(false)
277
278
        , npoint_gen()
        , output_npoint(false)
279
        , pop_data()
280
281
282
283
284
285
286
        , datasets()
        , covariables()
        , pleiotropic_traits_descr()
        , with_covariables()
        , with_traits()
        , with_lg()
        , covar_per_pop()
Damien Leroux's avatar
Damien Leroux committed
287
288
        , lod_support_inner(1.)
        , lod_support_outer(lod_support_inner * 2)
289
        , Rjoin(0.)
Damien Leroux's avatar
Damien Leroux committed
290
        , Rskip(1.)
291
292
        , report_format_string("text")
//         , report_format_string("excel")
Damien Leroux's avatar
Damien Leroux committed
293
294
//         , additional_table_output((ComputationResults) 0)
        , additional_table_output(ComputationResults::RSS | ComputationResults::Rank)
295
296
    {}

Damien Leroux's avatar
Damien Leroux committed
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
    std::vector<std::pair<const chromosome*, double>>
        selected_qtl() const
        {
            std::vector<std::pair<const chromosome*, double>> ret;
            size_t n = 0;
            for (auto& qc: working_set) {
                n += qc.qtl.size();
            }
            ret.reserve(n);
            for (auto& qc: working_set) {
                for (auto& q: qc.qtl) {
                    ret.emplace_back(qc.chr, q);
                }
            }
            return ret;
        }

    std::vector<double>
        selected_qtl(const std::string& name) const
        {
            for (auto& qc: working_set) {
                if (qc.chr->name == name) {
                    return qc.qtl;
                }
            }
            return {};
        }
324
325

    ~settings_t()
326
    {
327
//        if (pool) { delete pool; }
328
        /*if (marker_observation_specs) { delete marker_observation_specs; }*/
329
    }
330

331
332
333
334
335
336
337
338
339
340
341
    const chromosome*
        find_chromosome(const std::string& name) const
        {
            for (auto& c: map) {
                if (c.name == name) {
                    return &c;
                }
            }
            return NULL;
        }

342
    static settings_t* from_args(int argc, const char** argv);
343
344

    bool sanity_check() const;
345

346
    void finalize();
347

348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
//	template <class F, class... Args>
//    auto enqueue(int _Sync, F&& f, Args&&... args)
//        -> std::future<typename std::result_of<F(Args...)>::type>
//		{
//            return TaskPool::enqueue(std::forward<F>(f),
//                                     std::forward<Args>(args)...);
//
////			if (_Sync || parallel < 2 || std::this_thread::get_id() != main_thread) {
////                /*MSG_DEBUG(std::this_thread::get_id()  << " RUN SYNC");*/
////				return std::async(std::launch::deferred,
////								  std::forward<F>(f),
////								  std::forward<Args>(args)...);
////			}
////			return pool->enqueue(std::forward<F>(f),
////								 std::forward<Args>(args)...);
//		}
364

365
    std::string set_title(const std::string& t)
366
    {
367
        MSG_INFO(GREEN << "Current task: " << t);
368
369
370
371
372
373
374
375
376
        return {};
//        std::string ret;
//        if (pool) {
//            ret = pool->get_title();
//        }
//        if (pool && msg_handler_t::color()) {
//            pool->set_title(t);
//        }
//        return ret;
377
378
379
380
381
382
383
384
385
386
    }

    void set_parallel(const std::string& s)
    {
        if (s == "auto") {
            parallel = std::thread::hardware_concurrency();
        } else {
            std::stringstream ss(s);
            ss >> parallel;
        }
387
        TaskPool::set_slot_count(parallel);
388
    }
389

390
    double get_threshold(const std::string& trait, size_t y_block_size);
391
392
};

393
394
395

struct format_specification_t {
    std::string filename;
396
    ObservationDomain domain;
397
398
    std::map<std::string, marker_observation_spec> map;

399
400
401
    void operator () (std::map<std::string, marker_observation_spec>& settings)
    {
        settings.insert(map.begin(), map.end());
402
403
404
405
    }
};


406
407
std::ostream& operator << (std::ostream&, const settings_t&);

408
namespace read_data {
409
    settings_t* read_settings(std::istream& is);
410
411
412
}


Damien Leroux's avatar
Damien Leroux committed
413
414
extern settings_t* active_settings;

415
416
417
418
419

struct context_key_struc;
typedef std::shared_ptr<context_key_struc> context_key;

struct context_key_struc {
420
    const qtl_pop_type* pop;
421
    const geno_matrix* gen;
422
423
424
425
426
427
428
429
430
431
432
433
    const chromosome* chr;
    std::vector<double> loci;
    std::unordered_map<double, size_t> locus_indices;

    void _init_locus_indices()
    {
        size_t sz = loci.size();
        for (size_t i = 0; i < sz; ++i) {
            locus_indices[loci[i]] = i;
        }
    }

434
    context_key_struc(const qtl_pop_type* p, const chromosome* c, const std::vector<double>& l)
435
        : pop(p)
436
        , gen(p->gen.get())
437
438
439
440
441
442
        , chr(c)
        , loci(l)
    {
        _init_locus_indices();
    }

443
    context_key_struc(const geno_matrix* g)
Damien Leroux's avatar
Damien Leroux committed
444
445
446
447
448
449
450
        : pop(NULL)
        , gen(g)
        , chr(NULL)
        , loci()
    {
    }

451
    context_key_struc(const qtl_pop_type* p, const geno_matrix* g, const chromosome* c)
452
453
454
455
456
457
458
459
        : pop(p)
        , gen(g)
        , chr(c)
        , loci(compute_steps(chr->condensed.marker_locus, active_settings->step))
    {
        _init_locus_indices();
    }

460
461
#if 1
    context_key_struc(const qtl_pop_type* p, const chromosome* c)
462
        : pop(p)
463
        , gen(p->gen.get())
464
        , chr(c)
465
        , loci(compute_steps(chr->condensed.marker_locus, active_settings->step))
466
467
468
    {
        _init_locus_indices();
    }
469
#endif
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487

    bool operator == (const context_key_struc& ck) const
    {
        return pop == ck.pop && chr == ck.chr;
    }

    bool operator != (const context_key_struc& ck) const { return !(*this == ck); }

    friend
        bool operator == (const context_key& k1, const context_key& k2) { return *k1 == *k2; }

    friend
        bool operator != (const context_key& k1, const context_key& k2) { return *k1 != *k2; }

    friend
        std::ostream&
        operator << (std::ostream& os, const context_key& ck)
        {
Damien Leroux's avatar
Damien Leroux committed
488
489
490
491
492
            return os << "<context pop="
                      << (ck->pop ? ck->pop->name : "nil")
                      << " gen=" << (ck->gen ? ck->gen->name : "nil")
                      << " chr=" << (ck->chr ? ck->chr->name : "nil")
                      << ", " << ck->loci.size() << " loci>";
493
494
495
496
497
498
        }
};

namespace std {
    template <>
        struct hash<context_key> {
499
            size_t operator () (const context_key& ck) const {
500
                hash<string> h;
501
502
503
                size_t hn = 0;
                size_t hc = 0;
                size_t hs = 0;
504
                size_t hv = 0;
Damien Leroux's avatar
Damien Leroux committed
505
                if (ck->pop) {
506
507
                    hn = h(ck->pop->name);
                    hs = hash<size_t>()(ck->pop->size());
508
                    hv = hash<std::vector<std::string>>()(ck->pop->covariables.dim_names);
Damien Leroux's avatar
Damien Leroux committed
509
510
                }
                if (ck->chr) {
511
                    hc = h(ck->chr->name);
Damien Leroux's avatar
Damien Leroux committed
512
                }
513
                /*return (((hs * 0xdeadbe35) ^ hn) * 0xdeadbe35) ^ hc;*/
514
                return hv ^ hs ^ hn ^ hc;
515
516
517
518
            }
        };
}

519
520
521
522

int
get_n_parents(const context_key& ck, const locus_key& lk);

523
524
525
526
527

inline
std::vector<trait>
assemble_traits(const std::vector<trait>& src, const std::vector<std::string>& with_traits, const std::vector<pleiotropy_descr>& pleio) {
    std::vector<trait> ret(src.size() + pleio.size());
528
529
530
531
532
    MSG_DEBUG("with_traits.empty() ? " << with_traits.empty());
    auto it = with_traits.empty()
              ? std::copy(src.begin(), src.end(), ret.begin())
              : std::copy_if(src.begin(), src.end(), ret.begin(),
                             [&] (const trait& t) { return std::find(with_traits.begin(), with_traits.end(), t.name) != with_traits.end(); });
533
534
535
536
    ret.resize(it - ret.begin());
    for (const auto& d: pleio) {
        ret.emplace_back(src, d.traits, d.name, d.tolerance);
    }
537
    MSG_DEBUG("assembled traits " << ret);
538
539
540
541
    return ret;
}


542
543
544
545
inline
void
settings_t::finalize()
{
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
    {
        std::string log_file = MESSAGE(active_settings->work_directory << '/' << active_settings->name << ".spell-qtl.log");
        ensure_directories_exist(active_settings->work_directory);
        msg_handler_t::set_log_file(log_file);
        std::chrono::system_clock::time_point now = std::chrono::system_clock::now();
        std::time_t now_c = std::chrono::system_clock::to_time_t(now);
        (*msg_handler_t::instance().log_file)
                << "========================================================" << std::endl
                << "Log started on " <<  std::put_time(std::localtime(&now_c), "%F %T") << std::endl
                << "========================================================" << std::endl;
    }

    std::string pop_data_filename = work_directory + "/" + name + ".cache/" + name + ".spell-marker.data";
    if (!check_file(pop_data_filename, false, false, false)) {
        MSG_ERROR("The spell-marker data file was not found at " << pop_data_filename << ".", "Please run spell-pedigree and spell-marker prior to running spell-qtl.");
    }
562
563
564
565
566
567
568
569
570
571
572

    if (with_lg.size()) {
        std::vector<chromosome> with_map;
        for (const auto& chr: map) {
            if (std::find(with_lg.begin(), with_lg.end(), chr.name) != with_lg.end()) {
                with_map.push_back(chr);
            }
        }
        map.swap(with_map);
    }

573
574
575
    pop_data.reset(pop_data_type::load_from(pop_data_filename));
    pop_data->assemble_chromosomes(map);

576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
    auto all_variables = with_traits + with_covariables;
    std::vector<std::string> in_pleio;
    for (const auto& p: pleiotropic_traits_descr) {
        in_pleio = in_pleio + p.traits;
    }
    all_variables = all_variables + in_pleio;
    auto removed_traits = with_traits % in_pleio;
    auto removed_covar = with_covariables % in_pleio;
    with_traits = with_traits - in_pleio;
    with_covariables = with_covariables - in_pleio;
    if (removed_covar.size()) {
        MSG_ERROR("Some covariables are already used in pleiotropic traits: " << removed_covar, "Specify covariable names that are different from the trait names.");
        return;
    }
    if (removed_traits.size()) {
        MSG_ERROR("Some single traits are already used in pleiotropic traits: " << removed_traits, "Specify single trait names that are different from the trait names used in pleiotropic traits.");
        return;
    }
594
    std::set<std::string> actual_single_traits;
595
596
    for (const auto& gt: datasets) {
        ifile ifs(gt.second);
597
598
        auto all_traits = read_data::read_trait(ifs, all_variables);
        auto observed_traits = assemble_traits(all_traits, with_traits, pleiotropic_traits_descr);
599
600
601
602
603
        for (const auto& t: observed_traits) {
            if (t.dim_names.size() == 1) {
                actual_single_traits.insert(t.name);
            }
        }
604
        pop_data->set_qtl_generation(gt.first);
605
606
607
608
609
610
        MSG_DEBUG("with_covariables " << with_covariables);
        MSG_DEBUG("Creating covariables " << " from traits " << with_covariables << " in list " << observed_traits);
        covar_per_pop.emplace_back(all_traits, with_covariables, "Covariables", 0);
        MSG_DEBUG("COVAR struct " << covar_per_pop.back());
        pop_data->extract_subpops(populations, gt.second, observed_traits, linked_pops, covar_per_pop.back());
    }
Damien Leroux's avatar
Damien Leroux committed
611
612
613
614
615
616
617
    for (const auto& n_p: linked_pops) {
        for (const auto& vp: n_p.second) {
            for (const auto& p: vp) {
                MSG_DEBUG("trait " << n_p.first << " pop " << p->name << ' ' << MATRIX_SIZE(p->gen->inf_mat) << " labels " << p->gen->labels);
            }
        }
    }
618
619
620
621
622
    for (const auto& pleio: pleiotropic_traits_descr) {
        auto& pops = linked_pops[pleio.name];
        MatrixXd xtx = MatrixXd::Zero(pleio.traits.size(), pleio.traits.size());
        for (const auto& vp: pops) {
            for (const auto& p: vp) {
Damien Leroux's avatar
Damien Leroux committed
623
624
625
626
627
                if (!p->observed_traits.front().raw.size()) {
                    MSG_ERROR("Invalid dataset", "");
                } else {
                    xtx += p->observed_traits.front().XtX();
                }
628
629
630
631
632
633
634
635
636
637
638
            }
        }
        auto P = trait::do_PCA(xtx, pleio.tolerance);
        for (const auto& vp: pops) {
            for (const auto& p: vp) {
                p->observed_traits.front().setP(P);
            }
        }
        trait_metadata[pleio.name].P = P;
        trait_metadata[pleio.name].dim_names = pleio.traits;
    }
639
    for (const auto& t: actual_single_traits) {
640
641
        trait_metadata[t].P = MatrixXd::Ones(1, 1);
        trait_metadata[t].dim_names = {t};
642
643
    }

644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
    for (const auto& chr: map) {
        estimation_loci[&chr] = compute_steps(chr.condensed.marker_locus, step);
    }

    for (const auto& pptr: populations) {
        auto& parent_count_per_chr = parent_count_per_pop_per_chr[pptr->qtl_generation_name];
        auto ldi = ld_data.find(pptr->qtl_generation_name);
        if (ldi != ld_data.end()) {
            for (const auto& chr: map) {
                auto& parent_count = parent_count_per_chr[&chr];
                auto& ld = (*ldi).second[&chr];
                const auto& loci = estimation_loci[&chr];
                for (size_t i = 0; i < loci.size(); ++i) {
                    if (loci[i] != ld.loci[i]) {
                        MSG_ERROR("Discrepancy between LD locus (" << ld.loci[i] << ") and estimation locus (" << loci[i] << ") on chromosome " << chr.name, "Ensure the LD data is consistent with the map and step size you provided.");
                    }
                    parent_count[ld.loci[i]] = ld.ld[i].cols();
                }
            }
        } else {
Damien Leroux's avatar
Damien Leroux committed
664
665
666
667
            if (!pptr->gen) {
                MSG_ERROR("PROUT", "");
                continue;
            }
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
            context_key ck(new context_key_struc(pptr->gen.get()));
            locus_key lk;
            size_t n_par = get_n_parents(ck, lk);
            for (const auto& chr: map) {
                auto& parent_count = parent_count_per_chr[&chr];
                const auto& loci = estimation_loci[&chr];
                for (size_t i = 0; i < loci.size(); ++i) {
                    parent_count[loci[i]] = (int) n_par;
                }
            }
        }
    }
}


683
684
#endif