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

Question

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}}

Answer 1

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 = Compile`GetElement[a, i, 1];
        s = Compile`GetElement[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

Answer 2

deltav = 2.5; mw = 10.;

With[{lambdaAll = Range[2800, 3050]},
  Transpose[{lambdaAll, Lookup[Merge[#, Total], lambdaAll, 0]}] &[
       MapThread[Function[{v0, s}, With[{
            lambdas = Range[Max[Ceiling[v0 - mw], 2800], Min[Floor[v0 + mw], 3050]]},
            AssociationThread[lambdas -> s*PDF[NormalDistribution[v0, deltav], lambdas]]]],
       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}]]}]]