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 5820 : bool DRRAxis::isRealPeriodic() const {
37 5820 : if (periodic == true) {
38 17388 : if (std::abs(domainMax - max) < binWidth &&
39 5796 : std::abs(domainMin - min) < binWidth) {
40 : return true;
41 : } else {
42 0 : return false;
43 : }
44 : } else {
45 : return false;
46 : }
47 : }
48 :
49 2 : 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 2 : const double newWidth = d1.binWidth;
53 2 : const size_t newbins = size_t(std::nearbyint((newmax - newmin) / newWidth));
54 2 : 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 2 : return result;
59 : }
60 :
61 48 : vector<double> DRRAxis::getMiddlePoints() {
62 48 : vector<double> result(nbins, 0);
63 48 : const double width = binWidth;
64 48 : double temp = min - width / 2;
65 : std::generate(begin(result), end(result), [&temp, &width]() {
66 2556 : temp += width;
67 : return temp;
68 96 : });
69 48 : return result;
70 : }
71 :
72 2403142 : size_t DRRForceGrid::index1D(const DRRAxis &c, double x) {
73 2403142 : size_t idx = size_t(std::floor((x - c.min) / c.binWidth));
74 2403142 : idx = (idx == c.nbins) ? (c.nbins - 1) : idx;
75 2403142 : return idx;
76 : }
77 :
78 24 : void DRRForceGrid::fillTable(const vector<vector<double>> &in) {
79 48 : table.resize(ndims, vector<double>(sampleSize, 0));
80 84 : for (size_t i = 0; i < ndims; ++i) {
81 : size_t repeatAll = 1, repeatOne = 1;
82 48 : for (size_t j = i + 1; j < ndims; ++j)
83 12 : repeatOne *= in[j].size();
84 12 : for (size_t j = 0; j < i; ++j)
85 12 : repeatAll *= in[j].size();
86 : size_t in_i_sz = in[i].size();
87 2552 : for (size_t l = 0; l < in_i_sz; ++l)
88 2516 : std::fill_n(begin(table[i]) + l * repeatOne, repeatOne, in[i][l]);
89 768 : for (size_t k = 0; k < repeatAll - 1; ++k)
90 : std::copy_n(begin(table[i]), repeatOne * in_i_sz,
91 732 : begin(table[i]) + repeatOne * in_i_sz * (k + 1));
92 : }
93 24 : }
94 :
95 30 : DRRForceGrid::DRRForceGrid()
96 : : suffix(""), ndims(0), dimensions(0), sampleSize(0),
97 : headers(""), table(0), forces(0), samples(0), endpoints(0), shifts(0),
98 30 : outputunit(1.0) {}
99 :
100 20 : DRRForceGrid::DRRForceGrid(const vector<DRRAxis> &p_dimensions,
101 : const string &p_suffix, bool initializeTable)
102 40 : : suffix(p_suffix), ndims(p_dimensions.size()), dimensions(p_dimensions) {
103 20 : sampleSize = 1;
104 20 : vector<vector<double>> mp(ndims);
105 40 : std::stringstream ss;
106 20 : ss << "# " << ndims << '\n';
107 20 : shifts.resize(ndims, 0);
108 20 : shifts[0] = 1;
109 50 : for (size_t i = 0; i < ndims; ++i) {
110 30 : sampleSize = dimensions[i].nbins * sampleSize;
111 60 : mp[i] = dimensions[i].getMiddlePoints();
112 30 : if (i > 0) {
113 30 : shifts[i] = shifts[i - 1] * dimensions[i - 1].nbins;
114 : }
115 : ss.precision(std::numeric_limits<double>::max_digits10);
116 60 : ss << std::fixed << "# " << dimensions[i].min << ' '
117 30 : << dimensions[i].binWidth << ' ' << dimensions[i].nbins;
118 30 : if (dimensions[i].isPeriodic())
119 24 : ss << " 1" << '\n';
120 : else
121 6 : ss << " 0" << '\n';
122 : }
123 40 : headers = ss.str();
124 20 : if (initializeTable)
125 12 : fillTable(mp);
126 20 : forces.resize(sampleSize * ndims, 0.0);
127 20 : samples.resize(sampleSize, 0);
128 20 : outputunit = 1.0;
129 : // For 1D pmf
130 20 : if (ndims == 1) {
131 20 : endpoints.resize(dimensions[0].nbins + 1, 0);
132 10 : double ep = dimensions[0].min;
133 10 : double stride = dimensions[0].binWidth;
134 470 : for (auto &i : endpoints) {
135 450 : i = ep;
136 450 : ep += stride;
137 : }
138 20 : }
139 20 : }
140 :
141 391412 : bool DRRForceGrid::isInBoundary(const vector<double> &pos) const {
142 : bool result = true;
143 1170968 : for (size_t i = 0; i < ndims; ++i) {
144 1560712 : if (pos[i] < dimensions[i].min || pos[i] > dimensions[i].max)
145 : return false;
146 : }
147 : return result;
148 : }
149 :
150 1008926 : size_t DRRForceGrid::sampleAddress(const vector<double> &pos) const {
151 : size_t saddr = 0;
152 3022558 : for (size_t i = 0; i < ndims; ++i) {
153 6040896 : saddr += shifts[i] * index1D(dimensions[i], pos[i]);
154 : }
155 1008926 : return saddr;
156 : }
157 :
158 906 : bool DRRForceGrid::store(const vector<double> &pos, const vector<double> &f,
159 : unsigned long int nsamples) {
160 906 : if (isInBoundary(pos)) {
161 906 : if (nsamples == 0)
162 : return true;
163 506 : const size_t baseaddr = sampleAddress(pos) * ndims;
164 1012 : samples[baseaddr / ndims] += nsamples;
165 506 : auto it_fa = begin(forces) + baseaddr;
166 506 : std::transform(begin(f), end(f), it_fa, it_fa, std::plus<double>());
167 506 : return true;
168 : } else {
169 : return false;
170 : }
171 : }
172 :
173 2 : vector<DRRAxis> DRRForceGrid::merge(const vector<DRRAxis> &dA,
174 : const vector<DRRAxis> &dB) {
175 2 : vector<DRRAxis> dR(dA.size());
176 2 : std::transform(begin(dA), end(dA), begin(dB), begin(dR), DRRAxis::merge);
177 2 : return dR;
178 : }
179 :
180 : vector<double>
181 800 : DRRForceGrid::getAccumulatedForces(const vector<double> &pos) const {
182 800 : vector<double> result(ndims, 0);
183 800 : if (!isInBoundary(pos))
184 : return result;
185 400 : const size_t baseaddr = sampleAddress(pos) * ndims;
186 : std::copy(begin(forces) + baseaddr, begin(forces) + baseaddr + ndims,
187 : begin(result));
188 400 : return result;
189 : }
190 :
191 99065 : unsigned long int DRRForceGrid::getCount(const vector<double> &pos,
192 : bool SkipCheck) const {
193 99065 : if (!SkipCheck) {
194 800 : if (!isInBoundary(pos)) {
195 : return 0;
196 : }
197 : }
198 197330 : return samples[sampleAddress(pos)];
199 : }
200 :
201 195083 : vector<double> DRRForceGrid::getGradient(const vector<double> &pos,
202 : bool SkipCheck) const {
203 195083 : vector<double> result(ndims, 0);
204 195083 : if (!SkipCheck) {
205 162000 : if (!isInBoundary(pos)) {
206 : return result;
207 : }
208 : }
209 195083 : const size_t baseaddr = sampleAddress(pos);
210 : const unsigned long int &count = samples[baseaddr];
211 195083 : if (count == 0)
212 : return result;
213 194998 : auto it_fa = begin(forces) + baseaddr * ndims;
214 : std::transform(it_fa, it_fa + ndims, begin(result),
215 779384 : [&count](double fa) { return (-1.0) * fa / count; });
216 194998 : 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 129600 : 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 1440 : 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 : 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 1440 : 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 : 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 256320 : 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 195083 : DRRForceGrid::getCountsLogDerivative(const vector<double> &pos) const {
287 195083 : const size_t addr = sampleAddress(pos);
288 195083 : const unsigned long int count_this = samples[addr];
289 195083 : vector<double> result(ndims, 0);
290 584593 : for (size_t i = 0; i < ndims; ++i) {
291 389510 : const double binWidth = dimensions[i].binWidth;
292 : const size_t addr_first =
293 779020 : addr - shifts[i] * index1D(dimensions[i], pos[i]) + 0;
294 779020 : const size_t addr_last = addr_first + shifts[i] * (dimensions[i].nbins - 1);
295 389510 : if (addr == addr_first) {
296 2190 : 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 : result[i] =
301 4320 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
302 : } else {
303 12 : const unsigned long int &count_next = samples[addr + shifts[i]];
304 12 : const unsigned long int &count_next2 = samples[addr + shifts[i] * 2];
305 12 : if (count_next != 0 && count_this != 0 && count_next2 != 0)
306 : result[i] =
307 12 : (std::log(count_next2) * (-1.0) + std::log(count_next) * 4.0 -
308 8 : std::log(count_this) * 3.0) /
309 8 : (2.0 * binWidth);
310 : }
311 387320 : } else if (addr == addr_last) {
312 2190 : 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 : result[i] =
317 4320 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
318 : } else {
319 12 : const unsigned long int &count_prev = samples[addr - shifts[i]];
320 12 : const unsigned long int &count_prev2 = samples[addr - shifts[i] * 2];
321 12 : if (count_prev != 0 && count_this != 0 && count_prev2 != 0)
322 12 : result[i] = (std::log(count_this) * 3.0 - std::log(count_prev) * 4.0 +
323 8 : std::log(count_prev2)) /
324 8 : (2.0 * binWidth);
325 : }
326 : } else {
327 385130 : const unsigned long int &count_prev = samples[addr - shifts[i]];
328 385130 : const unsigned long int &count_next = samples[addr + shifts[i]];
329 385130 : if (count_next != 0 && count_prev != 0)
330 : result[i] =
331 770072 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
332 : }
333 : }
334 195083 : return result;
335 : }
336 :
337 12 : void DRRForceGrid::write1DPMF(string filename) const {
338 24 : filename += suffix + ".pmf";
339 : FILE *ppmf;
340 12 : ppmf = fopen(filename.c_str(), "w");
341 12 : const double w = dimensions[0].binWidth;
342 : double pmf = 0;
343 12 : fprintf(ppmf, "%.9f %.9f\n", endpoints[0], pmf);
344 680 : for (size_t i = 0; i < dimensions[0].nbins; ++i) {
345 656 : vector<double> pos(1, 0);
346 656 : pos[0] = table[0][i];
347 656 : const vector<double> f = getGradient(pos, true);
348 656 : pmf += f[0] * w / outputunit;
349 1312 : fprintf(ppmf, "%.9f %.9f\n", endpoints[i + 1], pmf);
350 : }
351 12 : fclose(ppmf);
352 12 : }
353 :
354 20 : void DRRForceGrid::writeAll(const string &filename) const {
355 40 : string countname = filename + suffix + ".count";
356 40 : string gradname = filename + suffix + ".grad";
357 20 : vector<double> pos(ndims, 0);
358 : FILE *pGrad, *pCount;
359 20 : pGrad = fopen(gradname.c_str(), "w");
360 20 : pCount = fopen(countname.c_str(), "w");
361 : char *buffer1, *buffer2;
362 20 : buffer1 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
363 20 : buffer2 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
364 20 : setvbuf(pGrad, buffer1, _IOFBF, (sizeof(double)) * sampleSize * ndims);
365 20 : setvbuf(pCount, buffer2, _IOFBF, (sizeof(double)) * sampleSize * ndims);
366 20 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
367 20 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
368 65510 : for (size_t i = 0; i < sampleSize; ++i) {
369 130364 : for (size_t j = 0; j < ndims; ++j) {
370 130364 : pos[j] = table[j][i];
371 130364 : fprintf(pGrad, " %.9f", table[j][i]);
372 130364 : fprintf(pCount, " %.9f", table[j][i]);
373 : }
374 65510 : fprintf(pCount, " %lu\n", getCount(pos, true));
375 65510 : vector<double> f = getGradient(pos, true);
376 130364 : for (size_t j = 0; j < ndims; ++j) {
377 130364 : fprintf(pGrad, " %.9f", (f[j] / outputunit));
378 : }
379 : fprintf(pGrad, "\n");
380 : }
381 20 : fclose(pGrad);
382 20 : fclose(pCount);
383 20 : free(buffer1);
384 20 : free(buffer2);
385 20 : if (ndims == 1) {
386 24 : write1DPMF(filename);
387 : }
388 20 : }
389 :
390 2 : void DRRForceGrid::writeDivergence(const string &filename) const {
391 4 : string divname = filename + suffix + ".div";
392 2 : vector<double> pos(ndims, 0);
393 : FILE *pDiv;
394 2 : pDiv = fopen(divname.c_str(), "w");
395 2 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pDiv);
396 64800 : for (size_t i = 0; i < sampleSize; ++i) {
397 129600 : for (size_t j = 0; j < ndims; ++j) {
398 129600 : pos[j] = table[j][i];
399 129600 : fprintf(pDiv, " %.9f", table[j][i]);
400 : }
401 64800 : const double divergence = getDivergence(pos);
402 64800 : fprintf(pDiv, " %.9f", (divergence / outputunit));
403 : fprintf(pDiv, "\n");
404 : }
405 2 : fclose(pDiv);
406 2 : }
407 :
408 106 : bool ABF::store_getbias(const vector<double> &pos, const vector<double> &f,
409 : vector<double> &fbias) {
410 106 : if (!isInBoundary(pos)) {
411 : std::fill(begin(fbias), end(fbias), 0.0);
412 : return false;
413 : }
414 106 : const size_t baseaddr = sampleAddress(pos);
415 : unsigned long int &count = samples[baseaddr];
416 106 : ++count;
417 106 : double factor = 2 * (static_cast<double>(count)) / mFullSamples - 1;
418 106 : auto it_fa = begin(forces) + baseaddr * ndims;
419 : auto it_fb = begin(fbias);
420 : auto it_f = begin(f);
421 : auto it_maxFactor = begin(mMaxFactors);
422 178 : do {
423 : // Clamp to [0,maxFactors]
424 178 : factor = factor < 0 ? 0 : factor > (*it_maxFactor) ? (*it_maxFactor) : factor;
425 178 : (*it_fa) += (*it_f); // Accumulate instantaneous force
426 178 : (*it_fb) = factor * (*it_fa) * (-1.0) /
427 178 : static_cast<double>(count); // Calculate bias force
428 : ++it_fa;
429 : ++it_fb;
430 : ++it_f;
431 : ++it_maxFactor;
432 : } while (it_f != end(f));
433 :
434 : return true;
435 : }
436 :
437 1 : ABF ABF::mergewindow(const ABF &aWA, const ABF &aWB) {
438 1 : const vector<DRRAxis> dR = merge(aWA.dimensions, aWB.dimensions);
439 1 : const string suffix = ".abf";
440 1 : ABF result(dR, suffix);
441 1 : const size_t nrows = result.sampleSize;
442 1 : const size_t ncols = result.ndims;
443 1 : vector<double> pos(ncols, 0);
444 201 : for (size_t i = 0; i < nrows; ++i) {
445 200 : for (size_t j = 0; j < ncols; ++j) {
446 200 : pos[j] = result.table[j][i];
447 : }
448 200 : const unsigned long int countA = aWA.getCount(pos);
449 200 : const unsigned long int countB = aWB.getCount(pos);
450 200 : const vector<double> aForceA = aWA.getAccumulatedForces(pos);
451 200 : const vector<double> aForceB = aWB.getAccumulatedForces(pos);
452 200 : result.store(pos, aForceA, countA);
453 200 : result.store(pos, aForceB, countB);
454 : }
455 1 : return result;
456 : }
457 :
458 195083 : vector<double> CZAR::getGradient(const vector<double> &pos,
459 : bool SkipCheck) const {
460 195083 : vector<double> result(ndims, 0);
461 195083 : if (!SkipCheck) {
462 162000 : if (!isInBoundary(pos)) {
463 : return result;
464 : }
465 : }
466 195083 : if (kbt <= std::numeric_limits<double>::epsilon()) {
467 0 : std::cerr << "ERROR! The kbt shouldn't be zero when use CZAR estimator. "
468 : << '\n';
469 0 : std::abort();
470 : }
471 195083 : const size_t baseaddr = sampleAddress(pos);
472 195083 : const vector<double> log_deriv(getCountsLogDerivative(pos));
473 : const unsigned long int &count = samples[baseaddr];
474 195083 : if (count == 0)
475 : return result;
476 195008 : auto it_fa = begin(forces) + baseaddr * ndims;
477 : std::transform(it_fa, it_fa + ndims, begin(log_deriv), begin(result),
478 : [&count, this](double fa, double ld) {
479 389408 : return fa * (-1.0) / count - kbt * ld;
480 779424 : });
481 195008 : return result;
482 : }
483 :
484 1 : CZAR CZAR::mergewindow(const CZAR &cWA, const CZAR &cWB) {
485 1 : const vector<DRRAxis> dR = merge(cWA.dimensions, cWB.dimensions);
486 1 : const double newkbt = cWA.kbt;
487 1 : const string suffix = ".czar";
488 : CZAR result(dR, suffix, newkbt);
489 1 : const size_t nrows = result.sampleSize;
490 1 : const size_t ncols = result.ndims;
491 1 : vector<double> pos(ncols, 0);
492 201 : for (size_t i = 0; i < nrows; ++i) {
493 200 : for (size_t j = 0; j < ncols; ++j) {
494 200 : pos[j] = result.table[j][i];
495 : }
496 200 : const unsigned long int countA = cWA.getCount(pos);
497 200 : const unsigned long int countB = cWB.getCount(pos);
498 200 : const vector<double> aForceA = cWA.getAccumulatedForces(pos);
499 200 : const vector<double> aForceB = cWB.getAccumulatedForces(pos);
500 200 : result.store(pos, aForceA, countA);
501 200 : result.store(pos, aForceB, countB);
502 : }
503 1 : return result;
504 : }
505 :
506 10 : void CZAR:: writeZCount(const string &filename) const {
507 10 : string countname = filename + ".zcount";
508 10 : vector<double> pos(ndims, 0);
509 : FILE *pCount;
510 10 : pCount = fopen(countname.c_str(), "w");
511 : char *buffer;
512 10 : buffer = (char *)malloc((sizeof(double)) * sampleSize * ndims);
513 10 : setvbuf(pCount, buffer, _IOFBF, (sizeof(double)) * sampleSize * ndims);
514 10 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
515 32755 : for (size_t i = 0; i < sampleSize; ++i) {
516 65182 : for (size_t j = 0; j < ndims; ++j) {
517 65182 : pos[j] = table[j][i];
518 65182 : fprintf(pCount, " %.9f", table[j][i]);
519 : }
520 32755 : fprintf(pCount, " %lu\n", getCount(pos, true));
521 : }
522 10 : fclose(pCount);
523 10 : free(buffer);
524 10 : }
525 :
526 : }
527 5874 : }
528 :
529 : #endif
|
