predict model and fit fail

[LydiaMak]

Hi

I have an issue I didn't have before with the same version (I believe i.e. 1.2.1.dev17+g0a5e85e). I build the observations like that

def build_obs(snr=10, ldist=10.0, **extras):
    """Build a dictionary of observational data.  

    :param snr:
        The S/N to assign to the photometry.

    :param ldist:
        The luminosity distance to assume for translating absolute magnitudes 
        into apparent magnitudes.

    :returns obs:
        A dictionary of observational data to use in the fit.
    """
    from prospect.utils.obsutils import fix_obs
    import sedpy

    # The obs dictionary, empty for now
    obs = {}

    # These are the names of the relevant filters, 
    # in the same order as the photometric data (see below)
    galex = ['galex_FUV', 'galex_NUV']
    #spitzer = ['spitzer_irac_ch'+n for n in ['1','2','3','4']]
    twomass = ['twomass_J', 'twomass_H', 'twomass_Ks']
    sdss = ['sdss_{0}0'.format(b) for b in ['u','g','r','i','z']]
    wise = ['wise_W1', 'wise_W2']
    filternames = galex + sdss +twomass +wise 
    # And here we instantiate the `Filter()` objects using methods in `sedpy`,
    # and put the resultinf list of Filter objects in the "filters" key of the `obs` dictionary
    obs["filters"] = sedpy.observate.load_filters(filternames)

    # Now we store the measured fluxes for a single object, **in the same order as "filters"**
    # In this example we use a row of absolute AB magnitudes from Johnson et al. 2013 (NGC4163)
    # We then turn them into apparent magnitudes based on the supplied `ldist` meta-parameter.
    # You could also, e.g. read from a catalog.
    # The units of the fluxes need to be maggies (Jy/3631) so we will do the conversion here too.
    M_AB = np.array([-13.05, -14.84, -17.71, -19.30, -19.98, -20.31, 
                     -20.56,-20.55,-20.71,-20.67, -19.81,-19.10]) #-
    dm = 25 + 5.0 * np.log10(ldist)
    mags = M_AB + dm
    obs["maggies"] = 10**(-0.4*mags)

    # And now we store the uncertainties (again in units of maggies)
    # In this example we are going to fudge the uncertainties based on the supplied `snr` meta-parameter.
    #obs["maggies_unc"] = (1./snr) * obs["maggies"]
    obs["maggies_unc"] = (np.log(10)*np.array([0.29,0.06,0.01,0.01,0.01,0.02,0.02,0.10,0.13,0.14,0.02,0.03])\
                          * obs["maggies"])/2.5

    # Now we need a mask, which says which flux values to consider in the likelihood.
    # IMPORTANT: the mask is *True* for values that you *want* to fit, 
    # and *False* for values you want to ignore.  Here we ignore the spitzer bands.
    #obs["phot_mask"] = np.array(['spitzer' not in f.name for f in obs["filters"]])

    # This is an array of effective wavelengths for each of the filters.  
    # It is not necessary, but it can be useful for plotting so we store it here as a convenience
    obs["phot_wave"] = np.array([f.wave_effective for f in obs["filters"]])

    # We do not have a spectrum, so we set some required elements of the obs dictionary to None.
    # (this would be a vector of vacuum wavelengths in angstroms)
    obs["wavelength"] = None
    # (this would be the spectrum in units of maggies)
    obs["spectrum"] = None
    # (spectral uncertainties are given here)
    obs['unc'] = None
    # (again, to ignore a particular wavelength set the value of the 
    #  corresponding elemnt of the mask to *False*)
    obs['mask'] = None

    obs['redshift'] = 0.0284

    # This function ensures all required keys are present in the obs dictionary,
    # adding default values if necessary
    obs = fix_obs(obs)

    return obs

run_params = {}
run_params["snr"] = 10.0
run_params["ldist"] = 124.28

obs = build_obs(**run_params)

Then I use the alpha model

from prospect.models.templates import TemplateLibrary
from prospect.models import SedModel
model_params = TemplateLibrary["alpha"]
#model_params.update(TemplateLibrary["nebular"])
model_params["zred"]["init"] = obs["redshift"]
model_params["total_mass"]["init"] = 1e10
#model_params["imf_type"]["init"] = 1

model = SedModel(model_params)
print(model.free_params)
print(model)

and

from prospect.sources import FastStepBasis
sps = FastStepBasis(zcontinuous=1)
print(sps.ssp.libraries)

Which eventually for

print(phot / obs["maggies"])

returns values of 1e-68.

I am not sure what happened. Any ideas?

[bd-j]

Hi @LydiaMak the commit referenced in the dev version number above (g0a5e85e) does not seem to exist in the propsector repo, so I'm not sure what version of the code you are working with, and I'm not sure what could have caused this new issue if the code did not change. Also v1.2, which appears to be the base version you have, is a bit older.

Given the demo(?) it looks like you are trying to run, I suggest installing the released v1.4.0 from PyPI, unless you have specific reasons not to.

python -m pip uninstall astro-prospector
python -m pip install astro-prospector==1.4
python -c "import prospect; print(prospect.__version__)"

Please let me know if issues persist after v1.4.0 is installed.

I also suggest switching from SedModel to SpecModel for the model object, as SedModel is deprecated, and will be removed in v2.0

from prospect.models import SpecModel
model = SpecModel(model_params)

[LydiaMak]

I changed both but still got the same values.

