# Very slow code for simulation of optical emission (Raman) spectra

Can anyone help me optimise my code? It currently takes up to 2 hours to run a calculation!

I have written a simple code to simulate a type of optical emission spectrum (Raman scattering) for different molecules. These Raman spectra are characterised by a wavelength of emission and a signal intensity.

I am dealing with a complex molecule with many spectral lines (3.5 million), which can make the code described below become very slow.

My problem occurrs when converting a "stick spectrum" (i.e. pairs of wavelengths and intensities) to a "simulated spectrum" that takes into account the spectral resolution of the spectrometer. By applying the spectrometer resolution in the simulation, each discrete spectral line at a specific wavelength becomes a distribution of light intensity over a range of wavelengths. The lineshape of each line is a normal (Gaussian) distribution. Once the lineshape is calculated, the sum of intensities of every spectral line at a respective wavelengths is calculated in order to determine an overall spectrum.

This set of equations to convert from a line to a stick spectrum is given below:

deltav = 2.5; "resolution in cm-1";
gauss[{v0_, s_}] := {lambda, s*PDF[NormalDistribution[v0, deltav],lambda]}
lineShape[{v0_, s_}] := gauss[{v0, s}]
lineCalc[{v0_, s_}] := Piecewise[{{lineShape[{v0, s}], Abs[v0 - lambda] <= (mw = 10.)}, {{lambda, 0.}, Abs[v0 - lambda] > mw}}]

calc[data_] := Table[{lambda,Total[Table[lineCalc[data[[i]]],{i,1,Length[data]}][[All, 2]]]},{lambda, 2800., 3050., 1.}];


Some random sample data that the function "calc" would be applied to is given here:

{{2916.48, 1.5}, {2988.05, 0.207}, {2965.87, 1.01}, {2803.22,
0.0265}, {2825.75, 0.00431}, {2999.06, 0.414}, {2826.21,
0.0275}, {2849.13, 0.0109}, {2916.5, 2.67}, {3019.95,
0.38}, {2864.11, 1.28*10^-7}, {2967.56, 0.24}, {3011.59,
1.29*10^-6}, {3013.13, 0.00018}, {2810.96, 0.0125}, {2824.37,
0.00177}, {2997.86, 0.208}, {3044.54, 0.0000482}, {2966.43,
0.403}, {3013.11, 0.0000933}, {2855.13, 0.0822}, {2946.9,
0.636}, {2992.88, 0.0019}, {2830.16, 0.021}, {2850.83,
0.0321}, {2916.54, 2.89}, {3018.59, 0.297}, {2843.21,
1.62*10^-6}, {2945.26, 0.367}, {2991.21, 0.0000448}, {2993.14,
0.000439}, {2841., 0.0261}, {2849.25, 0.0203}, {2870.29,
0.000138}, {2875.75, 0.00179}, {3039.61, 0.0934}, {3041.33,
0.538}, {2817.9, 0.02}, {2823.35, 0.000255}, {2987.22,
0.558}, {2988.93, 0.346}, {3034.83, 1.31*10^-6}, {2945.32,
0.345}, {2947.03, 0.135}, {2992.93, 0.000899}, {2818.96,
0.0468}, {2836., 0.0071}, {2847.76, 0.00219}, {2851.7,
0.0503}, {2916.54, 4.33}, {3017.54, 0.033}, {3020.47,
0.704}, {2804.57, 0.0038}, {2816.33, 0.0089}, {2820.27,
0.00287}, {2885.11, 1.72*10^-6}, {2986.1, 0.197}, {2989.03,
0.305}, {3033.18, 2.49*10^-9}, {3034.88, 6.97*10^-7}, {2843.21,
6.13*10^-7}, {2944.2, 0.75}, {2947.13, 0.0115}, {2991.28,
0.0000198}, {2992.97, 0.000326}, {2812.04, 0.0683}, {2976.63,
2.25}, {2813.51, 0.177}, {2855.26, 0.242}, {2916.62,
9.73}, {3020.91, 2.32}, {2845.98, 0.0346}, {2865.93,
0.0107}, {2885.1, 0.000134}, {2890.41, 0.00367}, {3048.17,
0.0383}, {2811.77, 0.0168}, {2817.08, 0.00122}, {2974.84,
0.279}, {2978.85, 0.312}, {3025.45, 0.000066}, {2922.47,
0.562}, {2926.49, 0.0386}, {2973.09, 0.000915}, {2868.15,
0.0306}, {2878.07, 0.0386}, {2888.07, 0.000507}, {2905.87,
0.0000254}, {2906.56, 0.000179}, {2968.99, 1.96*10^-7}, {2815.76,
0.0114}, {2825.68, 0.00785}, {2835.68, 0.042}, {2853.47,
0.0472}, {2854.17, 0.0615}, {2916.59, 5.84}, {3016.72,
0.179}, {3018.36, 0.314}, {3020.55, 0.0701}, {2811.57, 0.0107}}

