diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0d20b64 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.pyc diff --git a/doppler/__init__.py b/doppler/__init__.py index b81c6a8..25f191b 100644 --- a/doppler/__init__.py +++ b/doppler/__init__.py @@ -13,8 +13,8 @@ models = cannon.load_models() cannon.models = models -def read(filename=None,format=None): - return reader.read(filename=filename,format=None) +def read(filename=None,format=None,badval=None): + return reader.read(filename=filename,format=None,badval=None) def fit(*args,**kwargs): return rv.fit(*args,**kwargs) diff --git a/doppler/cannon.py b/doppler/cannon.py index 916318f..5643fc2 100644 --- a/doppler/cannon.py +++ b/doppler/cannon.py @@ -1146,7 +1146,10 @@ def prepare_cannon_model(model,spec,dointerp=False): # using 3 affects the final profile shape nbin = np.round(np.min(fwhmpix)//4).astype(int) if nbin==0: - raise Exception('Model has lower resolution than the observed spectrum') + import pdb + nbin=1 + print(spec.filename, 'Model has lower resolution than the observed spectrum',fwhmpix.min()) + #raise Exception('Model has lower resolution than the observed spectrum',spec.filename,fwhmpix.min()) if nbin>1: rmodel = rebin_cannon_model(tmodel,nbin) rmodel.rebin = True diff --git a/doppler/cannongrid_cubic_norm_training.py b/doppler/cannongrid_cubic_norm_training.py new file mode 100644 index 0000000..ffe5cc4 --- /dev/null +++ b/doppler/cannongrid_cubic_norm_training.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python + +# Imports +import os +import thecannon as tc +from astropy.io import fits +from astropy.table import Table +import numpy as np +import matplotlib.pyplot as plt +import pickle + +#tag = '3000_18000_allstars' +#tag = '3000_18000_rgb' +#tag = '3000_18000_cooldwarfs' +tag = '3000_18000_hotstars' +#tag = '3000_18000_hotdwarfs' +#tag = '3000_18000_coolstars' +#tag = '3000_18000_whitedwarfs' + + +print(tag) + +# Import the spectra and labels +normalized_flux = fits.getdata('cannongrid_'+tag+'_synth_data_norm.fits.gz') +labelled_set=Table.read('cannongrid_'+tag+'_synth_pars.fits') + +# Add wavelengths +nspec, npix = normalized_flux.shape +wave = np.arange(npix)*0.10+3000.0 +#model3.dispersion = wave +normalized_ivar = normalized_flux.copy()*0 + 1e4 +vec3 = tc.vectorizer.PolynomialVectorizer(labelled_set.colnames, 3) +model3 = tc.CannonModel(labelled_set, normalized_flux, normalized_ivar, vec3, wave) +model3.regularization = 0 # no regularization for now + +# Train the model +nr_theta, nr_s2, nr_metadata = model3.train() +if os.path.exists('cannongrid_'+tag+'_norm_cubic_model.pkl'): os.remove('cannongrid_'+tag+'_norm_cubic_model.pkl') +model3.write('cannongrid_'+tag+'_norm_cubic_model.pkl') + +# Check if there is a continuum model to add +contfile = 'cannongrid_'+tag+'_cont_cubic_model_logflux.pkl' +if os.path.exists(contfile): + print('Continuum model found. Adding it') + f = open(contfile,'rb') + cont = pickle.load(f) + f.close() +else: + cont = None + +# This adds my custom attributes +# model.write() doesn't same them +infile = open('cannongrid_'+tag+'_norm_cubic_model.pkl','rb') +temp = pickle.load(infile) +infile.close() +temp['fwhm'] = 0.001 +temp['wavevac'] = False +if cont is not None: + temp['continuum'] = cont +md = temp['metadata'] +ta = md['trained_attributes'] +if cont is not None: + ta = ta+('fwhm','wavevac','continuum') +else: + ta = ta+('fwhm','wavevac') +md['trained_attributes'] = ta +temp['metadata'] = md +os.remove('cannongrid_'+tag+'_norm_cubic_model.pkl') +outfile = open('cannongrid_'+tag+'_norm_cubic_model.pkl','wb') +pickle.dump(temp,outfile) +outfile.close() + +# Test/fit the spectra to get the labels +labels3, cov, meta = model3.test(normalized_flux, normalized_ivar) + + +# one-to-one comparisons +names = model3.vectorizer.label_names +fig = plt.figure(figsize=(15,10)) +for i in range(0, 3): + plt.subplot(2, 2, i+1) + plt.plot(labelled_set[names[i]],labels3[:,i],'+') + plt.plot(labelled_set[names[i]],labelled_set[names[i]]) + plt.xlabel('Input '+names[i]) + plt.ylabel('Recovered '+names[i]) +fig.savefig('cannongrid_'+tag+'_comparison.png') + + diff --git a/doppler/lsf.py b/doppler/lsf.py index 1b85859..2a178e6 100644 --- a/doppler/lsf.py +++ b/doppler/lsf.py @@ -118,7 +118,7 @@ def ghlsf(x,xcenter,params,nowings=False): wingparams2 = np.empty((nlsf*npix,params['nWpar']+1)) for i in range(params['nWpar']+1): wingparams2[:,i] = np.repeat(wingparams1[i,:],nlsf) - out += ghwingsbin(xlsf,wingparams2,params['binsize'],params['Wproftype']) + if not nowings : out += ghwingsbin(xlsf,wingparams2,params['binsize'],params['Wproftype']) # Reshape it to [Npix,Nlsf] out = out.reshape(npix,nlsf) @@ -528,7 +528,7 @@ def ghlsf_bovy(x,xcenter,params,nowings=False): # Calculate the GH part of the LSF out = gausshermitebin_bovy(x,ghparams,params['binsize']) # Calculate the Wing part of the LSF - out += ghwingsbin_bovy(x,wingparams,params['binsize'],params['Wproftype']) + if not nowings: out += ghwingsbin_bovy(x,wingparams,params['binsize'],params['Wproftype']) return out diff --git a/doppler/reader.py b/doppler/reader.py index eb80313..a84beff 100644 --- a/doppler/reader.py +++ b/doppler/reader.py @@ -17,6 +17,8 @@ from scipy.ndimage.filters import median_filter from dlnpyutils import utils as dln, bindata from .spec1d import Spec1D +import matplotlib.pyplot as plt +import pdb # Ignore these warnings, it's a bug warnings.filterwarnings("ignore", message="numpy.dtype size changed") @@ -24,7 +26,7 @@ # Load a spectrum -def read(filename=None,format=None): +def read(filename=None,format=None,badval=None): ''' This reads in a SDSS-IV MWM training set spectrum and returns an object that is guaranteed to have certain information. @@ -80,7 +82,7 @@ def read(filename=None,format=None): if format.lower() not in _readers.keys(): raise ValueError('reader '+format+' not found') # Use the requested reader/format - out = _readers[format](filename) + out = _readers[format](filename,badval=badval) if out is not None: return out # Loop over all readers until we get a spectrum out @@ -97,7 +99,7 @@ def read(filename=None,format=None): # Load APOGEE apVisit/asVisit spectra -def apvisit(filename): +def apvisit(filename,badval=20735): """ Read a SDSS APOGEE apVisit spectrum. @@ -119,7 +121,7 @@ def apvisit(filename): """ base, ext = os.path.splitext(os.path.basename(filename)) - + # APOGEE apVisit, visit-level spectrum if (base.find("apVisit") > -1) | (base.find("asVisit") > -1): # HISTORY AP1DVISIT: HDU0 = Header only @@ -165,17 +167,25 @@ def apvisit(filename): # 'LITTROW_GHOST','PERSIST_HIGH','PERSIST_MED','PERSIST_LOW','SIG_SKYLINE','SIG_TELLURIC','NOT_ENOUGH_PSF',''] # badflag = [1,1,1,1,1,1,1,1, # 0,0,0,0,0,0,1,0] - mask = (np.bitwise_and(spec.bitmask,16639)!=0) | (np.isfinite(spec.flux)==False) + #mask = (np.bitwise_and(spec.bitmask,16639)!=0) | (np.isfinite(spec.flux)==False) + mask = (np.bitwise_and(spec.bitmask,badval)!=0) | (np.isfinite(spec.flux)==False) # Extra masking for bright skylines - x = np.arange(spec.npix) - nsky = 4 - for i in range(spec.norder): - sky = spec.sky[:,i] - medsky = median_filter(sky,201,mode='reflect') - medcoef = dln.poly_fit(x,medsky/np.median(medsky),2) - medsky2 = dln.poly(x,medcoef)*np.median(medsky) - skymask1 = (sky>nsky*medsky2) # pixels Nsig above median sky - mask[:,i] = np.logical_or(mask[:,i],skymask1) # OR combine + # Commented out in favor of using SIG_SKYLINE in bitmask + # This can also mask too many pixels + #x = np.arange(spec.npix) + #nsky = 4 + ##plt.clf() + #for i in range(spec.norder): + # sky = spec.sky[:,i] + # medsky = median_filter(sky,201,mode='reflect') + # medcoef = dln.poly_fit(x,medsky/np.nanmedian(medsky),2) + # medsky2 = dln.poly(x,medcoef)*np.nanmedian(medsky) + # skymask1 = (sky>nsky*medsky2) # pixels Nsig above median sky + # #mask[:,i] = np.logical_or(mask[:,i],skymask1) # OR combine + # #plt.plot(spec.wave[:,i],sky) + # #plt.plot(spec.wave[:,i],nsky*medsky2) + # #plt.plot(spec.wave[:,i],spec.flux[:,i]) + ##plt.draw() spec.mask = mask # Fix NaN pixels for i in range(spec.norder): @@ -197,7 +207,7 @@ def apvisit(filename): # Load APOGEE apStar/asStar spectra -def apstar(filename): +def apstar(filename,badval=20735): """ Read an SDSS APOGEE apStar spectrum. @@ -281,7 +291,7 @@ def apstar(filename): # 'LITTROW_GHOST','PERSIST_HIGH','PERSIST_MED','PERSIST_LOW','SIG_SKYLINE','SIG_TELLURIC','NOT_ENOUGH_PSF',''] # badflag = [1,1,1,1,1,1,1,1, # 0,0,0,0,0,0,1,0] - mask = (np.bitwise_and(spec.bitmask,16639)!=0) | (np.isfinite(spec.flux)==False) + mask = (np.bitwise_and(spec.bitmask,badval)!=0) | (np.isfinite(spec.flux)==False) # Extra masking for bright skylines x = np.arange(spec.npix) nsky = 4 @@ -309,7 +319,7 @@ def apstar(filename): # Load SDSS BOSS spectra -def boss(filename): +def boss(filename,badval=0): """ Read a SDSS BOSS spectrum. @@ -377,7 +387,7 @@ def boss(filename): # Load SDSS MaStar spectra -def mastar(filename): +def mastar(filename,badval=0): """ Read a SDSS MaStar spectrum. @@ -441,7 +451,7 @@ def mastar(filename): # Load IRAF-style spectra -def iraf(filename): +def iraf(filename,badval=0): """ Read an IRAF-style spectrum. diff --git a/doppler/rv.py b/doppler/rv.py index afb3eec..e73a624 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -36,6 +36,7 @@ matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib.legend import Legend +import pdb # Ignore these warnings, it's a bug warnings.filterwarnings("ignore", message="numpy.dtype size changed") @@ -45,9 +46,9 @@ def xcorr_dtype(nlag): """Return the dtype for the xcorr structure""" - dtype = np.dtype([("xshift0",float),("ccp0",float),("xshift",float),("xshifterr",float), - ("xshift_interp",float),("ccf",(float,nlag)),("ccferr",(float,nlag)),("ccnlag",int), - ("cclag",(int,nlag)),("ccpeak",float),("ccpfwhm",float),("ccp_pars",(float,4)), + dtype = np.dtype([("xshift0",float),("ccp0",float),("vrel0",float),("xshift",float),("xshifterr",float), + ("xshift_interp",float), ("ccf",(float,nlag)),("ccferr",(float,nlag)),("ccnlag",int), + ("cclag",(int,nlag)),("ccvlag",(float,nlag)),("ccpeak",float),("ccpfwhm",float),("ccp_pars",(float,4)), ("ccp_perror",(float,4)),("ccp_polycoef",(float,4)),("vrel",float), ("vrelerr",float),("w0",float),("dw",float),("chisq",float)]) return dtype @@ -471,7 +472,7 @@ def ccorrelate(x, y, lag, yerr=None, covariance=False, double=None, nomean=False return cross -def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf=False,prior=None): +def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=[-200,200],gfilt=0,errccf=False,prior=None,plot=False): """This measures the radial velocity of a spectrum vs. a template using cross-correlation. This program measures the cross-correlation shift between @@ -526,10 +527,12 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf # Set up the cross-correlation parameters # this only gives +/-450 km/s with 2048 pixels, maybe use larger range - nlag = 2*np.round(np.abs(maxlag))+1 + #nlag = 2*np.round(np.abs(maxlag))+1 + nlag=maxlag[1]-maxlag[0]+1 if ((nlag % 2) == 0): nlag +=1 # make sure nlag is odd dlag = 1 - minlag = -np.int(np.ceil(nlag/2)) + #minlag = -np.int(np.ceil(nlag/2)) + minlag=maxlag[0] lag = np.arange(nlag)*dlag+minlag+1 # Initialize the output structure @@ -593,8 +596,8 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf ccf *= np.exp(-0.5*(((lag-prior[0])/prior[1])**4))*0.8+np.exp(-0.5*(((lag-prior[0])/150)**2))*0.2 else: # no good pixels - ccf = np.float(lag)*0.0 - if (errccf is True) | (nofit is False): ccferr=ccf + ccf = lag.astype(float)*0.0 + if (errccf is True) : ccferr=ccf # Remove the median ccf -= np.median(ccf) @@ -609,8 +612,22 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf gdmask = (np.isfinite(spec)==True) & (np.isfinite(temp)==True) & (spec>0.0) & (err>0.0) & (err < 1e5) ngdpix = np.sum(gdmask) if (ngdpix==0): - raise Exception('Bad spectrum') + raise RuntimeError('Bad spectrum') chisq = np.sqrt( np.sum( (spec[gdmask]-temp[gdmask])**2/err[gdmask]**2 )/ngdpix ) + + if plot : + print(len(spec),len(gdmask), chisq) + matplotlib.use('TkAgg') + from tools import plots + fig,ax=plots.multi(1,3,figsize=(8,6)) + for i in range(3) : + ax[0].plot(wobs[:,i],spec[:,i],color='k') + ax[0].plot(wobs[gdmask[:,i],i],spec[gdmask[:,i],i],color='g') + ax[0].plot(wobs[gdmask[:,i],i],err[gdmask[:,i],i],color='b') + ax[0].plot(wobs[gdmask[:,i],i],temp[gdmask[:,i],i],color='r') + ax[0].set_ylim(0,1.2) + ax[2].plot(wobs[gdmask[:,i],i],(spec[gdmask[:,i],i]-temp[gdmask[:,i],i])**2/err[gdmask[:,i],i]**2,color='r') + ax[2].set_ylim(0,100) outstr["chisq"] = chisq outstr["ccf"] = ccf @@ -619,8 +636,10 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf outstr["cclag"] = lag # Remove smooth background at large scales - cont = gaussian_filter1d(ccf,100) - ccf_diff = ccf-cont + if gfilt > 0 : + cont = gaussian_filter1d(ccf,filter) + ccf_diff = ccf-cont + else : ccf_diff = ccf # Get peak of CCF best_shiftind = np.argmax(ccf_diff) @@ -631,8 +650,8 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf # Some CCF peaks are SOOO wide that they span the whole width # do the first one without background subtraction estimates0 = [ccf_diff[best_shiftind0], best_xshift0, 4.0, 0.0] - lbounds0 = [1e-3, np.min(lag), 0.1, -np.inf] - ubounds0 = [np.inf, np.max(lag), np.max(lag), np.inf] + lbounds0 = [1e-6, np.min(lag), 0.1, -np.inf] + ubounds0 = [np.inf, np.max(lag), np.max(lag)-np.min(lag), np.inf] pars0, cov0 = dln.gaussfit(lag,ccf_diff,estimates0,ccferr,bounds=(lbounds0,ubounds0)) perror0 = np.sqrt(np.diag(cov0)) @@ -648,7 +667,7 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf yfit1 = dln.gaussian(lag[lo1:hi1],*pars1) perror1 = np.sqrt(np.diag(cov1)) - # Fefit and let constant vary more, keep width constrained + # Refit and let constant vary more, keep width constrained estimates2 = pars1 estimates2[1] = dln.limit(estimates2[1],np.min(lag),np.max(lag)) # must be in range estimates2[3] = np.median(ccf_diff[lo1:hi1]-yfit1) + pars1[3] @@ -683,7 +702,7 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf dpars3, dcov3 = dln.gaussfit(lag[lo3:hi3],ccf_diff[lo3:hi3],pars3,ccferr[lo3:hi3],bounds=(dlbounds3,dubounds3)) dyfit3 = dln.gaussian(lag[lo3:hi3],*pars3) perror3 = np.sqrt(np.diag(dcov3)) - + # Final parameters pars = pars3 perror = perror3 @@ -698,6 +717,7 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf #--------------------------------- # delta log(wave) = log(v/c+1) # v = (10^(delta log(wave))-1)*c + # why not v = delta ln(wave) * c? dwlog = np.median(dln.slope(np.log10(wave))) vrel = ( 10**(xshift*dwlog)-1 )*cspeed # Vrel uncertainty @@ -707,6 +727,7 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf # Make CCF structure and add to STR #------------------------------------ outstr["xshift0"] = best_xshift + outstr["vrel0"] = best_xshift*dwlog*np.log(10)*cspeed outstr["ccp0"] = np.max(ccf) outstr["xshift"] = xshift outstr["xshifterr"] = xshifterr @@ -720,6 +741,14 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,errccf outstr["vrelerr"] = vrelerr outstr["w0"] = np.min(wave) outstr["dw"] = dwlog + + tmp = np.median(dln.slope(np.log(wave)))*cspeed + if plot : + ax[1].plot(lag*tmp,ccf) + ax[1].plot(lag*tmp,ccf_diff) + ax[1].plot(lag[lo3:hi3]*tmp,yfit3) + plt.draw() + return outstr @@ -776,7 +805,7 @@ def normspec(spec=None,ncorder=6,fixbadpix=True,noerrcorr=False, # Can only do 1D or 2D arrays if spec.flux.ndim>2: - raise Exception("Flux can only be 1D or 2D arrays") + raise RuntimeError("Flux can only be 1D or 2D arrays") # Do special processing if the input is 2D # Loop over the shorter axis @@ -855,7 +884,7 @@ def normspec(spec=None,ncorder=6,fixbadpix=True,noerrcorr=False, gdbin = np.isfinite(ybin) ngdbin = np.sum(gdbin) if ngdbin<(ncorder+1): - raise Exception("Not enough good flux points to fit the continuum") + raise RuntimeError("Not enough good flux points to fit the continuum") # Fit with robust polynomial coef1 = dln.poly_fit(xbin[gdbin],ybin[gdbin],ncorder,robust=True) cont1 = dln.poly(x,coef1) @@ -881,7 +910,7 @@ def normspec(spec=None,ncorder=6,fixbadpix=True,noerrcorr=False, gdbin2 = np.isfinite(ybin2) ngdbin2 = np.sum(gdbin2) if ngdbin2<(ncorder+1): - raise Exception("Not enough good flux points to fit the continuum") + raise RuntimeError("Not enough good flux points to fit the continuum") # Fit with robust polynomial coef2 = dln.poly_fit(xbin2[gdbin2],ybin2[gdbin2],ncorder,robust=True) cont2 = dln.poly(x,coef2) @@ -1048,7 +1077,7 @@ def emcee_lnprob(theta, x, y, yerr, models, spec): return lp + emcee_lnlike(theta, x, y, yerr, models, spec) -def fit_xcorrgrid(spec,models=None,samples=None,verbose=False,maxvel=1000.0): +def fit_xcorrgrid(spec,models=None,samples=None,verbose=False,maxvel=[-1000.,1000],plot=False,usepeak=False): """ Fit spectrum using cross-correlation with models sampled in the parameter space. @@ -1111,8 +1140,9 @@ def fit_xcorrgrid(spec,models=None,samples=None,verbose=False,maxvel=1000.0): #------------------------------------------------------------------------------------------------ dwlog = np.median(dln.slope(np.log10(wavelog))) # vrel = ( 10**(xshift*dwlog)-1 )*cspeed - maxlag = np.int(np.ceil(np.log10(1+maxvel/cspeed)/dwlog)) - maxlag = np.maximum(maxlag,50) + maxlag = np.ceil(np.log10(1+np.array(maxvel)/cspeed)/dwlog).astype(int) + #maxlag = np.int(np.ceil(np.log10(1+np.array(maxvel)/cspeed)/dwlog)) + #maxlag = np.maximum(maxlag,50) if samples is None: #teff = [3500.0, 4000.0, 5000.0, 6000.0, 7500.0, 9000.0, 15000.0, 25000.0, 40000.0, 3500.0, 4300.0, 4700.0, 5200.0] #logg = [4.8, 4.8, 4.6, 4.4, 4.0, 4.0, 4.0, 4.0, 8.0, 0.5, 1.0, 2.0, 3.0] @@ -1125,34 +1155,41 @@ def fit_xcorrgrid(spec,models=None,samples=None,verbose=False,maxvel=1000.0): samples['teff'][:] = teff samples['logg'][:] = logg samples['feh'][:] = feh - outdtype = np.dtype([('xshift',np.float32),('vrel',np.float32),('vrelerr',np.float32),('ccpeak',np.float32),('ccpfwhm',np.float32), + outdtype = np.dtype([('xshift',np.float32),('vrel',np.float32),('vrelerr',np.float32),('ccpeak',np.float32),('ccp0',np.float32), + ('ccpfwhm',np.float32),('vrel0',np.float32), ('chisq',np.float32),('teff',np.float32),('logg',np.float32),('feh',np.float32)]) outstr = np.zeros(len(teff),dtype=outdtype) - if verbose is True: print('TEFF LOGG FEH VREL CCPEAK CHISQ') + if verbose is True: print('TEFF LOGG FEH VREL CCPEAK CCP0 CHISQ VREL0') for i in range(len(samples)): m = models([samples['teff'][i],samples['logg'][i],samples['feh'][i]],rv=0,wave=wavelog) - outstr1 = specxcorr(m.wave,m.flux,obs.flux,obs.err,maxlag) + mcont = polynorm(m.flux,obs.mask) + m.flux /= mcont + outstr1 = specxcorr(m.wave,m.flux,obs.flux,obs.err,maxlag,plot=plot) #if outstr1['chisq'] > 1000: # import pdb; pdb.set_trace() if verbose is True: - print('%-7.2f %5.2f %5.2f %5.2f %5.2f %5.2f' % (teff[i],logg[i],feh,outstr1['vrel'][0],outstr1['ccpeak'][0],outstr1['chisq'][0])) - for n in ['xshift','vrel','vrelerr','ccpeak','ccpfwhm','chisq']: outstr[n][i] = outstr1[n] + print('%-7.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f' % (teff[i],logg[i],feh,outstr1['vrel'][0],outstr1['ccpeak'][0],outstr1['ccp0'][0],outstr1['chisq'][0],outstr1['vrel0'][0])) + for n in ['xshift','vrel','vrelerr','ccpeak','ccpfwhm','chisq', 'ccp0', 'vrel0']: outstr[n][i] = outstr1[n] outstr['teff'][i] = teff[i] outstr['logg'][i] = logg[i] outstr['feh'][i] = feh + if plot : pdb.set_trace() # Get best fit bestind = np.argmin(outstr['chisq']) + bestind = np.argmax(outstr['ccp0']) beststr = outstr[bestind] - bestmodel = models(teff=beststr['teff'],logg=beststr['logg'],feh=beststr['feh'],rv=beststr['vrel']) + if usepeak : rv=beststr['vrel0'] + else : rv=beststr['vrel'] + bestmodel = models(teff=beststr['teff'],logg=beststr['logg'],feh=beststr['feh'],rv=rv) if verbose is True: print('Initial RV fit:') - printpars([beststr['teff'],beststr['logg'],beststr['feh'],beststr['vrel']],[None,None,None,beststr['vrelerr']]) + printpars([beststr['teff'],beststr['logg'],beststr['feh'],rv],[None,None,None,beststr['vrelerr']]) return beststr, bestmodel -def fit_lsq(spec,models=None,initpar=None,verbose=False): +def fit_lsq(spec,models=None,initpar=None,verbose=False,maxvel=[-1000,1000]): """ Least Squares fitting with forward modeling of the spectrum. @@ -1204,9 +1241,11 @@ def fit_lsq(spec,models=None,initpar=None,verbose=False): for p in models: lbounds[0:3] = np.minimum(lbounds[0:3],np.min(p.ranges,axis=1)) ubounds[0:3] = np.maximum(ubounds[0:3],np.max(p.ranges,axis=1)) - lbounds[3] = -1000 - ubounds[3] = 1000 + print('fit_lsq: ', maxvel) + lbounds[3] = maxvel[0] + ubounds[3] = maxvel[1] bounds = (lbounds, ubounds) + initpar = np.minimum(ubounds-0.01,np.maximum(lbounds+0.01,initpar)) # function to use with curve_fit def spec_interp(x,teff,logg,feh,rv): @@ -1215,11 +1254,14 @@ def spec_interp(x,teff,logg,feh,rv): m = models(teff=teff,logg=logg,feh=feh,rv=rv) if m is None: return np.zeros(spec.flux.shape,float).flatten()+1e30 - return m.flux.flatten() + else : + cont = polynorm(m.flux,spec.mask) + m.flux /= cont + return m.flux.flatten() # Use curve_fit lspars, lscov = curve_fit(spec_interp, spec.wave.flatten(), spec.flux.flatten(), sigma=spec.err.flatten(), - p0=initpar, bounds=bounds) + p0=initpar, bounds=bounds,diff_step=0.01) # If it hits a boundary then the solution won't change much compared to initpar # setting absolute_sigma=True gives crazy low lsperror values lsperror = np.sqrt(np.diag(lscov)) @@ -1342,7 +1384,7 @@ def fit_mcmc(spec,models=None,initpar=None,steps=100,cornername=None,verbose=Fal return out,mcmodel -def multifit_lsq(speclist,modlist,initpar=None,verbose=False): +def multifit_lsq(speclist,modlist,initpar=None,verbose=False,maxvel=[-1000,1000]): """ Least Squares fitting with forward modeling of multiple spectra simultaneously. @@ -1400,9 +1442,10 @@ def multifit_lsq(speclist,modlist,initpar=None,verbose=False): for p in modlist[0]: lbounds[0:3] = np.minimum(lbounds[0:3],np.min(p.ranges,axis=1)) ubounds[0:3] = np.maximum(ubounds[0:3],np.max(p.ranges,axis=1)) - lbounds[3:] = -1000 - ubounds[3:] = 1000 + lbounds[3:] = maxvel[0] + ubounds[3:] = maxvel[1] bounds = (lbounds, ubounds) + initpar = np.minimum(ubounds-0.01,np.maximum(lbounds+0.01,initpar)) # function to use with curve_fit def multispec_interp(x,*argv): @@ -1418,10 +1461,14 @@ def multispec_interp(x,*argv): flux = np.zeros(npix,float) cnt = 0 for i in range(nspec): - npx = speclist[i].npix*speclist[i].norder + if iorder< 0 : npx = speclist[i].npix*speclist[i].norder + else : npx = speclist[i].npix m = modlist[i]([teff,logg,feh],rv=vrel[i]) if m is not None: - flux[cnt:cnt+npx] = m.flux.T.flatten() + cont = polynorm(m.flux,speclist[i].mask) + m.flux /= cont + if iorder<0 : flux[cnt:cnt+npx] = m.flux.T.flatten() + else : flux[cnt:cnt+npx] = m.flux[:,iorder].T.flatten() else: flux[cnt:cnt+npx] = 1e30 cnt += npx @@ -1446,40 +1493,71 @@ def multispec_interp_jac(x,*argv): # Model at current values f0 = multispec_interp(x,*argv) # Compute full models for teff/logg/feh + step=np.array([5,0.01,0.01]) for i in range(3): pars = np.array(copy.deepcopy(argv)) - step = relstep*pars[i] - pars[i] += step + #step = relstep*pars[i] + #pars[i] += step + pars[i] += step[i] f1 = multispec_interp(x,*pars) # Hit an edge, try the negative value instead nbd = np.sum(f1>1000) if nbd>1000: pars = np.array(copy.deepcopy(argv)) - step = -relstep*pars[i] - pars[i] += step + #step = -relstep*pars[i] + #pars[i] += step + step[i] *= -1 + pars[i] += step[i] f1 = multispec_interp(x,*pars) - jac[:,i] = (f1-f0)/step + jac[:,i] = (f1-f0)/step[i] # Compute model for single spectra nspec = len(speclist) cnt = 0 + #step = 1.0 + step=0.1 for i in range(nspec): vrel1 = vrel[i] - step = 1.0 vrel1 += step - npx = speclist[i].npix*speclist[i].norder + if iorder<0: npx = speclist[i].npix*speclist[i].norder + else : npx = speclist[i].npix m = modlist[i]([teff,logg,feh],rv=vrel1) if m is not None: - jac[cnt:cnt+npx,i] = (m.flux.T.flatten()-f0[cnt:cnt+npx])/step + cont = polynorm(m.flux,speclist[i].mask) + m.flux /= cont + if iorder<0 : jac[cnt:cnt+npx,3+i] = (m.flux.T.flatten()-f0[cnt:cnt+npx])/step + else : jac[cnt:cnt+npx,3+i] = (m.flux[:,iorder].T.flatten()-f0[cnt:cnt+npx])/step else: - jac[cnt:cnt+npx,i] = 1e30 + jac[cnt:cnt+npx,3+i] = 1e30 cnt += npx - + return jac # We are fitting 3 stellar parameters and Nspec relative RVs # Put all of the spectra into a large 1D array + #for iorder in range(3) : + # ntotpix = 0 + # for s in speclist: + # ntotpix += s.npix + # wave = np.zeros(ntotpix) + # flux = np.zeros(ntotpix) + # err = np.zeros(ntotpix) + # cnt = 0 + # for i in range(nspec): + # sp = speclist[i] + # npx = sp.npix + # wave[cnt:cnt+npx] = sp.wave[:,iorder].T.flatten() + # flux[cnt:cnt+npx] = sp.flux[:,iorder].T.flatten() + # err[cnt:cnt+npx] = sp.err[:,iorder].T.flatten() + # cnt += npx + # + # # Use curve_fit + # pdb.set_trace() + # lspars, lscov = curve_fit(multispec_interp, wave, flux, sigma=err, p0=initpar, bounds=bounds, jac=multispec_interp_jac) + # print(iorder,lspars) + + iorder=-1 ntotpix = 0 for s in speclist: ntotpix += s.npix*s.norder @@ -1520,10 +1598,11 @@ def multispec_interp_jac(x,*argv): out['parcov'] = lscov out['chisq'] = lschisq + del wave, flux, err return out, lsmodel -def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=None,retpmodels=False): +def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=None,retpmodels=False,plot=False,tweak=True,usepeak=False,maxvel=[-1000,1000]) : """ Fit the spectrum. Find the best RV and stellar parameters using the Cannon models. @@ -1580,7 +1659,12 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No # Step 1: Prepare the spectrum #----------------------------- # Normalize and mask the spectrum - spec = utils.specprep(spec) + spec.normalized = True + spec = utils.specprep(spec) + spec.cont = polynorm(spec.flux,spec.mask) + spec.flux /= spec.cont + spec.err /= spec.cont + # Mask out any large positive outliers, e.g. badly subtracted sky lines specm = utils.maskoutliers(spec,verbose=verbose) @@ -1592,7 +1676,7 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No # Step 3: Get initial RV using cross-correlation with rough sampling of Teff/logg parameter space #------------------------------------------------------------------------------------------------ - beststr, xmodel = fit_xcorrgrid(specm,pmodels,verbose=verbose,maxvel=1000.0) + beststr, xmodel = fit_xcorrgrid(specm,pmodels,verbose=verbose,maxvel=maxvel,plot=plot) # Step 4: Get better Cannon stellar parameters using initial RV #-------------------------------------------------------------- @@ -1615,7 +1699,7 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No printpars(labels0) # Tweak the continuum normalization - specm = tweakcontinuum(specm,bestmodelspec0) + if tweak : specm = tweakcontinuum(specm,bestmodelspec0) # Mask out very discrepant pixels when compared to the best-fit model specm = utils.maskdiscrepant(specm,bestmodelspec0,verbose=verbose) @@ -1636,21 +1720,24 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No m = pmodels.get_best_model(labels).interp(wavelog)(labels,rv=0) dwlog = np.median(dln.slope(np.log10(wavelog))) # vrel = ( 10**(xshift*dwlog)-1 )*cspeed - maxlag = np.int(np.ceil(np.log10(1+1000.0/cspeed)/dwlog)) - maxlag = np.maximum(maxlag,50) + maxlag = np.ceil(np.log10(1+np.array(maxvel)/cspeed)/dwlog).astype(int) + #maxlag = np.int(np.ceil(np.log10(1+np.array(maxvel)/cspeed)/dwlog)) + #maxlag = np.maximum(maxlag,50) outstr2 = specxcorr(m.wave,m.flux,obs.flux,obs.err,maxlag) - outdtype = np.dtype([('xshift',np.float32),('vrel',np.float32),('vrelerr',np.float32),('ccpeak',np.float32),('ccpfwhm',np.float32), + outdtype = np.dtype([('xshift',np.float32),('vrel',np.float32),('vrel0',np.float32),('vrelerr',np.float32),('ccpeak',np.float32),('ccpfwhm',np.float32), ('chisq',np.float32),('teff',np.float32),('logg',np.float32),('feh',np.float32)]) beststr2= np.zeros(1,dtype=outdtype) - for n in ['xshift','vrel','vrelerr','ccpeak','ccpfwhm','chisq']: beststr2[n] = outstr2[n] + for n in ['xshift','vrel','vrelerr','ccpeak','ccpfwhm','chisq','vrel0']: beststr2[n] = outstr2[n] beststr2['teff'] = labels[0] beststr2['logg'] = labels[1] beststr2['feh'] = labels[2] + if usepeak : bestrv = beststr2['vrel0'] + else : bestrv = beststr['vrel'] # Step 6: Improved Cannon stellar parameters #------------------------------------------- - restwave = specm.wave*(1-beststr['vrel']/cspeed) + restwave = specm.wave*(1-bestrv/cspeed) bestmodel = pmodels.get_best_model([beststr2['teff'],beststr2['logg'],beststr2['feh']]) bestmodelinterp = bestmodel.interp(restwave) labels2, cov2, meta2 = bestmodelinterp.test(specm) @@ -1660,25 +1747,28 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No bestmodelspec2 = bestmodelinterp(labels2) if verbose is True: print('Improved RV and Cannon stellar parameters:') - printpars(np.concatenate((labels2,beststr2['vrel'])),[None,None,None,beststr2['vrelerr']]) + printpars(np.concatenate((labels2,bestrv)),[None,None,None,beststr2['vrelerr']]) # Step 7: Least Squares fitting with forward modeling #---------------------------------------------------- # Get best model so far - m = pmodels(teff=beststr2['teff'],logg=beststr2['logg'],feh=beststr2['feh'],rv=beststr2['vrel']) + m = pmodels(teff=beststr2['teff'],logg=beststr2['logg'],feh=beststr2['feh'],rv=bestrv) # Tweak the continuum - specm = tweakcontinuum(specm,m) + if tweak: specm = tweakcontinuum(specm,m) # Get initial estimates - initpar = [beststr2['teff'],beststr2['logg'],beststr2['feh'],beststr2['vrel']] + initpar = [beststr2['teff'],beststr2['logg'],beststr2['feh'],bestrv] initpar = np.array(initpar).flatten() - lsout, lsmodel = fit_lsq(specm,pmodels,initpar=initpar,verbose=verbose) + lsout, lsmodel = fit_lsq(specm,pmodels,initpar=initpar,verbose=verbose,maxvel=maxvel) lspars = lsout['pars'][0] lsperror = lsout['parerr'][0] # Step 8: Run fine grid in RV, forward modeling #---------------------------------------------- - maxv = np.maximum(beststr2['vrel'][0],20.0) + #maxv = np.maximum(bestrv,20.0) + # with maxv=bestrv that can be a big range! What is the point anyway, + # after doing the least squares? + maxv=20. vel = dln.scale_vector(np.arange(30),lspars[3]-maxv,lspars[3]+maxv) chisq = np.zeros(len(vel)) for i,v in enumerate(vel): @@ -1739,6 +1829,7 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No # Make diagnostic figure if figfile is not None: # Apply continuum tweak to original spectrum as well + if hasattr(spec,'cont') is False: spec.cont = spec.flux.copy()*0+1 cratio = specm.cont/spec.cont orig = spec.copy() orig.flux /= cratio @@ -1749,7 +1840,9 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No # How long did this take if verbose is True: print('dt = %5.2f sec.' % (time.time()-t0)) - + + del spec, wavelog, lsmodel, m + # Return the prpared models if retpmodels is True: return out, fmodel, specm, pmodels @@ -1758,28 +1851,39 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No -def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose=False,outdir=None): +def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose=False,outdir=None,plot=False,tweak=True,maxlag=400,maxvel=[-500,500],usepeak=True) : """This fits a Cannon model to multiple spectra of the same star.""" # speclist is list of Spec1D objects. nspec = len(speclist) t0 = time.time() - + # If list of filenames input, then load them # Creating catalog of info on each spectrum + nlag = 2*np.round(np.abs(maxlag))+1 dt = np.dtype([('filename',np.str,300),('snr',float),('vhelio',float),('vrel',float),('vrelerr',float), ('teff',float),('tefferr',float),('logg',float),('loggerr',float),('feh',float), - ('feherr',float),('chisq',float),('bc',float)]) + ('feherr',float),('chisq',float),('bc',float),('x_ccf',(float,nlag)),('ccf',(float,nlag)), + ('ccferr',(float,nlag)),('xcorr_vrel',float),('xcorr_vrelerr',float),('xcorr_vhelio',float), + ('ccpfwhm',float),('autofwhm',float)]) info = np.zeros(nspec,dtype=dt) for n in dt.names: info[n] = np.nan for i,s in enumerate(speclist): info['filename'][i] = s.filename info['snr'][i] = s.snr + # Create catalog of summary info + sumdt = np.dtype([('medsnr',float),('totsnr',float),('vhelio',float),('vscatter',float),('verr',float), + ('teff',float),('tefferr',float),('logg',float),('loggerr',float),('feh',float), + ('feherr',float),('chisq',float)]) + sumstr = np.zeros(1,dtype=sumdt) + # Make sure some spectra pass the S/N cut hisnr, nhisnr = dln.where(info['snr']>snrcut) - if nhisnr < np.ceil(0.25*nspec): + if nspec == 1 : + snrcut=0 + elif nhisnr < np.ceil(0.25*nspec): snr = np.flip(np.sort(info['snr'])) snrcut = snr[np.maximum(np.int(np.ceil(0.25*nspec)),np.minimum(4,nspec-1))] if verbose is True: @@ -1789,6 +1893,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if verbose is True: print('Step #1: Fitting the individual spectra') specmlist = [] modlist = [] + bdlist = [] for i in range(len(speclist)): spec = speclist[i].copy() if verbose is True: @@ -1804,7 +1909,31 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if outdir is not None: figfile = outdir+'/'+figfile if (outdir is None) & (fdir != ''): figfile = fdir+'/'+figfile # Fit the spectrum - out, model, specm, pmodels = fit(spec,verbose=verbose,mcmc=mcmc,figfile=figfile,retpmodels=True) + try : + out, model, specm, pmodels = \ + fit(spec,verbose=verbose,mcmc=mcmc,figfile=figfile,retpmodels=True, + plot=plot,tweak=tweak,usepeak=usepeak,maxvel=maxvel) + except RuntimeError as err : + print('Exception raised for: ', speclist[i].filename) + print("Runtime error: {0}".format(err)) + # if we had a failure in fit, treat it as lower S/N object and see if + # we can fit it in the multifit_lsq step + #print('removing from list ....') + #bdlist.append(i) + modlist.append(cannon.models.prepare(speclist[i]).copy()) + sp = speclist[i].copy() + sp.normalized = True + sp = utils.specprep(sp) # mask and normalize + sp.cont = polynorm(sp.flux,sp.mask) + sp.flux /= sp.cont + sp.err /= sp.cont + # Mask outliers + sp = utils.maskoutliers(sp) + specmlist.append(sp) + # at least need BC + info['bc'][i] = speclist[i].barycorr() + continue + modlist.append(pmodels.copy()) del pmodels specmlist.append(specm.copy()) @@ -1820,12 +1949,18 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= info['feherr'][i] = out['feherr'] info['chisq'][i] = out['chisq'] info['bc'][i] = out['bc'] + del out + del model else: if verbose is True: print('Skipping: S/N=%6.1f below threshold of %6.1f. Loading spectrum and preparing models.' % (spec.snr,snrcut)) modlist.append(cannon.models.prepare(speclist[i]).copy()) sp = speclist[i].copy() + sp.normalized = True sp = utils.specprep(sp) # mask and normalize + sp.cont = polynorm(sp.flux,sp.mask) + sp.flux /= sp.cont + sp.err /= sp.cont # Mask outliers sp = utils.maskoutliers(sp) specmlist.append(sp) @@ -1833,6 +1968,40 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= info['bc'][i] = speclist[i].barycorr() if verbose is True: print(' ') + # remove failed frames from list + info = np.delete(info,bdlist) + nspec -= len(bdlist) + + if len(speclist) == 1 : + # if we only have one image, we are done, create summary structure and return + sumstr['medsnr'] = info['snr'][0] + sumstr['totsnr'] = info['snr'][0] + sumstr['vhelio'] = info['vhelio'][0] + sumstr['verr'] = info['vrelerr'][0] + for key in ['teff','tefferr','logg','loggerr','feh','feherr','chisq'] : + sumstr[key] = info[key][0] + + pars1 = [info['teff'][0], info['logg'][0], info['feh'][0]] + vr1 = info['vrel'][0] + outstr,model_outstr = final_xcorr(specmlist[0],modlist[0],pars1,vr1,maxvel=maxvel,plot=plot) + m = modlist[0](pars1,rv=vr1) + cont = polynorm(m.flux,specmlist[0].mask) + m.flux /= cont + if usepeak : rv=outstr['vrel0'] + else : rv=outstr['vrel'] + if plot : pdb.set_trace() + final = info.copy() + nlag = len(outstr['ccf'][0]) + final['x_ccf'][0][0:nlag] = outstr['ccvlag']+info['bc'][i] + final['ccf'][0][0:nlag] = outstr['ccf'] + final['ccferr'][0][0:nlag] = outstr['ccferr'] + final['xcorr_vrel'][0] = rv+vr1 + final['xcorr_vrelerr'][0] = outstr['vrelerr'] + final['xcorr_vhelio'][0] = rv+vr1+info['bc'][i] + final['ccpfwhm'][0] = outstr['ccpfwhm'] + final['autofwhm'][0] = model_outstr['ccpfwhm'] + return sumstr, final, [m], specmlist, time.time()-t0 + # Step 2) find weighted stellar parameters if verbose is True: print('Step #2: Getting weighted stellar parameters') @@ -1853,6 +2022,8 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= # Unweighted else: wtpars[i] = np.mean(p) + # weighted by S/N + wtpars[i] = dln.wtmean(p,info['snr'][gd]) if verbose is True: print('Initial weighted parameters are:') printpars(wtpars) @@ -1867,7 +2038,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= print('No good fits. Using these as intial guesses:') printpars(wtpars) - # Make initial guesses for all the parameters, 3 stellar paramters and Nspec relative RVs + # Make initial guesses for all the parameters, 3 stellar parameters and Nspec relative RVs initpar1 = np.zeros(3+nspec,float) initpar1[0:3] = wtpars[0:3] # the default is to use mean vhelio + BC for all visit spectra @@ -1877,12 +2048,12 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if ngdinit>0: initpar1[gdinit+3] = info['vrel'][gdinit] - # Step 3) refit all spectra simultaneous fitting stellar parameters and RVs if verbose is True: print(' ') print('Step #3: Fitting all spectra simultaneously') - out1, fmodels1 = multifit_lsq(specmlist,modlist,initpar1) + print('initpar1: ', initpar1) + out1, fmodels1 = multifit_lsq(specmlist,modlist,initpar1,maxvel=maxvel) stelpars1 = out1['pars'][0,0:3] stelparerr1 = out1['parerr'][0,0:3] vrel1 = out1['pars'][0,3:] @@ -1904,12 +2075,19 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= print('Step #4: Tweaking continuum and masking outliers') for i,spm in enumerate(specmlist): bestm = modlist[i](stelpars1,rv=vrel1[i]) + if bestm is None: + raise RuntimeError('No valid model found for: ',stelpars1) + cont = polynorm(bestm.flux,specmlist[i].mask) + bestm.flux /= cont # Tweak the continuum normalization - spm = tweakcontinuum(spm,bestm) + if tweak : spm = tweakcontinuum(spm,bestm) # Mask out very discrepant pixels when compared to the best-fit model spm = utils.maskdiscrepant(spm,bestm,verbose=verbose) + if stelpars1[0] > 7500 : + bd = np.where((spm.wave[:,1]<15900) | (spm.wave[:,1]>16350) )[0] + spm.err[bd,1] = 1.e30 specmlist[i] = spm.copy() - + # Step 5) refit all spectra simultaneous fitting stellar parameters and RVs if verbose is True: print(' ') @@ -1917,14 +2095,16 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= # Initial guesses for all the parameters, 3 stellar paramters and Nspec relative RVs initpar2 = out1['pars'][0] - out2, fmodels2 = multifit_lsq(specmlist,modlist,initpar2) + del out1, fmodels1 + out2, fmodels2 = multifit_lsq(specmlist,modlist,initpar2,maxvel=maxvel) stelpars2 = out2['pars'][0,0:3] stelparerr2 = out2['parerr'][0,0:3] vrel2 = out2['pars'][0,3:] vrelerr2 = out2['parerr'][0,3:] vhelio2 = vrel2+info['bc'] medvhelio2 = np.median(vhelio2) - vscatter2 = dln.mad(vhelio2) + #vscatter2 = dln.mad(vhelio2) + vscatter2 = vhelio2.std(ddof=1) verr2 = vscatter2/np.sqrt(nspec) if verbose is True: print('Final parameters:') @@ -1932,7 +2112,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= print('Vhelio = %6.2f +/- %5.2f km/s' % (medvhelio2,verr2)) print('Vscatter = %6.3f km/s' % vscatter2) print(vhelio2) - + # Final output structure final = info.copy() final['teff'] = stelpars2[0] @@ -1950,20 +2130,42 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= for i in range(nspec): pars1 = [final['teff'][i], final['logg'][i], final['feh'][i]] vr1 = final['vrel'][i] - sp = specmlist[i] + sp = specmlist[i].copy() m = modlist[i](pars1,rv=vr1) + cont = polynorm(m.flux,specmlist[i].mask) + m.flux /= cont chisq = np.sqrt(np.sum(((sp.flux-m.flux)/sp.err)**2)/(sp.npix*sp.norder)) totchisq += np.sum(((sp.flux-m.flux)/sp.err)**2) totnpix += sp.npix*sp.norder final['chisq'][i] = chisq bmodel.append(m) + + # final cross-correlation + #outstr = final_xcorr(sp,modlist[i],pars1,vr1,maxvel=maxvel,plot=plot) + outstr,model_outstr = final_xcorr(sp,modlist[i],pars1,vr1,maxvel=[-1000,1000],plot=plot) + if usepeak : rv=outstr['vrel0'] + else : rv=outstr['vrel'] + if plot : pdb.set_trace() + #final['x_ccf'][i] = (np.arange(nlag)-maxlag)*(np.log(m.wave[1,0])-np.log(m.wave[0,0]))*3.e5+vr1 + nlag = len(outstr['ccf'][0]) + # put CCF on heliocentric scale + final['x_ccf'][i][0:nlag] = outstr['ccvlag']+info['bc'][i] + final['ccf'][i][0:nlag] = outstr['ccf'] + final['ccferr'][i][0:nlag] = outstr['ccferr'] + # for xcorr_vhelio, restrict to input velocity range + gd = np.where((final['x_ccf'][i] >= maxvel[0]+info['bc'][i]) & + (final['x_ccf'][i] <= maxvel[1]+info['bc'][i]) )[0] + imax = final['ccf'][i][gd].argmax() + xcorr_vhelio = final['x_ccf'][i][gd[imax]] + final['xcorr_vhelio'][i] = xcorr_vhelio + final['xcorr_vrel'][i] =xcorr_vhelio - info['bc'][i] + final['xcorr_vrelerr'][i] = outstr['vrelerr'] + final['ccpfwhm'][i] = outstr['ccpfwhm'] + final['autofwhm'][i] = model_outstr['ccpfwhm'] + totchisq = np.sqrt(totchisq/totnpix) # Average values - sumdt = np.dtype([('medsnr',float),('totsnr',float),('vhelio',float),('vscatter',float),('verr',float), - ('teff',float),('tefferr',float),('logg',float),('loggerr',float),('feh',float), - ('feherr',float),('chisq',float)]) - sumstr = np.zeros(1,dtype=sumdt) sumstr['medsnr'] = np.median(info['snr']) sumstr['totsnr'] = np.sqrt(np.sum(info['snr']**2)) sumstr['vhelio'] = medvhelio2 @@ -1980,4 +2182,42 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= # How long did this take if verbose is True: print('dt = %5.2f sec.' % (time.time()-t0)) - return sumstr, final, bmodel, specmlist + return sumstr, final, bmodel, specmlist, time.time()-t0 + +def final_xcorr(sp,model,pars,rv,maxvel=[-500,500],plot=False) : + """ Get a final cross correlation with best fit spectrum + """ + wavelog = utils.make_logwave_scale(sp.wave,vel=0.0) # get new wavelength solution + m = model(pars,rv=rv,wave=wavelog) + if m is None: + raise RuntimeError('No valid model found for: ',pars) + obs = sp.interp(wavelog) + mcont = polynorm(m.flux,obs.mask) + m.flux /= mcont + dwlog = np.median(dln.slope(np.log10(wavelog))) + maxlag = np.ceil(np.log10(1+np.array(np.array(maxvel)-rv)/cspeed)/dwlog).astype(int) + nlag = maxlag[1]-maxlag[0]+1 + if ((nlag % 2) == 0): nlag +=1 # make sure nlag is odd + outstr = specxcorr(m.wave,m.flux,obs.flux,obs.err,maxlag,plot=plot) + outstr['ccvlag'] = rv+outstr['cclag']*dwlog*cspeed*np.log(10) + auto_outstr = specxcorr(m.wave,m.flux,m.flux,obs.err,maxlag,plot=plot) + return outstr,auto_outstr + +def polynorm(flux,mask,order=4) : + """ simple polynomial continuum + """ + x = np.arange(flux.shape[0]) + cont = np.full_like(flux,1.) + if len(flux.shape) > 1 : + norder = flux.shape[1] + for iorder in range(norder) : + gd=np.where(~mask[:,iorder] & np.isfinite(flux[:,iorder]))[0] + if len(gd) > order : + coef = np.polyfit(x[gd],flux[gd,iorder],order) + cont[:,iorder] = np.polyval(coef,x) + else : + gd=np.where(~mask & np.isfinite(flux))[0] + if len(gd) > order : + coef = np.polyfit(x[gd],flux[gd],order) + cont = np.polyval(coef,x) + return cont diff --git a/doppler/spec1d.py b/doppler/spec1d.py index c629fef..9ca34b3 100644 --- a/doppler/spec1d.py +++ b/doppler/spec1d.py @@ -426,7 +426,7 @@ def interp(self,x=None,xtype='wave',order=None): smask1[bd] = 1 ind,nind = dln.where( (wave1>np.min(swave1)) & (wave1wr1[0]): ind,nind = dln.where( (wave1>np.min(swave1)) & (wave1