sg_proscode.py 文件源码

def lnprobfn(theta, model, obs, sps, verbose=False, spec_noise=None,
             phot_noise=None, residuals=False):
    """Define the likelihood function.
    Given a parameter vector and a dictionary of observational data and a model
    object, return the ln of the posterior. This requires that an sps object
    (and if using spectra and gaussian processes, a GP object) be instantiated.
    """
    from prospect.likelihood import lnlike_spec, lnlike_phot, chi_spec, chi_phot

    # Calculate the prior probability and exit if not within the prior
    lnp_prior = model.prior_product(theta)
    if not np.isfinite(lnp_prior):
        return -np.infty

    # Generate the mean model--
    t1 = time()
    try:
        model_spec, model_phot, model_extras = model.mean_model(theta, obs, sps=sps)
    except(ValueError):
        return -np.infty
    d1 = time() - t1

    # Return chi vectors for least-squares optimization--
    if residuals:
        chispec = chi_spec(model_spec, obs)
        chiphot = chi_phot(model_phot, obs)        
        return np.concatenate([chispec, chiphot])

    # Noise modeling--
    if spec_noise:
        spec_noise.update(**model.params)
    if phot_noise:
        phot_noise.update(**model.params)

    vectors = {
        'spec': model_spec,    # model spectrum
        'phot': model_phot,    # model photometry
        'sed': model._spec,    # object spectrum
        'cal': model._speccal, # object calibration spectrum
        }

    # Calculate likelihoods--
    t2 = time()
    lnp_spec = lnlike_spec(model_spec, obs=obs, spec_noise=spec_noise, **vectors)
    lnp_phot = lnlike_phot(model_phot, obs=obs, phot_noise=phot_noise, **vectors)
    d2 = time() - t2
    if False:
        from prospect.likelihood import write_log
        write_log(theta, lnp_prior, lnp_spec, lnp_phot, d1, d2)

    #if (lnp_prior + lnp_phot + lnp_spec) > 0:
    #    print('Total probability is positive!!', lnp_prior, lnp_phot)
    #    pdb.set_trace()

    return lnp_prior + lnp_phot + lnp_spec