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 numpy as np
import RichLucy as RL
%pylab inline
bdata = RL.LoadFromFile("Bandpass_Spaltbreite_2mm_620nm-645nm.csv")
Mdata = RL.LoadFromFile("Kamera_2mm_Spaltbreite.csv")
lb,b,lM,M = RL.fix_raw_inputs(bdata,Mdata)[:4]
delta = lb[1]-lb[0]
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")
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\)
figure(101,figsize=(10,6))
plot(bdata[:,0],bdata[:,2],"o-")
title("uncertainties assoc. with bandpass data")
xlabel("wavelength / nm")
xlim(625,640)
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")
figure(101,figsize=(10,6))
plot(Mdata[:,0],Mdata[:,2],"o-")
title("uncertainties assoc. with meas. spectrum")
xlabel("wavelength / nm")
xlim(630,670)
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.
Shat,RMS,CURV = RL.RichLucy(b,delta,M,maxiter=5000)
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.
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()
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()
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)
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()
figure(5,figsize=(12,6))
plot(raw_lM,np.sqrt(np.diag(UShat)),"o-")
title("Standard uncertainties assoc. with RL estimate")