From 98cbb20fc538af32d1f20cddada9732d98452cfb Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Wed, 3 Jun 2020 12:43:34 -0600 Subject: [PATCH 01/18] don't raise exception for nbin, assuming this is an APOGEE LSF bug issues ... --- doppler/cannon.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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 From ee49a85e9fc0e87f585a1398501a8543ce470be1 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Wed, 3 Jun 2020 12:44:03 -0600 Subject: [PATCH 02/18] add nowings processing --- doppler/lsf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From e0b93ddbcf30612dfc90ebead08474ce3d3e1093 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Wed, 3 Jun 2020 12:44:34 -0600 Subject: [PATCH 03/18] add use SIG_SKYLINE bit for masking turn off masking based on sky spectra only --- doppler/reader.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/doppler/reader.py b/doppler/reader.py index eb80313..f78e6bb 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") @@ -165,17 +167,24 @@ 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,20735)!=0) | (np.isfinite(spec.flux)==False) # Extra masking for bright skylines 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.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 + #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() + #pdb.set_trace() spec.mask = mask # Fix NaN pixels for i in range(spec.norder): From c856990d9c72a671deacb0e25b937143dfcab19d Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Wed, 3 Jun 2020 12:45:17 -0600 Subject: [PATCH 04/18] change interpolation for err to linear, to avoid wild swings near 1.e30 values --- doppler/spec1d.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 Date: Wed, 3 Jun 2020 12:46:20 -0600 Subject: [PATCH 05/18] add some keywords for options: plot, tweak change default to no gaussian smoothing for xcorr raise specific exceptions rather than general add raw ccpeak to xcorr output and use it for selecting initial model add check for initial pars within bounds fix major bug in multifit_lsq jacobian change step sizes for derivatives add some memory cleanup add xcorr of final model with each spectrum and put in output modify to allow jointfit to run with single spectrum, for uniform interface trap for failed initial RVs in jointfit, and remove these, but proceed with others change weighting for initial paramters to be by S/N change vscatter to use std and not mad --- doppler/rv.py | 160 +++++++++++++++++++++++++++++++++++++------------- 1 file changed, 120 insertions(+), 40 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index afb3eec..250649a 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") @@ -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,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 @@ -609,8 +610,20 @@ 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) + plt.close() + #plt.plot((spec[gdmask]-temp[gdmask])**2/err[gdmask]**2 ) + from tools import plots + fig,ax=plots.multi(1,2,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') + #plt.plot(wobs[gdmask],err[gdmask]) + for i in range(3) : ax[0].plot(wobs[gdmask[:,i],i],temp[gdmask[:,i],i],color='r') outstr["chisq"] = chisq outstr["ccf"] = ccf @@ -619,8 +632,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) @@ -683,6 +698,12 @@ 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)) + + if plot : + ax[1].plot(lag,ccf) + ax[1].plot(lag,ccf_diff) + ax[1].plot(lag[lo3:hi3],yfit3) + plt.draw() # Final parameters pars = pars3 @@ -776,7 +797,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 +876,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 +902,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 +1069,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.0,plot=False): """ Fit spectrum using cross-correlation with models sampled in the parameter space. @@ -1125,23 +1146,26 @@ 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), ('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') 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) + 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' % (teff[i],logg[i],feh,outstr1['vrel'][0],outstr1['ccpeak'][0],outstr1['ccp0'][0],outstr1['chisq'][0])) + for n in ['xshift','vrel','vrelerr','ccpeak','ccp0','ccpfwhm','chisq']: 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']) @@ -1207,6 +1231,7 @@ def fit_lsq(spec,models=None,initpar=None,verbose=False): lbounds[3] = -1000 ubounds[3] = 1000 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): @@ -1403,6 +1428,7 @@ def multifit_lsq(speclist,modlist,initpar=None,verbose=False): lbounds[3:] = -1000 ubounds[3:] = 1000 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): @@ -1446,34 +1472,39 @@ 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 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 + jac[cnt:cnt+npx,3+i] = (m.flux.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 @@ -1520,10 +1551,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) : """ Fit the spectrum. Find the best RV and stellar parameters using the Cannon models. @@ -1592,7 +1624,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=1000.0,plot=plot) # Step 4: Get better Cannon stellar parameters using initial RV #-------------------------------------------------------------- @@ -1615,7 +1647,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) @@ -1667,7 +1699,7 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No # Get best model so far m = pmodels(teff=beststr2['teff'],logg=beststr2['logg'],feh=beststr2['feh'],rv=beststr2['vrel']) # 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']] @@ -1739,6 +1771,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 +1782,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,25 +1793,33 @@ 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=200) : """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)), + ('xcorr_vrel',float),('xcorr_vrelerr',float),('xcorr_vhelio',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): @@ -1789,6 +1832,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 +1848,15 @@ 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) + except RuntimeError as err : + print('Exception raised for: ', speclist[i].filename) + print("Runtime error: {0}".format(err)) + print('removing from list ....') + bdlist.append(i) + continue + modlist.append(pmodels.copy()) del pmodels specmlist.append(specm.copy()) @@ -1820,6 +1872,8 @@ 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)) @@ -1833,6 +1887,25 @@ 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 + for bd in bdlist : + print('removing: ', bd) + del speclist[bd] + print(len(info)) + info = np.delete(info,bd) + print(len(info)) + nspec-=1 + + 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] + return 0, sumstr, info, modlist, specmlist, time.time()-t0 + # Step 2) find weighted stellar parameters if verbose is True: print('Step #2: Getting weighted stellar parameters') @@ -1853,6 +1926,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) @@ -1871,17 +1946,18 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= initpar1 = np.zeros(3+nspec,float) initpar1[0:3] = wtpars[0:3] # the default is to use mean vhelio + BC for all visit spectra + print(len(initpar1),len(info)) initpar1[3:] = wtpars[3]-info['bc'] # vhelio = vrel + BC # Use the Vrel values from the initial fitting if they are accurate enough gdinit,ngdinit = dln.where(np.isfinite(info['vrel']) & (info['snr']>5)) 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') + print('initpar1: ', initpar1) out1, fmodels1 = multifit_lsq(specmlist,modlist,initpar1) stelpars1 = out1['pars'][0,0:3] stelparerr1 = out1['parerr'][0,0:3] @@ -1905,11 +1981,11 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= for i,spm in enumerate(specmlist): bestm = modlist[i](stelpars1,rv=vrel1[i]) # 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) specmlist[i] = spm.copy() - + # Step 5) refit all spectra simultaneous fitting stellar parameters and RVs if verbose is True: print(' ') @@ -1917,6 +1993,7 @@ 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] + del out1, fmodels1 out2, fmodels2 = multifit_lsq(specmlist,modlist,initpar2) stelpars2 = out2['pars'][0,0:3] stelparerr2 = out2['parerr'][0,0:3] @@ -1924,7 +2001,8 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= 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:') @@ -1956,14 +2034,16 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= totchisq += np.sum(((sp.flux-m.flux)/sp.err)**2) totnpix += sp.npix*sp.norder final['chisq'][i] = chisq + outstr = specxcorr(m.wave,m.flux,sp.flux,sp.err,maxlag,plot=plot) + final['x_ccf'][i] = (np.arange(nlag)-maxlag)*(np.log(m.wave[1,0])-np.log(m.wave[0,0]))*3.e5+vr1 + final['ccf'][i] = outstr['ccf'] + final['xcorr_vrel'][i] = outstr['vrel']+vr1 + final['xcorr_vrelerr'][i] = outstr['vrelerr'] + final['xcorr_vhelio'][i] = outstr['vrel']+vr1+info['bc'][i] bmodel.append(m) 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 +2060,4 @@ 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 From e16202bc55cb9492e8ee54cf3855aebb23652ee4 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Wed, 3 Jun 2020 12:52:01 -0600 Subject: [PATCH 06/18] add .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0d20b64 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.pyc From 83c95c0c802f3802506c255ebc7f17f912597b4b Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Wed, 24 Jun 2020 22:04:37 -0600 Subject: [PATCH 07/18] option to use CCF peak instead of fit put in keywords for maxvel add new continuum normalization to normalize obs and models the same way --- doppler/rv.py | 202 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 141 insertions(+), 61 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index 250649a..4404055 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -46,7 +46,7 @@ def xcorr_dtype(nlag): """Return the dtype for the xcorr structure""" - dtype = np.dtype([("xshift0",float),("ccp0",float),("xshift",float),("xshifterr",float), + 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)),("ccpeak",float),("ccpfwhm",float),("ccp_pars",(float,4)), ("ccp_perror",(float,4)),("ccp_polycoef",(float,4)),("vrel",float), @@ -615,15 +615,17 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,gfilt= if plot : print(len(spec),len(gdmask), chisq) - plt.close() - #plt.plot((spec[gdmask]-temp[gdmask])**2/err[gdmask]**2 ) + matplotlib.use('TkAgg') from tools import plots - fig,ax=plots.multi(1,2,figsize=(8,6)) + 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') - #plt.plot(wobs[gdmask],err[gdmask]) - for i in range(3) : ax[0].plot(wobs[gdmask[:,i],i],temp[gdmask[:,i],i],color='r') + 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 @@ -699,12 +701,6 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,gfilt= dyfit3 = dln.gaussian(lag[lo3:hi3],*pars3) perror3 = np.sqrt(np.diag(dcov3)) - if plot : - ax[1].plot(lag,ccf) - ax[1].plot(lag,ccf_diff) - ax[1].plot(lag[lo3:hi3],yfit3) - plt.draw() - # Final parameters pars = pars3 perror = perror3 @@ -719,6 +715,7 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,gfilt= #--------------------------------- # 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 @@ -728,6 +725,7 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,gfilt= # 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 @@ -741,6 +739,14 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,gfilt= 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 @@ -1069,7 +1075,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,plot=False): +def fit_xcorrgrid(spec,models=None,samples=None,verbose=False,maxvel=1000.0,plot=False,usepeak=False): """ Fit spectrum using cross-correlation with models sampled in the parameter space. @@ -1147,36 +1153,40 @@ def fit_xcorrgrid(spec,models=None,samples=None,verbose=False,maxvel=1000.0,plot samples['logg'][:] = logg samples['feh'][:] = feh outdtype = np.dtype([('xshift',np.float32),('vrel',np.float32),('vrelerr',np.float32),('ccpeak',np.float32),('ccp0',np.float32), - ('ccpfwhm',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 CCP0 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) + 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 %5.2f' % (teff[i],logg[i],feh,outstr1['vrel'][0],outstr1['ccpeak'][0],outstr1['ccp0'][0],outstr1['chisq'][0])) - for n in ['xshift','vrel','vrelerr','ccpeak','ccp0','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() + 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): """ Least Squares fitting with forward modeling of the spectrum. @@ -1228,8 +1238,8 @@ 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 + lbounds[3] = -maxvel + ubounds[3] = maxvel bounds = (lbounds, ubounds) initpar = np.minimum(ubounds-0.01,np.maximum(lbounds+0.01,initpar)) @@ -1240,11 +1250,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)) @@ -1367,7 +1380,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): """ Least Squares fitting with forward modeling of multiple spectra simultaneously. @@ -1425,8 +1438,8 @@ 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 + ubounds[3:] = maxvel bounds = (lbounds, ubounds) initpar = np.minimum(ubounds-0.01,np.maximum(lbounds+0.01,initpar)) @@ -1444,10 +1457,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 @@ -1497,10 +1514,14 @@ def multispec_interp_jac(x,*argv): for i in range(nspec): vrel1 = vrel[i] 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,3+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,3+i] = 1e30 cnt += npx @@ -1511,6 +1532,28 @@ def multispec_interp_jac(x,*argv): # 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 @@ -1555,7 +1598,7 @@ def multispec_interp_jac(x,*argv): return out, lsmodel -def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=None,retpmodels=False,plot=False,tweak=True) : +def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=None,retpmodels=False,plot=False,tweak=True,usepeak=False,maxvel=1000) : """ Fit the spectrum. Find the best RV and stellar parameters using the Cannon models. @@ -1612,7 +1655,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) @@ -1624,7 +1672,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,plot=plot) + beststr, xmodel = fit_xcorrgrid(specm,pmodels,verbose=verbose,maxvel=maxvel,plot=plot) # Step 4: Get better Cannon stellar parameters using initial RV #-------------------------------------------------------------- @@ -1668,21 +1716,23 @@ 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.int(np.ceil(np.log10(1+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) @@ -1692,17 +1742,17 @@ 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 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) lspars = lsout['pars'][0] @@ -1710,7 +1760,7 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No # Step 8: Run fine grid in RV, forward modeling #---------------------------------------------- - maxv = np.maximum(beststr2['vrel'][0],20.0) + maxv = np.maximum(bestrv,20.0) vel = dln.scale_vector(np.arange(30),lspars[3]-maxv,lspars[3]+maxv) chisq = np.zeros(len(vel)) for i,v in enumerate(vel): @@ -1793,7 +1843,7 @@ 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,plot=False,tweak=True,maxlag=200) : +def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose=False,outdir=None,plot=False,tweak=True,maxlag=200,maxvel=500,usepeak=True) : """This fits a Cannon model to multiple spectra of the same star.""" # speclist is list of Spec1D objects. @@ -1807,7 +1857,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= 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),('x_ccf',(float,nlag)),('ccf',(float,nlag)), - ('xcorr_vrel',float),('xcorr_vrelerr',float),('xcorr_vhelio',float)]) + ('ccferr',(float,nlag)),('xcorr_vrel',float),('xcorr_vrelerr',float),('xcorr_vhelio',float)]) info = np.zeros(nspec,dtype=dt) for n in dt.names: info[n] = np.nan for i,s in enumerate(speclist): @@ -1849,7 +1899,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if (outdir is None) & (fdir != ''): figfile = fdir+'/'+figfile # Fit the spectrum try : - out, model, specm, pmodels = fit(spec,verbose=verbose,mcmc=mcmc,figfile=figfile,retpmodels=True,plot=plot,tweak=tweak) + 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)) @@ -1879,7 +1929,11 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= 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) @@ -1888,13 +1942,8 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if verbose is True: print(' ') # remove failed frames from list - for bd in bdlist : - print('removing: ', bd) - del speclist[bd] - print(len(info)) - info = np.delete(info,bd) - print(len(info)) - nspec-=1 + 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 @@ -1958,7 +2007,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= print(' ') print('Step #3: Fitting all spectra simultaneously') print('initpar1: ', initpar1) - out1, fmodels1 = multifit_lsq(specmlist,modlist,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:] @@ -1980,10 +2029,15 @@ 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]) + cont = polynorm(bestm.flux,specmlist[i].mask) + bestm.flux /= cont # Tweak the continuum normalization 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 @@ -1994,7 +2048,7 @@ 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] del out1, fmodels1 - out2, fmodels2 = multifit_lsq(specmlist,modlist,initpar2) + 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:] @@ -2010,7 +2064,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] @@ -2028,19 +2082,32 @@ 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 - outstr = specxcorr(m.wave,m.flux,sp.flux,sp.err,maxlag,plot=plot) + bmodel.append(m) + + # final cross-correlation + wavelog = utils.make_logwave_scale(sp.wave,vel=0.0) # get new wavelength solution + obs = sp.interp(wavelog) + m = modlist[i](pars1,rv=vr1,wave=wavelog) + mcont = polynorm(m.flux,obs.mask) + m.flux /= mcont + outstr = specxcorr(m.wave,m.flux,obs.flux,obs.err,maxlag,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 final['ccf'][i] = outstr['ccf'] - final['xcorr_vrel'][i] = outstr['vrel']+vr1 + final['ccferr'][i] = outstr['ccferr'] + final['xcorr_vrel'][i] = rv+vr1 final['xcorr_vrelerr'][i] = outstr['vrelerr'] - final['xcorr_vhelio'][i] = outstr['vrel']+vr1+info['bc'][i] - bmodel.append(m) + final['xcorr_vhelio'][i] = rv+vr1+info['bc'][i] totchisq = np.sqrt(totchisq/totnpix) # Average values @@ -2061,3 +2128,16 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if verbose is True: print('dt = %5.2f sec.' % (time.time()-t0)) return sumstr, final, bmodel, specmlist, time.time()-t0 + +def polynorm(flux,mask,order=4) : + """ simple polynomial continuum + """ + x = np.arange(flux.shape[0]) + norder = flux.shape[1] + cont = np.full_like(flux,1.) + for iorder in range(norder) : + gd=np.where(~mask[:,iorder] & np.isfinite(flux[:,iorder]))[0] + if len(gd) > norder : + coef = np.polyfit(x[gd],flux[gd,iorder],order) + cont[:,iorder] = np.polyval(coef,x) + return cont From 9317e62221add3884dc28cff4fea3c5d1ae5ae1d Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Tue, 30 Jun 2020 17:49:31 -0600 Subject: [PATCH 08/18] allow to specify velocity range to use mod for polynorm for 1D --- doppler/rv.py | 80 +++++++++++++++++++++++++++++++-------------------- 1 file changed, 49 insertions(+), 31 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index 4404055..c515c27 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -472,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,gfilt=0,errccf=False,prior=None,plot=False): +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 @@ -527,10 +527,12 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,gfilt= # 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 @@ -648,8 +650,8 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,gfilt= # 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)) @@ -665,7 +667,7 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=200,gfilt= 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] @@ -1075,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,plot=False,usepeak=False): +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. @@ -1138,8 +1140,9 @@ def fit_xcorrgrid(spec,models=None,samples=None,verbose=False,maxvel=1000.0,plot #------------------------------------------------------------------------------------------------ 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] @@ -1186,7 +1189,7 @@ def fit_xcorrgrid(spec,models=None,samples=None,verbose=False,maxvel=1000.0,plot return beststr, bestmodel -def fit_lsq(spec,models=None,initpar=None,verbose=False,maxvel=1000): +def fit_lsq(spec,models=None,initpar=None,verbose=False,maxvel=[-1000,1000]): """ Least Squares fitting with forward modeling of the spectrum. @@ -1238,8 +1241,8 @@ def fit_lsq(spec,models=None,initpar=None,verbose=False,maxvel=1000): 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] = -maxvel - ubounds[3] = maxvel + lbounds[3] = maxvel[0] + ubounds[3] = maxvel[1] bounds = (lbounds, ubounds) initpar = np.minimum(ubounds-0.01,np.maximum(lbounds+0.01,initpar)) @@ -1380,7 +1383,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,maxvel=1000): +def multifit_lsq(speclist,modlist,initpar=None,verbose=False,maxvel=[-1000,1000]): """ Least Squares fitting with forward modeling of multiple spectra simultaneously. @@ -1438,8 +1441,8 @@ def multifit_lsq(speclist,modlist,initpar=None,verbose=False,maxvel=1000): 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:] = -maxvel - ubounds[3:] = maxvel + lbounds[3:] = maxvel[0] + ubounds[3:] = maxvel[1] bounds = (lbounds, ubounds) initpar = np.minimum(ubounds-0.01,np.maximum(lbounds+0.01,initpar)) @@ -1598,7 +1601,7 @@ def multispec_interp_jac(x,*argv): return out, lsmodel -def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=None,retpmodels=False,plot=False,tweak=True,usepeak=False,maxvel=1000) : +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. @@ -1650,7 +1653,8 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No t0 = time.time() # Make internal copy - spec = spectrum.copy() + #spec = spectrum.copy() + spec = copy.deepcopy(spectrum) # Step 1: Prepare the spectrum #----------------------------- @@ -1716,8 +1720,9 @@ 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+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) outstr2 = specxcorr(m.wave,m.flux,obs.flux,obs.err,maxlag) 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)]) @@ -1843,10 +1848,11 @@ 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,plot=False,tweak=True,maxlag=200,maxvel=500,usepeak=True) : +def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose=False,outdir=None,plot=False,tweak=True,maxlag=200,maxvel=[-500,500],usepeak=True) : """This fits a Cannon model to multiple spectra of the same star.""" # speclist is list of Spec1D objects. + pdb.set_trace() nspec = len(speclist) t0 = time.time() @@ -1884,7 +1890,8 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= modlist = [] bdlist = [] for i in range(len(speclist)): - spec = speclist[i].copy() + #spec = speclist[i].copy() + spec = copy.deepcopy(speclist[i]) if verbose is True: print('Fitting spectrum '+str(i+1)) print(speclist[i].filename) @@ -1991,7 +1998,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 @@ -2098,13 +2105,18 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= m = modlist[i](pars1,rv=vr1,wave=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(maxvel)/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) 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 - final['ccf'][i] = outstr['ccf'] - final['ccferr'][i] = outstr['ccferr'] + #final['x_ccf'][i] = (np.arange(nlag)-maxlag)*(np.log(m.wave[1,0])-np.log(m.wave[0,0]))*3.e5+vr1 + final['x_ccf'][i][0:nlag] = maxlag[0]+np.arange(nlag)*dwlog*3.e5+vr1 + final['ccf'][i][0:nlag] = outstr['ccf'] + final['ccferr'][i][0:nlag] = outstr['ccferr'] final['xcorr_vrel'][i] = rv+vr1 final['xcorr_vrelerr'][i] = outstr['vrelerr'] final['xcorr_vhelio'][i] = rv+vr1+info['bc'][i] @@ -2133,11 +2145,17 @@ def polynorm(flux,mask,order=4) : """ simple polynomial continuum """ x = np.arange(flux.shape[0]) - norder = flux.shape[1] cont = np.full_like(flux,1.) - for iorder in range(norder) : - gd=np.where(~mask[:,iorder] & np.isfinite(flux[:,iorder]))[0] - if len(gd) > norder : - coef = np.polyfit(x[gd],flux[gd,iorder],order) - cont[:,iorder] = np.polyval(coef,x) + 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 From 1fbbc42ab6fcb8feeffba8cf318342e301d1b0b9 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Tue, 30 Jun 2020 17:51:53 -0600 Subject: [PATCH 09/18] switch deepcopy back to copy method of spec1D --- doppler/rv.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index c515c27..98e8046 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -1653,8 +1653,7 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No t0 = time.time() # Make internal copy - #spec = spectrum.copy() - spec = copy.deepcopy(spectrum) + spec = spectrum.copy() # Step 1: Prepare the spectrum #----------------------------- @@ -1852,7 +1851,6 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= """This fits a Cannon model to multiple spectra of the same star.""" # speclist is list of Spec1D objects. - pdb.set_trace() nspec = len(speclist) t0 = time.time() @@ -1890,8 +1888,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= modlist = [] bdlist = [] for i in range(len(speclist)): - #spec = speclist[i].copy() - spec = copy.deepcopy(speclist[i]) + spec = speclist[i].copy() if verbose is True: print('Fitting spectrum '+str(i+1)) print(speclist[i].filename) From d5ebcbb75706ca5929a4ec5881c15b545856f823 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Sun, 5 Jul 2020 12:04:22 -0600 Subject: [PATCH 10/18] add final_xcorr routine to use for both single and multiple visits fix velocity scale for xcorr --- doppler/rv.py | 52 +++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 38 insertions(+), 14 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index 98e8046..6ac32f6 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -47,8 +47,8 @@ def xcorr_dtype(nlag): """Return the dtype for the xcorr structure""" 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)),("ccpeak",float),("ccpfwhm",float),("ccp_pars",(float,4)), + ("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 @@ -1957,7 +1957,22 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= sumstr['verr'] = info['vrelerr'][0] for key in ['teff','tefferr','logg','loggerr','feh','feherr','chisq'] : sumstr[key] = info[key][0] - return 0, sumstr, info, modlist, specmlist, time.time()-t0 + + pars1 = [info['teff'][0], info['logg'][0], info['feh'][0]] + vr1 = info['vrel'][0] + outstr = final_xcorr(specmlist[0],modlist[0],pars1,vr1,maxvel=maxvel,plot=plot) + 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'] + 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] + return sumstr, final, modlist, specmlist, time.time()-t0 # Step 2) find weighted stellar parameters @@ -2097,26 +2112,19 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= bmodel.append(m) # final cross-correlation - wavelog = utils.make_logwave_scale(sp.wave,vel=0.0) # get new wavelength solution - obs = sp.interp(wavelog) - m = modlist[i](pars1,rv=vr1,wave=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(maxvel)/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 = final_xcorr(sp,modlist[i],pars1,vr1,maxvel=maxvel,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 - final['x_ccf'][i][0:nlag] = maxlag[0]+np.arange(nlag)*dwlog*3.e5+vr1 + nlag = len(outstr['ccf'][0]) + final['x_ccf'][i][0:nlag] = outstr['ccvlag'] final['ccf'][i][0:nlag] = outstr['ccf'] final['ccferr'][i][0:nlag] = outstr['ccferr'] final['xcorr_vrel'][i] = rv+vr1 final['xcorr_vrelerr'][i] = outstr['vrelerr'] final['xcorr_vhelio'][i] = rv+vr1+info['bc'][i] + totchisq = np.sqrt(totchisq/totnpix) # Average values @@ -2138,6 +2146,22 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= return sumstr, final, bmodel, specmlist, time.time()-t0 +def final_xcorr(sp,model,pars,rv,maxvel=[-500,500],usepeak=False,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) + 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(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'] = maxvel[0]+np.arange(nlag)*dwlog*cspeed*np.log(10) + return outstr + def polynorm(flux,mask,order=4) : """ simple polynomial continuum """ From cef8ac0ea28ab71a263694b1030c456544a79ab4 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Fri, 10 Jul 2020 12:06:12 -0600 Subject: [PATCH 11/18] add maxvel to fit_lsq call change fine grid range handling of single visit returns to make multi-visit handling of failed frames in single visit --- doppler/rv.py | 41 +++++++++++++++++++++++++++++++---------- 1 file changed, 31 insertions(+), 10 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index 6ac32f6..fa15c00 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -1241,6 +1241,7 @@ def fit_lsq(spec,models=None,initpar=None,verbose=False,maxvel=[-1000,1000]): 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)) + print('fit_lsq: ', maxvel) lbounds[3] = maxvel[0] ubounds[3] = maxvel[1] bounds = (lbounds, ubounds) @@ -1758,13 +1759,16 @@ def fit(spectrum,models=None,verbose=False,mcmc=False,figfile=None,cornername=No # Get initial estimates 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(bestrv,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): @@ -1847,7 +1851,7 @@ 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,plot=False,tweak=True,maxlag=200,maxvel=[-500,500],usepeak=True) : +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. @@ -1903,12 +1907,29 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if (outdir is None) & (fdir != ''): figfile = fdir+'/'+figfile # Fit the spectrum try : - out, model, specm, pmodels = fit(spec,verbose=verbose,mcmc=mcmc,figfile=figfile,retpmodels=True,plot=plot,tweak=tweak,usepeak=usepeak,maxvel=maxvel) + print('calling fit: ', maxvel) + 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)) - print('removing from list ....') - bdlist.append(i) + # 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()) @@ -1961,6 +1982,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= pars1 = [info['teff'][0], info['logg'][0], info['feh'][0]] vr1 = info['vrel'][0] outstr = final_xcorr(specmlist[0],modlist[0],pars1,vr1,maxvel=maxvel,plot=plot) + m = [modlist[0](pars1,rv=vr1)] if usepeak : rv=outstr['vrel0'] else : rv=outstr['vrel'] if plot : pdb.set_trace() @@ -1972,7 +1994,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= final['xcorr_vrel'][0] = rv+vr1 final['xcorr_vrelerr'][0] = outstr['vrelerr'] final['xcorr_vhelio'][0] = rv+vr1+info['bc'][i] - return sumstr, final, modlist, specmlist, time.time()-t0 + return sumstr, final, m, specmlist, time.time()-t0 # Step 2) find weighted stellar parameters @@ -2014,7 +2036,6 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= initpar1 = np.zeros(3+nspec,float) initpar1[0:3] = wtpars[0:3] # the default is to use mean vhelio + BC for all visit spectra - print(len(initpar1),len(info)) initpar1[3:] = wtpars[3]-info['bc'] # vhelio = vrel + BC # Use the Vrel values from the initial fitting if they are accurate enough gdinit,ngdinit = dln.where(np.isfinite(info['vrel']) & (info['snr']>5)) @@ -2146,7 +2167,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= return sumstr, final, bmodel, specmlist, time.time()-t0 -def final_xcorr(sp,model,pars,rv,maxvel=[-500,500],usepeak=False,plot=False) : +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 @@ -2159,7 +2180,7 @@ def final_xcorr(sp,model,pars,rv,maxvel=[-500,500],usepeak=False,plot=False) : 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'] = maxvel[0]+np.arange(nlag)*dwlog*cspeed*np.log(10) + outstr['ccvlag'] = rv+(maxlag[0]+np.arange(nlag))*dwlog*cspeed*np.log(10) return outstr def polynorm(flux,mask,order=4) : From ce8bdda63936682a753dd6fbb9aba34981e9742a Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Mon, 13 Jul 2020 15:49:15 -0600 Subject: [PATCH 12/18] fix to normalize returned model for single input spectrum --- doppler/rv.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index fa15c00..8fa340e 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -1982,7 +1982,9 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= pars1 = [info['teff'][0], info['logg'][0], info['feh'][0]] vr1 = info['vrel'][0] outstr = final_xcorr(specmlist[0],modlist[0],pars1,vr1,maxvel=maxvel,plot=plot) - m = [modlist[0](pars1,rv=vr1)] + 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() @@ -1994,7 +1996,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= final['xcorr_vrel'][0] = rv+vr1 final['xcorr_vrelerr'][0] = outstr['vrelerr'] final['xcorr_vhelio'][0] = rv+vr1+info['bc'][i] - return sumstr, final, m, specmlist, time.time()-t0 + return sumstr, final, [m], specmlist, time.time()-t0 # Step 2) find weighted stellar parameters From 0e08766bdc4e26ce2812df3afc69281646f0d35e Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Sat, 18 Jul 2020 18:28:28 -0600 Subject: [PATCH 13/18] add badval argument for reader --- doppler/__init__.py | 4 ++-- doppler/reader.py | 18 +++++++++--------- 2 files changed, 11 insertions(+), 11 deletions(-) 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/reader.py b/doppler/reader.py index f78e6bb..f2fd4a9 100644 --- a/doppler/reader.py +++ b/doppler/reader.py @@ -26,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. @@ -82,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 @@ -99,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. @@ -168,7 +168,7 @@ def apvisit(filename): # 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,20735)!=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 @@ -206,7 +206,7 @@ def apvisit(filename): # Load APOGEE apStar/asStar spectra -def apstar(filename): +def apstar(filename,badval=20735): """ Read an SDSS APOGEE apStar spectrum. @@ -290,7 +290,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 @@ -318,7 +318,7 @@ def apstar(filename): # Load SDSS BOSS spectra -def boss(filename): +def boss(filename,badval=0): """ Read a SDSS BOSS spectrum. @@ -386,7 +386,7 @@ def boss(filename): # Load SDSS MaStar spectra -def mastar(filename): +def mastar(filename,badval=0): """ Read a SDSS MaStar spectrum. @@ -450,7 +450,7 @@ def mastar(filename): # Load IRAF-style spectra -def iraf(filename): +def iraf(filename,badval=0): """ Read an IRAF-style spectrum. From faf6f90b11c0de68079c5e307cb12d89d5056045 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Sat, 18 Jul 2020 18:28:44 -0600 Subject: [PATCH 14/18] remove a print --- doppler/rv.py | 1 - 1 file changed, 1 deletion(-) diff --git a/doppler/rv.py b/doppler/rv.py index 8fa340e..9bbf84c 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -1907,7 +1907,6 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if (outdir is None) & (fdir != ''): figfile = fdir+'/'+figfile # Fit the spectrum try : - print('calling fit: ', maxvel) out, model, specm, pmodels = \ fit(spec,verbose=verbose,mcmc=mcmc,figfile=figfile,retpmodels=True, plot=plot,tweak=tweak,usepeak=usepeak,maxvel=maxvel) From ecbe9438791d13e93a90a51943effdde84f6b501 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Thu, 23 Jul 2020 14:15:55 -0600 Subject: [PATCH 15/18] small bug fixes add some exception trapping --- doppler/rv.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index 9bbf84c..693e7be 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -596,8 +596,8 @@ def specxcorr(wave=None,tempspec=None,obsspec=None,obserr=None,maxlag=[-200,200] 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) @@ -1880,7 +1880,9 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= # 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: @@ -2070,6 +2072,8 @@ 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 @@ -2173,6 +2177,8 @@ def final_xcorr(sp,model,pars,rv,maxvel=[-500,500],plot=False) : """ 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 From 3973ce0aaed762d8fc1efd74dca0116b7648d326 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Thu, 6 Aug 2020 12:38:11 -0600 Subject: [PATCH 16/18] comment out sky masking stuff in favor of bitmask --- doppler/reader.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/doppler/reader.py b/doppler/reader.py index f2fd4a9..a84beff 100644 --- a/doppler/reader.py +++ b/doppler/reader.py @@ -121,7 +121,7 @@ def apvisit(filename,badval=20735): """ 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 @@ -170,21 +170,22 @@ def apvisit(filename,badval=20735): #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 - #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.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 - #plt.plot(spec.wave[:,i],sky) - #plt.plot(spec.wave[:,i],nsky*medsky2) - #plt.plot(spec.wave[:,i],spec.flux[:,i]) - #plt.draw() - #pdb.set_trace() + # 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): From 2aebbdb6e63c82ab94ca3039ed50ec669cf8d01e Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Fri, 14 Aug 2020 18:28:30 -0600 Subject: [PATCH 17/18] driver routine to make cannon models (untested) --- doppler/cannongrid_cubic_norm_training.py | 88 +++++++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100644 doppler/cannongrid_cubic_norm_training.py 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') + + From 5fd76ed7bfa86ab5546f224f03c682daf9711ba4 Mon Sep 17 00:00:00 2001 From: Jon Holtzman Date: Wed, 4 Nov 2020 13:19:23 -0700 Subject: [PATCH 18/18] add returning CCF FWHM and autocorrelation FWHM fix velocity array for CCF add BC to CCF velocity array calculate full CCF, but just look for peak in input range --- doppler/rv.py | 33 +++++++++++++++++++++++---------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/doppler/rv.py b/doppler/rv.py index 693e7be..e73a624 100644 --- a/doppler/rv.py +++ b/doppler/rv.py @@ -1865,7 +1865,8 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= 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),('x_ccf',(float,nlag)),('ccf',(float,nlag)), - ('ccferr',(float,nlag)),('xcorr_vrel',float),('xcorr_vrelerr',float),('xcorr_vhelio',float)]) + ('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): @@ -1982,7 +1983,7 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= pars1 = [info['teff'][0], info['logg'][0], info['feh'][0]] vr1 = info['vrel'][0] - outstr = final_xcorr(specmlist[0],modlist[0],pars1,vr1,maxvel=maxvel,plot=plot) + 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 @@ -1991,12 +1992,14 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= if plot : pdb.set_trace() final = info.copy() nlag = len(outstr['ccf'][0]) - final['x_ccf'][0][0:nlag] = outstr['ccvlag'] + 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 @@ -2138,18 +2141,27 @@ def jointfit(speclist,models=None,mcmc=False,snrcut=10.0,saveplot=False,verbose= bmodel.append(m) # final cross-correlation - outstr = final_xcorr(sp,modlist[i],pars1,vr1,maxvel=maxvel,plot=plot) + #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]) - final['x_ccf'][i][0:nlag] = outstr['ccvlag'] + # 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'] - final['xcorr_vrel'][i] = rv+vr1 + # 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['xcorr_vhelio'][i] = rv+vr1+info['bc'][i] + final['ccpfwhm'][i] = outstr['ccpfwhm'] + final['autofwhm'][i] = model_outstr['ccpfwhm'] totchisq = np.sqrt(totchisq/totnpix) @@ -2183,12 +2195,13 @@ def final_xcorr(sp,model,pars,rv,maxvel=[-500,500],plot=False) : 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(maxvel-rv)/cspeed)/dwlog).astype(int) + 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+(maxlag[0]+np.arange(nlag))*dwlog*cspeed*np.log(10) - return outstr + 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