VTDA-Group / superphot-plus

Superphot+, upgraded version of Superphot that uses nested sampling
MIT License
10 stars 0 forks source link

Lightcurve.find_max_flux could output negative values #159

Open hombit opened 1 year ago

hombit commented 1 year ago

Testing new samplers I've found that in some cases (probably when S/N < 1, so not scientifically important cases) Lightcurve.find_max_flux() could output negative maximum fluxes. It would cause unphysical scaling of parameters, such as amplitude and extra error.

jeremykubica commented 1 year ago

Happy to make the code change for this, but Kaylee should make the call on what we want the behavior to be.

kdesoto-astro commented 1 year ago

There's probably a more nuanced way to handle this, but the possible options I can think of right now are: (1) Change max_flux to include all band info (2) Either have an exception where it uses the flux itself when S/N < 1 (instead of F - Ferr), or always use the flux value itself. The downside of the latter is sometimes there's outlier points with abnormally high flux but also high uncertainty. (3) Exit the sampling with an exception if S/N < 1 in the reference band.

I also want to note that for LSST the SNR values tend to be much lower than ZTF, so option 3 may skip lots of light curves. Thoughts on best option?

hombit commented 1 year ago

I think option 3) has a really good sense from physical point of view. However, it would require fixing tests/benchmarks using randomly-generated data to handle this new exception

camposandro commented 1 year ago

I've stumbled upon an error generating the posteriors for the MOSFiT data and it seems to be related to this issue. I was testing generation for the LiCu sampler (algorithm "mcmc-ceres").

$ PYTHONUNBUFFERED=1;RUST_BACKTRACE=full python scripts/generate_mosfit.py --sampler licu-mcmc-ceres 2> out.txt

  0%|          | 0/3332 [00:00<?, ?it/s]
Worker 0:   0%|          | 0/3332 [00:00<?, ?it/s]
Worker 0:   0%|          | 1/3332 [00:00<16:58,  3.27it/s]
Worker 0:   0%|          | 2/3332 [00:00<17:34,  3.16it/s]
Worker 0:   0%|          | 3/3332 [00:00<16:01,  3.46it/s]
Worker 0:   0%|          | 4/3332 [00:01<15:23,  3.60it/s]
Worker 0:   0%|          | 5/3332 [00:01<14:20,  3.87it/s]thread '<unnamed>' panicked at 'Left boundary is larger than right one: -1.2309942 > -30.921175', /Users/hombit/.cargo/registry/src/index.crates.io-6f17d22bba15001f/light-curve-feature-0.6.0/src/nl_fit/mcmc.rs:192:21
stack backtrace:
   0:        0x151cfb0f4 - <std::sys_common::backtrace::_print::DisplayBacktrace as core::fmt::Display>::fmt::h56d7672c82815b65
   1:        0x151c99338 - core::fmt::write::h632b3cc8e66b4f04
   2:        0x151cdf5bc - std::io::Write::write_fmt::hc28441b249a4971d
   3:        0x151cfeb08 - std::sys_common::backtrace::print::h8057ced4b0f9fdc1
   4:        0x151cfe72c - std::panicking::default_hook::{{closure}}::h8157fa8f0f7934b5
   5:        0x151cff62c - std::panicking::rust_panic_with_hook::h21091a3c79c5da9c
   6:        0x151cff1a4 - std::panicking::begin_panic_handler::{{closure}}::hd2f65b958d3068b8
   7:        0x151cff118 - std::sys_common::backtrace::__rust_end_short_backtrace::h53ec33e49ec66621
   8:        0x151cff10c - _rust_begin_unwind
   9:        0x151d0f230 - core::panicking::panic_fmt::h3bbf9265d206434c
  10:        0x151af8ae0 - light_curve_feature::nl_fit::mcmc::generate_initial_guesses::h92be647a63223a83
  11:        0x151ad7910 - <light_curve_feature::nl_fit::curve_fit::CurveFitAlgorithm as light_curve_feature::nl_fit::curve_fit::CurveFitTrait>::curve_fit::hdfb6ef9e4f836412
  12:        0x151ad3168 - <light_curve_feature::features::villar_fit::VillarFit as light_curve_feature::evaluator::FeatureEvaluator<T>>::eval::ha6788e0f7343e014
  13:        0x151c45160 - light_curve::features::PyFeatureEvaluator::call_impl::h0c092e83c5bfaf1d
  14:        0x151c5bbc8 - light_curve::features::PyFeatureEvaluator::__call__::h5f252bc2229f0a18
  15:        0x151c7a124 - light_curve::features::_::_::__INVENTORY::trampoline::ha55682a2637ead45
  16:        0x10033baec - __PyObject_MakeTpCall
  17:        0x10042c6c8 - _call_function
  18:        0x1004284a4 - __PyEval_EvalFrameDefault
  19:        0x10042150c - __PyEval_Vector
  20:        0x10033c330 - _PyVectorcall_Call
  21:        0x100428740 - __PyEval_EvalFrameDefault
  22:        0x10042150c - __PyEval_Vector
  23:        0x10033f1ac - _method_vectorcall
  24:        0x10033c330 - _PyVectorcall_Call
  25:        0x100428740 - __PyEval_EvalFrameDefault
  26:        0x10042150c - __PyEval_Vector
  27:        0x10033f1ac - _method_vectorcall
  28:        0x10033c330 - _PyVectorcall_Call
  29:        0x100428740 - __PyEval_EvalFrameDefault
  30:        0x10042150c - __PyEval_Vector
  31:        0x100428740 - __PyEval_EvalFrameDefault
  32:        0x10042150c - __PyEval_Vector
  33:        0x10042c630 - _call_function
  34:        0x100427f28 - __PyEval_EvalFrameDefault
  35:        0x10042150c - __PyEval_Vector
  36:        0x10042c630 - _call_function
  37:        0x100427f28 - __PyEval_EvalFrameDefault
  38:        0x10042150c - __PyEval_Vector
  39:        0x10042c630 - _call_function
  40:        0x1004284a4 - __PyEval_EvalFrameDefault
  41:        0x10042150c - __PyEval_Vector
  42:        0x10042c630 - _call_function
  43:        0x100428518 - __PyEval_EvalFrameDefault
  44:        0x10042150c - __PyEval_Vector
  45:        0x10047c444 - _run_mod
  46:        0x10047f018 - _PyRun_SimpleStringFlags
  47:        0x10049f6d8 - _Py_RunMain
  48:        0x1004a0b88 - _pymain_main
  49:        0x1002e5e4c - _main
hombit commented 1 year ago

Thank you @camposandro. This part of the message shows that priors are ill-defined, most probably because of the negative max flux issue we are talking about here

Left boundary is larger than right one: -1.2309942 > -30.921175