pythonscipyconvolutionastronomyspectra

Efficient method for generating combined spectra from N-body simulation data using Python (convolution)

I'm working with some N-body simulation data, and attempting to generate a combined spectrum from particles within a specific region (i.e. inside one 'pixel' or voronoi bin). Each particle has associated line-of-sight (LOS) velocities. My current approach involves Doppler shifting a stellar template spectrum (see top panel of figure below) for each particle and interpolating it back to original wavelenght points to sum up the fluxes. This method works, but it is a bit slow when applied to the full dataset of about 8 million particles.

Example of stellar template (obtained from Miles ). The template includes 4300 wavelength points. The lower panel represent the combined spectrum of 112,866 particles.

Example stellar template and combined spectrum

Current method :


# Variables used in this method
template_wavelength_ln # Wavelengths of the template spectra with even intervals in ln scale
template_flux_ln # Fluxes of this template spectra

LOSv  # Line-Of-Sight velocities of particles corresponding to certain 'pixel'
# --------------
# The coding method itself


wavelengths_shift = DPF.doppler_shift(np.e**template_wavelength_ln, LOSv.reshape(-1,1)) 
# This calculates the doppler shift of the template wavelengths for each particle. 
# Returns 2d array. Due to high particle number, I have used some chunking method here but I took it off here as its not relevant to the problem itself. (i think)

df_tot = pd.DataFrame({'lambda': np.e**template_wavelength_ln, 'Flux_tot':np.zeros_like(template_flux_ln), 'Flux_1':np.zeros_like(template_flux_ln),
'Flux_2':np.zeros_like(template_flux_ln)})

# The above is just data frame which would include the spectrum infos. Currently I'm doing this for two different stellar templates, thus Flux_1 and Flux_2, but for the sake of example we use only one template.

for wave in wavelenght_shift:
    flux_1 = np.interp(df_tot['lambda'], wave, template_flux_ln) #interpolate the shifted wavelenghts back to the original positions
    df_tot['Flux_1'] += flux_1
    # df_tot['Flux_tot'] += (flux_1 + flux_2)

# -----------------

And this generated the nice combined spectrum of the particles. However, its quite slow when the doing 8 million particles, mainly due to np.interp() I guess.

where DPF.doppler_shift is

def doppler_shift(wavelenghts, velocity, c=299792.458):
    ...
    lambda_shifted = wavelenghts * (1 + velocity/c)
    return lambda_shifted

Method 2 This method would use the advantage of convolution. Here is a link to a paper which describes how one could generate the spectrum of galaxy using a template and LOSVD (line of sight velocity distribution) Improving the full spectrum fitting method: accurate convolution with Gauss-Hermite functions by Cappellari equation 1.

So, what I have done so far: Obtained LOSVD of the particles in form of a Gauss-Hermite function. First I do a histogram of the LOS velocity particles, and then fit the Gauss-Hermite function to that histogram.

from scipy.optimize import curve_fit
# Same initial variables are used here as in the code above

# Initial guess for the Gauss-hermite parameters 
initial_guess = [LOSv.mean(), LOSv.std(), 0, 0] # V_mean, sigma, h3, h4

# Histogram binning
bin_width = 10 # km/s
bin_edges = np.arange(LOSv.min() , LOSv.max() + bin_width, bin_width)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

hist_density = np.histogram(LOSv, bins=bin_edges, density=True)[0]

params_fit_total, error_fit_total = curve_fit(DPF.gauss_hermite, bin_centers, hist_density, p0=initial_guess, maxfev=10000)

# ---------------
# where DPF.gauss_hermite is 
from scipy.special import hermite
def gauss_hermite(v, v_mean, sigma, h3, h4):
   ...
   w = (v - v_mean) / sigma
   H3 = hermite(3, monic=True)(w)
   H4 = hermite(4, monic=True)(w)
   # One could also use these, I think
   # H3 = (2*np.sqrt(2)*w**3 - 3*np.sqrt(2)*w)/np.sqrt(6)
   # H4 = (4*w**4 - 12*w**2+3)/np.sqrt(24)
   L = 1/(sigma*np.sqrt(2*np.pi)) * np.exp(-0.5*w**2) * (1+h3*H4 + h4*H4)
   return L

