Triumvirate C++ API 0.5.0
Three-point clustering measurements in large-scale structure analyses.
Loading...
Searching...
No Matches
parameters.cpp
Go to the documentation of this file.
1// Copyright (C) [GPLv3 Licence]
2//
3// This file is part of the Triumvirate program. See the COPYRIGHT
4// and LICENCE files at the top-level directory of this distribution
5// for details of copyright and licensing.
6//
7// This program is free software: you can redistribute it and/or modify it
8// under the terms of the GNU General Public License as published by
9// the Free Software Foundation, either version 3 of the License, or
10// (at your option) any later version.
11//
12// This program is distributed in the hope that it will be useful,
13// but WITHOUT ANY WARRANTY; without even the implied warranty of
14// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
15// See the GNU General Public License for more details.
16//
17// You should have received a copy of the GNU General Public License
18// along with this program. If not, see <https://www.gnu.org/licenses/>.
19
27#include "parameters.hpp"
28
29namespace trvs = trv::sys;
30
31namespace trv {
32
34 // Copy I/O parameters.
35 this->catalogue_dir = other.catalogue_dir;
36 this->measurement_dir = other.measurement_dir;
40 this->output_tag = other.output_tag;
41
42 // Copy mesh sampling parameters.
43 for (int i = 0; i < 3; i++) {
44 this->boxsize[i] = other.boxsize[i];
45 this->ngrid[i] = other.ngrid[i];
46 }
47 this->alignment = other.alignment;
48 this->padscale = other.padscale;
49 this->padfactor = other.padfactor;
50 this->assignment = other.assignment;
51 this->interlace = other.interlace;
52 this->volume = other.volume;
53 this->nmesh = other.nmesh;
55
56 // Copy measurement parameters.
57 this->catalogue_type = other.catalogue_type;
58 this->statistic_type = other.statistic_type;
59 this->npoint = other.npoint;
60 this->space = other.space;
61 this->ell1 = other.ell1;
62 this->ell2 = other.ell2;
63 this->ELL = other.ELL;
64 this->i_wa = other.i_wa;
65 this->j_wa = other.j_wa;
66 this->form = other.form;
67 this->norm_convention = other.norm_convention;
68 this->binning = other.binning;
69 this->bin_min = other.bin_min;
70 this->bin_max = other.bin_max;
71 this->num_bins = other.num_bins;
72 this->idx_bin = other.idx_bin;
73
74 // Copy misc parameters.
75 this->fftw_scheme = other.fftw_scheme;
77 this->use_fftw_wisdom = other.use_fftw_wisdom;
81 this->verbose = other.verbose;
82}
83
84int ParameterSet::read_from_file(char* parameter_filepath) {
85 // ---------------------------------------------------------------------
86 // Initialisation
87 // ---------------------------------------------------------------------
88
89 // Load parameter file.
90 std::string param_filepath = parameter_filepath;
91
92 std::ifstream fin(param_filepath.c_str());
93
94 // Initialise temporary variables to hold the extracted parameters.
95 char catalogue_dir_[1024] = "";
96 char measurement_dir_[1024] = "";
97 char data_catalogue_file_[1024] = "";
98 char rand_catalogue_file_[1024] = "";
99 char catalogue_columns_[1024] = "";
100 char output_tag_[1024] = "";
101
102 double boxsize_x, boxsize_y, boxsize_z;
103 int ngrid_x, ngrid_y, ngrid_z;
104
105 char alignment_[16] = "";
106 char padscale_[16] = "";
107 char assignment_[16] = "";
108 char interlace_[16] = "";
109
110 char catalogue_type_[16] = "";
111 char statistic_type_[16] = "";
112 char form_[16] = "";
113 char norm_convention_[16] = "";
114 char binning_[16] = "";
115
116 char fftw_scheme_[16] = "";
117 char use_fftw_wisdom_[1024] = "";
118 char save_binned_vectors_[1024] = "";
119
120 // ---------------------------------------------------------------------
121 // Extraction
122 // ---------------------------------------------------------------------
123
124 std::string line_str;
125 char dummy_str[1024], dummy_equal[1024];
126 while (std::getline(fin, line_str)) {
127 // Check if the line is a parameter assignment.
128 if (line_str.find("#") == 0) {
129 continue;
130 } // skip comment lines
131 if (
132 std::sscanf(
133 line_str.data(), "%1023s %1023s %1023s",
134 dummy_str, dummy_equal, dummy_str
135 ) != 3
136 ) {
137 continue;
138 } // skip non-relation lines
139 if (std::strcmp(dummy_equal, "=") != 0) {
140 continue;
141 } // skip non-assignment lines
142
143 // Define convenience function for scanning string parameters.
144 auto scan_par_str = [line_str, dummy_str, dummy_equal](
145 const char* par_name, const char* fmt, const char* par_value
146 ) {
147 if (line_str.find(par_name) != std::string::npos) {
148 std::sscanf(
149 line_str.data(), fmt, dummy_str, dummy_equal, par_value
150 );
151 }
152 };
153
154 // -- I/O ------------------------------------------------------------
155
156 scan_par_str("catalogue_dir", "%s %s %s", catalogue_dir_);
157 scan_par_str("measurement_dir", "%s %s %s", measurement_dir_);
158 scan_par_str("data_catalogue_file", "%s %s %s", data_catalogue_file_);
159 scan_par_str("rand_catalogue_file", "%s %s %s", rand_catalogue_file_);
160 scan_par_str("catalogue_columns", "%s %s %s", catalogue_columns_);
161 scan_par_str("output_tag", "%s %s %s", output_tag_);
162
163 // -- Mesh sampling --------------------------------------------------
164
165 if (line_str.find("boxsize_x") != std::string::npos) {
166 std::sscanf(
167 line_str.data(), "%1023s %1023s %lg",
168 dummy_str, dummy_equal, &boxsize_x
169 );
170 }
171 if (line_str.find("boxsize_y") != std::string::npos) {
172 std::sscanf(
173 line_str.data(), "%1023s %1023s %lg",
174 dummy_str, dummy_equal, &boxsize_y
175 );
176 }
177 if (line_str.find("boxsize_z") != std::string::npos) {
178 std::sscanf(
179 line_str.data(), "%1023s %1023s %lg",
180 dummy_str, dummy_equal, &boxsize_z
181 );
182 }
183
184 if (line_str.find("ngrid_x") != std::string::npos) {
185 std::sscanf(
186 line_str.data(), "%1023s %1023s %d",
187 dummy_str, dummy_equal, &ngrid_x
188 );
189 }
190 if (line_str.find("ngrid_y") != std::string::npos) {
191 std::sscanf(
192 line_str.data(), "%1023s %1023s %d",
193 dummy_str, dummy_equal, &ngrid_y
194 );
195 }
196 if (line_str.find("ngrid_z") != std::string::npos) {
197 std::sscanf(
198 line_str.data(), "%1023s %1023s %d", dummy_str, dummy_equal, &ngrid_z
199 );
200 }
201
202 scan_par_str("alignment", "%1023s %1023s %1023s", alignment_);
203 scan_par_str("padscale", "%1023s %1023s %1023s", padscale_);
204
205 if (line_str.find("padfactor") != std::string::npos) {
206 std::sscanf(
207 line_str.data(), "%1023s %1023s %lg",
208 dummy_str, dummy_equal, &this->padfactor
209 );
210 }
211
212 scan_par_str("assignment", "%1023s %1023s %1023s", assignment_);
213 scan_par_str("interlace", "%1023s %1023s %1023s", interlace_);
214
215 // -- Measurement ----------------------------------------------------
216
217 scan_par_str("catalogue_type", "%1023s %1023s %1023s", catalogue_type_);
218 scan_par_str("statistic_type", "%1023s %1023s %1023s", statistic_type_);
219 scan_par_str("form", "%1023s %1023s %1023s", form_);
220 scan_par_str("norm_convention", "%1023s %1023s %1023s", norm_convention_);
221 scan_par_str("binning", "%1023s %1023s %1023s", binning_);
222
223 if (line_str.find("ell1") != std::string::npos) {
224 std::sscanf(
225 line_str.data(), "%1023s %1023s %d",
226 dummy_str, dummy_equal, &this->ell1
227 );
228 }
229 if (line_str.find("ell2") != std::string::npos) {
230 std::sscanf(
231 line_str.data(), "%1023s %1023s %d",
232 dummy_str, dummy_equal, &this->ell2
233 );
234 }
235 if (line_str.find("ELL") != std::string::npos) {
236 std::sscanf(
237 line_str.data(), "%1023s %1023s %d",
238 dummy_str, dummy_equal, &this->ELL
239 );
240 }
241
242 if (line_str.find("i_wa") != std::string::npos) {
243 std::sscanf(
244 line_str.data(), "%1023s %1023s %d",
245 dummy_str, dummy_equal, &this->i_wa
246 );
247 }
248 if (line_str.find("j_wa") != std::string::npos) {
249 std::sscanf(
250 line_str.data(), "%1023s %1023s %d",
251 dummy_str, dummy_equal, &this->j_wa
252 );
253 }
254
255 if (line_str.find("bin_min") != std::string::npos) {
256 std::sscanf(
257 line_str.data(), "%1023s %1023s %lg",
258 dummy_str, dummy_equal, &this->bin_min
259 );
260 }
261 if (line_str.find("bin_max") != std::string::npos) {
262 std::sscanf(
263 line_str.data(), "%1023s %1023s %lg",
264 dummy_str, dummy_equal, &this->bin_max
265 );
266 }
267
268 if (line_str.find("num_bins") != std::string::npos) {
269 std::sscanf(
270 line_str.data(), "%1023s %1023s %d",
271 dummy_str, dummy_equal, &this->num_bins
272 );
273 }
274 if (line_str.find("idx_bin") != std::string::npos) {
275 std::sscanf(
276 line_str.data(), "%1023s %1023s %d",
277 dummy_str, dummy_equal, &this->idx_bin
278 );
279 }
280
281 // -- Misc -------------------------------------------------------------
282
283 scan_par_str("fftw_scheme", "%1023s %1023s %1023s", fftw_scheme_);
284 scan_par_str("use_fftw_wisdom", "%1023s %1023s %1023s", use_fftw_wisdom_);
285 scan_par_str(
286 "save_binned_vectors", "%1023s %1023s %1023s", save_binned_vectors_
287 );
288
289 if (line_str.find("verbose") != std::string::npos) {
290 std::sscanf(
291 line_str.data(), "%1023s %1023s %d",
292 dummy_str, dummy_equal, &this->verbose
293 );
294 }
295 }
296
297 // ---------------------------------------------------------------------
298 // Attribution
299 // ---------------------------------------------------------------------
300
301 // Attribute numerical parameters (directly extracted above).
302
303 // Attribute string parameters.
304 this->catalogue_dir = catalogue_dir_;
305 this->measurement_dir = measurement_dir_;
306 this->data_catalogue_file = data_catalogue_file_;
307 this->rand_catalogue_file = rand_catalogue_file_;
308 this->catalogue_columns = catalogue_columns_;
309 this->output_tag = output_tag_;
310
311 this->alignment = alignment_;
312 this->padscale = padscale_;
313 this->assignment = assignment_;
314 this->interlace = interlace_;
315
316 this->catalogue_type = catalogue_type_;
317 this->statistic_type = statistic_type_;
318 this->form = form_;
319 this->norm_convention = norm_convention_;
320 this->binning = binning_;
321
322 this->fftw_scheme = fftw_scheme_;
323 this->use_fftw_wisdom = use_fftw_wisdom_;
324 this->save_binned_vectors = save_binned_vectors_;
325
326 // Attribute derived parameters.
327 this->boxsize[0] = boxsize_x;
328 this->boxsize[1] = boxsize_y;
329 this->boxsize[2] = boxsize_z;
330
331 this->ngrid[0] = ngrid_x;
332 this->ngrid[1] = ngrid_y;
333 this->ngrid[2] = ngrid_z;
334
335 this->volume = boxsize_x * boxsize_y * boxsize_z;
336 this->nmesh = ngrid_x * ngrid_y * ngrid_z;
337
338 // ---------------------------------------------------------------------
339 // Debugging mode
340 // ---------------------------------------------------------------------
341
342#ifdef DBG_PARS
343 // Define convenience function for displaying debugged parameters.
344 auto debug_par_str = [](const std::string& name, const std::string& value) {
345 std::cout << name << ": " << value << std::endl;
346 };
347 auto debug_par_int = [](const std::string& name, int value) {
348 std::cout << name << ": " << value << std::endl;
349 };
350 auto debug_par_longlong = [](const std::string& name, long long value) {
351 std::cout << name << ": " << value << std::endl;
352 };
353 auto debug_par_double = [](const std::string& name, double value) {
354 std::cout << name << ": " << value << std::endl;
355 };
356
357 // Display debugged parameters.
358 debug_par_str("catalogue_dir", this->catalogue_dir);
359 debug_par_str("measurement_dir", this->measurement_dir);
360 debug_par_str("data_catalogue_file", this->data_catalogue_file);
361 debug_par_str("rand_catalogue_file", this->rand_catalogue_file);
362 debug_par_str("catalogue_columns", this->catalogue_columns);
363 debug_par_str("output_tag", this->output_tag);
364
365 debug_par_str("alignment", this->alignment);
366 debug_par_str("padscale", this->padscale);
367 debug_par_str("assignment", this->assignment);
368 debug_par_str("interlace", this->interlace);
369
370 debug_par_str("catalogue_type", this->catalogue_type);
371 debug_par_str("statistic_type", this->statistic_type);
372 debug_par_str("form", this->form);
373 debug_par_str("norm_convention", this->norm_convention);
374 debug_par_str("binning", this->binning);
375
376 debug_par_str("fftw_scheme", this->fftw_scheme);
377 debug_par_str("use_fftw_wisdom", this->use_fftw_wisdom);
378 debug_par_str("save_binned_vectors", this->save_binned_vectors);
379
380 debug_par_int("ngrid[0]", this->ngrid[0]);
381 debug_par_int("ngrid[1]", this->ngrid[1]);
382 debug_par_int("ngrid[2]", this->ngrid[2]);
383
384 debug_par_longlong("nmesh", this->nmesh);
385
386 debug_par_int("ell1", this->ell1);
387 debug_par_int("ell2", this->ell2);
388 debug_par_int("ELL", this->ELL);
389 debug_par_int("i_wa", this->i_wa);
390 debug_par_int("j_wa", this->j_wa);
391
392 debug_par_int("num_bins", this->num_bins);
393 debug_par_int("idx_bin", this->idx_bin);
394
395 debug_par_double("boxsize[0]", this->boxsize[0]);
396 debug_par_double("boxsize[1]", this->boxsize[1]);
397 debug_par_double("boxsize[2]", this->boxsize[2]);
398 debug_par_double("volume", this->volume);
399 debug_par_double("padfactor", this->padfactor);
400 debug_par_double("bin_min", this->bin_min);
401 debug_par_double("bin_max", this->bin_max);
402#endif // DBG_PARS
403
404 return this->validate();
405}
406
409
410 // Validate and derive string parameters.
411 // Any duplicate '/' has no effect.
412 if (
413 this->catalogue_dir.find_first_not_of(" \t\n\r\v\f") != std::string::npos
414 ) {
415 this->catalogue_dir += "/"; // transmutation
416 }
417 if (
418 this->measurement_dir.find_first_not_of(" \t\n\r\v\f") == std::string::npos
419 ) {
420 this->measurement_dir = "./"; // transmutation
421 } else {
422 this->measurement_dir += "/"; // transmutation
423 }
424 if (this->catalogue_type == "survey") {
425 if (this->data_catalogue_file != "") {
426 if (this->data_catalogue_file.rfind("/", 0) != 0) {
428 + this->data_catalogue_file;
429 } // transmutation
430 }
431 if (this->rand_catalogue_file != "") {
432 if (this->rand_catalogue_file.rfind("/", 0) != 0) {
434 + this->rand_catalogue_file;
435 } // transmutation
436 }
437 } else
438 if (this->catalogue_type == "random") {
439 this->data_catalogue_file = ""; // transmutation
440 if (this->rand_catalogue_file != "") {
441 if (this->rand_catalogue_file.rfind("/", 0) != 0) {
443 + this->rand_catalogue_file;
444 } // transmutation
445 }
446 } else
447 if (this->catalogue_type == "sim") {
448 if (this->data_catalogue_file != "") {
449 if (this->data_catalogue_file.rfind("/", 0) != 0) {
451 + this->data_catalogue_file;
452 } // transmutation
453 }
454 this->rand_catalogue_file = ""; // transmutation
455 } else
456 if (this->catalogue_type == "none") {
457 // Nothing should happen.
458 } else {
459#ifndef TRV_EXTCALL
460 if (trvs::currTask == 0) {
462 "Catalogue type must be 'survey', 'random', 'sim' or 'none: "
463 "`catalogue_type` = '%s'.",
464 this->catalogue_type.c_str()
465 );
466 }
468 "Catalogue type must be 'survey', 'random', 'sim' or 'none': "
469 "`catalogue_type` = '%s'.\n",
470 this->catalogue_type.c_str()
471 );
472#endif // !TRV_EXTCALL
473 }
474
475 if (!(this->alignment == "centre" || this->alignment == "pad")) {
476 if (trvs::currTask == 0) {
478 "Box alignment must be 'centre' or 'pad': `alignment` = '%s'.",
479 this->alignment.c_str()
480 );
481 }
483 "Box alignment must be 'centre' or 'pad': `alignment` = '%s'.\n",
484 this->alignment.c_str()
485 );
486 }
487 if (!(this->padscale == "box" || this->padscale == "grid")) {
488 if (trvs::currTask == 0) {
490 "Pad scale must be 'box' or 'grid': `padscale` = '%s'.",
491 this->padscale.c_str()
492 );
493 }
495 "Pad scale must be 'box' or 'grid': `padscale` = '%s'.\n",
496 this->padscale.c_str()
497 );
498 }
499
500 if (this->assignment == "ngp") {
501 this->assignment_order = 1;
502 } else
503 if (this->assignment == "cic") {
504 this->assignment_order = 2;
505 } else
506 if (this->assignment == "tsc") {
507 this->assignment_order = 3;
508 } else
509 if (this->assignment == "pcs") {
510 this->assignment_order = 4;
511 } else {
512 if (trvs::currTask == 0) {
514 "Mesh assignment scheme must be "
515 "'ngp', 'cic', 'tsc' or 'pcs': `assignment` = '%s'.",
516 this->assignment.c_str()
517 );
518 }
520 "Mesh assignment scheme must be "
521 "'ngp', 'cic', 'tsc' or 'pcs': `assignment` = '%s'.\n",
522 this->assignment.c_str()
523 );
524 }
525 if (this->interlace == "true" || this->interlace == "on") {
526 this->interlace = "true"; // transmutation
527 } else
528 if (this->interlace == "false" || this->interlace == "off") {
529 this->interlace = "false"; // transmutation
530 } else {
531 if (trvs::currTask == 0) {
533 "Interlacing must be 'true'/'on' or 'false'/'off': "
534 "`interlace` = '%s'.",
535 this->interlace.c_str()
536 );
537 }
539 "Interlacing must be 'true'/'on' or 'false'/'off': "
540 "`interlace` = '%s'.\n",
541 this->interlace.c_str()
542 );
543 }
544
545 if (this->statistic_type == "powspec") {
546 this->npoint = "2pt"; this->space = "fourier"; // derivation
547 if (this->catalogue_type == "random" || this->catalogue_type == "none") {
548 if (trvs::currTask == 0) {
550 "Power spectrum requires 'sim' or 'survey' catalogue(s): "
551 "`catalogue_type` = '%s'.",
552 this->catalogue_type.c_str()
553 );
554 }
556 "Power spectrum requires 'sim' or 'survey' catalogue(s): "
557 "`catalogue_type` = '%s'.\n",
558 this->catalogue_type.c_str()
559 );
560 }
561 } else
562 if (this->statistic_type == "2pcf") {
563 this->npoint = "2pt"; this->space = "config"; // derivation
564 if (this->catalogue_type == "random" || this->catalogue_type == "none") {
565 if (trvs::currTask == 0) {
567 "Two-point correlation function requires 'sim' or 'survey' "
568 "catalogue(s): `catalogue_type` = '%s'.",
569 this->catalogue_type.c_str()
570 );
571 }
573 "Two-point correlation function requires 'sim' or 'survey' "
574 "catalogue(s): `catalogue_type` = '%s'.\n",
575 this->catalogue_type.c_str()
576 );
577 }
578 } else
579 if (this->statistic_type == "2pcf-win") {
580 this->npoint = "2pt"; this->space = "config"; // derivation
581 if (this->catalogue_type != "random") {
582 if (trvs::currTask == 0) {
584 "Two-point correlation function window requires 'random' catalogue: "
585 "`catalogue_type` = '%s'.",
586 this->catalogue_type.c_str()
587 );
588 }
590 "Two-point correlation function window requires 'random' catalogue: "
591 "`catalogue_type` = '%s'.\n",
592 this->catalogue_type.c_str()
593 );
594 }
595 } else
596 if (this->statistic_type == "bispec") {
597 this->npoint = "3pt"; this->space = "fourier"; // derivation
598 if (this->catalogue_type == "random" || this->catalogue_type == "none") {
599 if (trvs::currTask == 0) {
601 "Bispectrum requires 'sim' or 'survey' catalogue(s): "
602 "`catalogue_type` = '%s'.",
603 this->catalogue_type.c_str()
604 );
605 }
607 "Bispectrum requires 'sim' or 'survey' catalogue(s): "
608 "`catalogue_type` = '%s'.\n",
609 this->catalogue_type.c_str()
610 );
611 }
612 } else
613 if (this->statistic_type == "3pcf") {
614 this->npoint = "3pt"; this->space = "config"; // derivation
615 if (this->catalogue_type == "random" || this->catalogue_type == "none") {
616 if (trvs::currTask == 0) {
618 "Three-point correlation function requires 'sim' or 'survey' "
619 "catalogue(s): `catalogue_type` = '%s'.",
620 this->catalogue_type.c_str()
621 );
622 }
624 "Three-point correlation function requires 'sim' or 'survey' "
625 "catalogue(s): `catalogue_type` = '%s'.\n",
626 this->catalogue_type.c_str()
627 );
628 }
629 } else
630 if (
631 this->statistic_type == "3pcf-win"
632 || this->statistic_type == "3pcf-win-wa"
633 ) {
634 this->npoint = "3pt"; this->space = "config"; // derivation
635 if (this->catalogue_type != "random") {
636 if (trvs::currTask == 0) {
638 "Three-point correlation function window requires 'random' "
639 "catalogue: `catalogue_type` = '%s'.",
640 this->catalogue_type.c_str()
641 );
642 }
644 "Three-point correlation function window requires 'random' "
645 "catalogue: `catalogue_type` = '%s'.\n",
646 this->catalogue_type.c_str()
647 );
648 }
649 } else
650 if (this->statistic_type == "modes") {
651 this->npoint = "none"; this->space = "fourier"; // derivation
652 } else
653 if (this->statistic_type == "pairs") {
654 this->npoint = "none"; this->space = "config"; // derivation
655 } else {
656#ifndef TRV_EXTCALL
657 if (trvs::currTask == 0) {
659 "Statistic type is not recognised: `statistic_type` = '%s'.",
660 this->statistic_type.c_str()
661 );
662 }
664 "Statistic type is not recognised: `statistic_type` = '%s'.\n",
665 this->statistic_type.c_str()
666 );
667#endif // !TRV_EXTCALL
668 }
669 if (!(
670 this->form == "full"
671 || this->form == "diag"
672 || this->form == "off-diag"
673 || this->form == "row"
674 )) {
675 if (trvs::currTask == 0) {
677 "Three-point statistic form is not recognised: `form` = '%s'.",
678 this->form.c_str()
679 );
680 }
682 "Three-point statistic form is not recognised: `form` = '%s'.\n",
683 this->form.c_str()
684 );
685 } else {
686 if (this->form == "full" && this->ell1 == this->ell2) {
687 this->shape = "triu"; // derivation
688 } else {
689 this->shape = this->form; // derivation
690 }
691 }
692 if (!(
693 this->norm_convention == "none"
694 || this->norm_convention == "particle"
695 || this->norm_convention == "mesh"
696 || this->norm_convention == "mesh-mixed"
697 )) {
698 if (trvs::currTask == 0) {
700 "Normalisation convention must be 'mesh', 'particle', "
701 "'mesh' or 'mesh-mixed': `norm_convention` = '%s'.",
702 this->norm_convention.c_str()
703 );
704 }
706 "Normalisation convention must be 'mesh', 'particle', "
707 "'mesh' or 'mesh-mixed': `norm_convention` = '%s'.\n",
708 this->norm_convention.c_str()
709 );
710 }
711 if (this->norm_convention == "mesh-mixed" && this->npoint != "2pt") {
712 if (trvs::currTask == 0) {
714 "Normalisation convention 'mesh-mixed' only applies to "
715 "two-point statistics: `npoint` = '%s'.",
716 this->npoint.c_str()
717 );
718 }
720 "Normalisation convention 'mesh-mixed' only applies to "
721 "two-point statistics: `npoint` = '%s'.\n",
722 this->npoint.c_str()
723 );
724 }
725 if (!(
726 this->binning == "lin"
727 || this->binning == "log"
728 || this->binning == "linpad"
729 || this->binning == "logpad"
730 || this->binning == "custom"
731 )) {
732 if (trvs::currTask == 0) {
734 "Binning scheme is unrecognised: `binning` = '%s'.",
735 this->binning.c_str()
736 );
737 }
739 "Binning scheme is unrecognised: `binning` = '%s'.\n",
740 this->binning.c_str()
741 );
742 }
743
744 if (this->fftw_scheme == "estimate") {
745 this->fftw_planner_flag = FFTW_ESTIMATE; // derivation
746 } else
747 if (this->fftw_scheme == "measure" || this->fftw_scheme == "") {
748 this->fftw_scheme = "measure"; // transmutation
749 this->fftw_planner_flag = FFTW_MEASURE; // derivation
750 } else
751 if (this->fftw_scheme == "patient") {
752 this->fftw_planner_flag = FFTW_PATIENT; // derivation
753 } else {
754 if (trvs::currTask == 0) {
756 "FFTW planner scheme is not supported: `fftw_scheme` = '%s'.",
757 this->fftw_scheme.c_str()
758 );
759 }
761 "FFTW planner scheme is not supported: `fftw_scheme` = '%s'.\n",
762 this->fftw_scheme.c_str()
763 );
764 }
765
766 if (this->use_fftw_wisdom == "false" || this->use_fftw_wisdom == "") {
767 this->use_fftw_wisdom = ""; // transmutation
768 } else {
769 this->use_fftw_wisdom += "/"; // transmutation
770 }
771
772 if (this->use_fftw_wisdom != "") {
773 if (this->fftw_scheme != "measure" && this->fftw_scheme != "patient") {
774 if (trvs::currTask == 0) {
776 "FFTW wisdom is enabled but the planner scheme "
777 "is not 'measure' or 'patient': "
778 "`fftw_scheme` = '%s'.",
779 this->fftw_scheme.c_str()
780 );
781 }
783 "FFTW wisdom is enabled but the planner scheme "
784 "is not 'measure' or 'patient': "
785 "`fftw_scheme` = '%s'.\n",
786 this->fftw_scheme.c_str()
787 );
788 }
789
790 char fftw_wisdom_file_f_[1024];
791 char fftw_wisdom_file_b_[1024];
792
793#if defined(TRV_USE_OMP) && defined(TRV_USE_FFTWOMP)
794 std::snprintf(
795 fftw_wisdom_file_f_, sizeof(fftw_wisdom_file_f_),
796 "%sfftw_omp_cif_%dx%dx%d.wisdom",
797 this->use_fftw_wisdom.c_str(),
798 this->ngrid[0], this->ngrid[1], this->ngrid[2]
799 );
800 std::snprintf(
801 fftw_wisdom_file_b_, sizeof(fftw_wisdom_file_b_),
802 "%sfftw_omp_cib_%dx%dx%d.wisdom",
803 this->use_fftw_wisdom.c_str(),
804 this->ngrid[0], this->ngrid[1], this->ngrid[2]
805 );
806#else // !TRV_USE_OMP || !TRV_USE_FFTWOMP
807 std::snprintf(
808 fftw_wisdom_file_f_, sizeof(fftw_wisdom_file_f_),
809 "%sfftw_cif_%dx%dx%d.wisdom",
810 this->use_fftw_wisdom.c_str(),
811 this->ngrid[0], this->ngrid[1], this->ngrid[2]
812 );
813 std::snprintf(
814 fftw_wisdom_file_b_, sizeof(fftw_wisdom_file_b_),
815 "%sfftw_cib_%dx%dx%d.wisdom",
816 this->use_fftw_wisdom.c_str(),
817 this->ngrid[0], this->ngrid[1], this->ngrid[2]
818 );
819#endif // TRV_USE_OMP && TRV_USE_FFTWOMP
820
821 this->fftw_wisdom_file_f = fftw_wisdom_file_f_;
822 this->fftw_wisdom_file_b = fftw_wisdom_file_b_;
823 }
824
825 char default_bvec_sfilepath[1024];
826 std::snprintf(
827 default_bvec_sfilepath, sizeof(default_bvec_sfilepath),
828 "%s/binned_vectors%s",
829 this->measurement_dir.c_str(), this->output_tag.c_str()
830 );
831 if (this->save_binned_vectors == "false") {
832 this->save_binned_vectors = ""; // transmutation
833 } else
834 if (this->save_binned_vectors == "true") {
835 this->save_binned_vectors = default_bvec_sfilepath; // transmutation
836 } else
837 if (this->save_binned_vectors != "") {
838 // Check whether path is absolute.
839 if (this->save_binned_vectors.rfind("/", 0) != 0) {
841 + this->save_binned_vectors;
842 } // transmutation
843 }
844
845 // Validate and derive numerical parameters.
846 this->volume =
847 this->boxsize[0] * this->boxsize[1] * this->boxsize[2]; // derivation
848 this->nmesh = static_cast<long long>(this->ngrid[0])
849 * this->ngrid[1] * this->ngrid[2]; // derivation
850
851 if (this->volume <= 0.) {
852 if (trvs::currTask == 0) {
854 "Derived total box volume is non-positive: `volume` = '%.6e'. "
855 "Possible numerical overflow due to large `boxsize`, "
856 "or `boxsize` is unset.",
857 this->volume
858 );
859 }
861 "Derived total box volume is non-positive: `volume` = '%.6e'. "
862 "Possible numerical overflow due to large `boxsize`, "
863 "or `boxsize` is unset.\n",
864 this->volume
865 );
866 }
867 if (this->nmesh <= 0) {
868 if (trvs::currTask == 0) {
870 "Derived total mesh grid number is non-positive: `nmesh` = '%lld'. "
871 "Possible numerical overflow due to large `ngrid`, "
872 "or `ngrid` is unset.",
873 this->nmesh
874 );
875 }
877 "Derived total mesh grid number is non-positive: `nmesh` = '%lld'. "
878 "Possible numerical overflow due to large `ngrid`, "
879 "or `ngrid` is unset.\n",
880 this->nmesh
881 );
882 }
883
884 if (this->alignment == "pad") {
885 if (this->padfactor < 0.) {
886 if (trvs::currTask == 0) {
888 "Padding is enabled but the padding factor is negative: "
889 "`padfactor` = '%lg'.",
890 this->padfactor
891 );
892 }
894 "Padding is enabled but the padding factor is negative: "
895 "`padfactor` = '%lg'.\n",
896 this->padfactor
897 );
898 }
899 if (this->padscale == "box" && this->padfactor >= 1.) {
900 if (trvs::currTask == 0) {
902 "Padding is enabled but the %s padding factor is too large "
903 "for the box size: `padfactor` = '%lg'.",
904 this->padscale.c_str(), this->padfactor
905 );
906 }
908 "Padding is enabled but the %s padding factor is too large "
909 "for the box size: `padfactor` = '%lg'.\n",
910 this->padscale.c_str(), this->padfactor
911 );
912 }
913 if (this->padscale == "grid" && (
914 this->padfactor >= this->ngrid[0]
915 || this->padfactor >= this->ngrid[1]
916 || this->padfactor >= this->ngrid[2]
917 )) {
918 if (trvs::currTask == 0) {
920 "Padding is enabled but the %s padding factor is too large "
921 "for the mesh grid numbers: `padfactor` = '%lg'.",
922 this->padscale.c_str(), this->padfactor
923 );
924 }
926 "Padding is enabled but the %s padding factor is too large "
927 "for the mesh grid numbers: `padfactor` = '%lg'.\n",
928 this->padscale.c_str(), this->padfactor
929 );
930 }
931 }
932
933 if (this->bin_min < 0.) {
934 if (trvs::currTask == 0) {
935 trvs::logger.error("Lower bin edge must be non-negative.");
936 }
938 "Lower bin edge must be non-negative.\n"
939 );
940 }
941 if (this->bin_min >= this->bin_max) {
942 if (trvs::currTask == 0) {
944 "Lower bin edge must be less than the upper bin edge."
945 );
946 }
948 "Lower bin edge must be less than the upper bin edge.\n"
949 );
950 }
951 if (this->space == "fourier") {
952 double wavenum_nyquist = M_PI
953 * *std::min_element(this->ngrid, this->ngrid + 3)
954 / *std::max_element(this->boxsize, this->boxsize + 3);
955 if (this->bin_min > wavenum_nyquist) {
956 if (trvs::currTask == 0) {
958 "Lower wavenumber limit exceeds the Nyquist wavenumber %.4f.",
959 wavenum_nyquist
960 );
961 }
962 }
963 } else
964 if (this->space == "config") {
965 double separation_nyquist = 2
966 * *std::max_element(this->boxsize, this->boxsize + 3)
967 / *std::min_element(this->ngrid, this->ngrid + 3);
968 if (this->bin_max < separation_nyquist) {
969 if (trvs::currTask == 0) {
971 "Upper separation limit undershoots the Nyquist scale %.4f.",
972 separation_nyquist
973 );
974 }
975 }
976 }
977
978 if (this->num_bins < 2) {
979 if (trvs::currTask == 0) {
980 trvs::logger.error("Number of bins `num_bins` must be >= 2.");
981 }
983 "Number of bins `num_bins` must be >= 2.\n"
984 );
985 }
986
987 if (
988 this->idx_bin < 0
989 && this->npoint == "3pt"
990 && this->form == "row"
991 ) {
992 if (trvs::currTask == 0) {
993 trvs::logger.error("Fixed row bin index `idx_bin` must be >= 0.");
994 }
996 "Fixed row bin index `idx_bin` must be >= 0.\n"
997 );
998 }
999
1000 // Check for parameter conflicts.
1001 if (this->binning == "linpad" || this->binning == "logpad") {
1002 // CAVEAT: See @ref trv::Binning.
1003 int nbin_pad = 5;
1004
1005 if (this->num_bins < nbin_pad + 2) {
1006 if (trvs::currTask == 0) {
1008 "Binning scheme '%s' requires `num_bins` >= %d.",
1009 this->binning.c_str(), nbin_pad + 2
1010 );
1011 }
1013 "Binning scheme '%s' requires `num_bins` >= %d.\n",
1014 this->binning.c_str(), nbin_pad + 2
1015 );
1016 }
1017 }
1018
1019 if (std::abs(this->idx_bin) >= this->num_bins) {
1020 if (trvs::currTask == 0) {
1022 "Bin index `idx_bin` must be < `num_bins` in absolute value."
1023 );
1024 }
1026 "Bin index `idx_bin` must be < `num_bins` in absolute value.\n"
1027 );
1028 }
1029
1030 if (this->npoint == "3pt" && this->interlace == "true") {
1031 this->interlace = "false"; // transmutation
1032
1033 if (trvs::currTask == 0) {
1035 "Interlacing is unsupported for three-point measurements. "
1036 "`interlace` is set to 'false'."
1037 );
1038 }
1039 }
1040
1041 if (
1042 (this->statistic_type == "modes" || this->statistic_type == "pairs")
1043 && this->save_binned_vectors == ""
1044 ) {
1045 this->save_binned_vectors = default_bvec_sfilepath; // transmutation
1046 if (trvs::currTask == 0) {
1048 "`save_binned_vectors` is overridden, as `statistic_type` is '%s' "
1049 "so binned vectors are saved as the output to the default path.",
1050 this->statistic_type.c_str()
1051 );
1052 }
1053 }
1054
1055 if (trvs::currTask == 0) {
1056 trvs::logger.stat("Parameters validated.");
1057 }
1058
1059 return 0;
1060}
1061
1062int ParameterSet::print_to_file(char* out_parameter_filepath) {
1063 // Create output file.
1064 std::FILE* ofileptr;
1065 if (!(ofileptr = std::fopen(out_parameter_filepath, "w"))) {
1066 if (trvs::currTask == 0) {
1068 "Non-existent or unwritable output directory: %s",
1069 this->measurement_dir.c_str()
1070 );
1071 }
1072 throw trvs::IOError(
1073 "Non-existent or unwritable output directory: %s\n",
1074 this->measurement_dir.c_str()
1075 );
1076 }
1077
1078 // Define convenience function for printing parameters.
1079 auto print_par_str = [ofileptr](
1080 const char* fmt, const std::string& par_val
1081 ) {
1082 std::fprintf(ofileptr, fmt, par_val.c_str());
1083 };
1084 auto print_par_int = [ofileptr](const char* fmt, int par_val) {
1085 std::fprintf(ofileptr, fmt, par_val);
1086 };
1087 auto print_par_double = [ofileptr](const char* fmt, double par_val) {
1088 std::fprintf(ofileptr, fmt, par_val);
1089 };
1090
1091 // Print parameters to file.
1092 print_par_str("catalogue_dir = %s\n", this->catalogue_dir);
1093 print_par_str("measurement_dir = %s\n", this->measurement_dir);
1094 print_par_str("data_catalogue_file = %s\n", this->data_catalogue_file);
1095 print_par_str("rand_catalogue_file = %s\n", this->rand_catalogue_file);
1096 print_par_str("catalogue_columns = %s\n", this->catalogue_columns);
1097 print_par_str("output_tag = %s\n", this->output_tag);
1098
1099 print_par_double("boxsize_x = %.3f\n", this->boxsize[0]);
1100 print_par_double("boxsize_y = %.3f\n", this->boxsize[1]);
1101 print_par_double("boxsize_z = %.3f\n", this->boxsize[2]);
1102 print_par_int("ngrid_x = %d\n", this->ngrid[0]);
1103 print_par_int("ngrid_y = %d\n", this->ngrid[1]);
1104 print_par_int("ngrid_z = %d\n", this->ngrid[2]);
1105
1106 print_par_double("volume = %.6e\n", this->volume);
1107 print_par_int("nmesh = %lld\n", this->nmesh);
1108
1109 print_par_str("alignment = %s\n", this->alignment);
1110 print_par_str("padscale = %s\n", this->padscale);
1111 print_par_double("padfactor = %.4f\n", this->padfactor);
1112
1113 print_par_str("assignment = %s\n", this->assignment);
1114 print_par_str("interlace = %s\n", this->interlace);
1115 print_par_int("assignment_order = %d\n", this->assignment_order);
1116
1117 print_par_str("catalogue_type = %s\n", this->catalogue_type);
1118 print_par_str("statistic_type = %s\n", this->statistic_type);
1119 print_par_str("npoint = %s\n", this->npoint);
1120 print_par_str("space = %s\n", this->space);
1121
1122 print_par_int("ell1 = %d\n", this->ell1);
1123 print_par_int("ell2 = %d\n", this->ell2);
1124 print_par_int("ELL = %d\n", this->ELL);
1125
1126 print_par_int("i_wa = %d\n", this->i_wa);
1127 print_par_int("j_wa = %d\n", this->j_wa);
1128
1129 print_par_str("form = %s\n", this->form);
1130 print_par_str("norm_convention = %s\n", this->norm_convention);
1131 print_par_str("binning = %s\n", this->binning);
1132 print_par_str("shape = %s\n", this->shape);
1133
1134 print_par_double("bin_min = %.4f\n", this->bin_min);
1135 print_par_double("bin_max = %.4f\n", this->bin_max);
1136 print_par_int("num_bins = %d\n", this->num_bins);
1137 print_par_int("idx_bin = %d\n", this->idx_bin);
1138
1139 print_par_str("fftw_scheme = %s\n", this->fftw_scheme);
1140 print_par_str("use_fftw_wisdom = %s\n", this->use_fftw_wisdom.c_str());
1141 print_par_str("fftw_wisdom_file_f = %s\n", this->fftw_wisdom_file_f.c_str());
1142 print_par_str("fftw_wisdom_file_b = %s\n", this->fftw_wisdom_file_b.c_str());
1143 print_par_str("save_binned_vectors = %s\n", this->save_binned_vectors);
1144 print_par_int("verbose = %d\n", this->verbose);
1145 print_par_int("fftw_planner_flag = %d\n", this->fftw_planner_flag);
1146
1147 std::fclose(ofileptr);
1148
1149 if (trvs::currTask == 0) {
1151 "Check used-parameter file for reference: %s", out_parameter_filepath
1152 );
1153 }
1154
1155 return 0;
1156}
1157
1159 // Set output file path to default.
1160 char ofilepath[1024];
1161 std::snprintf(
1162 ofilepath, sizeof(ofilepath), "%s/parameters_used%s",
1163 this->measurement_dir.c_str(), this->output_tag.c_str()
1164 );
1165
1166 return ParameterSet::print_to_file(ofilepath);
1167}
1168
1169} // namespace trv
Parameter set.
int validate()
Validate parameters.
std::string alignment
box alignment: {"centre" (default), "pad"}
std::string shape
std::string catalogue_columns
catalogue data columns (comma-separated without space)
long long nmesh
number of mesh grid cells
double bin_min
measurement range minimum (in Mpc/h or h/Mpc)
std::string output_tag
output tag
std::string assignment
mesh assignment scheme: {"ngp", "cic", "tsc" (default), "pcs"}
std::string save_binned_vectors
int ELL
spherical degree associated with the line of sight
unsigned fftw_planner_flag
derived FFTW planner flag
int i_wa
first order of the wide-angle correction term
int print_to_file()
Print out extracted parameters to the default file path in the output measurement directory.
std::string catalogue_dir
catalogue directory
double padfactor
padding factor
ParameterSet()=default
Construct a parameter set.
std::string fftw_scheme
FFTW scheme: {"estimate", "measure" (default), "patient"}.
double volume
box volume (in Mpc^3/h^3)
int read_from_file(char *parameter_filepath)
Read parameters from a file.
std::string norm_convention
std::string statistic_type
std::string measurement_dir
measurement/output directory
std::string space
coordinate space: {"fourier", "config"}
std::string catalogue_type
catalogue type: {"survey", "random", "sim", "none"}
int ell1
spherical degree associated with the first wavevector
std::string fftw_wisdom_file_b
backward-transform wisdom file path
int num_bins
number of measurement bins
std::string padscale
padding scale (if alignment is "pad"): {"box" (default), "grid"}
std::string interlace
interlacing switch: {"true"/"on", "false"/"off" (default)}
int assignment_order
order of the assignment scheme
std::string binning
std::string fftw_wisdom_file_f
derived FFTW wisdom file paths
int ell2
spherical degree associated with the second wavevector
int ngrid[3]
grid number in each dimension
std::string data_catalogue_file
data catalogue file
std::string npoint
N-point case: {"2pt", "3pt", "none"}
std::string rand_catalogue_file
random catalogue file
int idx_bin
fixed bin index in "off-fiag"/"row" form three-point measurements
int j_wa
second order of the wide-angle correction term
double boxsize[3]
box size (in Mpc/h) in each dimension
std::string use_fftw_wisdom
use FFTW wisdom: {"false" (default), <path-to-dir>}
double bin_max
measurement range maximum (in Mpc/h or h/Mpc)
Exception raised when an input/output operation fails.
Definition monitor.hpp:408
Exception raised when parameters are invalid.
Definition monitor.hpp:432
void info(const char *fmt_string,...)
Emit a information-level message.
Definition monitor.cpp:264
void error(const char *fmt_string,...)
Emit a warning-level message.
Definition monitor.cpp:282
void warn(const char *fmt_string,...)
Emit a warning-level message.
Definition monitor.cpp:273
void stat(const char *fmt_string,...)
Emit a status-level message.
Definition monitor.cpp:255
void reset_level(LogLevel level)
Reset the logger threshold level.
Definition monitor.cpp:156
Logger logger
default logger (at NSET logging level)
int currTask
current task
Definition monitor.cpp:44
Program parameter configuration.