Return the most likely distance between two contigs and the number * of pairs that support that estimate. */
| 100 | /** Return the most likely distance between two contigs and the number |
| 101 | * of pairs that support that estimate. */ |
| 102 | static pair<int, unsigned> |
| 103 | maximumLikelihoodEstimate(int first, int last, |
| 104 | const Histogram& samples, |
| 105 | const PMF& pmf, |
| 106 | unsigned len0, unsigned len1) |
| 107 | { |
| 108 | int filterSize = 2 * (int)(0.05 * pmf.mean()) + 3; // want an odd filter size |
| 109 | first = max(first, (int)pmf.minValue() - samples.maximum()) - filterSize/2; |
| 110 | last = min(last, (int)pmf.maxValue() - samples.minimum()) + filterSize/2 + 1; |
| 111 | |
| 112 | /* When randomly selecting fragments that span a given point, |
| 113 | * longer fragments are more likely to be selected than |
| 114 | * shorter fragments. |
| 115 | */ |
| 116 | WindowFunction window(len0, len1); |
| 117 | |
| 118 | unsigned nsamples = samples.size(); |
| 119 | double bestLikelihood = -numeric_limits<double>::max(); |
| 120 | int bestTheta = first; |
| 121 | unsigned bestn = 0; |
| 122 | vector<double> le; |
| 123 | vector<unsigned> le_n; |
| 124 | vector<int> le_theta; |
| 125 | for (int theta = first; theta <= last; theta++) { |
| 126 | // Calculate the normalizing constant of the PMF, f_theta(x). |
| 127 | double c = 0; |
| 128 | for (int i = pmf.minValue(); i <= (int)pmf.maxValue(); ++i) |
| 129 | c += pmf[i] * window(i - theta); |
| 130 | |
| 131 | double likelihood; |
| 132 | unsigned n; |
| 133 | tie(likelihood, n) = computeLikelihood(theta, samples, pmf); |
| 134 | likelihood -= nsamples * log(c); |
| 135 | le.push_back(likelihood); |
| 136 | le_n.push_back(n); |
| 137 | le_theta.push_back(theta); |
| 138 | } |
| 139 | |
| 140 | HannWindow filter(filterSize); |
| 141 | for (int i = filterSize / 2; i < (int)le.size()-(filterSize / 2); i++) { |
| 142 | double likelihood = 0; |
| 143 | for (int j = -filterSize / 2; j <= filterSize / 2; j++) { |
| 144 | assert((unsigned)(i + j) < le.size() && i + j >= 0); |
| 145 | likelihood += filter(j) * le[i + j]; |
| 146 | } |
| 147 | |
| 148 | if (le_n[i] > 0 && likelihood > bestLikelihood) { |
| 149 | bestLikelihood = likelihood; |
| 150 | bestTheta = le_theta[i]; |
| 151 | bestn = le_n[i]; |
| 152 | } |
| 153 | } |
| 154 | return make_pair(bestTheta, bestn); |
| 155 | } |
| 156 | |
| 157 | /** Return the most likely distance between two contigs and the number |
| 158 | * of pairs that support that distance estimate. |
no test coverage detected