# and then the LOSVD function i.e.

LOSVD_func = lambda v : DPF.gauss_hermite(v, *params_fit_total)

Now this params_fit_total includes the v0, sigma, h3 and h4 parameters of the LOSVD. Below is plot of the histogram and the LOSVD, including the v0 sigma h3 and h4 values in the legend.

LOSVD histogram and Gauss-Hermite function

The challenge now is to generate the 'total' spectrum using the initial stellar template and this LOSVD function via convolution. However I have some difficulties on understanding how to do this, as I have not used convolution very much. I believe that the units of 'parameters' within the convolution factors have to match, so one might need to convert everything to either wavelengths or velocities. I remember that I have heard/read somewhere that "uniform sampling in the logarithm of the wavelength corresponds to a uniform sampling in velocity", which suggests a method for converting the units, but I'm not entirely sure how to implement this in practice.

Here are some equation algebra how I tried to justify the statement: Doppler shift

$\lambda_0 = \lambda_e (1+\frac{v}{c})$

$\rightarrow \lambda_0 - \lambda_e = \lambda_e \frac{v}{c}$

$\rightarrow \frac{\lambda_0-\lambda_e}{\lambda_e} = \frac{\Delta \lambda}{\lambda_e} = \frac{v}{c}$

which for small changes in $\lambda$, thus small change in $v$ yields into

$\Delta \ln(\lambda) = \frac{v}{c}$

meaning that small and uniform change in the ln scale indeed relates to uniform sampling in velocity.

Yet... I'm not completely sure how to implement this information to the convolution. I feel like this is close to the answer, but I have not got the convolution working...

I have tried something like

from scipy.signal import fftconvolve

step_h_log = wave_log[1] - wave_log[0]
c = 299792.458 # km/s

# step_h_log = 0.0001717861122...
# so this would mean velocity bin of 51.5 km/s (step_h_log * c)

vel_space = np.arange( np.min(LOSv.min()), np.min(LOSv.max())+ step_h_log*c, step_h_log*c  )

test_convo = fftconvolve( DPF.doppler_shift(np.e**wave_log, vel_space), LOSVD_func(vel_space), mode='same' )

but indeed this gives error as the shapes differs... I have tried different ways but none gave the correct output.

If one is interested where to obtain the stellar templates, one can get them i.e. from here http://research.iac.es/proyecto/miles/pages/webtools/tune-ssp-models.php. I am not sure if I could paste the .txt files here which contain the wavelengths and fluxes

Any ideas? I would greatly appreciate the help for both coding aspect and perhaps understanding how the convolution actually works in this context (unit conversion).

