Do not output results for bad inputs (incl zero uncertainty)

Author: cschreib

src/fast++-fitter.cpp

@@ -227,6 +227,12 @@ fitter_t::fitter_t(const options_t& opt, const input_state_t& inp, const gridder
 
             for (uint_t is : range(input.id)) {
                 chi2_filename[is] = odir+file::get_basename(opts.catalog)+"_"+input.id[is]+".chi2.grid";
+
+                if (!input.good[is]) {
+                    file::remove(chi2_filename[is]);
+                    continue;
+                }
+
                 std::ofstream out;
                 out.open(chi2_filename[is], std::ios::binary | std::ios::out | std::ios::trunc);
                 if (!out.is_open()) {
@@ -352,6 +358,8 @@ void fitter_t::write_chi2(uint_t igrid, const vec1f& chi2, const vec2f& props, u
         for (uint_t cis : range(chi2)) {
             uint_t is = cis + i0;
 
+            if (!input.good.safe[is]) continue;
+
             if (chi2.safe[cis] > best_chi2.safe[is] + opts.save_bestchi) continue;
 
             if (chi2.safe[cis] < best_chi2.safe[is]) {
@@ -477,6 +485,11 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
     for (uint_t i : range(i1-i0)) {
         uint_t is = i + i0;
 
+        if (!input.good.safe[is]) {
+            // Skip this galaxy
+            continue;
+        }
+
         // Apply constraints on redshift
         bool dofit = (idzl.safe[is] <= iz && iz <= idzu.safe[is]);
         bool keepfit = true;
@@ -729,6 +742,11 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
 
         for (uint_t i : range(i1-i0)) {
             uint_t is = i + i0;
+            if (!input.good.safe[is]) {
+                // Skip this galaxy
+                continue;
+            }
+
             ++output.num_models[is];
             if (output.best_chi2.safe[is]  > wsp.chi2[i]) {
                 output.best_chi2.safe[is]  = wsp.chi2[i];
@@ -843,16 +861,34 @@ void fitter_t::find_best_fits() {
 
     if (opts.verbose) note("finding best fits...");
 
+    std::string best_fits_odir = opts.output_dir+"best_fits/";
+    bool save_sim = opts.save_sim;
+    if (save_sim) {
+        if (!file::mkdir(best_fits_odir)) {
+            warning("could not save simulations");
+            warning("the output directory '", best_fits_odir, "' could not be created");
+            save_sim = false;
+        }
+    }
+
     for (uint_t is : range(input.id)) {
+        std::string best_fits_output_file =
+            best_fits_odir+file::get_basename(opts.catalog)+"_"+input.id[is]+".sims.fits";
+
         if (!is_finite(output.best_chi2[is])) {
-            if (!silence_invalid_chi2) {
+            if (input.good[is] && !silence_invalid_chi2) {
                 warning("galaxy ", input.id[is], " has no best fit solution");
                 warning("there is probably a problem with the models, please re-run FAST++ with DEBUG=1");
                 warning("(further occurences of this warning for other galaxies will be suppressed)");
                 silence_invalid_chi2 = true;
             }
 
             output.best_params(is,_,_) = dnan;
+
+            if (save_sim) {
+                file::remove(best_fits_output_file);
+            }
+
             continue;
         }
 
@@ -881,17 +917,11 @@ void fitter_t::find_best_fits() {
                 bparams.safe(gridder.nparam+ip,_) = output.mc_best_props.safe(is,ip,_);
             }
 
-            if (opts.save_sim) {
-                std::string odir = opts.output_dir+"best_fits/";
-                if (!file::mkdir(odir)) {
-                    warning("could not save simulations");
-                    warning("the output directory '", odir, "' could not be created");
-                } else {
-                    fits::write_table(odir+file::get_basename(opts.catalog)+"_"+input.id[is]+".sims.fits",
-                        "result", bparams, "params", output.param_names,
-                        "chi2", output.mc_best_chi2(is,_)
-                    );
-                }
+            if (save_sim) {
+                fits::write_table(best_fits_output_file,
+                    "result", bparams, "params", output.param_names,
+                    "chi2", output.mc_best_chi2(is,_)
+                );
             }
 
             if (!opts.interval_from_chi2) {

src/fast++-read_input.cpp

@@ -980,6 +980,7 @@ bool read_fluxes(const options_t& opts, input_state_t& state) {
     state.zspec = replicate(fnan, ngal);
     state.flux = replicate(fnan, ngal, state.no_filt.size());
     state.eflux = replicate(fnan, ngal, state.no_filt.size());
+    state.good = replicate(true, ngal);
 
     // Now read the catalog itself, only keeping the columns we are interested in
     uint_t l = 0;
@@ -1053,12 +1054,20 @@ bool read_fluxes(const options_t& opts, input_state_t& state) {
             err *= totcor;
         }
 
+        // Check for zero uncertainty (ambiguous; could mean no data or no uncertainty)
+        for (uint_t idz : where(err == 0)) {
+            state.good.safe[gid] = false;
+            warning("object ", state.id.safe[gid], " (l.", l, ") has zero uncertainty for ",
+                header[col_eflux[idz]], "; will not fit");
+        }
+
         // Flag bad values
         vec1u idb = where(err < 0 || !is_finite(flx) || !is_finite(err));
         err.safe[idb] = finf; flx.safe[idb] = 0;
 
         if (idb.size() == flx.size()) {
-            warning("object ", state.id.safe[gid], " (l.", l, ") has no valid photometry");
+            state.good.safe[gid] = false;
+            warning("object ", state.id.safe[gid], " (l.", l, ") has no valid photometry; will not fit");
         }
 
         // Save flux and uncertainties in the input state

src/fast++-write_output.cpp

@@ -251,6 +251,16 @@ void write_best_fits(const options_t& opts, const input_state_t& input, const gr
     vec1f lam, sed, sed_nodust, sed_lsf, flx;
     auto pg = progress_start(input.id.size());
     for (uint_t is : range(input.id)) {
+        std::string output_file_tpl = odir+file::get_basename(opts.catalog)+"_"+input.id[is]+".fit";
+        std::string output_file_flx = odir+file::get_basename(opts.catalog)+"_"+input.id[is]+".input_res.fit";
+
+        if (!input.good[is]) {
+            file::remove(output_file_tpl);
+            file::remove(output_file_flx);
+            if (opts.verbose) progress(pg, 13);
+            continue;
+        }
+
         float scale = finf;
         uint_t model = npos;
 
@@ -303,7 +313,7 @@ void write_best_fits(const options_t& opts, const input_state_t& input, const gr
         }
 
         // Save model
-        std::ofstream fout(odir+file::get_basename(opts.catalog)+"_"+input.id[is]+".fit");
+        std::ofstream fout(output_file_tpl);
         fout << "# wl fl";
         if (opts.intrinsic_best_fit) {
             fout << " fl_nodust";
@@ -328,7 +338,7 @@ void write_best_fits(const options_t& opts, const input_state_t& input, const gr
         fout.close();
 
         // Save fluxes
-        fout.open(odir+file::get_basename(opts.catalog)+"_"+input.id[is]+".input_res.fit");
+        fout.open(output_file_flx);
         fout << "# wl fl_model fl_obs unc_obs (x 10^-19 ergs s^-1 cm^-2 Angstrom^-1)\n";
         for (uint_t il : range(input.lambda)) {
             fout << std::setw(13) << float(input.lambda[il])
@@ -386,6 +396,14 @@ void write_sfhs(const options_t& opts, const input_state_t& input, const gridder
 
     auto pg = progress_start(input.id.size());
     for (uint_t is : range(input.id)) {
+        std::string output_file = odir+opts.catalog+"_"+input.id[is]+".sfh";
+
+        if (!input.good[is]) {
+            file::remove(output_file);
+            if (opts.verbose) progress(pg, 13);
+            continue;
+        }
+
         vec1u idm = gridder.grid_ids(output.best_model[is]);
         vec1d t = rgen_step(1e6, e10(gridder.auniv[idm[grid_id::z]]), opts.sfh_output_step);
         vec2f sfh = replicate(fnan, t.size(), 1+nconf);
@@ -459,7 +477,7 @@ void write_sfhs(const options_t& opts, const input_state_t& input, const gridder
         }
 
         // Save SFH
-        std::ofstream fout(odir+opts.catalog+"_"+input.id[is]+".sfh");
+        std::ofstream fout(output_file);
         fout << header;
 
         for (uint_t it : range(t)) {

src/fast++.hpp

@@ -220,6 +220,7 @@ struct input_state_t {
     vec1b lir_log;                      // [ngal]
     vec1u lir_comp;                     // [ngal]
     vec2f flux, eflux;                  // [ngal,nfilt+nspec]
+    vec1b good;                         // [ngal]
     vec1f conf_interval;                // [nconf]
     vec1f delta_chi2;                   // [nconf]