can't get cigale to work with host data in repo

[caseyjlaw]

I'm playing with the cigale modeling of FRB hosts. My goal is to integrate future DSA-110 info with the repo. I was happy to see that you have already worked on integrating cigale with the FRBs repo! I have cigale installed and can import pcigale to a python session. However, it seems that many of the bands are not recognized. For example:

> f0 = frb.FRB.by_name("FRB20190714A")
> h0 = f0.grab_host()
> h0.run_cigale()
No AGN module found. Options are: fritz2006, skirtor2016.
No radio module found. Options are: radio.
Initialising the analysis module... 
The out/ directory was renamed to 2022-09-16_12:37:41_out/
Warning: Pan-STARRS_g in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_r in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_i in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_z in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_y in the input file but not to be taken into account in the fit.
Warning: VLT_FORS2_g in the input file but not to be taken into account in the fit.
Warning: VLT_FORS2_I in the input file but not to be taken into account in the fit.
Warning: WFC3_F300X in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_g_err in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_r_err in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_i_err in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_z_err in the input file but not to be taken into account in the fit.
Warning: Pan-STARRS_y_err in the input file but not to be taken into account in the fit.
Warning: VLT_FORS2_g_err in the input file but not to be taken into account in the fit.
Warning: VLT_FORS2_I_err in the input file but not to be taken into account in the fit.
Warning: WFC3_F300X_err in the input file but not to be taken into account in the fit.

Processing block 1/1...

Computing models ...
16800/16800 in 00:00:21.8 (769.7/s)

Estimating the physical properties ...
1/1 in 00:00:08.2 (0.1/s)

Block processed.

Estimating physical properties on all blocks

Computing the best fit spectra
1/1 in 00:00:01.0 (1.0/s)

Sanity check of the analysis results...
0.0% of the objects have χ²_red~0 and 0.0% χ²_red<0.5

Saving the analysis results...
Run completed!
1/1 in 00:00:01.3 (0.8/s)
No SED found for HG20190714A. No plot created.

And ultimately, an output directory and files made, but no plot is generated. I see that the pcigale-plots sed script exists, but it doesn't properly parse the output of the run_cigale method. I guess I was wondering if there was anything obviously wrong I'm doing here? I could work on adding new bands by following the cigale instructions, but the lack of an output make me wonder if something else isn't right. Any tips would be appreciated.

[SunilSimha]

Hi Casey, Sorry I didn't see this sooner. It appears you need to install the custom filters that we have in the frb/data/analysis/CIGALE folder. You can find the script to do so in frb/docs/installing.rst. Let me know if you if it works.

[caseyjlaw]

Ok, thanks. This is making sense now. I can also see that running the parse_cigale fills the "derived" property in the FRBGalaxy instance. Nice!

I'm still seeing this after running the modeling of either yours or my own cigale runs (even if they produce useful output): No SED found for HG20190714A. No plot created.

Any idea why it does not find an SED or make a plot?

[SunilSimha]

I think you need to pass the argument save_sed=True to run_cigale.

[caseyjlaw]

Ok, thanks. I see things working as expected now. I'm seeing another issue with plotting the SED, but I think that's a matplotlib versioning issue. Feel free to close this issue.

[caseyjlaw]

Actually, this versioning issue is turning into a real blocker for me. Could you confirm that the following works for you?

in python:
> f0 = frb.FRB.by_name("FRB20121102A")
> h0 = f0.grab_host()
> h0.run_cigale(save_sed=True, save_best_sed=True)
then in terminal:
> pcigale-plots sed

When I try this, it fails in a few different ways. It seems that the path to the files isn't right, but if I fix that, the configuration parser in cigale fails, etc.. I suspect I'm mixing versions, but it isn't clear what combo works for others. I'm using Python 3.8, cigale 2020.0, and matplotlib 3.3.4.

[SunilSimha]

Can you share the error message? Those versions should work.

[caseyjlaw]

Sure. Thanks. In Python (3.8):

> f0 = frb.FRB.by_name("FRB20121102A")
> h0 = f0.grab_host()
> h0.run_cigale(save_sed=True, save_best_sed=True, compare_obs_model=True)

Estimating the physical properties ...
1/1 in 00:00:09.0 (0.1/s)

Block processed.

Estimating physical properties on all blocks

Computing the best fit spectra
1/1 in 00:00:01.1 (0.9/s)

Sanity check of the analysis results...
0.0% of the objects have χ²_red~0 and 0.0% χ²_red<0.5

Saving the analysis results...
Run completed!
1/1 in 00:00:01.2 (0.8/s)

Then at the command line:

