Python sncosmo.fit_lc函数代码示例

    def test_wrong_param_names(self):
        """Supplying parameter names that are not part of the model should
        raise an error."""

        # a parameter not in the model
        with pytest.raises(ValueError):
            res, fitmodel = sncosmo.fit_lc(self.data, self.model, ["t0", "not_a_param"])

        # no parameters
        with pytest.raises(ValueError):
            res, fitmodel = sncosmo.fit_lc(self.data, self.model, [])

	def fit_lc(self, taxis, band):	
		
		data = simlc().reals(taxis, band)[0]
		mod = self.model
		res, fitted_model = sncosmo.fit_lc(data, mod, ['t0'])
		
		return res, fitted_model

def fit_curve(hml_dat, source, hml):
    # Retrieving model from source, (requires internet access)
    model = sncosmo.Model(source=source)
    # Fitting model: model parameter z bound.
    res, fitted_model = sncosmo.fit_lc(hml_dat, model, ['z', 't0', 'amplitude'], bounds={'z': (0.005, 0.35)})

    # The following excludes data sets with fewer than 4 data points and not enough distribution..

    # print 'peak', float(res.parameters[1])
    j = 0
    k = 0
    for i in hml[:, 1]:
        if float(i) > float(res.parameters[1]):
            j += 1
        else:
            k += 1

    if j >= 2 and k >= 2:
        print len(hml[:, 1])
        if len(hml[:, 1]) > 3:
            # print res.errors
            # print model

            sncosmo.plot_lc(hml_dat, model=fitted_model, errors=res.errors, color='blue', figtext=str(
                hml[:, 0][0] + ' Type' + types[hml[:, 0][0]] + '\n' + 'Model name: ' + source.name), xfigsize=10)
            plt.show()
            print 'Done:', hml[:, 0][0], 'z:', float(
                res.parameters[0]), 'Reduced chi^2:', res.chisq / res.ndof, 'Data points:', len(
                hml[:, 1]), '\n'  # 'Dof:', res.ndof

        else:
            pass

    def fit_all(self, modelcov=False):
        """
        """
        if modelcov:
            # TODO: Implement fitting with modelcov
            # from hacked_sncosmo_fitting import fit_lc_hacked
            raise NotImplementedError("Model covariance coming soon.")

        param0 = self._get_current_model_param_()
        self._derived_properties["raw_fit"] = []
        self._derived_properties["idx_good"] = []

        for k, lc in enumerate(self.lightcurves):
            for pname in [name for name in self.model.param_names if name not in self.fit_param and name != "mwr_v"]:
                if pname not in lc.meta.keys():
                    raise KeyError("Parameter '%s' not in lightcurve meta dict" % pname)
                self.model.set(**{pname: lc.meta[pname]})

            try:
                if modelcov:
                    res, _ = fit_lc_hacked(lc, self.model, self.fit_param)
                else:
                    res, _ = sncosmo.fit_lc(lc, self.model, self.fit_param)

                if res["covariance"] is not None:
                    self._derived_properties["raw_fit"].append(res)
                    self._derived_properties["idx_good"].append(k)
                else:
                    print "Light curve fit #%i failed" % k
            except sncosmo.fitting.DataQualityError:
                print "Light curve fit #%i failed because of data quality" % k
            except RuntimeError:
                print "Light curve fit #%i failed to converge" % k

        self.model.set(**param0)

