WriteResults.cpp 4.93 KB
Newer Older
Iakov Davydov's avatar
Iakov Davydov committed
1
2
3
4
5
6
7
8

#include <iostream>
#include <iomanip>
#include <fstream>
#include <limits>
#include "WriteResults.h"
#include <string>

Iakov Davydov's avatar
Iakov Davydov committed
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
void WriteResults::outputResults(void) {
  // If no file set, then do nothing
  if (!mFilename)
    return;

  // Range of branches to be output (for H0 and H1)
  std::map<size_t, double>::const_iterator im;
  std::map<size_t, std::string>::const_iterator ims;

  size_t min_branch = std::numeric_limits<size_t>::max();
  size_t max_branch = 0;
  for (im = mLnL[0].begin(); im != mLnL[0].end(); ++im) {
    size_t v = im->first;
    if (v < min_branch)
      min_branch = v;
    if (v > max_branch)
      max_branch = v;
  }
  for (im = mLnL[1].begin(); im != mLnL[1].end(); ++im) {
    size_t v = im->first;
    if (v < min_branch)
      min_branch = v;
    if (v > max_branch)
      max_branch = v;
  }

  // No entries, so do nothing
  if (min_branch == std::numeric_limits<size_t>::max())
    return;

  // Open the output file
  std::ofstream out(mFilename, std::ios_base::trunc | std::ios_base::out);
  if (!out.good()) {
    std::cout << "Cannot create results file <" << mFilename << ">"
              << std::endl;
    return;
  }

  // Write the log-likelihood values (if a value is not present, write NA)

  for (size_t branch = min_branch; branch <= max_branch; ++branch) {
    out << "Branch: " << std::setw(4) << branch << std::endl
        << std::endl
        << "  LnL0: ";

    // Prints LnL for H0 if present
    im = mLnL[0].find(branch);
    if (im == mLnL[0].end()) {
      out << std::setw(22) << "NA" << std::endl
          << std::endl;
    } else {
      out << std::setw(22) << std::setprecision(15) << std::fixed << im->second;

      ims = mParamStr[0].find(branch);
      // This output does not work on BG/Q.
64
65
      out << std::endl
          << std::fixed << "Branch lengths:" << ims->second << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    }
    out << "  LnL1: ";

    // Prints LnL for H1 if present
    im = mLnL[1].find(branch);
    if (im == mLnL[1].end()) {
      out << std::setw(22) << "NA";
    } else {
      out << std::setw(22) << std::setprecision(15) << std::fixed << im->second;

      ims = mParamStr[1].find(branch);
      out << std::endl
          << std::fixed << "Branch lengths:" << ims->second << std::endl;
    }
    out << std::endl;
  }

  // Write the positive selection sites (adding one to the site number because
  // they start from 1 and not zero)
  for (size_t branch = min_branch; branch <= max_branch; ++branch) {
    std::map<size_t, std::pair<std::vector<unsigned int>,
87
                               std::vector<double> > >::const_iterator ipss;
Iakov Davydov's avatar
Iakov Davydov committed
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
    ipss = mPositiveSelSites.find(branch);
    if (ipss != mPositiveSelSites.end()) {
      const std::vector<unsigned int> &site = ipss->second.first;
      const std::vector<double> &prob = ipss->second.second;

      const std::vector<size_t> &idx = orderSites(site);

      size_t ns = site.size();
      for (size_t s = 0; s < ns; ++s) {
        size_t i = idx[s];
        out << "PositiveSelectionSite for branch: " << std::setw(4) << branch;
        out << "  Site: " << std::setw(6) << site[i] + 1
            << "  Prob: " << std::setw(9) << std::setprecision(6) << std::fixed
            << prob[i] << std::endl;
      }
    }
  }

  out.close();
Iakov Davydov's avatar
Iakov Davydov committed
107
108
}

Iakov Davydov's avatar
Iakov Davydov committed
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
const std::vector<size_t> &
WriteResults::orderSites(const std::vector<unsigned int> &aSites) const {
  // Fill the map with pairs(site, its index)
  std::multimap<unsigned int, size_t> ordered_map;
  size_t end(aSites.size());
  for (size_t i = 0; i < end; ++i)
    ordered_map.insert(std::pair<unsigned int, size_t>(aSites[i], i));

  // Take the list of indices after ordering the list of sites
  mSiteOrder.clear();
  std::multimap<unsigned int, size_t>::const_iterator im(ordered_map.begin());
  std::multimap<unsigned int, size_t>::const_iterator endm(ordered_map.end());
  for (; im != endm; ++im)
    mSiteOrder.push_back(im->second);

  return mSiteOrder;
Iakov Davydov's avatar
Iakov Davydov committed
125
126
}

Iakov Davydov's avatar
Iakov Davydov committed
127
128
129
130
131
void WriteResults::saveLnL(size_t aFgBranch, double aLnL,
                           unsigned int aHypothesis) {
  // If no file set, then do nothing
  if (!mFilename)
    return;
Iakov Davydov's avatar
Iakov Davydov committed
132

Iakov Davydov's avatar
Iakov Davydov committed
133
134
135
  // Sanity check
  if (aHypothesis > 1)
    return;
Iakov Davydov's avatar
Iakov Davydov committed
136

Iakov Davydov's avatar
Iakov Davydov committed
137
138
  // Save the likelihood for later printing
  mLnL[aHypothesis][aFgBranch] = aLnL;
Iakov Davydov's avatar
Iakov Davydov committed
139
140
}

Iakov Davydov's avatar
Iakov Davydov committed
141
142
143
144
145
146
147
148
149
150
151
void WriteResults::savePositiveSelSites(
    size_t aFgBranch, const std::vector<unsigned int> &aPositiveSelSites,
    const std::vector<double> &aPositiveSelSitesProb) {
  // If no file set, then do nothing
  if (!mFilename)
    return;

  // Save the positive selection sites and corresponding probabilities for later
  // printing
  mPositiveSelSites[aFgBranch] =
      std::make_pair(aPositiveSelSites, aPositiveSelSitesProb);
Iakov Davydov's avatar
Iakov Davydov committed
152
153
}

Iakov Davydov's avatar
Iakov Davydov committed
154
155
156
157
158
void WriteResults::saveParameters(size_t aFgBranch, std::string &aParamStr,
                                  unsigned int aHypothesis) {
  // If no file set, then do nothing
  if (!mFilename)
    return;
Iakov Davydov's avatar
Iakov Davydov committed
159

Iakov Davydov's avatar
Iakov Davydov committed
160
161
162
  // Sanity check
  if (aHypothesis > 1)
    return;
Iakov Davydov's avatar
Iakov Davydov committed
163

Iakov Davydov's avatar
Iakov Davydov committed
164
165
  // Save the likelihood for later printing
  mParamStr[aHypothesis][aFgBranch] = aParamStr;
Iakov Davydov's avatar
Iakov Davydov committed
166
}