请求优化拉曼光谱模拟代码：计算耗时过长需大幅提速

Hey there! Let's break down why your Raman spectrum simulation is dragging its feet and fix it—3.5 million spectral lines is a heavy load, but we can slash that 15-minute runtime down to something way closer to Mathematica's 10-second mark.

The Core Problem #

Your current approach is probably looping through each of the 3.5 million lines, generating a full Gaussian curve for every single one, then adding them all together. That's an O(N*M) operation (N = number of peaks, M = number of wavelength points), which gets brutally slow for large N. Mathematica uses optimized vectorized operations and FFT-based convolution under the hood—we can replicate that in Python.

Optimized Solutions (From Fastest to Most Flexible) #

1. Use FFT-Based Convolution (Best for Speed)

The entire process of turning stick spectra into resolution-broadened spectra is mathematically equivalent to convolving your stick spectrum with a Gaussian kernel. This cuts the time complexity to O((N+M)log(N+M)), which is night-and-day faster for large datasets.

Here's how to implement it:

import numpy as np
from scipy.signal import fftconvolve

# Your original parameters
vlow = 2800.0  # cm⁻¹
vhigh = 3050.0  # cm⁻¹
dv = 1  # cm⁻¹
resolution = 2.0  # cm⁻¹ (assuming this is FWHM)

# Target wavelength grid
wavelength = np.arange(vlow, vhigh, dv)
n_wave_points = len(wavelength)

# Assume you have these two arrays (3.5M entries each)
peak_wavenumbers = np.random.uniform(vlow, vhigh, 3_500_000)  # Example peaks
peak_intensities = np.random.rand(3_500_000)  # Example intensities

# Step 1: Build a stick spectrum array (map peaks to wavelength grid)
stick_spectrum = np.zeros(n_wave_points)
# Find which grid index each peak belongs to
peak_indices = np.searchsorted(wavelength, peak_wavenumbers)
# Clip indices to stay within valid range
peak_indices = np.clip(peak_indices, 0, n_wave_points - 1)
# Accumulate intensities into the stick spectrum
np.add.at(stick_spectrum, peak_indices, peak_intensities)

# Step 2: Create the Gaussian kernel (match spectrometer resolution)
# Convert FWHM to sigma: FWHM = 2*sqrt(2*ln2)*sigma ≈ 2.3548*sigma
sigma = resolution / (2 * np.sqrt(2 * np.log(2)))
# Build kernel with ±3σ coverage (captures ~99.7% of Gaussian area)
kernel_half_width = int(np.ceil(3 * sigma / dv))
kernel_wavenumbers = np.arange(-kernel_half_width, kernel_half_width + 1) * dv
gaussian_kernel = np.exp(-(kernel_wavenumbers ** 2) / (2 * sigma ** 2))
# Normalize kernel to preserve total intensity
gaussian_kernel /= np.sum(gaussian_kernel)

# Step 3: Run FFT convolution
simulated_spectrum = fftconvolve(stick_spectrum, gaussian_kernel, mode='same')

2. Vectorized Peak Calculation (No Loops)

If you need more control over individual peak shapes (e.g., non-uniform resolution), you can avoid Python loops entirely by vectorizing the Gaussian calculations:

# Vectorized Gaussian calculation for all peaks at once
# Calculate the distance from every peak to every wavelength point (broadcasted)
delta_v = wavelength[np.newaxis, :] - peak_wavenumbers[:, np.newaxis]
# Compute Gaussian intensities for all peaks
gaussians = peak_intensities[:, np.newaxis] * np.exp(-(delta_v ** 2) / (2 * sigma ** 2))
# Sum across all peaks to get the final spectrum
simulated_spectrum = np.sum(gaussians, axis=0)

Note: This uses more memory than convolution (it creates a 3.5M x 250 array), so only use it if you need per-peak customization.

3. Numba-Accelerated Loops (For Custom Logic)

If you have custom peak shapes that can't be easily vectorized, use Numba to compile your loop to machine code (with parallelization):

from numba import njit, prange

@njit(parallel=True)
def simulate_with_numba(wavelength, peak_wavenumbers, peak_intensities, sigma):
    n_wave = len(wavelength)
    n_peaks = len(peak_wavenumbers)
    spectrum = np.zeros(n_wave)
    
    # Parallelize over peaks
    for i in prange(n_peaks):
        v0 = peak_wavenumbers[i]
        I0 = peak_intensities[i]
        # Only calculate near the peak (skip irrelevant wavelength points)
        mask = np.abs(wavelength - v0) < 3 * sigma
        spectrum[mask] += I0 * np.exp(-((wavelength[mask] - v0) ** 2) / (2 * sigma ** 2))
    
    return spectrum

# Run the accelerated function
simulated_spectrum = simulate_with_numba(wavelength, peak_wavenumbers, peak_intensities, sigma)

Quick Additional Tips #

• Use float32 instead of float64: If your precision allows, downcast your arrays to np.float32 to cut memory usage in half and speed up calculations.
• Precompute the Gaussian kernel: Don't regenerate it every time you run a simulation—calculate it once and reuse it.
• Batch process if needed: If even FFT convolution uses too much memory, split your peak list into smaller batches and sum the results.

These changes should get your runtime down to seconds, matching Mathematica's performance. The key is ditching slow Python loops and leveraging optimized numerical operations!

内容的提问来源于stack exchange，提问作者Butters