def fitcurve(x, source_name, hml):
    try:
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 
        
        #adding zpsys and filter columns
        ab=np.zeros(len(hml[:,1]), dtype='|S2')
        for i in range(len(ab)):
            ab[i]='ab'
        hml=np.column_stack((hml,ab))   
        band=np.zeros(len(hml[:,1]), dtype='|S6')
        for i in range(len(band)):
            band[i]='ptf48r'
        hml=np.column_stack((hml,band))
        #fit to model
        z0 = float(spec_z(x[:-7]))
        hml_dat=astropy.table.Table(data=hml, names=('ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys', 'filter'), dtype=('str','float','float','float','float','float','float','float','float','float','str', 'str'))
        res, fitted_model=sncosmo.fit_lc(hml_dat, model, ['z','t0','amplitude'], bounds={'z':(z0,z0 + 0.0001)}, nburn=10000, nsamples=50000)
        
        #The following excludes data points not in the range of the model and data sets with fewer than 4 data points
        limit = modellimit(source_name, x[:-7], res.parameters[1])
        hml2 = []
        for j in range(len(hml[:,1])):
            datapoint = hml [:,1][j]
            if (res.parameters[1]- limit[0])< float(datapoint) < (res.parameters[1]+limit[1]):
                hml2.append(hml[j]) 
        hml2 = np.array(hml2)
        if len(hml2)>3:
            
            return finalfitcurve(x, source_name, hml2) 
               
        
    except ValueError:
        print 'error'  

def get_summary_stats(z, lc, method="emcee", convert_x0_to_mb=True):
    model = sncosmo.Model(source="salt2-extended")
    model.set(z=z)
    if method == "emcee":
        res, fitted_model = sncosmo.mcmc_lc(lc, model, ["t0", "x0", "x1", "c"])
    elif method == "minuit":
        res, fitted_model = sncosmo.fit_lc(lc, model, ["t0", "x0", "x1", "c"])
    else:
        raise ValueError("Method %s not recognised" % method)
    parameters = res.parameters[2:]

    if convert_x0_to_mb:
        determined_parameters = {k: v for k, v in zip(res.param_names, res.parameters)}
        model.set(**determined_parameters)
        mb = fitted_model.source.peakmag("bessellb", "ab")
        parameters = np.array([mb, parameters[1], parameters[2]])
        x0, x1, c = 1, 2, 3
        sigma_mb2 = 5 * np.sqrt(res.covariance[x0, x0]) / (2 * x0 * np.log(10))
        sigma_mbx1 = -5 * res.covariance[x0, x1] / (2 * x0 * np.log(10))
        sigma_mbc = -5 * res.covariance[x0, c] / (2 * x0 * np.log(10))
        cov = res.covariance
        cov = np.array(
            [
                [sigma_mb2, sigma_mbx1, sigma_mbc],
                [sigma_mbx1, cov[x1, x1], cov[x1, c]],
                [sigma_mbc, cov[x1, c], cov[c, c]],
            ]
        )
    else:
        cov = res.covariance[1:, :][:, 1:]
    return parameters, cov

def get_gaussian_fit(z, t0, x0, x1, c, lc, seed, temp_dir, interped, type="iminuit"):
    model = sncosmo.Model(source='salt2-extended')
    p = {'z': z, 't0': t0, 'x0': x0, 'x1': x1, 'c': c}
    model.set(**p)

    correct_model = sncosmo.Model(source='salt2-extended')
    correct_model.set(**p)
    if type == "iminuit":
        res, fitted_model = sncosmo.fit_lc(lc, model, ['t0', 'x0', 'x1', 'c'],
                                           guess_amplitude=False, guess_t0=False)
        chain = np.random.multivariate_normal(res.parameters[1:], res.covariance, size=int(1e5))
        fig = sncosmo.plot_lc(lc, model=[fitted_model, correct_model], errors=res.errors)
        fig.savefig(temp_dir + os.sep + "lc_%d.png" % seed, bbox_inches="tight", dpi=300)
    elif type == "mcmc":
        res, fitted_model = sncosmo.mcmc_lc(lc, model, ['t0', 'x0', 'x1', 'c'], nburn=500, nwalkers=20,
                                            nsamples=1500, guess_amplitude=False, guess_t0=False)
        chain = np.random.multivariate_normal(res.parameters[1:], res.covariance, size=int(1e5))
    elif type == "nestle":
        bounds = {"t0": [980, 1020], "x0": [0.1e-6, 9e-3], "x1": [-10, 10], "c": [-1, 1]}
        res, fitted_model = sncosmo.nest_lc(lc, model, ['t0', 'x0', 'x1', 'c'], bounds,
                                            guess_amplitude=False, guess_t0=False)
        chain = np.random.multivariate_normal(res.parameters[1:], res.covariance, size=int(1e5))
    else:
        raise ValueError("Type %s not recognised" % type)

    map = {"x0": "$x_0$", "x1": "$x_1$", "c": "$c$", "t0": "$t_0$"}
    parameters = [map[a] for a in res.vparam_names]

    chain, parameters = add_mu_to_chain(interped, chain, parameters)
    return chain, parameters, res.parameters[1:], res.covariance

