Notebook Version 2: This notebook has been updated to reflect changes made to the package DifferentialMobilityAnalyzers.jl to work with the Julia v1 series (tested with Julia 1.1.0). To read the original supplement published with the paper please switch to v1.0.0 of the package DifferentialMobilityAnalyzers.jl and/or download the virtual machine on zenodo.org which contains a complete installation that works with Julia 0.6.4
This notebook demonstrates how to invert a size distribution from a measured noisy response function. The notebook is a supplement to the manuscript
Petters, M. D. (2018) A language to simplify computation of differential mobility analyzer response functions, Aerosol Science & Technology.
Let R be a measured response vector from a DMA. What is the measured true size distribution?
Figure 1. Convolution of true size distribution resulting in the response vector.For situations where instrument noise is neglible, $\epsilon_i \approx 0$ and
Routines for regularization are contained in regularization.jl. The method reproduced below is taken from Hansen (2000). The inverse of 𝕣 is found as follows. The residual norm and size of the regularized system are defined as
using Plots, Plots.PlotMeasures, Distributions, DifferentialMobilityAnalyzers, Random, LinearAlgebra, Printf
plotlyjs();
t,p = 295.15, 1e5 # Temperature [K], Pressure [Pa]
qsa,β = 1.66e-5, 1/5 # Qsample [m3 s-1], Sample-to-sheath ratio,
r₁,r₂,l = 9.37e-3,1.961e-2,0.44369 # DMA geometry [m]
leff = 13.0 # DMA effective diffusion length [m]
m = 3 # Upper number of charges
Λ = DMAconfig(t,p,qsa,qsa/β,r₁,r₂,l,leff,:-,m,:cylindrical) # Specify DMA with negative polarity
bins,z₁,z₂ = 128, dtoz(Λ,1000e-9), dtoz(Λ,10e-9) # bins, upper, lower mobility limit
δ = setupDMA(Λ, z₁, z₂, bins); # Compute matrices
The example provided here is for an assumed size distribution. From that distribution, the response function is computed and Poisson counting noise superimposed. The inversion algorithm is then applied and the result compared to the original input distribution.
figure("Nimbus Sans L", 2, 3, 2, 8)
𝕟 = DMALognormalDistribution([[400, 30, 1.2],[500, 110, 1.7]], δ)
p = plot(𝕟.Dp, 𝕟.S, xaxis = :log10, xticks = [10, 100, 1000], lt = :steppre, right_margin = 35px,
legend = :none, ylabel = "dN/dlnD (cm⁻³)", xlim = (8,1000), color = :black,
xlabel = "Mobility Diameter (nm)", fmt = :svg)
Figure 2. Lognormal aerosol size distribution 𝕟, represented as dN/dlnD. This is the assumed mobility size distribution from which the response vector is computed.
The response size distribution 𝕣 is computed as
Random.seed!(703) # random number seed; fix for repeatable runs
tscan = 120 # SMPS scan time [s]
Qcpc = 16.66 # CPC flow rate [cm³ s⁻¹]. 16.66 cm³ s⁻¹ = 1 L min⁻¹
t = tscan./bins # time in each bin
𝕣 = δ.𝐀*𝕟; # DMA response function
c = 𝕣.N*Qcpc*t; # number of counts in each bin
R = Float64[] # Construct a noisy response function 𝕣+ϵ, where ϵ is counting uncertainty
for i = c
f = rand(Poisson(i),1) # draw Poisson random number with mean = counts
push!(R,f[1]/(Qcpc*t)) # convert back to concentration
end
p = plot(𝕣.Dp, R, xaxis = :log10, xticks = [10, 100, 1000], lt = :steppre, label = "𝕣 = 𝐀*𝕟+ϵ",
ylabel = "Concentration (cm⁻³)", xlim = (8,1000), color = color = RGBA(0,0,0.8,1),
xlabel = "Apparent +1 Mobility Diam. (nm)", left_margin = 25px)
p = plot!(𝕣.Dp, 𝕣.N, label = "𝕣 = 𝐀*𝕟", color = :black, ls = :dash, fmt = :svg)
Figure 2. Response function for system with no noise (𝕣 = A𝕟) and system with counting uncertainty (𝕣 = A𝕟+ϵ) represented as measured concentration. The counting uncertainty depends on the total number concentration, the CPC flow rate and the time spent in each bin, which is determined by the scan rate.
The matrices $\bf{A}$ and $\bf{S}$ are defined in the module DifferentialMobilityAnalyzer.jl. The identity matrix is $\mathbf{I}$ = eye(n), where n is the number of bins. Calculations can be performed on either the number concentration or dN/dlnD vectors. Here calculations are performed on the number concentration vector R and the dN/dlnD and resulting 𝕟ᵢₙᵥ is assembled manually at the end of the computation. This is to avoid unnecessary complexity of the code in regularization.jl.
setupRegularization(A,I,R,Ni) defines the system of equations in terms of the matrices A, and I, the response R and the initial guess Ni.
lcurve(λ₁,λ₂;n=200) computes the L₁ and L₂ over the interval [λ₁,λ₂] for n=200 points. The function discretizes [λ₁,λ₂] in logspace and returns the array of λ-values (λs) and index where the curvature is maximum (ii). The function is only used for smooth visulation of the L-curve and not needed in the actual inversion.
lcorner(λ₁,λ₂;n=10,r=3) computes the L₁ and L₂ over the interval [λ₁,λ₂] for n=10 points calling lcurve(). The grid is recursively adapted three times to zoom into the optimal value, resulting in 30 total inversions. The optimal lambda value (λopt) is returned
(reginv(λopt, r = :Nλ))[1] return the inverted number vector for λopt. reginv() returns a vector of inversions so the first element is extracted. For low number concentration, the inversion can return negative numbers. The clean function zeros negative values in the vector.
Finally, the regularized size distribution 𝕟ᵢₙᵥ is reconstructed and actual vs. inverted number concentration are compared. The time required to compute lcurve and lcorner is reported. For n = 128 bins, λopt can be obtained within ~0.5s
λ₁,λ₂ = 1e-3, 1e1 # bounds [λ₁,λ₂] within which the optimal distribution lies
eyeM = Matrix{Float64}(I, bins, bins)
setupRegularization(δ.𝐀,eyeM,R,inv(δ.𝐒)*R) # setup the system of equations to regularize
@time L1,L2,λs,ii = lcurve(λ₁,λ₂;n=200) # compute the L-curve for plotting only
@time λopt = lcorner(λ₁,λ₂;n=10,r=3) # compute the optimal λ using recursive algorithn
N = clean((reginv(λopt, r = :Nλ))[1]) # find the inverted size distribution
𝕟ᵢₙᵥ= SizeDistribution([],𝕟.De,𝕟.Dp,𝕟.ΔlnD,N./𝕟.ΔlnD,N,:regularized) # store as AerosolSizeDistribution
@printf("Input number concentration %i\n", sum(𝕟.N))
@printf("Inverted number concentration %i\n", sum(𝕟ᵢₙᵥ.N))
8.032454 seconds (4.28 M allocations: 2.109 GiB, 4.14% gc time) 1.515819 seconds (113.09 k allocations: 394.505 MiB, 3.58% gc time) Input number concentration 900 Inverted number concentration 894
figure("Nimbus Sans L", 2, 6, 4.5, 10)
p1 = plot(L1, L2, xaxis = :log10, yaxis = :log10, xlabel = "𝐀*(𝕟-𝕣)", ylabel = "𝐈*(𝕟-𝕟ᵢ)", color = :black,
label = "L-curve"*(@sprintf(" λ ∈ [%.0e %.0e]", λ₁,λ₂)), bottom_margin = 50px)
p1 = plot!([L1[ii], L1[ii]], [L2[ii], L2[ii]], marker = :square, color = RGBA(0.8,0,0,0.5),
label = "Optimum λ ="*(@sprintf("%.1e", λopt)),lw = 0, ms = 4)
p2 = plot(𝕟.Dp, 𝕟.S, xaxis = :log10, xticks = [10, 100, 1000],left_margin = 65px, ylabel = "dN/dlnD (cm⁻³)",
label = "𝕟 = 𝓛𝓝 (N,Dₘ,σ)", xlabel = "Mobility Diameter (nm)", xlim = (8,1000), ls = :dash,
lw = 3, color = :black)
p2 = plot!(𝕟ᵢₙᵥ.Dp, 𝕟ᵢₙᵥ.S, color = RGBA(0.8,0,0,1), lt = :steppre,
label = "N<sub>inv</sub>=(𝐀ᵀ𝐀+λ²𝐈)⁻¹(𝐀ᵀR-λ²𝐒⁻¹R)", xlim = (8,1000))
p3 = plot(𝕣.Dp, R, xaxis = :log10, xticks = [10, 100, 1000], lt = :steppre,
label = "𝕣 = 𝐀*𝕟+ϵ", ylabel = "Concentration (cm⁻³)", xlim = (8,1000),
color = RGBA(0,0,0.8,1), xlabel = "Apparent +1 Mobility Diam. (nm)", left_margin = 65px)
plot(p1,p2,p3, layout = (l = @layout [a; b c]), right_margin = 0px, legend=:right,
top_margin = 15px, fmt = :svg)
Figure 3. Top: L-curve of $L_1$ vs. $L_2$ as defined in Block 4. The L-curve is computed from the interval [λ₁,λ₂] using the lcurve() function. The optimum $\lambda_{opt}$ is determined from the interative algorithm lcorner(). Bottom left: dashed line is the input size distribution (as in Figure 1). Red line is the regularized inverse of the noisy response function shown in the bottom right (as in Figure 2).
Figure 3 shows that the algorithm is able to correctly invert noisy data. This notebook can be used to test the algorithm over a wide range of conditions, including DMA flow ratios, input aerosol size distribution, scan rates, and bin density. Semi-systematic testing with this notebook can be used to convince oneself that the method is robust. The optimum regularization parameter $\lambda_{opt}$ (and thus solution) can be found within ~0.5-1s for a 128 bin grid (depending on computer speed and number of processes running).
Hansen, P. C. (2000) The L-Curve and its Use in the Numerical Treatment of Inverse Problems, in Computational Inverse Problems in Electrocardiology, ed. P. Johnston, Advances in Computational Bioengineering, 119-142, WIT Press.
Kandlikar, M., & Ramachandran, G. (1999) Inverse methods for analysing aerosol spectrometer measurements: A critical review. Journal of Aerosol Science, 30(4), 413-437, DOI: 10.1016/S0021-8502(98)00066-4.
Reineking, A. & J. Porstendörfer (1986) Measurements of Particle Loss Functions in a Differential Mobility Analyzer (TSI, Model 3071) for Different Flow Rates, Aerosol Science and Technology, 5:4, 483-486, DOI: 10.1080/02786828608959112.
Talukdar, Suddha S. & Mark T. Swihart (2010) An Improved Data Inversion Program for Obtaining Aerosol Size Distributions from Scanning Differential Mobility Analyzer Data, Aerosol Science and Technology, 37:2, 145-161, DOI: 10.1080/02786820300952.
Twomey, S. (1963) On the numerical solution of Fredholm integral equations of the first kind by inversion of the linear system produced by quadrature, Journal of the ACM, 19(1963), 97–101, DOI:0.1145/321150.321157.