SED fitting of a transient neighbourhood¶
In this jupyter notebook we go over how to integrate the SED of the transient neighbourhood and remove the contribution from the transient spectrum before fitting with BPASS SEDs using
For a kilonova the region around the transient is unlikely to be the region where the system was born so it is not all that relevant but you could do the same thing for a core collapse supernova in which case you transient is a lot more likely to be near its birth environment (although there are many caveats to that such as the type of CCSN and how distant the galaxy etc…).
# standard python libraries import matplotlib.pyplot as plt import numpy as np import pandas as pd import pickle # ppxf stuff for fitting from ppxf.ppxf import ppxf import ppxf.ppxf_util as util # hoki things for fitting from hoki import load from hoki import sedfitting as hsed plt.style.use('../tuto.mplstyle') c = 299792.458
1 - Separating Transient and Neighbourhood stars SEDs¶
1.1 - Neighbourhood SED¶
First let’s plot our galaxy and the situate our transient. I’ve already marked it and given you the pixel coordinate of the “centroid” (or as close as you can get it when the center is a 2 by 2 pixel section).
with open('data/cropped_NGC4993.pkl', 'rb') as pickle_file: cube = pickle.load(pickle_file)
f,ax = plt.subplots(nrows=1, ncols=1) ima = cube.sum(axis=0) ima.plot(scale='arcsinh', colorbar='v', ax=ax) ax.scatter(, , color='w', marker='*', label='AT 2017gfo') ax.legend()
<matplotlib.legend.Legend at 0x7f3e1c2019a0>
A very important thing to note is that the # NOTE THAT THE X AND Y ARE REVERESED IN THE DATA COMPARED TO MATPLOTLIB. So when plotting and handling data we need to keep in mind the axes are transposed and all is well. If you forget it gets confusing quick.
Here we’re going to integrate the SEDs for the pixels that are between 1.5 and 2.5 arcseconds away from the centroid of the transient. That follows Levant et al. 2017 for GW170817 and we didn’t see a need to change it. Play with it for your transient and pick the best option for your particular cases.
# NOTE THAT THE X AND Y ARE REVERESED IN THE DATA COMPARED TO MATPLOTLIB!!! cube2 = cube.subcube_circle_aperture(center=[119,49], # centroid y,x radius=2.5, # this is for the OUTER BOUNDARY unit_center=None, # the center is given in Pixels so the Units are None unit_radius='arcsec' # the radii are in arcsec which we need to state here ) cube15 = cube.subcube_circle_aperture(center=[119,49], radius=1.5, # this will be for the INNER BOUNDARY unit_center=None, unit_radius='arcsec')
f,ax = plt.subplots(nrows=1, ncols=1, figsize=(4,4)) ima = cube2.sum(axis=0) maplot = ima.plot(scale='arcsinh', colorbar='v', ax=ax) ax.set_title('2.5 Arcsec annulus')
Text(0.5, 1.0, '2.5 Arcsec annulus')
f,ax = plt.subplots(nrows=1, ncols=1, figsize=(4,4)) ima = cube15.sum(axis=0) maplot = ima.plot(scale='arcsinh', colorbar='v', ax=ax) ax.set_title('1.5 Arcsec annulus')
Text(0.5, 1.0, '1.5 Arcsec annulus')
Now all we need to do is integrate the SED in the 1.5 Arsec radius and subtract it from the integrated 2.5 arcsec flux and we’ll ahve recovered the flux from the pixels between the inner and the outer boundary.
Book keeping: We need a wavelength scale¶
Like in the voronoi binning notebook we create our 1D wavelength vector from the information we have in the headers.
# The first thing we need to do is create the wavelength scale. We've only got a starting # wavelength, a Delta lambda and a number of bins in the spectrum (3861 - see above) # The code below creates the wavelength array # minimum and maximum wavelength from the header wl_min = cube.data_header['CRVAL3'] wl_max = cube.data_header['CRVAL3']+cube.data_header['CD3_3']*cube. data.shape # creating the wavelength array WL = np.arange(cube.data_header['CRVAL3'], wl_max, cube.data_header['CD3_3'])
Now we go fetch our data¶
In order to go get our data and extract the flux and the noise to integrate we need to do a bit of reshapeing.
# we need the shape information shape_2 = cube2.data.data.shape shape_15 = cube15.data.data.shape shape_2, shape_15
((3681, 25, 25), (3681, 15, 15))
# then we reshape the data to be 2D data_2arsec = cube2.data.data.reshape(shape_2, shape_2*shape_2) data_15arsec = cube15.data.data.reshape(shape_15, shape_15*shape_15)
# axis = 0 is the spectrum; axis = 1 is the pixels all in a row data_2arsec.shape, data_15arsec.shape
((3681, 625), (3681, 225))
# now we sum our spectra over all pixels and we subtract the 1.5 arcsec region from the 2.5 arcsec region flux_annulus=data_2arsec.sum(axis=1) - data_15arsec.sum(axis=1)
# for the noise I did the same thing but in one line that's harder to look at # take your pick on how you want to do that stuff in your own work # also note I sqrt the variance to get the standard deviation noise_annulus=np.sqrt(cube2.var.data.sum(axis=1).sum(axis=1) - cube15.var.data.sum(axis=1).sum(axis=1))
# Now we need to know what area on the sky (in pixels) we covered so we can # turn our sum SED into an average SED arrea_annulus=shape_2*shape_2-shape_15*shape_15
Let’s plot it!¶
plt.figure(figsize=(5,3)) plt.step(WL, flux_annulus/np.nanmedian(flux_annulus), zorder=10, c='k', where='mid' ) plt.errorbar(WL, flux_annulus/np.nanmedian(flux_annulus), yerr=noise_annulus/np.nanmedian(flux_annulus) ) plt.title("Annulus 1.5' to 2.5'")
Text(0.5, 1.0, "Annulus 1.5' to 2.5'")
Let’s store this in a neat dataframe.
data_annulus = pd.DataFrame(np.array([WL, flux_annulus/np.nanmedian(flux_annulus), noise_annulus/np.nanmedian(flux_annulus)]).T, columns = ['WL', 'F', 'ERR']) data_annulus.drop(data_annulus.tail(1).index,inplace=True) # because the last one was messed up
3680 rows × 3 columns
1.2 Extracting the Transient data¶
We’re not done! We can now use this to extract our transient data and remove the galaxy background. This is actually not important for the galaxy SED fitting we’re doing in these tutorials but it is likely to be something that you will want to do when you are studying your shiny new transient.
Here we’re going to focus on the inner 1 arcsecond so a radius of 0.5 arcseconds. This value depends on your seeing when you’re observing as this is what will affect the full width half maximum.
# same as we did before but smaller cube1 = cube.subcube_circle_aperture(center=[119,49], radius=0.5, unit_center=None, unit_radius='arcsec') plt.figure(figsize=(3,3)) ima = cube1.sum(axis=0) ima.plot(scale='arcsinh', colorbar='v')
<matplotlib.image.AxesImage at 0x7f3e12b5e6a0>
ima.unmask() # we need to do this to access the data - because. # same as earlier we do some reshaping to get all our ducks in a row, i mean pixels shape_1 = cube1.data.data.shape data_1arsec = cube1.data.data.reshape(shape_1, shape_1*shape_1) data_1arsec.shape
# let's see what they all look like! plt.figure(figsize=(5,3)) for spec in data_1arsec.T: plt.plot(WL, spec, c='k', alpha=0.1) plt.xlim([4800,7500])
# let's sum our caterpillars transient_flux = data_1arsec.sum(axis=1) # spectrum noise_transient=np.sqrt(cube1.var.data.sum(axis=1).sum(axis=1)) # standard deviation
Let’s plot the transient (uncorrected)¶
plt.figure(figsize=(10,5)) plt.plot(WL, transient_flux/np.nanmedian(transient_flux), zorder=12) plt.errorbar(WL, transient_flux/np.nanmedian(transient_flux), yerr = noise_transient/np.nanmedian(transient_flux)) plt.xlim([4800,7500]) plt.ylim([0.8,1.2]) plt.title('Within 0.5 arsecond of the transient')
Text(0.5, 1.0, 'Within 0.5 arsecond of the transient')
Now we need to remove the galaxy background¶
To do this we need to first figure out how much “weight” the galactic SED has in our transient SED given the number of pixels we have summed together.
The weight is the fraction of the light in the transient spectra that comes from the galaxy based on the surface brightness of the neighbours in the annulus. We need to compare the areas of we’ve been looking at so far:
cube2.shape, cube15.shape, cube1.shape
((3681, 25, 25), (3681, 15, 15), (3681, 5, 5))
area_cube1 = (cube1.shape*cube1.shape) area_cube2 = (cube2.shape*cube2.shape) area_cube15= (cube15.shape*cube15.shape) # fraction we need to multiply to our galactic SED weight_gal_spec = area_cube1/(area_cube2-area_cube15)
# correcting the transient flux by removing (subtracting) the weighted galactic SED flux_gal_removed = (transient_flux/np.nanmedian(transient_flux))[:-1] - data_annulus.F*weight_gal_spec
Let’s plot the whole thing!¶
as it so happens here it didn’t make much of a difference to our transient but it might have! Also note the scales between the transient and the galactic neighrbourhood because everything is normalised to be around 1 (which we need for fitting later).
f, ax = plt.subplots(nrows=2, figsize=(10,7)) ax.plot(WL[:-1],flux_gal_removed, label='Galaxy removed') ax.plot(WL, transient_flux/np.nanmedian(transient_flux)-0.12, label='raw - 0.12') ax.set_xlim([4800,9000]) ax.set_ylim([0.7,1.2]) ax.legend() ax.set_title('AT 2017gfo') ax.plot(WL[:-1],data_annulus.F, label='Galaxy aroudn transient', c='k') ax.set_xlim([4800,9000]) ax.set_ylim([0.6,1.2]) ax.set_title('NGC 4993 around AT 2017gfo') ax.set_xlabel(r'Wavelength ($\AA$)')
Text(0.5, 0, 'Wavelength ($\\AA$)')
Let’s save our extracted spectra to a data frame¶
And we can then pickle the data frame (there are other ways to save dataframe data I just used to pickle them and I’ve left it as such here - make your prefered choice)
data = pd.DataFrame(columns=['WL', 'transient', 't_err', 'transient_raw', 'neighbours', 'n_err']) data.WL = WL[:-1] data.transient = flux_gal_removed data.transient_raw = transient_flux[:-1]/np.nanmedian(transient_flux) data.t_err=noise_transient[:-1]/np.nanmedian(transient_flux) data.neighbours=data_annulus.F data.n_err=data_annulus.ERR
data.metadata="WL: Wavelength in angstrom // Transient: Normalised transiant flux GALAXY REMOVED // Trasnient_raw: Normalised transient flux - galaxy NOT removed // neighbours: galaxy flux in anulus 1.5-2 arcsec around transient"
2 - Fitting the neighbourhood SED with BPASS!¶
Now that we have extracted our host environment spectral information let’s fit it to see what it is made of.
2.1 - Create ppxf templates¶
I recommend checking the ppxf manual for extra detail but I am going to walk you through the essentials
# just to make some of our code less long we'll store some columns in short variable names WL = data.WL.values F = data.neighbours.values noise = data.n_err.values
# Now we need to de-redshift z = 0.009783 recessional_vel = z*c WL = (WL - WL*z)
2.1.1 - BPASS models¶
To make templates we’ll need the BPASS models. To download the models you can check the BPASSv2.2.1 starter kit on Zenodo. Then we’ll need to be able to locate them so you’ll need to change the code below
BPASS_MODEL_PATH = '../../BPASS_hoki_dev/bpass_v2.2.1_imf135_300/'
2.1.3 - Log (natural not log10) rebin the wavelength¶
This is required by the
ppxf magic. Feel free to check out Capellari 2003 and 2017 if you want to read more on how
ppxf works and why the log rebin is a really neat mathematical trick ;)
# ppxf.util gives us the log_rebin function flux, wl_log_bin, velscale = util.log_rebin([WL, WL[-1]], F) # later we'll need the linear scale of our log rebinned spectra to plot our fits # this step ensures we have a linear wavelength scale for those fits with the # right number (and location) of wavelenght bins! wl_fits = np.exp(wl_log_bin)
2.1.4 - Creating templates for ppxf with KVN¶
KVN is a helper class in
hoki.sedfitting that does the dirty work of creating templates that can be fed to
ppxf - it also has utilities to decipher what
ppxf has found.
kvn=hsed.KVN() kvn.make_templates(BPASS_MODEL_PATH, fwhm_obs=1.25, wl_obs=WL, # Wavelength bins of the observations (not log rebinned) wl_range_obs=[WL, WL[-1]], # wavelength range of the observations velscale_obs=velscale, wl_range_padding=[-50,50], # number of extra Angstroms befrore and after the range z_list=['z010','z020','z030'], verbose=True, )
[---INFO---] TemplateMaker Starting [--RUNNING-] Initial Checks [--RUNNING-] Loading model spectrum [--RUNNING-] Calulating obs. velocity scale -- No dispersion [--RUNNING-] Calculating template wavelength (log rebin) and velocity scale [--RUNNING-] Calculating sigma [--RUNNING-] Instanciating container arrays [---INFO---] Using all ages from 6.0 to 10.2 (included) - log_age_cols now set [--RUNNING-] Compiling your templates Progress: |██████████████████████████████████████████████████| 100.00% Complete [-COMPLETE-] Templates compiled successfully
Padding: When creating the templates some rebinning has to occur and if we only load the exact wavelength range this can cause problems with the templates and the fitting later. Then you might ask why we crop our BPASS SEDs at all then, and the answer is simple: size! The BPASS SEDs cover everything from 1 A to 100,000 A. We’d be creating unnecessary work for our computer so it’s better to crop.
Where and how are the templates stored¶
The SED templates are BPASS SEDs that have been appropriately processed to be used by
ppxf for use on our specific data set. Unfortunately processing needs to happen on the templates for different observations with different properties like resolution element and wavelength range. That is why the helper class is so necessary and I can’t just give you a library of template to slap into your fit.
Let’s look at where to find our newly created templates:
(3766, 42, 3)
The dimensions are
[Wavelength, age, metallicity]. But obviously age = 42 and metallicity = 1 is not useful physically, we need to know what that corresponds to. The template properties are stored in the attribute
template_properties which has 2D
[N templates, properties] where the properties are numpy arrays containing [metallicity, age].
array([0.01, 6. ])
Tempalte 0 has a metallicity Z=0.010 and an age of 10^6 or 1 million years.
Note that here metallicity comes before the age dimension (in practice you should never need to think about it). The order the templates come in and the order the template properties are arranged is according to what is needed to talk to
ppxf: feed it the templates in the order it requries and reading its solutions in the order it provides.
This ordering thing is why we have a helper class so we can abstract it away.
2.2 - Seting up fit with PPXF¶
2.2.1 - Start guesses for the kinematics¶
Below we define two start guesses for the recession velocity and the dispersion velocity of the galaxy.
ppxf will fit the kinematics at the same time as looking for a best SED template combination so we need to give it a fighting chance.
Here we have already corrected for the average recession velocity so our initial guess is 0. It’s likely the final fitted guess will differ since galaxies rotated and our transient is on the outskirts so likely spining towards or away from us.
In terms of dispersion velocity you can make guesses using breadth of lines you know to be narrow, or you can eyeball it given the galaxy type and then tweak the number if the fits come out garbage. It doens’t need to be spot on to be decent usually.
start = [0, 180] #180km/s from previous fits
2.2.2- Good pixels¶
ppxf you can provide information on which pixels to take into account and which to ignore. Let’s see an example here (some more artificial than others but here for the pedagogy):
CUTTING BASED ON FLUX
Imagine if we have an emission line we don’t want to fit (NOTE: ppxf has some stuff to take into account emission features which you can read into in their manual). A stupid but often effective way to handle this is to crop anything that rises to high! Here we don’t have have any to worry about but when fitting all voronoi bins we do because the center of NGC 4993 is a LIER region (see Levan et al. 2017). Say we want to crop everything that rises above 2.5 (determined through visual anlysis in other spectra)
Our flux should be positive - if it’s not something is probably messed up so let’s cut everything that drops below 0.
CUTTING BASED ON WAVELENGTH
Let’s say we want to hide specific line because we thing our templates won’t be able to fit it well; here we pick the first absorption region; we can also do this by telling ppxf the pixels between that range are not good.
# 1 FLUXMAX =2.5 # 2 FLUXMIN = 0.0 # 3 WLMIN = 4875 WLMAX = 4975 # This looks convoluted but this is the format ppxf wants # A list of which wavelenght pixels are good (THEIR INDEX not their value) goodpixels = np.argwhere((flux < FLUXMAX) & (flux > FLUXMIN) # Flux condition & ((WL < WLMIN) | (WL > WLMAX)) # Wavelength condition ).T # Make it a 1D horizontal vector
2.2.3 - Fitting with ppxf¶
The method was originally described in 2004 but the most extensive paper on ppxf is Capellari 2017. When you download
ppxf there are also exmaple scripts provided with the code and what I use is adapted from those scripts. You can also consult the online PyPi manual.
A few quick start things to know:
moments: if set to 2, will only fit the kinematic parameters [recession velocity, dispersion]. We set it to 4 because we are also fitting for SFH and metallcity
degree: adding a polynomial can help if your spectrum is not super well calibrated or if there is a little bit of reddening are not fitting for with the reddening utilities (which I’m not here because there is not a lot of dust and I don’t care about the nature of the dust).
clean: this does sigma clipping and can help or hinder. See the Supplementary Information in Stevance et al. 2023 (UPDATE) to see how that played a role in our fitting of NGC 4993.
# From the ppxf manual: # Reference velocity in km/s (default 0). The input initial guess and the output velocities # are measured with respect to this velocity. This keyword can be used to account for # the difference in the starting wavelength of the templates and the galaxy spectrum. dv = c*np.log(kvn.wl_tem/kvn.wl_obs) ### FEEDING PPXF AND FITTING pp = ppxf(kvn.templates, # these are the BPASS SED templates flux, # observed flux vector noise, # the standard deviation vector velscale, # the velocity scale calculated by the log rebin ppxf utility start, # the start guesses for recession vel and dispersion goodpixels=goodpixels, # our goodpixels plot=False, # I set this to false because I made my own in hoki.sedfitting moments=4, # order of the Gauss-Hermit moment to fit (see above) degree=2, # order of the polynomial to fit on top of the SED mixture vsyst=dv, # see above clean=False, # whether to perform sigma clipping and itterative fitting lam=WL # the observed wavelength vector ) plt.ylim([0.4,1.2]) # FITTING THE RESULTS! res = hsed.ppxf_plot_fixed(pp, wl_fits=wl_fits, WL=WL, F=F)
Best Fit: Vel sigma h3 h4 comp. 0: -36 164 0.005 0.064 chi2/DOF: 23.05; degree = 2; mdegree = 0 method = capfit; Jac calls: 2; Func calls: 12; Status: 2 linear_method = lsq_box; Nonzero Templates (>0.1%): 5 / 126
2.3 - Saving and exploring the results with KVN¶
2.3.1 ``make_results`` and SED weights¶
KVN is not only a template creator, it also deciphers the results from ppxf to give them to us in a neat table and (as we’ll see later) we can save KVN as a pickle file so that our templates and fitting results are stored in one place and can be called back into other codes or shared with colleagues without needed to rerun the fitting.
First we need to ‘make the results’ by putting the
ppxf.ppxf.ppxf object we created (here called
pp) during fitting into our KVN
The results are stored in a
met is the metallicity (absolute, not with respect to solar Z=0.010 is 1% metals) *
age is the log10 age. That’s the way ages are binned in BPASS there is not other option *
weights are the weights ppxf gives those respective SEDs. They don’t necessarily add up to 1. If you want percentages or fractions you’ll need to rescale.
The only SED properties shown in the table are the SEDs ppxf actually uses for the fit.
2.3.2 - The
matching_ spectra and polynomials¶
When we make the results KVN also goes to fetch the SEDs used to match the observed SED so we can have a convenient access to them: *
kvn.matching_raw_spectra: the raw templates pre-fit (NOT THE RAW BPASS SEDs - they are processed for ppxf) *
kvn.matching_spectra: the spectra processed with the kinematics found by ppxf *
kvn.matching_apolynomial: the additive polynomial (if you included one) *
kvn.matching_mpolynomial: the multiplicative polynomial (if you included one)
MAJOR NOTE: The
kvn.matching_spectra need to be plotted with
wl_fits NOT with
kvn.wl_obs and it can’t be either any of the logarithming wavelength scales we have flowting around otherwise it will have weird offsets compared to NGC4993 (or your galaxy). There are a lot of different wavelength scales we deal with in this process, remembering which goes with what is one of the main beginner pit falls so if your spectra don’t line up this should be your
first port of call - what wavelength scale am I using and is it the right one?
CORRESPONDING WAVELENGTH SCALE FOR PLOTTING
Let’s have a look at what these various templates look like. We’ll also plot the fit.
# TEMPLATE (RAW POST-FIT) plt.plot(kvn.wl_tem, kvn.matching_raw_spectra, color='orange', alpha=0.8, label='A RAW template' ) # TEMPLATE (WITH KINEMATICS) plt.plot(wl_fits, kvn.matching_spectra, color='red', lw=1.5, label='A template kinematically processed' ) # TEMPLATE (WITH KINEMATICS) plt.plot(wl_fits, kvn.matching_apolynomial, color='maroon', lw=1.5, label='The added polynomial' ) # OBSERVED NGC 4993 plt.plot(wl_fits, flux, color='k', lw=1, alpha=0.9, label='NGC 4993' ) # FIT TO NGC 4993 plt.plot(wl_fits, kvn.ppxf.bestfit, color='blue', lw=0.9, label='The fit' ) plt.legend(fontsize=8, ncol=2)
<matplotlib.legend.Legend at 0x7f3e0aada8e0>
Note that to get the full fit that has put together the SEDs with the right weights and included the kinematics and the polynomial we call
kvn.ppxf.bestfit. Here we are directly looking into the ppxf object which contains all the results. The kinematically process tempaltes and polynomial can also be found in there alongside a lot more info (consult the manual for full details) such as the chi2 and a record of all the parameters we gave like the velscale.
Now a final point we need to address is What’s the age and metallicity of that template we just plotted?. I don’t know that it will necessarily match the order in which things are reported in the
results table so you need a way to link your matching templates to their template properties - for that we have
kvn.match_indices. As the name suggestes they are the indices identifying the template spectra in our full long list of teemplate spectra.
array([ 88, 92, 105, 111, 117])
So the properties are:
array([[0.02, 8.9 ], [0.03, 9. ], [0.01, 9.5 ], [0.01, 9.7 ], [0.01, 9.9 ]])
Which matches what’s in the result table at the minute, but if you had processed that table and re-ordered things (which you might do in some circumstances) you’d be in trouble so this is the most sure way to get your matching template properties.
2.3.3 Pickling KVN¶
The final thing we’re going to do is pickle our KVN object so we can save our work and then load it again next time we want to plot our fit or look at the star formation weights found by ppxf. This way we have everything in one place: the templates, the full fit solutions (with the ppxf object contained within KVN), and some convenient summaries of our results like the
To load the pickle file we need create a KVN object and then load
pwiz2 = hsed.KVN() pwiz2.load('data/pwiz_at2017gfo.pkl')
# just to show it's the same sum(pwiz2.ppxf.bestfit-kvn.ppxf.bestfit)
BONUS - Example of a bad fit¶
A common source of bad fit with ppxf is something wrong with the wavelenghts or velocities, either the templates weren’t created right because a log wavelength was given where a linear wavelength vector was expected (or vice versa),or the dv are wrong.
The result will be a spectrum that CLEARLY doesn’t have the right kinematics: it’s too broad. Let’s see an example below so you can identify that problem quickly and double check your steps.
pp = ppxf(kvn.templates, # these are the BPASS SED templates flux, # observed flux vector noise, # the standard deviation vector velscale, # the velocity scale calculated by the log rebin ppxf utility start, # the start guesses for recession vel and dispersion goodpixels=goodpixels, # our goodpixels plot=False, # I set this to false because I made my own in hoki.sedfitting moments=4, # order of the Gauss-Hermit moment to fit (see above) degree=2, # order of the polynomial to fit on top of the SED mixture vsyst=0, # CHANGED THIS TO 0 clean=False, # whether to perform sigma clipping and itterative fitting lam=WL # the observed wavelength vector ) plt.ylim([0.4,1.2]) # PLOTTING THE RESULTS! hsed.ppxf_plot_fixed(pp, wl_fits=wl_fits, WL=WL, F=F)
Best Fit: Vel sigma h3 h4 comp. 0: -2000 1000 -0.300 -0.081 chi2/DOF: 45.83; degree = 2; mdegree = 0 method = capfit; Jac calls: 10; Func calls: 52; Status: 2 linear_method = lsq_box; Nonzero Templates (>0.1%): 3 / 126
array([-0.11513668, -0.14540687, -0.12683344, ..., -0.11569206, -0.04607372, -0.03827896])