def showcurve(sn, source_name):
    try:
        hml=np.load('/home/fcm1g13/Documents/Supernova/CorecollapseLC/ccdata/' + sn)    
    
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 
    # print len(hml[:,1]), 'data points'
    
        # Adding zpsys and filter columns. zp system used is ab and filter used is ptf48r
        ab = np.zeros(len(hml[:, 1]), dtype='|S2')
        band = np.zeros(len(hml[:, 1]), dtype='|S6')
        for i in range(len(ab)):
            ab[i] = 'ab'
            band[i] = 'ptf48r'
        hml = np.column_stack((hml, ab, band))

        hml_dat=astropy.table.Table(data=hml, names=('ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys', 'filter'), dtype=('str','float','float','float','float','float','float','float','float','float','str', 'str'))
        
        ###fitting model
        res, fitted_model=sncosmo.fit_lc(hml_dat, model, ['z','t0','amplitude'], bounds={'z':(0.005,0.35)}, nburn=10000, nsamples=50000)
            #print 'peak', float(res.parameters[1])
        print res                              
        sncosmo.plot_lc(hml_dat, model=fitted_model, errors=res.errors, color='blue', figtext=str(hml[:,0][0]+' Type'+types[hml[:,0][0]]+'\n'+'Model name: '+ source.name), xfigsize=10)
        #plt.clfig()
        plt.show()

        #print #'### Parameters ###'
        #print 'SN',str(hml[:,0][0]), 'z:',float(res.parameters[0])# float(res.errors['z']), float(res.parameters[1]), float(res.errors['t0']),float(res.parameters[2]), float(res.errors['x0']),  float(res.parameters[3]), float(res.errors['x1']), float(res.parameters[4]), float(res.errors['c']), float(hml[:,8][0]), float(hml[:,9][0])
        print 'Done:', hml[:,0][0], 'z:',float(res.parameters[0]), 'Reduced chi^2:', res.chisq/res.ndof,#'Data points:', len(hml[:,1]),' Type'+types[hml[:,0][0]]+'Model name: '+ source.name,'\n' #'Dof:', res.ndof

    except ValueError:
        sn, source_name, 'cannot be plotted'

    def get_supernova(self, z, num_obs, tmin, tmax, cosmology, alpha, beta,
                      x0_mean=-19.3, x0_sigma=0.1):
        t0 = np.random.uniform(tmin, tmax)
        ts = np.linspace(t0 - 60, t0 + 60, num_obs)

        times = np.array([[t, t + 0.1, t + 0.2] for t in ts]).flatten()
        bands = [b for t in ts for b in ['desg', 'desr', 'desi']]
        gains = np.ones(times.shape)
        skynoise = 20 * np.ones(times.shape)
        zp = 30 * np.ones(times.shape)
        zpsys = ['ab'] * times.size

        obs = Table({'time': times,
                     'band': bands,
                     'gain': gains,
                     'skynoise': skynoise,
                     'zp': zp,
                     'zpsys': zpsys})

        model = sncosmo.Model(source='salt2')

        mabs = np.random.normal(x0_mean, x0_sigma)
        model.set(z=z)
        x1 = np.random.normal(0., 1.)
        c = np.random.normal(0., 0.2)
        mabs = mabs - alpha * x1 + beta * c
        model.set_source_peakabsmag(mabs, 'bessellb', 'ab', cosmo=cosmology)

        x0 = model.get('x0')

        p = {'z': z,
             't0': t0,
             'x0': x0,
             'x1': x1,
             'c': c
             }

        lcs = sncosmo.realize_lcs(obs, model, [p])
        res, fitted_model = sncosmo.fit_lc(lcs[0], model, ['t0', 'x0', 'x1', 'c'])
        determined_parameters = {k: v for k, v in zip(res.param_names, res.parameters)}
        model.set(**determined_parameters)
        mb = fitted_model.source.peakmag("bessellb", "ab")
        x0_ind = res.vparam_names.index('x0')
        x1_ind = res.vparam_names.index('x1')
        c_ind = res.vparam_names.index('c')

        x0 = res.parameters[res.param_names.index('x0')]
        x1 = res.parameters[res.param_names.index('x1')]
        c = res.parameters[res.param_names.index('c')]
        sigma_mb2 = 5 * np.sqrt(res.covariance[x0_ind, x0_ind]) / (2 * x0 * np.log(10))
        sigma_mbx1 = -5 * res.covariance[x0_ind, x1_ind] / (2 * x0 * np.log(10))
        sigma_mbc = -5 * res.covariance[x0_ind, c_ind] / (2 * x0 * np.log(10))
        covariance = np.array([[sigma_mb2, sigma_mbx1, sigma_mbc],
                              [sigma_mbx1, res.covariance[x1_ind, x1_ind],
                               res.covariance[x1_ind, c_ind]],
                              [sigma_mbc, res.covariance[x1_ind, c_ind],
                               res.covariance[c_ind, c_ind]]])
        icov = np.linalg.inv(covariance)
        return [mb, x1, c, covariance, icov]

    def test_fit_lc(self):
        res, fitmodel = sncosmo.fit_lc(self.data, self.model,
                                       ['z', 't0', 'amplitude'],
                                       bounds={'z': (0., 1.0)})

        # set model to true parameters
        self.model.set(**self.params)
        assert_allclose(res.parameters, self.model.parameters, rtol=1.e-3)

 def fits(self):
     """
     `sncosmo.fit_lc` return which is a tuple to fit results, and sncosmo
     model set to best fit model.
     """
     return sncosmo.fit_lc(self.SNCosmoLC(scattered=True,
                                         nightlyCoadd=True, seed=0), model=self.sncosmoModel,
                     vparam_names=['t0', 'x0', 'x1', 'c'], minsnr=0.01)

    def fitmodel(self, modeldir='/Users/rodney/Dropbox/WFIRST/SALT2IR',
                 salt2subdir='salt2ir', fitbands=None):
        """ fit a new model to the data
        :param modeldir:
        :return:
        """
        salt2irmodeldir = os.path.join(modeldir, salt2subdir)
        salt2irsource = sncosmo.models.SALT2Source(
            modeldir=salt2irmodeldir, name=salt2subdir)
        salt2irmodel = sncosmo.Model(source=salt2irsource)

        # limit the fit to specific bands if requested
        if fitbands is None:
            datatofit = self.lcdata
        else:
            if isinstance(fitbands, basestring):
                fitbands = [fitbands]
            itofit = np.array([], dtype=int)
            for fitband in fitbands:
                ithisband = np.where(self.band == fitband)[0]
                itofit = np.append(itofit, ithisband)
            datatofit = self.lcdata[itofit]

        if 'salt2x1' not in self.__dict__:
            print("No SALT2 fitres data available for this SN."
                  " No predefined model")
            self.model_fixedx1c = None
            salt2irmodel.set(z=self.z, x1=0, c=0.1)
        else:
            # evaluate the SALT2 model for the given redshift, x1, c
            salt2irmodel.set(z=self.z, t0=self.TBmax,
                             x1=self.salt2x1, c=self.salt2c)
            salt2irmodel.set_source_peakabsmag(-19.3, 'bessellb','AB')

            # fit the model without allowing x1 and c to vary
            res, model_fixedx1c = sncosmo.fit_lc(
                datatofit, salt2irmodel, ['t0', 'x0'], bounds={})
            self.model_fixedx1c = model_fixedx1c

        # fit the model allowing all parameters to vary
        res, model_fitted = sncosmo.fit_lc(
            datatofit, salt2irmodel, ['t0', 'x0', 'x1', 'c'],
            bounds={'x1':[-3,3], 'c':[-0.5,0.5]})
        self.model_fitted = model_fitted

    def test_fit_lc(self):
        """Ensure that fit results match input model parameters (data are
        noise-free).

        Pass in parameter names in order different from that stored in
        model; tests parameter re-ordering."""
        res, fitmodel = sncosmo.fit_lc(self.data, self.model, ["amplitude", "z", "t0"], bounds={"z": (0.0, 1.0)})

        # set model to true parameters and compare to fit results.
        self.model.set(**self.params)
        assert_allclose(res.parameters, self.model.parameters, rtol=1.0e-3)