[LydiaMak]

I believe that the issue is in this lines

current_parameters = ",".join([f"{p}={v}" for p, v in zip(model.free_params, model.theta)])

because probably of the five bins for the fraction but for some reason it's not assigned correctly now.

logzsol=-0.5,dust2=0.6,z_fraction=0.8333333333333334,total_mass=0.8,duste_umin=0.75,duste_qpah=0.6666666666666666,duste_gamma=0.5,fagn=10000000000.0,agn_tau=1.0,dust_ratio=4.0,dust_index=0.001

fagn cannot be that large this should be the mass I gave.

Am I right?

What's the proper way to do it?

[LydiaMak]

Ah no. Actually I don't think this plays any role since this should be fine

spec, phot, mfrac = model.predict(model.theta, obs=obs, sps=sps)

[LydiaMak]

An additional comment is that for continuity_sfh I get a good fit. Not sure what changed for alpha.

[bd-j]

Hi @LydiaMak, have you compared the output of print(model) for both continuity_sfh (which works) and alpha (which doesn't)? Assuming the rest of the code has not changed, I would guess there's some difference in the free or fixed parameters that is causing this issue.

[LydiaMak]

I cannot really see a difference causing this issue

For alpha

<class 'prospect.models.sedmodel.SpecModel'>

Free Parameters: (name: prior) 
-----------
  logzsol: <class 'prospect.models.priors.TopHat'>(mini=-2,maxi=0.19)
  dust2: <class 'prospect.models.priors.TopHat'>(mini=0.0,maxi=4.0)
  z_fraction: <class 'prospect.models.priors.Beta'>(mini=0.0,maxi=1.0,alpha=[5 4 3 2 1],beta=[1 1 1 1 1])
  total_mass: <class 'prospect.models.priors.LogUniform'>(mini=100000000.0,maxi=1000000000000.0)
  duste_umin: <class 'prospect.models.priors.TopHat'>(mini=0.1,maxi=25)
  duste_qpah: <class 'prospect.models.priors.TopHat'>(mini=0.5,maxi=7.0)
  duste_gamma: <class 'prospect.models.priors.LogUniform'>(mini=0.001,maxi=0.15)
  fagn: <class 'prospect.models.priors.LogUniform'>(mini=1e-05,maxi=3.0)
  agn_tau: <class 'prospect.models.priors.LogUniform'>(mini=5.0,maxi=150.0)
  dust_ratio: <class 'prospect.models.priors.ClippedNormal'>(mean=1.0,sigma=0.3,mini=0.0,maxi=2.0)
  dust_index: <class 'prospect.models.priors.TopHat'>(mini=-2.0,maxi=0.5)

Fixed Parameters: (name: value [, depends_on]) 
-----------
  zred: [0.0284] 
  mass: [1.] <function zfrac_to_masses at 0x177c02820>
  sfh: [3] 
  imf_type: [2] 
  dust_type: [4] 
  agebins: [[ 0.          8.        ]
 [ 8.          8.47712125]
 [ 8.47712125  9.        ]
 [ 9.          9.47712125]
 [ 9.47712125  9.77815125]
 [ 9.77815125 10.13353891]] 
  add_dust_emission: [ True] 
  add_neb_emission: [ True] 
  add_neb_continuum: [ True] 
  nebemlineinspec: [ True] 
  gas_logz: [0.] <function stellar_logzsol at 0x141ff3a60>
  gas_logu: [-2.] 
  add_agn_dust: [False] 
  dust1: [0.] <function dustratio_to_dust1 at 0x177c023a0>

and for continuity_sfh

Free Parameters: (name: prior) 
-----------
  logzsol: <class 'prospect.models.priors.TopHat'>(mini=-2,maxi=0.19)
  dust2: <class 'prospect.models.priors.TopHat'>(mini=0.0,maxi=2.0)
  logmass: <class 'prospect.models.priors.TopHat'>(mini=7,maxi=12)
  logsfr_ratios: <class 'prospect.models.priors.StudentT'>(mean=[0. 0. 0. 0. 0.],scale=[0.3 0.3 0.3 0.3 0.3],df=[2 2 2 2 2])

Fixed Parameters: (name: value [, depends_on]) 
-----------
  zred: [0.0284] 
  mass: [1000000.] <function logsfr_ratios_to_masses at 0x17863db80>
  sfh: [3] 
  imf_type: [2] 
  dust_type: [0] 
  agebins: [[ 0.  8.]
 [ 8.  9.]
 [ 9. 10.]]

[bd-j]

Thanks for these details. The models are pretty different, with the latter model not including nenbular emission, etc.

Can you post a gist of the full code leading to the unexpected result? Also, can you let me know your fsps version (python -c "import fsps; print(fsps.__version__)")? Thanks!

[LydiaMak]

The fsps version is this fsps: 0.4.3.dev7+g8ae9b1f (same used for continuity_sfh).

The code I run after building the obs and model as described in the initial post is:

current_parameters = ",".join([f"{p}={v}" for p, v in zip(model.free_params, model.theta)])
print(current_parameters)
spec, phot, mfrac = model.predict(model.theta, obs=obs, sps=sps)
print(phot / obs["maggies"])

which returns very small values for phot. Same happens when running dynesty properly.

[LydiaMak]

Removing nebular emission fixed the issue.

[bd-j]

ok, I'm a little surprised by that, but glad it's working. I suggest upgrading to fsps v0.4.7 when you get a chance.