Solution

  • EDIT : This is my own reply/answer to my problem.

    I completely forgot to post this answer here once I got the convolution method to work. Thanks to a friend of mine for reminding me!

    Recap

    First, a quick recap: As mentioned in the original question, I have N-body simulation data of a simulated galaxy. I want to generate mock spectra using two different stellar templates from the MILES library (see the original question for details). Each particle is weighted to contribute to each of the two templates, with the total weight being unity (W_1 + W_2 = 1.0).

    Background: Gauss-Hermite function/parameters

    Figure 1 below shows the velocity distribution of the particles within a single Voronoi bin of the galaxy:

    Velocity distribution of particles within Voronoi Bin 296

    The fitted GH parameters in this case are:

    1. Mean velocity V_0
    2. Velocity dispersion \sigma
    3. Gauss-Hermite parameters h_3 and h_4, which can be calculated using the equations provided in various references

    One can calculate Gauss-Hermite parameter values using existing Python libraries (e.g., NumPy), but one can also find the equations on Wikipedia. Below are a few Python functions representing how I have implemented the Gauss-Hermite function:

    def fhn(x,n):
    """
    Physicists Hermite polynomials
    
    Parameters:
    x : array-like, the input values.
    n : int, the order of the Hermite polynomial.
    
    Returns:
    array-like, the values of the Hermite polynomial.
    Source:
        - Michele Cappellari: MNRAS 466, 798–811 (2017) Eq. 14. 
        - https://en.wikipedia.org/wiki/Hermite_polynomials
        - B. Neureiter : arXiv:2212.06173
        - K. Mehrgan : The mass composition of massive early-type galaxies
        - 
    """
    n_fak = math.factorial(n)
    norm = n_fak * 2**n
    k_up = int(n/2)
    fhn = 0
    for k in range(0,k_up+1):
        beta_kn = n-2*k
        k_fak = math.factorial(k)
        n_2k_fak = math.factorial(n-2*k)
        vz = (-1) ** k
        a_kn = vz * n_fak/k_fak/n_2k_fak * 2**beta_kn
        fhn = fhn + a_kn*x**beta_kn
    fhn = fhn/np.sqrt(norm)
    return fhn

