Skip to content

Commit 353fe24

Browse files
cschreibcschreib
committed
Added support for chi2 grid for SFH confidence intervals
1 parent 31850b1 commit 353fe24

File tree

6 files changed

+216
-154
lines changed

6 files changed

+216
-154
lines changed

‎src/fast++-fitter.cpp‎

Lines changed: 23 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -778,18 +778,8 @@ void fitter_t::find_best_fits() {
778778
}
779779
}
780780

781-
vec1f delta_chi2;
782-
if (opts.interval_from_chi2) {
783-
for (auto& c : input.conf_interval) {
784-
if (c < 0.5) {
785-
delta_chi2.push_back(-get_chi2_from_conf_interval(1.0 - 2*c));
786-
} else {
787-
delta_chi2.push_back(+get_chi2_from_conf_interval(2*c - 1.0));
788-
}
789-
}
790-
}
791-
792781
if (opts.verbose) note("finding best fits...");
782+
793783
for (uint_t is : range(input.id)) {
794784
if (!silence_invalid_chi2 && !is_finite(output.best_chi2[is])) {
795785
warning("galaxy ", input.id[is], " has no best fit solution");
@@ -897,41 +887,34 @@ void fitter_t::find_best_fits() {
897887
// If asked, obtain confidence intervals from chi2 grid saved on disk
898888

899889
if (opts.interval_from_chi2) {
900-
std::ifstream in(chi2_filename[is], std::ios::binary);
901-
in.seekg(obchi2.hpos);
902-
903-
while (in) {
904-
uint32_t id;
905-
float chi2;
906-
vec1f p(gridder.nprop);
907-
908-
if (file::read(in, id) && file::read(in, chi2) && file::read(in, p)) {
909-
vec1u ids = gridder.grid_ids(id);
910-
for (uint_t ip : range(gridder.nparam+gridder.nprop)) {
911-
double v;
912-
if (ip < gridder.nparam) {
913-
v = output.grid[ip][ids[ip]];
914-
} else {
915-
v = p[ip - gridder.nparam];
916-
}
890+
iterate_best_chi2(is, [&](uint_t id, const vec1f& p, float chi2) {
891+
if (chi2 - best_chi2.safe[is] > input.delta_chi2.back()) return;
892+
893+
vec1u ids = gridder.grid_ids(id);
894+
for (uint_t ip : range(gridder.nparam+gridder.nprop)) {
895+
double v;
896+
if (ip < gridder.nparam) {
897+
v = output.grid[ip][ids[ip]];
898+
} else {
899+
v = p[ip - gridder.nparam];
900+
}
917901

918-
for (uint_t ic : range(input.conf_interval)) {
919-
if (chi2 - best_chi2.safe[is] < abs(delta_chi2[ic])) {
920-
float& saved = output.best_params(is,ip,1+ic);
921-
if (delta_chi2[ic] < 0.0) {
922-
if (saved > v || !is_finite(saved)) {
923-
saved = v;
924-
}
925-
} else {
926-
if (saved < v || !is_finite(saved)) {
927-
saved = v;
928-
}
902+
for (uint_t ic : range(input.conf_interval)) {
903+
if (chi2 - best_chi2.safe[is] < abs(input.delta_chi2[ic])) {
904+
float& saved = output.best_params(is,ip,1+ic);
905+
if (input.delta_chi2[ic] < 0.0) {
906+
if (saved > v || !is_finite(saved)) {
907+
saved = v;
908+
}
909+
} else {
910+
if (saved < v || !is_finite(saved)) {
911+
saved = v;
929912
}
930913
}
931914
}
932915
}
933916
}
934-
}
917+
});
935918
}
936919
}
937920
}

‎src/fast++-io.hpp‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
#ifndef FASTPP_IO_HPP
2+
#define FASTPP_IO_HPP
3+
4+
// Helper functions
5+
// ----------------
6+
7+
namespace vif {
8+
namespace file {
9+
template<typename S, typename T>
10+
bool write(S& s, const T& t) {
11+
s.write(reinterpret_cast<const char*>(&t), sizeof(T));
12+
return !s.fail();
13+
}
14+
15+
template<typename R, typename S, typename T>
16+
bool write_as(S& s, const T& t) {
17+
R r = t;
18+
s.write(reinterpret_cast<const char*>(&r), sizeof(R));
19+
return !s.fail();
20+
}
21+
22+
template<typename S>
23+
bool write(S& s, const std::string& t) {
24+
write_as<std::uint8_t>(s, t.size());
25+
s.write(t.c_str(), t.size());
26+
return !s.fail();
27+
}
28+
29+
template<typename S, std::size_t D, typename T>
30+
bool write(S& s, const vec<D,T>& t) {
31+
s.write(reinterpret_cast<const char*>(t.data.data()), sizeof(T)*t.size());
32+
return !s.fail();
33+
}
34+
35+
template<typename S, typename T>
36+
bool read(S& s, T& t) {
37+
s.read(reinterpret_cast<char*>(&t), sizeof(T));
38+
return !s.fail();
39+
}
40+
41+
template<typename S>
42+
bool read(S& s, std::string& t) {
43+
std::uint8_t n = 0;
44+
if (read(s, n)) {
45+
t.resize(n);
46+
s.read(&t[0], t.size());
47+
return !s.fail();
48+
} else {
49+
return false;
50+
}
51+
}
52+
53+
template<typename R, typename S, typename T>
54+
bool read_as(S& s, T& t) {
55+
R r;
56+
if (s.read(reinterpret_cast<char*>(&r), sizeof(R))) {
57+
t = r;
58+
}
59+
60+
return !s.fail();
61+
}
62+
63+
template<typename S, std::size_t D, typename T>
64+
bool read(S& s, vec<D,T>& t) {
65+
s.read(reinterpret_cast<char*>(t.data.data()), sizeof(T)*t.size());
66+
return !s.fail();
67+
}
68+
}
69+
}
70+
71+
#endif

‎src/fast++-read_input.cpp‎

Lines changed: 12 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -372,14 +372,19 @@ bool read_params(options_t& opts, input_state_t& state, const std::string& filen
372372
}
373373

374374
if (opts.n_sim != 0 || opts.interval_from_chi2) {
375-
state.conf_interval = 0.5*(1.0 - opts.c_interval/100.0);
376-
inplace_sort(opts.c_interval);
377-
vec1f cint = state.conf_interval;
378-
state.conf_interval.clear();
379-
for (float c : cint) {
380-
state.conf_interval.push_back(c);
381-
state.conf_interval.push_back(1.0 - c);
375+
vec1f cint = 0.5*(1.0 - opts.c_interval/100.0);
376+
inplace_sort(cint);
377+
for (uint_t ic : range(cint)) {
378+
state.conf_interval.push_back(cint[ic]);
379+
state.conf_interval.push_back(1.0 - cint[ic]);
380+
381+
double chi2 = get_chi2_from_conf_interval(opts.c_interval[ic]/100.0);
382+
state.delta_chi2.push_back(-chi2);
383+
state.delta_chi2.push_back(+chi2);
382384
}
385+
386+
print(state.conf_interval);
387+
print(state.delta_chi2);
383388
} else {
384389
opts.c_interval.clear();
385390
}

‎src/fast++-write_output.cpp‎

Lines changed: 67 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -274,7 +274,7 @@ void write_best_fits(const options_t& opts, const input_state_t& input, const gr
274274
}
275275

276276
void write_sfhs(const options_t& opts, const input_state_t& input, const gridder_t& gridder,
277-
const output_state_t& output) {
277+
const fitter_t& fitter, const output_state_t& output) {
278278

279279
if (opts.verbose) note("saving star formation histories");
280280

@@ -292,9 +292,20 @@ void write_sfhs(const options_t& opts, const input_state_t& input, const gridder
292292
unit = "Msol";
293293
}
294294

295+
uint_t nconf = input.conf_interval.size();
296+
bool has_median = false;
297+
if (!input.conf_interval.empty() && !opts.interval_from_chi2) {
298+
// Also add the median SFH for MC simulations
299+
++nconf;
300+
has_median = true;
301+
}
302+
295303
std::string header = "# t "+quantity+"(t) ("+unit+")";
296304
if (!opts.c_interval.empty()) {
297-
header += " med_"+quantity;
305+
if (has_median) {
306+
header += " med_"+quantity;
307+
}
308+
298309
for (float c : input.conf_interval) {
299310
float cc = 100*(1-2*c);
300311
std::string is = (cc < 0.0 ? "u" : "l")+to_string(round(abs(cc)));
@@ -308,9 +319,10 @@ void write_sfhs(const options_t& opts, const input_state_t& input, const gridder
308319
for (uint_t is : range(input.id)) {
309320
vec1u idm = gridder.grid_ids(output.best_model[is]);
310321
vec1d t = rgen_step(1e6, e10(gridder.auniv[idm[grid_id::z]]), opts.sfh_output_step);
322+
vec2f sfh = replicate(fnan, t.size(), 1+nconf);
311323

312-
// Get best fit SFH
313-
vec2f sfh(t.size(), 1+input.conf_interval.size()+(input.conf_interval.empty() ? 0 : 1)); {
324+
// Best fit
325+
{
314326
vec1d tsfh;
315327
if (!gridder.get_sfh(output.best_model[is], t, tsfh)) {
316328
return;
@@ -322,23 +334,58 @@ void write_sfhs(const options_t& opts, const input_state_t& input, const gridder
322334

323335
// Get confidence intervals
324336
if (!opts.c_interval.empty()) {
325-
vec2f sim_sfh(t.size(), opts.n_sim);
326-
for (uint_t ir : range(opts.n_sim)) {
327-
vec1d tsfh;
328-
if (!gridder.get_sfh(output.mc_best_model(is,ir), t, tsfh)) {
329-
return;
330-
}
337+
if (opts.n_sim != 0) {
338+
// From Monte Carlo simulations
339+
vec2f sim_sfh(t.size(), opts.n_sim);
340+
for (uint_t ir : range(opts.n_sim)) {
341+
vec1d tsfh;
342+
if (!gridder.get_sfh(output.mc_best_model(is,ir), t, tsfh)) {
343+
return;
344+
}
331345

332-
float mass = output.mc_best_props(is,prop_id::mass,ir);
333-
sim_sfh(_,ir) = tsfh*mass;
334-
}
346+
float mass = output.mc_best_props(is,prop_id::mass,ir);
347+
sim_sfh(_,ir) = tsfh*mass;
348+
}
335349

336-
for (uint_t it : range(t)) {
337-
vec1f tsfh = sim_sfh(it,_);
338-
sfh(it,1) = inplace_median(tsfh);
339-
for (uint_t ic : range(input.conf_interval)) {
340-
sfh(it,2+ic) = inplace_percentile(tsfh, input.conf_interval[ic]);
350+
for (uint_t it : range(t)) {
351+
vec1f tsfh = sim_sfh.safe(it,_);
352+
sfh.safe(it,1) = inplace_median(tsfh);
353+
for (uint_t ic : range(input.conf_interval)) {
354+
sfh.safe(it,2+ic) = inplace_percentile(tsfh, input.conf_interval[ic]);
355+
}
341356
}
357+
} else {
358+
// From chi2 grid
359+
360+
fitter.iterate_best_chi2(is, [&](uint_t id, const vec1f& p, float chi2) {
361+
if (chi2 - output.best_chi2.safe[is] > input.delta_chi2.back()) return;
362+
363+
vec1d tsfh;
364+
if (!gridder.get_sfh(id, t, tsfh)) {
365+
return;
366+
}
367+
368+
float mass = p[prop_id::mass];
369+
tsfh *= mass;
370+
371+
for (uint_t ic : range(input.conf_interval)) {
372+
if (chi2 - output.best_chi2.safe[is] < abs(input.delta_chi2[ic])) {
373+
for (uint_t it : range(t)) {
374+
double v = tsfh.safe[it];
375+
float& saved = sfh.safe(it,1+ic);
376+
if (input.delta_chi2[ic] < 0.0) {
377+
if (saved > v || !is_finite(saved)) {
378+
saved = v;
379+
}
380+
} else {
381+
if (saved < v || !is_finite(saved)) {
382+
saved = v;
383+
}
384+
}
385+
}
386+
}
387+
}
388+
});
342389
}
343390
}
344391

@@ -363,7 +410,7 @@ void write_sfhs(const options_t& opts, const input_state_t& input, const gridder
363410
}
364411

365412
void write_output(const options_t& opts, const input_state_t& input, const gridder_t& gridder,
366-
const output_state_t& output) {
413+
const fitter_t& fitter, const output_state_t& output) {
367414

368415
write_catalog(opts, input, gridder, output);
369416

@@ -372,7 +419,7 @@ void write_output(const options_t& opts, const input_state_t& input, const gridd
372419
}
373420

374421
if (opts.best_sfhs) {
375-
write_sfhs(opts, input, gridder, output);
422+
write_sfhs(opts, input, gridder, fitter, output);
376423
}
377424

378425
if (opts.verbose) {

‎src/fast++.cpp‎

Lines changed: 23 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -41,33 +41,34 @@ int vif_main(int argc, char* argv[]) {
4141
return 1;
4242
}
4343

44-
if (opts.make_seds.empty()) {
45-
if (gridder.read_from_cache && opts.parallel == parallel_choice::generators &&
46-
opts.n_thread > 0) {
47-
if (opts.verbose) {
48-
note("using cache, switched parallel execution from 'generators' to 'models'");
49-
}
50-
51-
opts.parallel = parallel_choice::models;
52-
}
53-
54-
// Initizalize the fitter
55-
fitter_t fitter(opts, input, gridder, output);
44+
if (!opts.make_seds.empty()) {
45+
// Write SEDs if asked, then terminate
46+
gridder.write_seds();
47+
return 0;
48+
}
5649

57-
// Build/read the grid and fit galaxies
58-
if (!gridder.build_and_send(fitter)) {
59-
return 1;
50+
if (gridder.read_from_cache && opts.parallel == parallel_choice::generators &&
51+
opts.n_thread > 0) {
52+
if (opts.verbose) {
53+
note("using cache, switched parallel execution from 'generators' to 'models'");
6054
}
6155

62-
// Compile results
63-
fitter.find_best_fits();
56+
opts.parallel = parallel_choice::models;
57+
}
58+
59+
// Initizalize the fitter
60+
fitter_t fitter(opts, input, gridder, output);
6461

65-
// Write output to disk
66-
write_output(opts, input, gridder, output);
67-
} else {
68-
// Write SEDs if asked
69-
gridder.write_seds();
62+
// Build/read the grid and fit galaxies
63+
if (!gridder.build_and_send(fitter)) {
64+
return 1;
7065
}
7166

67+
// Compile results
68+
fitter.find_best_fits();
69+
70+
// Write output to disk
71+
write_output(opts, input, gridder, fitter, output);
72+
7273
return 0;
7374
}

0 commit comments

Comments
 (0)