> pcigale-plots sed
1/1 in 00:00:01.1 (0.9/s)
multiprocessing.pool.RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/multiprocessing/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/multiprocessing/pool.py", line 51, in starmapstar
    return list(itertools.starmap(args[0], args[1]))
  File "/Users/claw/code/cigale/pcigale_plots/plot_types/sed.py", line 172, in _sed_worker
    ax1.loglog(wavelength_spec[wsed],
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/site-packages/matplotlib/axes/_axes.py", line 1802, in loglog
    self.set_yscale('log', **dy)
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/site-packages/matplotlib/axes/_base.py", line 73, in wrapper
    return get_method(self)(*args, **kwargs)
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/site-packages/matplotlib/axis.py", line 816, in _set_axes_scale
    ax._axis_map[name]._set_scale(value, **kwargs)
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/site-packages/matplotlib/axis.py", line 773, in _set_scale
    self._scale = mscale.scale_factory(value, self, **kwargs)
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/site-packages/matplotlib/scale.py", line 714, in scale_factory
    return scale_cls(axis, **kwargs)
TypeError: __init__() got an unexpected keyword argument 'nonposy'
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/Users/claw/anaconda3arm/envs/py38/bin/pcigale-plots", line 33, in <module>
    sys.exit(load_entry_point('pcigale', 'console_scripts', 'pcigale-plots')())
  File "/Users/claw/code/cigale/pcigale_plots/__init__.py", line 126, in main
    sed_action(config, args.type, args.nologo, args.xrange, args.yrange,
  File "/Users/claw/code/cigale/pcigale_plots/plot_types/sed.py", line 74, in sed
    pool.starmap(_sed_worker, zip(obs, mod, repeat(filters),
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/multiprocessing/pool.py", line 372, in starmap
    return self._map_async(func, iterable, starmapstar, chunksize).get()
  File "/Users/claw/anaconda3arm/envs/py38/lib/python3.8/multiprocessing/pool.py", line 771, in get
    raise self._value
TypeError: __init__() got an unexpected keyword argument 'nonposy'

This error comes up in older matplotlib issues, because they renamed this argument to the plot method. That suggests some kind of version mismatch.

[SunilSimha]

Aha! This is my fault. So a couple of months back, I noticed this issue with CIGALE. It is because that nonposy is not a valid argument for pyplot.loglog and line 180 of cigale.pcigale_plots.plot_types.sed in CIGALE v.2020.0has this. Instead of redoing my entire wrapper to accommodate the latest version of the wrapper, I took the easy way out (because barely anyone was using this code) and edited the CIGALE file to use the correct argument (nonposy->nonpositive). Now that more people are using it, I guess it is time to rewrite the wrapper. Sorry, Casey. So that's the quick fix for now. I will work on this once I have some free energy, so I'll leave this issue open until then.

[caseyjlaw]

So I did consider that. After making the edit, I ran into other problems. The first of which is (for the same plot command on same data): ValueError: 'yerr' must not contain negative values

[SunilSimha]

Of the three instances of yerr in cigale.pcigale_plots.plot_types.sed, the only place I think you'd trip is line 283. I think this can only happen if you have some filter in your photometry table that has some very low negative value for the error. I would recommend setting all error columns to -99.0 for filters where you want to specify flux upper limits.

[caseyjlaw]

I did have a few values like that, so I can remedy it. However, that should not be an issue for the example posted above (using the FRB20121102A data in the FRB repo). There is only one field in photom with negative value and it is set to -99. Adding print statements localizes the error to this line in sed.py (removed the yerr arg):

ax1.errorbar(filters_wl[mask_uplim], obs_fluxes[mask_uplim],
                             ls='', marker='v',
                             label='Observed upper limits',
                             markerfacecolor='None', markersize=6,
                             markeredgecolor='g', capsize=0.)

After removing the yerr argument there, there is one more below at this line:

            ax2.errorbar(filters_wl[mask],
                         (obs_fluxes[mask]-mod_fluxes[mask])/obs_fluxes[mask],
#                         yerr=obs_fluxes_err[mask]/obs_fluxes[mask],
                         marker='_', label="(Obs-Mod)/Obs", color='k',
                         capsize=0., ls='None', lw=1)

In both cases, it works if I remove the yerr line. It really seems like something isn't right with how it parses the photometry within the repo. I don't see why it fails, because the mask should select the positive values. Does this test pass for you?

[SunilSimha]

It does. Which branch are you on?

[caseyjlaw]

I forked and merged your latest into my main branch as of Sept 16 (commit 24f3f0b5d8b11685c9d2ae627092f262be79d3d4). I wonder if it is a pandas or astropy version issue? At some point, I saw an error that implied it was parsing ints as strings in the configuration file. All that said, I may just move on to using prospector for this work, unless you think cigale is the best way to go.

[SunilSimha]

I am baffled as to why it's failing on the main branch. It might be an issue with some dependency, as you suggest. I switched to the latest commit on the main branch (ed407d6), and the test works for me. @profxj, can you try running this on your system to see if it's a dependency issue?

I think prospector is better, especially if 1) you're only analyzing a single galaxy, 2) if you're running a non-parametric SFH mode and 3) are simultaneously fitting a spectrum. CIGALE is much faster, however, and I think is better suited for batches of galaxies. @muethela is our prospector expert and can comment on this.