def LOSVD_gauss_hermite(v, amp, v_mean, sigma, *coeff):
    """
    Compute the Gauss-Hermite series for given parameters.
    
    Parameters:
    v : array-like, velocity values.
    amp : float, amplitude.
    v_mean : float, mean velocity.
    sigma : float, velocity dispersion.
    coeff : tuple, Gauss-Hermite coefficients.
    
    Returns:
    array-like, the values of the Gauss-Hermite series.
    """
    y = (v - v_mean) / sigma
    gauss = amp/np.sqrt(2*np.pi) * np.exp(-y**2 /2)
    len_h_coeff = len(coeff)
    sum = 1
    for i in range(len_h_coeff):
        sum += coeff[i] * fhn(y, i+3)
    L = gauss/sigma * sum 
    return L

    These functions can be fitted to the data using the scipy.optimize.curve_fit function. The reason for introducing the Gauss-Hermite parameters is that they provide a way to double-check that the results are correct.

    Creating spectra

    Method 1: I explained this method in the original question, but briefly: We iterate through each particle, shift the stellar templates according to the line-of-sight velocity using the Doppler shift, and then sum the shifted spectra together. As a side note, the code for doing this has improved since the example given, and I can share that if needed.

    Method 2: This is the main part of the answer: Convolution method

    I used the idea from Michele Cappellari's 2016 paper, "Improving the Full Spectrum Fitting Method: Accurate Convolution with Gauss-Hermite Functions," specifically chapters 2.1 - 2.3 and 3.2.

    This is the function I came up with. It is quite lengthy, as it includes a few additional features that are handy for me, but let me explain the important steps verbally:

    1. Step 6.1 : Compute the GH parameter values for the distribution. This is then used to 'functionalize' the LOSVD of the particles (one could use, for example, an instrumental broadening function instead of this).

    2. Step 6.2 : Here, we compute the LOSVD of the particles and use it as the convolution kernel.

    3. Step 6.3 : The key step: convolution using scipy.signal.convolve(). The convolution is done using the fluxes of the stellar templates and the computed LOSVD kernel.

    That’s pretty much it; all the other steps are for bookkeeping and solving other values. Here is the function:

    def compute_spectra_convolution_method(LOSv, weight_list, wavelengths_list, fluxes_list,
                                       silence=False,
                                       return_GH=False, len_h_coeff=2, bin_width='auto',
                                       interpolate_to_linear=False):
    """
    USER COMMENT: This also assumes that each stellar template uses same wavelenght 'spacing' and 'range'.
                 - Usually case if one uses i.e MILES templates.
    
    Computes the convolved spectra for multiple templates using LOSVD modeled by Gauss-Hermite polynomials.

    Parameters:
    - LOSv : array-like, Line-Of-Sight velocities. (of particles/stars)
    - weight_list : list of arrays, weights for each template's contribution.
        -- Each element in the list corresponds to a template and contains weights that indicate how much each particle 
           contributes to that template.

    - wavelengths_list : list of arrays, The wavelengths corresponding to each template's spectrum.
    - fluxes_list : list of arrays, The flux (intensity) values for each template at the given wavelengths.
         
    - bin_width : float or 'auto'.  default='auto'
        The bin width for the velocity histogram and the LOSVD calculation. If 'auto', it is calculated from the 
        wavelength difference assuming a linear wavelength scale.

    - return_GH : bool, optional, default=False
        If True, the function also returns the computed Gauss-Hermite parameters along with their errors.
    
    - len_h_coeff : int, optional, default=2
        The number of Gauss-Hermite coefficients to consider beyond the Gaussian (default=2 corresponds to h3 and h4).

    - interpolate_to_linear : bool, optional, default=False
        If True, the computed spectra are interpolated onto a linear wavelenght scale. This is useful to keep scales constant in my case.
        NOTE: If interpolation to linear scale is desired, the first wavelengths_list should be
            in logarithmic (ln) form and the second should be in linear form which is used to interpolation.
    
    - Silence : bool, optional, default=False
        If True, the function prints information about the number of particles and templates.
            
    Returns:
    - df_tot : DataFrame
        A DataFrame containing the wavelengths ('lambda'), the total convolved flux ('Flux_tot'), and the convolved flux 
        for each template.

    - (Optionally) GH_params : DataFrame or dictionary
        If return_GH is True, this also returns a DataFrame containing the Gauss-Hermite parameters and their errors for each template.
   """
    # Step 1: Verify and prepare the wavelengths for correct velocity calculation.
    # Ensure the first wavelenghts_list is ln for correct velocity calculation.
    wavelengths = np.exp(wavelengths_list[0])
    differences = np.diff(wavelengths_list[0])
    if not np.allclose(differences, differences[0]):
        print('Warning: Wavelengths are not equally spaced in ln scale.')

    # Step 2: Calculate the bin width for the velocity histogram.
    # This step is crucial for accurate LOSVD computation and convolution.
    if bin_width == 'auto':
        bin_width = np.diff(wavelengths_list[0])[0] * 299792.458 # Converting to km/s ; (51.5 km/s) ; M. Capellari 2016 eq. 7 & 8
        print('Auto bin width:', bin_width)
        
    # Step 3: Print basic information if not silenced.
    # Determining the number of templates and particles for printing info
    NN_templates = len(fluxes_list)
    NN_particles = len(LOSv)
    if not silence: print(f'CONVOLUTION METHOD : Number of particles: {NN_particles}, NN_templates: {NN_templates}')
    
    # Step 4: Compute the Gauss-Hermite parameters, and save to dataframe if return_GH is True.
    if return_GH == True:
        GH_params, GH_error = compute_gauss_hermite_parameters(LOSv=LOSv,
                                            len_h_coeff=len_h_coeff, density=True,
                                            bin_width=bin_width)
        GH_params_data = []
        GH_params_data.append(['Total', *GH_params, *GH_error])
    
    # Step 5: Define the velocity range for the LOSVD kernel and prepare the output DataFrame.
    # The velocity range determines the extent of the convolution and should cover most of the velocity distribution.
    s4 = LOSv.std()*4 # 4-sigma range, should be more than enough 
    vel_range = np.arange(-s4, s4+bin_width, bin_width)  
    df_tot = pd.DataFrame({'lambda': wavelengths, 'Flux_tot': np.zeros_like(wavelengths_list[0])})

    # Step 6: Loop through each template to compute and convolve their spectra.
    # This is the core of the function, where each template's contribution is individually convolved with the LOSVD.
    for i_template in range(NN_templates):
        df_tot[f'Flux_{i_template+1}'] = np.zeros_like(wavelengths_list[0])

        # Step 6.1 : Compute GH parameters for each template, using the LOSVD and weights.
        GH_params, GH_error = compute_gauss_hermite_parameters(LOSv=LOSv,
                                                 len_h_coeff=len_h_coeff,
                                                 density=True,
                                                 bin_width=bin_width,
                                                 weights=weight_list[i_template] )

        # Step 6.2 :  Create the LOSVD kernel for convolution using the fitted velocity distribution.
        LOSVD_kernel = LOSVD_gauss_hermite(vel_range, *GH_params)
        LOSVD_kernel /= np.sum(LOSVD_kernel) # Normalize the kernel
        
        template_weight = np.sum(weight_list[i_template])

        # Step 6.3 : Convolute the template flux with the LOSVD kernel and add contributions to dataframe
        convoluted_flux = convolve(fluxes_list[i_template], LOSVD_kernel, mode='same', method='auto')*template_weight
        df_tot[f'Flux_{i_template+1}'] = convoluted_flux
        df_tot['Flux_tot'] += convoluted_flux

        # Store GH paramters if return_GH is True
        if return_GH:
            GH_params_data.append([f'T{i_template+1}', *GH_params, *GH_error])

    # Step 7: (Optional) Return Gauss-Hermite parameters in a DataFrame.
    if return_GH:
        columns = ['Label','amp', 'v0', 'sigma']
        error_columns = ['error_amp','error_v0', 'error_sigma']
        h_columns = [f'h{i}' for i in range(3, 3+len_h_coeff)]
        h_columns_error = [f'error_h{i}' for i in range(3, 3+len_h_coeff)]
        all_columns = columns + h_columns + error_columns + h_columns_error
        
        GH_params_df = pd.DataFrame(GH_params_data, columns=all_columns)
        
    # Step 8: (Optional) Interpolate to a linear wavelength scale if needed.
    # This is useful for ensuring that the output spectra are on a consistent wavelength grid/spacing.
    if interpolate_to_linear==True:
        df_tot = df_tot.rename(columns={'lambda':'lambda_ln'})
        df_tot.rename(columns=lambda x: x + '_ln' if 'Flux' in x else x, inplace=True)
        wavelengths_linear = wavelengths_list[1]
        
        for col in df_tot.columns[1:]:
            interpolated = np.interp(wavelengths_linear, df_tot['lambda_ln'], df_tot[col].values   )
            new_col_name = col[:-3]
            df_tot[new_col_name] = interpolated
        df_tot['lambda'] = wavelengths_linear

        # Just ordering the columns to nicer order
        column_order = ['lambda']
        non_ln_cols = column_order + [col for col in df_tot.columns if not col.endswith('_ln') and 'Flux' in col ]
        ln_cols = [col + '_ln' for col in non_ln_cols]
        column_order = non_ln_cols + ln_cols
        df_tot = df_tot[column_order]

    # Step 9: Return the final DataFrame, optionally including the GH parameters.
    if return_GH:
        return df_tot, GH_params_df
    else:
        return df_tot

    Results

    Below, in Figure 2, I show the spectra generated using both methods, interpolation and convolution. In terms of time (N_particles= 11,850, N_templates= 2):

    Comparison between interpolation and convolution spectras

    As can be seen, the results are remarkably similar (the largest difference I observed was on the order of 1%). However, as a disclaimer, the convolution method is not perfect. After producing the spectra with both methods and then giving them as input files into programs that solve LOSVD fittings (Gauss-Hermite), both methods return similar sigma and Hermite parameter values (h3, h4, etc.) but the line-of-sight velocity might vary. Based on my tests, I have noticed that the width of the convolution kernel plays a crucial role, as small changes in that are visible in the obtained line-of-sight velocities. This is still a mystery for me regarding how to accurately solve and improve this.