• The code returns errors because the conditions in the piecewise function don't cover the reals. Which of the inequalities should be non-strict? Commented Jan 9, 2018 at 14:27
• I hadn't noticed the inequality problem. The non-strict inequality can be the first one, I will update the code in the original post. Commented Jan 9, 2018 at 14:31

Piecewise is very, very slow but it can easily compiled. Have a look at this function:

cLineCalc = Quiet[Block[{a, lambda, v0, s},
With[{code = N[Piecewise[{{lineShape[{v0, s}][[2]],Abs[v0 - lambda] <= 10.}}, 0.]]},
Compile[{{lambda, _Real}, {a, _Real, 2}},
Block[{sum = 0., v0, s},
Do[
v0 = CompileGetElement[a, i, 1];
s = CompileGetElement[a, i, 2];
sum += code,
{i, 1, Length[a]}];
{lambda, sum}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
]
]]


With your example dataset stored in data, I obtain these timings:

lambdalist = Table[lambda, {lambda, 2800., 3050., 1.}];
a = calc[data]; // AbsoluteTiming // First
b = cLineCalc[lambdalist, data]; // AbsoluteTiming // First
a == b


0.161786

0.00013

True

• Having tested both of the suggested methods from Coolwater and Henrik I found different results on the machine I ran the code on. On my MacBook Coolwater's runs faster, but is still relatively slow (20+ mins for 30,000 lines) for very large data sets. Using my work PC, Henrik's solution is very fast (about 2.5 seconds for 3.7 million lines) where as Coolwater's is similarly slow as running on the mac. Thanks for all the help! Both proposed solutions are very helpful, both for this problem and in helping me to further understand Mathematica. Commented Jan 10, 2018 at 11:24
• Butters, that weird, I ran my tests on a MacBook Pro, so the OS is not the problem. Have you installed a C compiler on the MacBook? If not, you can, for example, install XCode. Alternatively, you can change CompilationTarget -> "C" to CompilationTarget -> "WVM". That's only a bit slower. Commented Jan 10, 2018 at 11:29
deltav = 2.5; mw = 10.;

With[{lambdaAll = Range[2800, 3050]},
Transpose[{lambdaAll, Lookup[Merge[#, Total], lambdaAll, 0]}] &[
lambdas = Range[Max[Ceiling[v0 - mw], 2800], Min[Floor[v0 + mw], 3050]]},
Transpose[data]]]]


Using MapThread an association is created for each data point. The keys are those lambdas between 2800 and 3050 that are within radius 10 from the data point x-value, i.e.
lambdas = Range[Max[Ceiling[v0 - mw], 2800], Min[Floor[v0 + mw], 3050]]

The values are s*PDF[NormalDistribution[v0, deltav], lambdas], and after MapThread the associations are merged to sum for each lambda.

Edit: Using Nearest runs about 4x faster:

deltav = 2.5; mw = 10.;

With[{lambdaAll = Range[2800, 3050]}, Transpose[{lambdaAll,
With[{near = Nearest[data[[All, 1]] -> {"Distance", "Index"}, lambdaAll, {∞, 10}]},
Total[TakeList[PDF[NormalDistribution[0, deltav], Join @@ near[[All, All, 1]]]
data[[Join @@ near[[All, All, 2]], 2]], Length /@ near], {2}]]}]]