search

kevin218 / Eureka

Eureka! is a data reduction and analysis pipeline intended for time-series observations with JWST and HST.
https://eurekadocs.readthedocs.io/
MIT License
62 stars 48 forks source link

Plotting problem in S5 using Polynomial and Step #703

Closed sstamer23 closed 1 month ago

sstamer23 commented 2 months ago

Instrument

Light curve fitting (Stages 4-6)

What happened?

I'm using the WASP-39b ERS NIRSpec dataset, and I'm trying to fit a polynomial and step to it. I put values into the S5 fit_par file, and tried a couple different combinations of clipping the unbinned and binned time series in S4, as well as fitting for just c0 versus fitting for c0, c1, step0, and steptime0. I've tried: both clips false and fit for only c0, clip unbinned false and binned true and fit for only c0, both clips false and fit for all the parameters, clip unbinned false and binned true and fit for all parameters.

I also printed out hist_values, and it looks completely empty for some reason.

Error traceback output

Traceback (most recent call last):
  File "/Users/sarahstamer/Desktop/Grad/2024-2025/Research/WASP-39b_Eureka_Analysis/DataAnalysis/JWST/WASP-39b_Code_From_S3/run_eureka.py", line 32, in <module>
    meta = s5.fitlc(eventlabel, ecf_path=ecf_path)
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/eureka/S5_lightcurve_fitting/s5_fit.py", line 432, in fitlc
    meta, params = fit_channel(meta, time, flux, 0, flux_err,
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/eureka/S5_lightcurve_fitting/s5_fit.py", line 850, in fit_channel
    lc_model.fit(model, meta, log, fitter='emcee')
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/eureka/S5_lightcurve_fitting/lightcurve.py", line 179, in fit
    fit_model = self.fitter_func(self, model, meta, log, **kwargs)
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/eureka/S5_lightcurve_fitting/fitters.py", line 448, in emceefitter
    plots.plot_res_distr(lc, model, meta, fitter='emcee')
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/eureka/S5_lightcurve_fitting/plots_s5.py", line 662, in plot_res_distr
    n, bins, patches = plt.hist(hist_vals, alpha=0.5, color='b',
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/matplotlib/pyplot.py", line 3440, in hist
    return gca().hist(
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/matplotlib/__init__.py", line 1473, in inner
    return func(
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/matplotlib/axes/_axes.py", line 7001, in hist
    m, bins = np.histogram(x[i], bins, weights=w[i], **hist_kwargs)
  File "<__array_function__ internals>", line 200, in histogram
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/numpy/lib/histograms.py", line 780, in histogram
    bin_edges, uniform_bins = _get_bin_edges(a, bins, range, weights)
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/numpy/lib/histograms.py", line 426, in _get_bin_edges
    first_edge, last_edge = _get_outer_edges(a, range)
  File "/opt/miniconda3/envs/eureka/lib/python3.10/site-packages/numpy/lib/histograms.py", line 323, in _get_outer_edges
    raise ValueError(
ValueError: autodetected range of [nan, nan] is not finite

What operating system are you using?

MacOS Sonoma 14.6.1

What version of Python are you running?

Python 3.10.14

What Python packages do you have installed?

# packages in environment at /opt/miniconda3/envs/eureka:
#
# Name                    Version                   Build  Channel
alabaster                 1.0.0                    pypi_0    pypi
asciitree                 0.3.3                    pypi_0    pypi
asdf                      3.4.0                    pypi_0    pypi
asdf-astropy              0.6.1                    pypi_0    pypi
asdf-coordinates-schemas  0.3.0                    pypi_0    pypi
asdf-standard             1.1.1                    pypi_0    pypi
asdf-transform-schemas    0.5.0                    pypi_0    pypi
asdf-wcs-schemas          0.4.0                    pypi_0    pypi
asteval                   1.0.2                    pypi_0    pypi
astraeus                  0.3                      pypi_0    pypi
astropy                   6.1.2                    pypi_0    pypi
astropy-healpix           1.0.3                    pypi_0    pypi
astropy-iers-data         0.2024.8.19.0.32.16          pypi_0    pypi
astroquery                0.4.7                    pypi_0    pypi
astroscrappy              1.2.0                    pypi_0    pypi
asttokens                 2.4.1                    pypi_0    pypi
attrs                     24.2.0                   pypi_0    pypi
babel                     2.16.0                   pypi_0    pypi
backports-tarfile         1.2.0                    pypi_0    pypi
batman-package            2.4.9                    pypi_0    pypi
bayesicfitting            3.2.1                    pypi_0    pypi
beautifulsoup4            4.12.3                   pypi_0    pypi
bokeh                     2.4.3                    pypi_0    pypi
bottleneck                1.4.0                    pypi_0    pypi
bzip2                     1.0.8                h6c40b1e_6  
ca-certificates           2024.7.2             hecd8cb5_0  
ccdproc                   2.4.2                    pypi_0    pypi
celerite                  0.4.3                    pypi_0    pypi
certifi                   2024.7.4                 pypi_0    pypi
cftime                    1.6.4                    pypi_0    pypi
charset-normalizer        3.3.2                    pypi_0    pypi
click                     8.1.7                    pypi_0    pypi
cloudpickle               3.0.0                    pypi_0    pypi
contourpy                 1.2.1                    pypi_0    pypi
corner                    2.2.2                    pypi_0    pypi
crds                      11.18.1                  pypi_0    pypi
cycler                    0.12.1                   pypi_0    pypi
dask                      2024.8.1                 pypi_0    pypi
decorator                 5.1.1                    pypi_0    pypi
dill                      0.3.8                    pypi_0    pypi
docutils                  0.21.2                   pypi_0    pypi
drizzle                   1.15.3                   pypi_0    pypi
dynesty                   2.1.4                    pypi_0    pypi
emcee                     3.1.6                    pypi_0    pypi
eureka                    0.10                     pypi_0    pypi
exceptiongroup            1.2.2                    pypi_0    pypi
executing                 2.0.1                    pypi_0    pypi
exotic-ld                 3.2.0                    pypi_0    pypi
fasteners                 0.19                     pypi_0    pypi
filelock                  3.15.4                   pypi_0    pypi
fonttools                 4.53.1                   pypi_0    pypi
fsspec                    2024.6.1                 pypi_0    pypi
future                    1.0.0                    pypi_0    pypi
george                    0.4.3                    pypi_0    pypi
gwcs                      0.21.0                   pypi_0    pypi
h5netcdf                  1.3.0                    pypi_0    pypi
h5py                      3.11.0                   pypi_0    pypi
html5lib                  1.1                      pypi_0    pypi
idna                      3.7                      pypi_0    pypi
imageio                   2.35.1                   pypi_0    pypi
imagesize                 1.4.1                    pypi_0    pypi
importlib-metadata        8.4.0                    pypi_0    pypi
iniconfig                 2.0.0                    pypi_0    pypi
ipython                   8.26.0                   pypi_0    pypi
jaraco-classes            3.4.0                    pypi_0    pypi
jaraco-context            6.0.1                    pypi_0    pypi
jaraco-functools          4.0.2                    pypi_0    pypi
jedi                      0.19.1                   pypi_0    pypi
jinja2                    3.1.4                    pypi_0    pypi
jmespath                  1.0.1                    pypi_0    pypi
jsonschema                4.23.0                   pypi_0    pypi
jsonschema-specifications 2023.12.1                pypi_0    pypi
jwst                      1.11.4                   pypi_0    pypi
keyring                   25.3.0                   pypi_0    pypi
kiwisolver                1.4.5                    pypi_0    pypi
lazy-loader               0.4                      pypi_0    pypi
libffi                    3.4.4                hecd8cb5_1  
lmfit                     1.3.2                    pypi_0    pypi
locket                    1.0.0                    pypi_0    pypi
markupsafe                2.1.5                    pypi_0    pypi
matplotlib                3.9.2                    pypi_0    pypi
matplotlib-inline         0.1.7                    pypi_0    pypi
more-itertools            10.4.0                   pypi_0    pypi
ncurses                   6.4                  hcec6c5f_0  
netcdf4                   1.6.5                    pypi_0    pypi
networkx                  3.3                      pypi_0    pypi
numcodecs                 0.13.0                   pypi_0    pypi
numpy                     1.24.4                   pypi_0    pypi
numpydoc                  1.8.0                    pypi_0    pypi
opencv-python-headless    4.10.0.84                pypi_0    pypi
openssl                   3.0.14               h46256e1_0  
packaging                 24.1                     pypi_0    pypi
pandas                    2.2.2                    pypi_0    pypi
parsley                   1.3                      pypi_0    pypi
parso                     0.8.4                    pypi_0    pypi
partd                     1.4.2                    pypi_0    pypi
pexpect                   4.9.0                    pypi_0    pypi
photutils                 1.13.0                   pypi_0    pypi
pillow                    10.4.0                   pypi_0    pypi
pip                       24.2            py310hecd8cb5_0  
pluggy                    1.5.0                    pypi_0    pypi
poppy                     1.1.1                    pypi_0    pypi
prompt-toolkit            3.0.47                   pypi_0    pypi
psutil                    6.0.0                    pypi_0    pypi
ptyprocess                0.7.0                    pypi_0    pypi
pure-eval                 0.2.3                    pypi_0    pypi
pyerfa                    2.0.1.4                  pypi_0    pypi
pygments                  2.18.0                   pypi_0    pypi
pyparsing                 3.1.2                    pypi_0    pypi
pysynphot                 2.0.0                    pypi_0    pypi
pytest                    8.3.2                    pypi_0    pypi
python                    3.10.14              h5ee71fb_1  
python-dateutil           2.9.0.post0              pypi_0    pypi
pytz                      2024.1                   pypi_0    pypi
pyvo                      1.5.2                    pypi_0    pypi
pyyaml                    6.0.2                    pypi_0    pypi
readline                  8.2                  hca72f7f_0  
referencing               0.35.1                   pypi_0    pypi
reproject                 0.14.0                   pypi_0    pypi
requests                  2.32.3                   pypi_0    pypi
rpds-py                   0.20.0                   pypi_0    pypi
scikit-image              0.24.0                   pypi_0    pypi
scipy                     1.9.3                    pypi_0    pypi
semantic-version          2.10.0                   pypi_0    pypi
setuptools                72.1.0          py310hecd8cb5_0  
setuptools-scm            8.1.0                    pypi_0    pypi
six                       1.16.0                   pypi_0    pypi
snowballstemmer           2.2.0                    pypi_0    pypi
soupsieve                 2.6                      pypi_0    pypi
spherical-geometry        1.3.2                    pypi_0    pypi
sphinx                    8.0.2                    pypi_0    pypi
sphinxcontrib-applehelp   2.0.0                    pypi_0    pypi
sphinxcontrib-devhelp     2.0.0                    pypi_0    pypi
sphinxcontrib-htmlhelp    2.1.0                    pypi_0    pypi
sphinxcontrib-jsmath      1.0.1                    pypi_0    pypi
sphinxcontrib-qthelp      2.0.0                    pypi_0    pypi
sphinxcontrib-serializinghtml 2.0.0                    pypi_0    pypi
sqlite                    3.45.3               h6c40b1e_0  
stack-data                0.6.3                    pypi_0    pypi
stcal                     1.4.4                    pypi_0    pypi
stdatamodels              1.7.2                    pypi_0    pypi
stpipe                    0.7.0                    pypi_0    pypi
stsci-image               2.3.9                    pypi_0    pypi
stsci-imagestats          1.8.3                    pypi_0    pypi
stsci-stimage             0.2.9                    pypi_0    pypi
svo-filters               0.4.4                    pypi_0    pypi
tabulate                  0.9.0                    pypi_0    pypi
tifffile                  2024.8.10                pypi_0    pypi
tk                        8.6.14               h4d00af3_0  
tomli                     2.0.1                    pypi_0    pypi
toolz                     0.12.1                   pypi_0    pypi
tornado                   6.4.1                    pypi_0    pypi
tqdm                      4.66.5                   pypi_0    pypi
traitlets                 5.14.3                   pypi_0    pypi
tweakwcs                  0.8.8                    pypi_0    pypi
typing-extensions         4.12.2                   pypi_0    pypi
tzdata                    2024.1                   pypi_0    pypi
uncertainties             3.2.2                    pypi_0    pypi
urllib3                   2.2.2                    pypi_0    pypi
wcwidth                   0.2.13                   pypi_0    pypi
webencodings              0.5.1                    pypi_0    pypi
wheel                     0.43.0          py310hecd8cb5_0  
wiimatch                  0.3.2                    pypi_0    pypi
xarray                    2024.7.0                 pypi_0    pypi
xz                        5.4.6                h6c40b1e_1  
zarr                      2.18.2                   pypi_0    pypi
zipp                      3.20.0                   pypi_0    pypi
zlib                      1.2.13               h4b97444_1

Code of Conduct

taylorbell57 commented 2 months ago

Hi Sarah, at present I don't have enough information to troubleshoot this with you.

First, please carefully look at all of the log statements and output plots made in Stages 3-5 which might help you to troubleshoot the problem yourself. For example, are all of the error bars equal to zero for some reason, was the Stage 3 or 4 MAD value NaN for some reason, do the priors for all your Stage 5 fitted parameters make sense, etc.

If that isn't sufficient, please next look at the Solution Manual to Hands-on Session II (SSW2023_Workshop2_NIRSpec_Solutons.ipynb) which I made for the 2023 Sagan Summer Workshop. I tried to explain all the relevant steps and settings for the participants of that workshop, and I used the same NIRSpec G395H transit of WASP-39b dataset for that tutorial.

Finally, if you're convinced that you're doing everything right and that there is a bug with Eureka! itself, please attach your Stages 1-5 ECFs and EPF so that I can try to reproduce your issue myself.

sstamer23 commented 2 months ago

So I just updated Eureka and made the changes in the ecf's first with the edits that my data needed for SSW, and then went through the newest control files and added anything that I didn't have before. After the masking of the bad columns in S4, it gives this error: 

Traceback (most recent call last):
  File "/Users/sarahstamer/Desktop/Grad/2024-2025/Research/WASP-39b_Eureka_Analysis/DataAnalysis/JWST/WASP-39b_Code_From_S3/run_eureka.py", line 30, in <module>
    spec, lc, meta = s4.genlc(eventlabel, ecf_path=ecf_path)
  File "/Users/sarahstamer/Desktop/Grad/2024-2025/Research/WASP-39b_Eureka_Analysis/Eureka/src/eureka/S4_generate_lightcurves/s4_genLC.py", line 304, in genlc
    index = np.where(spec.optmask.x == w)[0][0]
IndexError: index 0 is out of bounds for axis 0 with size 0

I'm not sure how I was able to get past this before, but I haven't seen anyone get this error.

taylorbell57 commented 2 months ago

Hi @sstamer23, I still don't have much to go off of to help troubleshoot this. If you have carefully looked through all your Stage 3 figures and your Stage 4 ECF and everything seems good as far as you can tell, then please copy-paste the contents of you Stage 3 and 4 ECFs into this thread so that I can help figure out what is going on.

sstamer23 commented 2 months ago 
# Eureka! Control File for Stage 3: Data Reduction

ncpu        1           # Number of CPUs
nfiles      1           # The number of data files to analyze simultaneously
max_memory  0.5         # The maximum fraction of memory you want utilized by read-in frames (this will reduce nfiles if need be)
indep_batches   False       # Independently treat each batch of files? Strongly recommended to leave this as False unless you have a clear reason to set it to True.
suffix      calints     # Data file suffix

calibrated_spectra  False   # Set True to generate flux-calibrated spectra/photometry in mJy
                            # Set False to convert to electrons

# Subarray region of interest
ywindow         [0,28]      # Vertical axis as seen in DS9
xwindow         [0,1272]   # Horizontal axis as seen in DS9
src_pos_type    gaussian    # Determine source position when not given in header (Options: header, gaussian, weighted, max, or hst)
record_ypos     True        # Option to record the y position and width for each integration (only records if src_pos_type is gaussian)
dqmask          True        # Mask pixels with an odd entry in the DQ array
expand          1           # Super-sampling factor along cross-dispersion direction

# Outlier rejection along time axis
ff_outlier      False       # Set False to use only background region (recommended for deep transits)
                            # Set True to use full frame (works well for shallow transits/eclipses)
bg_thresh       [10,10]     # Double-iteration X-sigma threshold for outlier rejection along time axis

# Background parameters
bg_hw           7           # Half-width of exclusion region for BG subtraction (relative to source position)
bg_deg          0           # Polynomial order for column-by-column background subtraction, -1 for median of entire frame
bg_method       mean        # Options: std (Standard Deviation), median (Median Absolute Deviation), mean (Mean Absolute Deviation)
p3thresh        2.5         # X-sigma threshold for outlier rejection during background subtraction

# Spectral extraction parameters
spec_hw         5           # Half-width of aperture region for spectral extraction (relative to source position)
fittype         meddata     # Method for constructing spatial profile (Options: smooth, meddata, poly, gauss, wavelet, or wavelet2D)
median_thresh   5           # Sigma threshold when flagging outliers in median frame, when fittype=meddata and window_len > 1
window_len      1           # Smoothing window length, for median frame or when fittype = smooth or meddata (when computing median frame). Can set to 1 for no smoothing when computing median frame for fittype=meddata.
prof_deg        3           # Polynomial degree, when fittype = poly
p5thresh        10          # X-sigma threshold for outlier rejection while constructing spatial profile
p7thresh        15          # X-sigma threshold for outlier rejection during optimal spectral extraction

# G395H curvature treatment
curvature       correct     # How to manage the curved trace on the detector (Options: None, correct). Use correct for NIRSpec/G395.

# Diagnostics
isplots_S3      4           # Generate few (1), some (3), or many (5) figures (Options: 1 - 5)
nplots          5           # How many of each type of figure do you want to make per file?
vmin            0.98        # Sets the vmin of the color bar for Figure 3101.
vmax            1.02        # Sets the vmax of the color bar for Figure 3101.
time_axis       'y'         # Determines whether the time axis in Figure 3101 is along the y-axis ('y') or the x-axis ('x')
testing_S3      False       # Boolean, set True to only use last file and generate select figures
hide_plots      True       # If True, plots will automatically be closed rather than popping up
save_output     True        # Save outputs for use in S4
save_fluxdata   False       # Save FluxData outputs for debugging or use with other tools (can be quite large files)
verbose         True        # If True, more details will be printed about steps

# Project directory
topdir     /Users/sarahstamer/Desktop/Grad/2024-2025/Research/WASP-39b_Eureka_Analysis/Data/JWST-Sim/NIRSpec/WASP-39b_Data_From_S3

# Directories relative to project dir
inputdir     /calints
outputdir    /Stage3
# Eureka! Control File for Stage 4: Generate Lightcurves
# SSW Code for Spectroscopic Lightcurves

# Number of spectroscopic channels spread evenly over given wavelength range
nspecchan       50      # Number of spectroscopic channels
compute_white   True    # Also compute the white-light lightcurve
wave_min        3.9     # Minimum wavelength
wave_max        5.1     # Maximum wavelength
wave_input      None        # Full path to a txt file with pre-defined wavelength bins in two columns separated by whitespace.
allapers        False   # Run S4 on all of the apertures considered in S3? Otherwise will use newest output in the inputdir

# *BONUS* step - mask bad columns
mask_columns    [298, 363, 680, 704, 1103, 1121, 1141, 1142, 1607,] # Manually masking pixel columns by index number

# Parameters for drift correction of 1D spectra
recordDrift     False    # Set True to record drift/jitter in 1D spectra (always recorded if correctDrift is True)
correctDrift    False   # Set True to correct drift/jitter in 1D spectra (not recommended for simulated data)

# Parameters for sigma clipping
clip_unbinned   False   # Whether or not sigma clipping should be performed on the 1D time series
clip_binned     True    # Whether or not sigma clipping should be performed on the 1D time series
sigma           4       # The number of sigmas a point must be from the rolling median to be considered an outlier
box_width       5       # The width of the box-car filter (used to calculated the rolling median) in units of number of data points
maxiters        20      # The number of iterations of sigma clipping that should be performed.
boundary        fill    # Use 'fill' to extend the boundary values by the median of all data points (recommended), 'wrap' to use a periodic boundary, or 'extend' to use the first/last data points
fill_value      mask    # Either the string 'mask' to mask the outlier values (recommended), 'boxcar' to replace data with the mean from the box-car filter, or a constant float-type fill value.

# Limb-darkening parameters needed to compute exotic-ld
compute_ld      False
# Options for ExoTiC-LD
metallicity     0.1     # Metallicity of the star
teff            6000    # Effective temperature of the star in K
logg            4.0     # Surface gravity in log g
exotic_ld_direc /home/User/exotic-ld_data/ # Directory for ancillary files for exotic-ld, download from: https://zenodo.org/doi/10.5281/zenodo.6047317
exotic_ld_grid  stagger # You can choose from kurucz (or 1D), stagger (or 3D), mps1, or mps2 model grids, or custom (to use custom_si_grid below)
# custom_si_grid  /home/User/path/to/custom/stellar/intensity/profile  #Custom Stellar Intensity profile. For examples see Eureka/demos/JWST/Custom_Stellar_Intensity_Files
# exotic_ld_file  /home/User/exotic-ld_data/Custom_throughput_file.csv # Custom throughput file, for examples see Eureka/demos/JWST/Custom_throughput_files
# Options for SPAM
spam_file       /home/User/spam-ldcoeff-u1-u2.txt # Pre-generated file of SPAM limb-darkening values at high resolution

# Diagnostics
isplots_S4      3       # Generate few (1), some (3), or many (5) figures (Options: 1 - 5)
vmin            0.98
vmax            1.02
time_axis       'y'
hide_plots      True   # If True, plots will automatically be closed rather than popping up
verbose         True    # If True, more details will be printed about steps

# Project directory
topdir     /Users/sarahstamer/Desktop/Grad/2024-2025/Research/WASP-39b_Eureka_Analysis/Data/JWST-Sim/NIRSpec/WASP-39b_Data_From_S3

# Directories relative to project dir
inputdir     /Stage3
outputdir    /Stage4

The one thing I haven't been able to figure out yet is the xwindow and ywindow in S3, as the ones I have are the only ones that won't give me an IndexError.

taylorbell57 commented 2 months ago

Oh, one thing I forgot to ask earlier is whether you're working on the NIRSpec PRISM or the NIRSpec G395H observations of WASP-39b. If you're working on the PRISM data, are you only working on the relevant NRS1 files or did you accidentally include any of the NRS2 files? And if you're working on the G395H observations, did you separately run Stages 3+ on NRS1 and separately run Stages 3+ on NRS2; the G395H data falls on two separate detectors which must be treated independently for Stages 3-6. If you've accidentally been including both detectors at the same time, then you can adjust the suffix ECF parameter in Stage 3 to be nrs1_calints and output those into a NRS1-specific folder by setting the outputdir parameter to something like nrs1/Stage3 and then do the equivalent thing for NRS2.

If you're confident you've done all that correctly, can you please also upload the Fig 3101 image that would've been output during your Stage 3 analysis?

sstamer23 commented 2 months ago

I split them up into NRS1 and NRS2, here are the two output figures. fig3101-2D_LC_nrs1 fig3101-2D_LC_nrs2

kevin218 commented 2 months ago

It looks like there's an issue with the curvature correction. Can you share that plot too?

sstamer23 commented 1 month ago

fig3107_file0_Curvature_nrs2 fig3107_file1_Curvature_nrs2 fig3107_file2_Curvature_nrs2 fig3107_file0_Curvature_nrs1 fig3107_file1_Curvature_nrs1 fig3107_file2_Curvature_nrs1

kevin218 commented 1 month ago

I suggest trimming the edges where there is no flux, but that probably won't help elsewhere.

sstamer23 commented 1 month ago

I tried some window adjustments along with some other minor changes, and it seems to work now. Thanks for all the help!