Richardson-Lucy example

This file contains a complete example for the application of the Richardson-Lucy method for an actual data set. It starts with reading and pre-processing data files, calculating the RichLucy estimate and associated uncertainties.

Import of modules

In [16]:
import numpy as np
import RichLucy as RL
%pylab inline
Populating the interactive namespace from numpy and matplotlib

Loading and pre-processing data

In [2]:
bdata = RL.LoadFromFile("Bandpass_Spaltbreite_2mm_620nm-645nm.csv")
Mdata = RL.LoadFromFile("Kamera_2mm_Spaltbreite.csv")
In [3]:
lb,b,lM,M = RL.fix_raw_inputs(bdata,Mdata)[:4]
Result of bandpass check:
wavelength scale values are strictly monotonically increasing?	'True'
wavelength scale values are spaced equidistantly?		'False'
bandpass is normalized?						'False (integral = 2.7486e-09)'
wavelength scale axis (minimum) spacing is:			 0.0485
Result of checkSpectrum:
wavelength scale values are strictly monotonically increasing?	'True'
wavelength scale values are spaced equidistantly?		'False'
wavelength scale (minimal) spacing: 					 1.9918
Fixing of input data:


   Check for equidistant wavelength scale axis 
Wavelength scale is not equidistant.New wavelength scale step size is 0.0485. Interpolating...
Interpolation resulted in 0 NaN(s) and 0 negative value(s).


   Check for equidistant wavelength scale axis
Wavelength step size in measured spectrum requires interpolation.New step size is 0.0485. Interpolating...
Interpolation resulted in 0 NaN(s) and 0 negative value(s).


Fixing of input data finished


In [4]:
delta = lb[1]-lb[0]
In [46]:
figure(0,figsize=(10,6))
plot(bdata[:,0],bdata[:,1],"ro",label="raw data")
xlabel("wavelength / nm"); xlim(625,640)
ylabel("raw bandpass data amplitude",color="r",fontsize=16)
gca().twinx().plot(lb,b,"b",label="interpolated and scaled data")
gca().set_ylabel("interpolated and scaled bandpass amplitude",color="b",fontsize=16)
gca().set_xlim(625,640)
title("bandpass data")
Out[46]:
<matplotlib.text.Text at 0x13f80b10>

Note that the blue curve (interpolated bandpass) does not match the red dots (raw data) due to normalization of the bandpass function such that \(\int b(\lambda) d\lambda = 1\)

In [52]:
figure(101,figsize=(10,6))
plot(bdata[:,0],bdata[:,2],"o-")
title("uncertainties assoc. with bandpass data")
xlabel("wavelength / nm")
xlim(625,640)
Out[52]:
(625, 640)
In [48]:
figure(0,figsize=(10,6))
plot(Mdata[:,0],Mdata[:,1],"ro",label="raw data")
xlabel("wavelength / nm"); xlim(625,640)
ylabel("raw meas. spectrum amplitude",color="r",fontsize=16)
gca().twinx().plot(lM,M,"b",label="interpolated data")
gca().set_ylabel("interpolated spectrum amplitude",color="b",fontsize=16)
gca().set_xlim(630,670)
title("meas. spectrum data")
Out[48]:
<matplotlib.text.Text at 0x1422c370>
In [53]:
figure(101,figsize=(10,6))
plot(Mdata[:,0],Mdata[:,2],"o-")
title("uncertainties assoc. with meas. spectrum")
xlabel("wavelength / nm")
xlim(630,670)
Out[53]:
(630, 670)

Running the RichLucy method

Before calculating uncertainties it is advocated to run the estimation once for the best estimate of the input data to get an idea of suitable maximum number of iterations.

In [22]:
Shat,RMS,CURV = RL.RichLucy(b,delta,M,maxiter=5000)

 -------------- Richardson-Lucy algorithm version 1.95 --------------
Richardson-Lucy method calculating optimal stopping criterion using calcOptCurv
Step 1: Carry out 5000 iterations.
.................... Done.
Step 2: Calculating optimal stopping using calcOptCurv
Calculating optimal stopping iteration using curvature information
Calculating root-mean-squared value of estimation progress....................................................................................done.
Calculating curvature of rms progress.
potential stopping points are:
 [1, 16, 97, 757]

Optimal number of iterations = 97

It turns out that the chosen maximum number of iterations (5000) is much too high. Hence, for the Monte Carlo uncertainty calculations a much smaller value can be choosen which reduces computational burden significantly.

In [14]:
figure(1,figsize=(10,6))
plot(lM,M,label="measured spectrum")
plot(lM,Shat,label="estimated spectrum")
xlim(630,670)
xlabel("wavelength / nm"); ylabel("spectrum amplitude / au")
legend()
Out[14]:
<matplotlib.legend.Legend at 0xeac4830>
In [21]:
figure(2,figsize=(10,6))
loglog(range(1000),RMS)
xlabel("iteration"); ylabel("rms of estimation progress")
ax2 = gca().twinx()
ax2.semilogx(10**np.linspace(0,3,1000),CURV,'r')
ax2.set_ylabel("curvature of rms of progress",color="r")
draw()

Evaluation of uncertainties

In [26]:
raw_lb = bdata[:,0]
raw_b  = bdata[:,1]
raw_Ub = bdata[:,2]

raw_lM = Mdata[:,0]
raw_M  = Mdata[:,1]
raw_UM = Mdata[:,2]

Shat, UShat = RL.RichLucyUncertainty(raw_b,raw_Ub,raw_lb,raw_M,raw_UM,raw_lM,runs=1000,maxiter=200)


Starting 1000 Monte Carlo trials with 6 parallel workers
-----------------------------------------------------------
All done. Monte Carlo with 1000 runs took about 328.00 seconds
-----------------------------------------------------------

In [34]:
figure(4,figsize=(12,6))
plot(raw_lM,raw_M,"o-",label="measured spectrum")
plot(raw_lM,Shat,"o-",label="Monte Carlo mean of RL estimates")
legend()
xlim(630,670)
draw()
In [37]:
figure(5,figsize=(12,6))
plot(raw_lM,np.sqrt(np.diag(UShat)),"o-")
title("Standard uncertainties assoc. with RL estimate")
Out[37]:
<matplotlib.text.Text at 0x121965f0>