Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017 of Haochuan Chen (excluding colvar_UIestimator.h)
3 : Copyright (c) 2017 of Haohao Fu (colvar_UIestimator.h)
4 :
5 : This program is free software: you can redistribute it and/or modify
6 : it under the terms of the GNU Lesser General Public License as published
7 : by the Free Software Foundation, either version 3 of the License, or
8 : (at your option) any later version.
9 :
10 : This program is distributed in the hope that it will be useful,
11 : but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 : GNU Lesser General Public License for more details.
14 :
15 : You should have received a copy of the GNU Lesser General Public License
16 : along with this program. If not, see <http://www.gnu.org/licenses/>.
17 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
18 : #ifdef __PLUMED_HAS_BOOST_SERIALIZATION
19 : #include "DRR.h"
20 :
21 : namespace PLMD {
22 : namespace drr {
23 :
24 : using std::vector;
25 : using std::string;
26 : using std::begin;
27 : using std::end;
28 :
29 0 : bool DRRAxis::isInBoundary(double x) const {
30 0 : if (x < min || x > max)
31 : return false;
32 : else
33 0 : return true;
34 : }
35 :
36 5860 : bool DRRAxis::isRealPeriodic() const {
37 5860 : if (periodic == true) {
38 5808 : if (std::abs(domainMax - max) < binWidth &&
39 5796 : std::abs(domainMin - min) < binWidth) {
40 : return true;
41 : } else {
42 12 : return false;
43 : }
44 : } else {
45 : return false;
46 : }
47 : }
48 :
49 4 : DRRAxis DRRAxis::merge(const DRRAxis &d1, const DRRAxis &d2) {
50 4 : const double newmin = std::min(d1.min, d2.min);
51 4 : const double newmax = std::max(d1.max, d2.max);
52 4 : const double newWidth = d1.binWidth;
53 4 : const size_t newbins = size_t(std::nearbyint((newmax - newmin) / newWidth));
54 4 : const bool newpbc = d1.periodic;
55 4 : const double newdmin = std::min(d1.domainMin, d2.domainMin);
56 4 : const double newdmax = std::max(d1.domainMax, d2.domainMax);
57 : DRRAxis result(newmin, newmax, newbins, newpbc, newdmin, newdmax);
58 4 : return result;
59 : }
60 :
61 60 : vector<double> DRRAxis::getMiddlePoints() {
62 60 : vector<double> result(nbins, 0);
63 60 : const double width = binWidth;
64 60 : double temp = min - width / 2;
65 60 : std::generate(begin(result), end(result), [&temp, &width]() {
66 3382 : temp += width;
67 : return temp;
68 : });
69 60 : return result;
70 : }
71 :
72 2407170 : size_t DRRForceGrid::index1D(const DRRAxis &c, double x) {
73 2407170 : size_t idx = size_t(std::floor((x - c.min) / c.binWidth));
74 2407170 : idx = (idx == c.nbins) ? (c.nbins - 1) : idx;
75 2407170 : return idx;
76 : }
77 :
78 34 : void DRRForceGrid::fillTable(const vector<vector<double>> &in) {
79 34 : table.resize(ndims, vector<double>(sampleSize, 0));
80 82 : for (size_t i = 0; i < ndims; ++i) {
81 : size_t repeatAll = 1, repeatOne = 1;
82 62 : for (size_t j = i + 1; j < ndims; ++j)
83 14 : repeatOne *= in[j].size();
84 62 : for (size_t j = 0; j < i; ++j)
85 14 : repeatAll *= in[j].size();
86 : size_t in_i_sz = in[i].size();
87 3390 : for (size_t l = 0; l < in_i_sz; ++l)
88 3342 : std::fill_n(begin(table[i]) + l * repeatOne, repeatOne, in[i][l]);
89 784 : for (size_t k = 0; k < repeatAll - 1; ++k)
90 736 : std::copy_n(begin(table[i]), repeatOne * in_i_sz,
91 736 : begin(table[i]) + repeatOne * in_i_sz * (k + 1));
92 : }
93 34 : }
94 :
95 38 : DRRForceGrid::DRRForceGrid()
96 38 : : suffix(""), ndims(0), dimensions(0), sampleSize(0),
97 38 : headers(""), table(0), forces(0), samples(0), endpoints(0), shifts(0),
98 76 : outputunit(1.0) {}
99 :
100 26 : DRRForceGrid::DRRForceGrid(const vector<DRRAxis> &p_dimensions,
101 26 : const string &p_suffix, bool initializeTable)
102 26 : : suffix(p_suffix), ndims(p_dimensions.size()), dimensions(p_dimensions) {
103 26 : sampleSize = 1;
104 26 : vector<vector<double>> mp(ndims);
105 26 : std::stringstream ss;
106 26 : ss << "# " << ndims << '\n';
107 26 : shifts.resize(ndims, 0);
108 26 : shifts[0] = 1;
109 64 : for (size_t i = 0; i < ndims; ++i) {
110 38 : sampleSize = dimensions[i].nbins * sampleSize;
111 38 : mp[i] = dimensions[i].getMiddlePoints();
112 38 : if (i > 0) {
113 12 : shifts[i] = shifts[i - 1] * dimensions[i - 1].nbins;
114 : }
115 : ss.precision(std::numeric_limits<double>::max_digits10);
116 76 : ss << std::fixed << "# " << dimensions[i].min << ' '
117 76 : << dimensions[i].binWidth << ' ' << dimensions[i].nbins;
118 38 : if (dimensions[i].isPeriodic())
119 56 : ss << " 1" << '\n';
120 : else
121 20 : ss << " 0" << '\n';
122 : }
123 26 : headers = ss.str();
124 26 : if (initializeTable)
125 18 : fillTable(mp);
126 26 : forces.resize(sampleSize * ndims, 0.0);
127 26 : samples.resize(sampleSize, 0);
128 26 : outputunit = 1.0;
129 : // For 1D pmf
130 26 : if (ndims == 1) {
131 14 : endpoints.resize(dimensions[0].nbins + 1, 0);
132 14 : double ep = dimensions[0].min;
133 14 : double stride = dimensions[0].binWidth;
134 882 : for (auto &i : endpoints) {
135 868 : i = ep;
136 868 : ep += stride;
137 : }
138 : }
139 26 : }
140 :
141 393846 : bool DRRForceGrid::isInBoundary(const vector<double> &pos) const {
142 : bool result = true;
143 1175042 : for (size_t i = 0; i < ndims; ++i) {
144 782805 : if (pos[i] < dimensions[i].min || pos[i] > dimensions[i].max)
145 : return false;
146 : }
147 : return result;
148 : }
149 :
150 1012383 : size_t DRRForceGrid::sampleAddress(const vector<double> &pos) const {
151 : size_t saddr = 0;
152 3029541 : for (size_t i = 0; i < ndims; ++i) {
153 2017158 : saddr += shifts[i] * index1D(dimensions[i], pos[i]);
154 : }
155 1012383 : return saddr;
156 : }
157 :
158 1723 : bool DRRForceGrid::store(const vector<double> &pos, const vector<double> &f,
159 : unsigned long int nsamples) {
160 1723 : if (isInBoundary(pos)) {
161 1714 : if (nsamples == 0)
162 : return true;
163 914 : const size_t baseaddr = sampleAddress(pos) * ndims;
164 914 : samples[baseaddr / ndims] += nsamples;
165 : auto it_fa = begin(forces) + baseaddr;
166 914 : std::transform(begin(f), end(f), it_fa, it_fa, std::plus<double>());
167 914 : return true;
168 : } else {
169 : return false;
170 : }
171 : }
172 :
173 4 : vector<DRRAxis> DRRForceGrid::merge(const vector<DRRAxis> &dA,
174 : const vector<DRRAxis> &dB) {
175 4 : vector<DRRAxis> dR(dA.size());
176 4 : std::transform(begin(dA), end(dA), begin(dB), begin(dR), DRRAxis::merge);
177 4 : return dR;
178 : }
179 :
180 : vector<double>
181 1600 : DRRForceGrid::getAccumulatedForces(const vector<double> &pos) const {
182 1600 : vector<double> result(ndims, 0);
183 1600 : if (!isInBoundary(pos))
184 : return result;
185 800 : const size_t baseaddr = sampleAddress(pos) * ndims;
186 800 : std::copy(begin(forces) + baseaddr, begin(forces) + baseaddr + ndims,
187 : begin(result));
188 800 : return result;
189 : }
190 :
191 100618 : unsigned long int DRRForceGrid::getCount(const vector<double> &pos,
192 : bool SkipCheck) const {
193 100618 : if (!SkipCheck) {
194 1600 : if (!isInBoundary(pos)) {
195 : return 0;
196 : }
197 : }
198 99818 : return samples[sampleAddress(pos)];
199 : }
200 :
201 195576 : vector<double> DRRForceGrid::getGradient(const vector<double> &pos,
202 : bool SkipCheck) const {
203 195576 : vector<double> result(ndims, 0);
204 195576 : if (!SkipCheck) {
205 162000 : if (!isInBoundary(pos)) {
206 : return result;
207 : }
208 : }
209 195576 : const size_t baseaddr = sampleAddress(pos);
210 : const unsigned long int &count = samples[baseaddr];
211 195576 : if (count == 0)
212 : return result;
213 195411 : auto it_fa = begin(forces) + baseaddr * ndims;
214 195411 : std::transform(it_fa, it_fa + ndims, begin(result),
215 389802 : [&count](double fa) { return (-1.0) * fa / count; });
216 195411 : return result;
217 : }
218 :
219 64800 : double DRRForceGrid::getDivergence(const vector<double> &pos) const {
220 : double div = 0.0;
221 64800 : vector<double> grad_deriv(ndims, 0.0);
222 64800 : if (!isInBoundary(pos)) {
223 : return div;
224 : }
225 64800 : const size_t force_addr = sampleAddress(pos) * ndims;
226 64800 : vector<double> grad = getGradient(pos);
227 194400 : for (size_t i = 0; i < ndims; ++i) {
228 129600 : const double binWidth = dimensions[i].binWidth;
229 129600 : vector<double> first = pos;
230 129600 : first[i] = dimensions[i].min + binWidth * 0.5;
231 129600 : vector<double> last = pos;
232 129600 : last[i] = dimensions[i].max - binWidth * 0.5;
233 129600 : const size_t force_addr_first = sampleAddress(first) * ndims;
234 129600 : const size_t force_addr_last = sampleAddress(last) * ndims;
235 129600 : if (force_addr == force_addr_first) {
236 720 : if (dimensions[i].isRealPeriodic() == true) {
237 720 : vector<double> next = pos;
238 720 : next[i] += binWidth;
239 720 : const vector<double> grad_next = getGradient(next);
240 720 : const vector<double> grad_prev = getGradient(last);
241 720 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
242 : } else {
243 0 : vector<double> next = pos;
244 0 : next[i] += binWidth;
245 0 : vector<double> next_2 = next;
246 0 : next_2[i] += binWidth;
247 0 : const vector<double> grad_next = getGradient(next);
248 0 : const vector<double> grad_next_2 = getGradient(next_2);
249 0 : grad_deriv[i] =
250 0 : (grad_next_2[i] * -1.0 + grad_next[i] * 4.0 - grad[i] * 3.0) /
251 0 : (2.0 * binWidth);
252 : }
253 128880 : } else if (force_addr == force_addr_last) {
254 720 : if (dimensions[i].isRealPeriodic() == true) {
255 720 : vector<double> prev = pos;
256 720 : prev[i] -= binWidth;
257 720 : const vector<double> grad_next = getGradient(first);
258 720 : const vector<double> grad_prev = getGradient(prev);
259 720 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
260 : } else {
261 0 : vector<double> prev = pos;
262 0 : prev[i] -= binWidth;
263 0 : vector<double> prev_2 = prev;
264 0 : prev_2[i] -= binWidth;
265 0 : const vector<double> grad_prev = getGradient(prev);
266 0 : const vector<double> grad_prev_2 = getGradient(prev_2);
267 0 : grad_deriv[i] =
268 0 : (grad[i] * 3.0 - grad_prev[i] * 4.0 + grad_prev_2[i] * 1.0) /
269 0 : (2.0 * binWidth);
270 : }
271 : } else {
272 128160 : vector<double> prev = pos;
273 128160 : prev[i] -= binWidth;
274 128160 : vector<double> next = pos;
275 128160 : next[i] += binWidth;
276 128160 : const vector<double> grad_next = getGradient(next);
277 128160 : const vector<double> grad_prev = getGradient(prev);
278 128160 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
279 : }
280 : }
281 : div = std::accumulate(begin(grad_deriv), end(grad_deriv), 0.0);
282 : return div;
283 : }
284 :
285 : vector<double>
286 195576 : DRRForceGrid::getCountsLogDerivative(const vector<double> &pos) const {
287 195576 : const size_t addr = sampleAddress(pos);
288 195576 : const unsigned long int count_this = samples[addr];
289 195576 : vector<double> result(ndims, 0);
290 585588 : for (size_t i = 0; i < ndims; ++i) {
291 390012 : const double binWidth = dimensions[i].binWidth;
292 : const size_t addr_first =
293 390012 : addr - shifts[i] * index1D(dimensions[i], pos[i]) + 0;
294 390012 : const size_t addr_last = addr_first + shifts[i] * (dimensions[i].nbins - 1);
295 390012 : if (addr == addr_first) {
296 2210 : if (dimensions[i].isRealPeriodic() == true) {
297 2178 : const unsigned long int &count_next = samples[addr + shifts[i]];
298 : const unsigned long int &count_prev = samples[addr_last];
299 2178 : if (count_next != 0 && count_prev != 0)
300 2160 : result[i] =
301 2160 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
302 : } else {
303 32 : const unsigned long int &count_next = samples[addr + shifts[i]];
304 32 : const unsigned long int &count_next2 = samples[addr + shifts[i] * 2];
305 32 : if (count_next != 0 && count_this != 0 && count_next2 != 0)
306 6 : result[i] =
307 6 : (std::log(count_next2) * (-1.0) + std::log(count_next) * 4.0 -
308 6 : std::log(count_this) * 3.0) /
309 6 : (2.0 * binWidth);
310 : }
311 387802 : } else if (addr == addr_last) {
312 2210 : if (dimensions[i].isRealPeriodic() == true) {
313 2178 : const unsigned long int &count_prev = samples[addr - shifts[i]];
314 : const unsigned long int &count_next = samples[addr_first];
315 2178 : if (count_next != 0 && count_prev != 0)
316 2160 : result[i] =
317 2160 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
318 : } else {
319 32 : const unsigned long int &count_prev = samples[addr - shifts[i]];
320 32 : const unsigned long int &count_prev2 = samples[addr - shifts[i] * 2];
321 32 : if (count_prev != 0 && count_this != 0 && count_prev2 != 0)
322 6 : result[i] = (std::log(count_this) * 3.0 - std::log(count_prev) * 4.0 +
323 6 : std::log(count_prev2)) /
324 6 : (2.0 * binWidth);
325 : }
326 : } else {
327 385592 : const unsigned long int &count_prev = samples[addr - shifts[i]];
328 385592 : const unsigned long int &count_next = samples[addr + shifts[i]];
329 385592 : if (count_next != 0 && count_prev != 0)
330 385432 : result[i] =
331 385432 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
332 : }
333 : }
334 195576 : return result;
335 : }
336 :
337 26 : void DRRForceGrid::write1DPMF(string filename) const {
338 52 : filename += suffix + ".pmf";
339 : FILE *ppmf;
340 26 : ppmf = fopen(filename.c_str(), "w");
341 26 : const double w = dimensions[0].binWidth;
342 : double pmf = 0;
343 26 : fprintf(ppmf, "%.9f %.9f\n", endpoints[0], pmf);
344 1166 : for (size_t i = 0; i < dimensions[0].nbins; ++i) {
345 1140 : vector<double> pos(1, 0);
346 1140 : pos[0] = table[0][i];
347 1140 : const vector<double> f = getGradient(pos, true);
348 1140 : pmf += f[0] * w / outputunit;
349 1140 : fprintf(ppmf, "%.9f %.9f\n", endpoints[i + 1], pmf);
350 : }
351 26 : fclose(ppmf);
352 26 : }
353 :
354 36 : void DRRForceGrid::writeAll(const string &filename, bool addition) const {
355 36 : string countname = filename + suffix + ".count";
356 36 : string gradname = filename + suffix + ".grad";
357 36 : vector<double> pos(ndims, 0);
358 : FILE *pGrad, *pCount;
359 36 : if (addition) {
360 8 : pGrad = fopen(gradname.c_str(), "a");
361 8 : pCount = fopen(countname.c_str(), "a");
362 : } else {
363 28 : pGrad = fopen(gradname.c_str(), "w");
364 28 : pCount = fopen(countname.c_str(), "w");
365 : }
366 :
367 : char *buffer1, *buffer2;
368 36 : buffer1 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
369 36 : buffer2 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
370 36 : setvbuf(pGrad, buffer1, _IOFBF, (sizeof(double)) * sampleSize * ndims);
371 36 : setvbuf(pCount, buffer2, _IOFBF, (sizeof(double)) * sampleSize * ndims);
372 36 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
373 36 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
374 66048 : for (size_t i = 0; i < sampleSize; ++i) {
375 196896 : for (size_t j = 0; j < ndims; ++j) {
376 130884 : pos[j] = table[j][i];
377 130884 : fprintf(pGrad, " %.9f", table[j][i]);
378 130884 : fprintf(pCount, " %.9f", table[j][i]);
379 : }
380 66012 : fprintf(pCount, " %lu\n", getCount(pos, true));
381 66012 : vector<double> f = getGradient(pos, true);
382 196896 : for (size_t j = 0; j < ndims; ++j) {
383 130884 : fprintf(pGrad, " %.9f", (f[j] / outputunit));
384 : }
385 : fprintf(pGrad, "\n");
386 : }
387 36 : fclose(pGrad);
388 36 : fclose(pCount);
389 36 : free(buffer1);
390 36 : free(buffer2);
391 36 : if (ndims == 1) {
392 52 : write1DPMF(filename);
393 : }
394 36 : }
395 :
396 2 : void DRRForceGrid::writeDivergence(const string &filename) const {
397 2 : string divname = filename + suffix + ".div";
398 2 : vector<double> pos(ndims, 0);
399 : FILE *pDiv;
400 2 : pDiv = fopen(divname.c_str(), "w");
401 2 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pDiv);
402 64802 : for (size_t i = 0; i < sampleSize; ++i) {
403 194400 : for (size_t j = 0; j < ndims; ++j) {
404 129600 : pos[j] = table[j][i];
405 129600 : fprintf(pDiv, " %.9f", table[j][i]);
406 : }
407 64800 : const double divergence = getDivergence(pos);
408 64800 : fprintf(pDiv, " %.9f", (divergence / outputunit));
409 : fprintf(pDiv, "\n");
410 : }
411 2 : fclose(pDiv);
412 2 : }
413 :
414 123 : bool ABF::store_getbias(const vector<double> &pos, const vector<double> &f,
415 : vector<double> &fbias) {
416 123 : if (!isInBoundary(pos)) {
417 : std::fill(begin(fbias), end(fbias), 0.0);
418 : return false;
419 : }
420 123 : const size_t baseaddr = sampleAddress(pos);
421 : unsigned long int &count = samples[baseaddr];
422 123 : ++count;
423 123 : double factor = 2 * (static_cast<double>(count)) / mFullSamples - 1;
424 123 : auto it_fa = begin(forces) + baseaddr * ndims;
425 : auto it_fb = begin(fbias);
426 : auto it_f = begin(f);
427 : auto it_maxFactor = begin(mMaxFactors);
428 : do {
429 : // Clamp to [0,maxFactors]
430 207 : factor = factor < 0 ? 0 : factor > (*it_maxFactor) ? (*it_maxFactor) : factor;
431 207 : (*it_fa) += (*it_f); // Accumulate instantaneous force
432 207 : (*it_fb) = factor * (*it_fa) * (-1.0) /
433 207 : static_cast<double>(count); // Calculate bias force
434 : ++it_fa;
435 : ++it_fb;
436 : ++it_f;
437 : ++it_maxFactor;
438 207 : } while (it_f != end(f));
439 :
440 : return true;
441 : }
442 :
443 2 : ABF ABF::mergewindow(const ABF &aWA, const ABF &aWB) {
444 2 : const vector<DRRAxis> dR = merge(aWA.dimensions, aWB.dimensions);
445 2 : const string suffix = ".abf";
446 2 : ABF result(dR, suffix);
447 2 : const size_t nrows = result.sampleSize;
448 2 : const size_t ncols = result.ndims;
449 2 : vector<double> pos(ncols, 0);
450 402 : for (size_t i = 0; i < nrows; ++i) {
451 800 : for (size_t j = 0; j < ncols; ++j) {
452 400 : pos[j] = result.table[j][i];
453 : }
454 400 : const unsigned long int countA = aWA.getCount(pos);
455 400 : const unsigned long int countB = aWB.getCount(pos);
456 400 : const vector<double> aForceA = aWA.getAccumulatedForces(pos);
457 400 : const vector<double> aForceB = aWB.getAccumulatedForces(pos);
458 400 : result.store(pos, aForceA, countA);
459 400 : result.store(pos, aForceB, countB);
460 : }
461 2 : return result;
462 0 : }
463 :
464 195576 : vector<double> CZAR::getGradient(const vector<double> &pos,
465 : bool SkipCheck) const {
466 195576 : vector<double> result(ndims, 0);
467 195576 : if (!SkipCheck) {
468 162000 : if (!isInBoundary(pos)) {
469 : return result;
470 : }
471 : }
472 195576 : if (kbt <= std::numeric_limits<double>::epsilon()) {
473 : std::cerr << "ERROR! The kbt shouldn't be zero when use CZAR estimator. "
474 0 : << '\n';
475 0 : std::abort();
476 : }
477 195576 : const size_t baseaddr = sampleAddress(pos);
478 195576 : const vector<double> log_deriv(getCountsLogDerivative(pos));
479 : const unsigned long int &count = samples[baseaddr];
480 195576 : if (count == 0)
481 : return result;
482 195421 : auto it_fa = begin(forces) + baseaddr * ndims;
483 195421 : std::transform(it_fa, it_fa + ndims, begin(log_deriv), begin(result),
484 389822 : [&count, this](double fa, double ld) {
485 389822 : return fa * (-1.0) / count - kbt * ld;
486 : });
487 195421 : return result;
488 : }
489 :
490 2 : CZAR CZAR::mergewindow(const CZAR &cWA, const CZAR &cWB) {
491 2 : const vector<DRRAxis> dR = merge(cWA.dimensions, cWB.dimensions);
492 2 : const double newkbt = cWA.kbt;
493 2 : const string suffix = ".czar";
494 : CZAR result(dR, suffix, newkbt);
495 2 : const size_t nrows = result.sampleSize;
496 2 : const size_t ncols = result.ndims;
497 2 : vector<double> pos(ncols, 0);
498 402 : for (size_t i = 0; i < nrows; ++i) {
499 800 : for (size_t j = 0; j < ncols; ++j) {
500 400 : pos[j] = result.table[j][i];
501 : }
502 400 : const unsigned long int countA = cWA.getCount(pos);
503 400 : const unsigned long int countB = cWB.getCount(pos);
504 400 : const vector<double> aForceA = cWA.getAccumulatedForces(pos);
505 400 : const vector<double> aForceB = cWB.getAccumulatedForces(pos);
506 400 : result.store(pos, aForceA, countA);
507 400 : result.store(pos, aForceB, countB);
508 : }
509 2 : return result;
510 : }
511 :
512 18 : void CZAR:: writeZCount(const string &filename, bool addition) const {
513 18 : string countname = filename + ".zcount";
514 18 : vector<double> pos(ndims, 0);
515 : FILE *pCount;
516 18 : if (addition) {
517 4 : pCount = fopen(countname.c_str(), "a");
518 : } else {
519 14 : pCount = fopen(countname.c_str(), "w");
520 : }
521 : char *buffer;
522 18 : buffer = (char *)malloc((sizeof(double)) * sampleSize * ndims);
523 18 : setvbuf(pCount, buffer, _IOFBF, (sizeof(double)) * sampleSize * ndims);
524 18 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
525 33024 : for (size_t i = 0; i < sampleSize; ++i) {
526 98448 : for (size_t j = 0; j < ndims; ++j) {
527 65442 : pos[j] = table[j][i];
528 65442 : fprintf(pCount, " %.9f", table[j][i]);
529 : }
530 33006 : fprintf(pCount, " %lu\n", getCount(pos, true));
531 : }
532 18 : fclose(pCount);
533 18 : free(buffer);
534 18 : }
535 :
536 : }
537 : }
538 :
539 : #endif
|