def fitcurve(x, source_name):
    # Get data set needed
    hml=np.load('/home/fcm1g13/Documents/Supernova/CorecollapseLC/ccdata/' + x)    
    try:
        # Import model
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 
       # print len(hml[:,1]), 'data points'
       
        # Adding zpsys and filter columns. zp system used is ab and filter used is ptf48r
        ab = np.zeros(len(hml[:, 1]), dtype='|S2')
        band = np.zeros(len(hml[:, 1]), dtype='|S6')
        for i in range(len(ab)):
            ab[i] = 'ab'
            band[i] = 'ptf48r'
        hml = np.column_stack((hml, ab, band))

        # Converting into a table using astropy with titles: ptfname, time,
        # magnitude, mag_err, flux, flux_err, zp_new, zp, ra, dec, zpsys and filter
        hml_dat = astropy.table.Table(data=hml, names=(
             'ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys',
               'filter'),
                                      dtype=(
                                          'str', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float',
                                            'float', 'str', 'str'))

        # Fitting model: model parameter z bound.
        res, fitted_model = sncosmo.fit_lc(hml_dat, model, ['z', 't0', 'amplitude'], bounds={'z': (0.005, 0.35)})

        #The following excludes data with fewer than 4 data points and not enough distribution..
        j=0; k=0
        for i in hml[:,1]:
            if float(i)>float(res.parameters[1]):
                j+=1
            else:
                k+=1
        # print hml[:,0][0],k,j
        if j>=2 and k>=2:            
            if len(hml[:,1])>3:                                
                #sncosmo.plot_lc(hml_dat, model=fitted_model, errors=res.errors, color='blue', figtext=str(hml[:,0][0]+' Type'+types[hml[:,0][0]]+'\n'+'Model name: '+ source.name), xfigsize=10)
                #plt.close()                
                
                return np.array([float(res.chisq/res.ndof), hml[:,0][0], source_name]) #returns reduced chi squared, sn name and model name. 
            else:
                pass
        else:
            pass
                
    except ValueError:
        pass

