convert_to_mapmaker.cc 2.38 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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




void convert_to_mapmaker(const pop_data_type& popdata, const std::vector<chromosome>& gm, const std::string& gen_name, const marker_observation_spec& spec, char unknown, const std::string& data_type, const std::string& filename)
{
    std::vector<char> obs_symbols;
    const auto& labels = popdata.generations[gen_name]->unique_labels();
    MatrixXd obs_to_symbol = MatrixXd::Zero(spec.scores.size(), labels.size());

    obs_symbols.resize(spec.scores.size());

    for (size_t c = 0; c < spec.scores.size(); ++c) {
        std::vector<std::pair<char, char>> geno_vec;
        std::tie(obs_symbols[c], geno_vec) = scores[c];
        double coef = 1. / geno_vec.size();
        for (size_t g = 0; g < geno_vec.size(); ++g) {
            allele_pair l(geno_vec[g].first, geno_vec[g].second);
            for (size_t r = 0; r < obs_to_symbol.rows(); ++r) {
                if (labels[r] == l) {
                    obs_to_symbol(r, c) = coef;
                }
            }
        }
    }

    std::ofstream output(filename);

    size_t n_ind = popdata.LV.data.begin()->second.find(gen_name)->second.size();
    size_t n_mark = 0;
    for (const auto& chrom: gm) {
        n_mark += chrom.raw.marker_names.size();
    }

    output << "data type " << data_type << std::endl;
    output << n_ind << ' ' << n_mark << " 0 0" << std::endl;
    
    for (const auto& chrom: gm) {
        size_t m = 0;
        for (const auto& mark: chrom.raw.marker_names) {
            output << '*' << mark << ' ';
            for (size_t i = 0; i < n_ind; ++i) {
                const auto& lv = popdata.LV(chrom.name, gen_name, i);
                const VectorXd obs = lv.col(i);
                VectorXd obs_like = obs_to_symbol * obs;
                int max_i = 0;
                double max = 0;
                bool fail = true;
                for (int i = 0; i < obs_like.size(); ++i) {
                    if (obs_like(i) > max) {
                        fail = false;
                        max_i = i;
                        max = obs_like(i);
                    } else if (obs_like(i) == max) {
                        fail = true;
                    }
                    if (fail) {
                        output << unknown;
                    } else {
                        output << obs_symbols[max_i];
                    }
                }
            }
            output << std::endl;
            ++m;
        }
    }
}