def showcurve(sn, source_name, hml):
    try:
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 

        #fit to model
        z0 = float(spec_z(sn[:-7]))
        hml_dat=astropy.table.Table(data=hml, names=('ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys', 'filter'), dtype=('str','float','float','float','float','float','float','float','float','float','str', 'str'))
        res, fitted_model=sncosmo.fit_lc(hml_dat, model, ['z','t0','amplitude'], bounds={'z':(z0, z0+0.001)}, nburn=10000, nsamples=50000)

        #plot model
        sncosmo.plot_lc(hml_dat, model=fitted_model, errors=res.errors, color='blue', figtext=str(hml[:,0][0]+' Type'+types[hml[:,0][0]]+'\n'+'Model name: '+ source.name + '\n'+'Reduced Chi Squared: ')+ str(res.chisq/res.ndof), xfigsize=10)
        plt.show()

        #print 'Parameters:''z:',float(res.parameters[0]), float(res.errors['z']), float(res.parameters[1]), float(res.errors['t0']),float(res.parameters[2]), float(res.errors['x0']),  float(res.parameters[3]), float(res.errors['x1']), float(res.parameters[4]), float(res.errors['c']), float(hml[:,8][0]), float(hml[:,9][0])'
        print 'Done:', hml[:,0][0], 'z:',float(res.parameters[0]), 'Reduced chi^2:', res.chisq/res.ndof,#'Data points:', len(hml[:,1]),' Type'+types[hml[:,0][0]]+'Model name: '+ source.name,'\n' #'Dof:', res.ndof

    except ValueError:
        print sn, source_name, 'cannot be plotted'