In [1]:
%matplotlib inline
import infer_structcol as ifs
import time
import os
from matplotlib import pyplot as plt
import seaborn as sns
The goal of this package is to infer the volume fraction, particle radius, and sample thickness of a colloidal structural color sample from its reflectance and/or transmittance spectrum. The generative model must be able to reproduce experimental reflectance and/or transmittance spectra from a set of parameters defining the sample. The foundation of the generative model is a stochastic multiple scattering algorithm included in a separate python package at https://github.com/manoharan-lab/structural-color. Here, we describe the multiple scattering model, our generative model (and its uncertainties), and our likelihood function.
These calculations are a part of the structural-color package (see above).
This model predicts the reflectance and transmittance spectra of structurally colored samples of thin films made from nanoscale spheres in a disordered packing. Past models have accounted for only a single scattering event in the sample [1], but this model uses a Monte Carlo approach to account for many scattering events, leading to more accurate results.
The multiple scattering model simulates the trajectories of individual photons as random walkers that move through a material by sampling both the step size and the direction from distributions that are calculated using the material parameters.
Each simulation sends in a user-specified number of photon packets at the surface of the sample film. To calculate the reflectance, the number of photon packets that exit the sample through the backscattering hemisphere is divided by the total number of photon packets sent into the sample. Analogously, transmittance is calculated by dividing the number of photon packets that exit the sample through the forward hemisphere by the total number of packets sent into the sample.
A model sample consists of an infinitely wide film of a finite thickness. The code also accounts for reflections due to the index contrast at the interface using the Fresnel equations. Refraction of the backscattered light exiting the sample is calculated using Snell's law, so that the angle of photons exiting the sample can be tracked in case the user is only interested in collecting a certain angle range of light.
The generative model involves corrections on the reflectance and/or transmittance spectrum produced by the calculations in the structural-color package (see above).
The goal of the generative model is to calculate reflectance and transmittance spectra from structural color samples defined by a set of parameters. The multiple scattering model described above generates theoretical spectra but these spectra do not correspond to experimental measurements. The full generative model incorporates the multiple scattering calculations as well as additional factors to produce reflectance and transmittance spectra that resemble experimental data.
The arguments taken by the multiple scattering model can be reliably measured for typical samples, except for particle volume fraction $\phi$, particle radius, and sample thickness, whose measurements can often be inaccurate or imprecise. These three parameters are included in the generative model as inference parameters. The multiple scattering model returns theoretical reflectance and transmittance calculations, $R^t(\lambda, \phi)$ and $T^t(\lambda, \phi)$, where $\lambda$ is the light wavelength at which reflectance and transmittance are calculated. A spectrum $\left \{ R^t(\phi) \right \}$ or $\left \{ T^t(\phi) \right \}$ consists of a set of $R^t_i(\lambda_i, \phi)$ or $T^t_i(\lambda_i, \phi)$ values, where $\left \{ \lambda_i \right \}$ represents all wavelengths used to produce the spectrum.
The theoretical reflectance $R^t$ and transmittance $T^t$ must be transformed into a corrected reflectance $R^c$ and transmittance $T^c$ before it can be compared with experimental data. Frequently, not all of the incident light on a non-absorbing sample is recorded as being reflected or transmitted since some light is lost in the experimental system. Since the exact reasons for loss are not known, they cannot be accurately measured. Instead, we parameterize a constant loss level $l_0$ with a linear modulation by wavelength with strength $l_1$, so that total losses $L_i = l_0 + l_1\lambda_i$. Therefore corrected spectra $\left \{ R^c(\phi) \right \}$ and $\left \{ T^c(\phi) \right \}$ are given by the elementwise product of $\left \{ R^t\left(\phi\right) \right \}$ or $\left \{ T^t\left(\phi\right) \right \}$ and $\left \{ 1-L_i \right \}$, or:
$$ R^c_i(\lambda_i,\phi, l_{0,r}, l_{1,r})=(1-l_{0,r}-l_{1,r}\lambda_i)\,R^t_i\left(\lambda_i, \phi\right) $$$$ T^c_i(\lambda_i,\phi, l_{0,t}, l_{1,t})=(1-l_{0,t}-l_{1,t}\lambda_i)\,T^t_i\left(\lambda_i, \phi\right) $$The $l_{0,r}$, $l_{1,r}$, $l_{0,t}$ and $l_{1,t}$ values are not generally known (or physically interesting), so they are marginalized over in the inference calculation (performed with MCMC so marginalization is trivial). Also, note that the loss parameters are different for reflectance and transmittance.
There are ways that the model could be extended in future work. A potential extension involves increasing the amount of information that goes into the multiple scattering calculation by allowing parameters to vary that are currently held fixed. Examples include sphere polydispersity and incident beam divergence, both of which are currently fixed at zero. Both values could be measured experimentally or sampled as parameters and marginalized over in an inference calculation. This would help to make the generative model more physical by accounting for phenomena that are currently neglected.
Due to the stochastic nature of the multiple scattering model, a calculated spectrum has uncertainties that must be included and propagated in the generative model. In general, these uncertainties will depend on the number of trajectories and scattering events we choose for a single calculation, and on the total number of runs we perform.
The choice of the number of scattering events is determined by the "equilibration" time for the scattering calculations. This means that we need enough scattering events such that most of the photons exit the sample, which would resemble more closely the behavior of light in an experimental measurement. The default setting is 200 scattering events, which was found to result in equilibration for the specific set of sample parameters that we used.
The number of trajectories will affect the noise level and reproducibility of the resulting spectra. The more trajectories we use, the less noisy and more reproducible the data. From previous calculations, we found that choosing 10,000 trajectories leads to very reproducible results; however, we also want to reduce the computation time of each run such that the later inference calculations are computationally accessible. The default value is 600 trajectories per run. The uncertainty range (set to 1 standard deviation) across the spectrum of visible wavelengths calculated from 100 of these runs is
$$ 0.00598 < \sigma_i < 0.02027, $$where $\sigma_i$ is the uncertainty in reflectance and transmittance at wavelength $\lambda_i$. The model uncertainties should be smaller than the measured uncertainties from experimental data for the model to have good predictive power.
Note, however, that the user can choose their own values for the number of trajectories. If the data uncertainties are large, then the user can reduce the number of trajectories to decrease computation time, as long as the model uncertainties are still smaller than the data's. The user can also change the number of events, though we advise setting its value to at least 200 scattering events.
Our basic likelihood function accepts theoretical spectra for a given volume fraction, particle radius, and thickness $\left \{ R^t(\phi, vf, thickness) \right \}$ and $\left \{ T^t(\phi, vf, thickness) \right \}$ calculated from the separate multiple scattering package. There are uncertainties associated with both the stochastic multiple scattering model and the experimental data. The model and experimental uncertainties are convoluted into a Gaussian with variance equal to the sum of the two independent noise source variances \cite{gregory}.
In the approximation of a linear model, likelihood is proportional to $\exp{\left(-\chi^2/2\right)}$. The likelihood is a product over all wavelengths $\lambda_i$, and depends on the experimental spectrum values $\{D_i\}$ with uncertainties $\{\sigma_{di}\}$, theoretical spectrum values $\{S_i^t\}$ with uncertainties $\{\sigma_{si}\}$, and loss parameters $l_0$ and $l_1$: $$ p\left(D|M,\phi,l_0,l_1,I\right)=\left(2\pi\right)^{-N/2}\left[\prod_{i=1} ^N \left(\sigma_{di}^2+\sigma_{si}^2\right)^{-1/2}\right]\exp{\left[\sum_{i=1}^N\frac{-\left(D_i-(1-l_0-l_1\lambda_i)S_i^t(\lambda_i,\phi)\right)^2}{2\left(\sigma_{di}^2+\sigma_{si}^2\right)}\right]} $$
Since the transmittance and reflectance spectra are recorded as separate, independent measurements, the likelihood of obtaining both spectra is just given by their product:
$$ \mathcal{L}\left(M,\phi,l_{0,r},l_{1,r},l_{0,t},l_{1,t}\right)=p\left(R|M,\phi,l_{0,r},l_{1,r},I\right)p\left(T|M,\phi,l_{0,t},l_{1,t},I\right) $$The prior for the volume fraction is given by previous knowledge of the experimental system. Structurally-colored samples are made of packings of spherical particles. The highest theoretical packing for spheres is 74%, which corresponds to face-centered cubic (fcc) or hexagonal closed-packed (hcp) crystal structures. Since the samples of interest in this model are amorphous, their volume fraction will always be less than 0.74 [5]. Therefore, the default prior for the volume fraction is given by a uniform probability over the interval:
$$ 0.1 < \mathrm{volume \, fraction} < 0.74 $$Note, however, that the user can set a narrower prior range if they have more information about the system.
The default prior for the radius is given by a uniform probability over the interval:
$$ 10 \,nm< \mathrm{particle \, radius} < 1000 \,nm $$The default prior for the sample thickness is also uniform over the range:
$$ 1 \,um< \mathrm{thickness} < 1000 \,um $$Realistically, the user will likely have better prior knowledge of the radius and thickness. For example, the prior for the particle radius will be dictated by the fact that the radius needs to be within the range that would produce a structural color peak in the visible part of the spectrum, in which case the range would be much smaller. The ranges for the priors should be adjusted on a sample-to-sample basis.
The prior for the losses is chosen based on the fact that the losses must be within the same range of possible values as the reflectance. A negative loss or a loss that is larger than 100% is not physically feasible. Thus, the prior is given by a uniform probability over the interval:
$$ 0 < l_0 < 1 \\ 0 < l_0 + l_1 < 1 $$
In [2]:
def calc_likelihood(spect1, spect2):
'''
Returns likelihood of obtaining an experimental dataset from a given theoretical spectrum
Parameters
-------
spect1: Spectrum object
experimental dataset
spect2: Spectrum object
calculated dataset
'''
resid_spect = calc_resid_spect(spect1, spect2)
chi_square = 0.
prefactor = 1.
if 'reflectance' in spect1.keys():
chi_square += np.sum(resid_spect.reflectance**2/resid_spect.sigma_r**2)
prefactor *= 1/np.prod(resid_spect.sigma_r * np.sqrt(2*np.pi))
if 'transmittance' in spect1.keys():
chi_square += np.sum(resid_spect.transmittance**2/resid_spect.sigma_t**2)
prefactor *= 1/np.prod(resid_spect.sigma_t * np.sqrt(2*np.pi))
return prefactor * np.exp(-chi_square/2)
Inference calculations with our generative model can take a long time due to the many calculations required in the multiple scattering algorithm. Here we present the steps to load in a simulated dataset and perform inference with a small number of walkers and steps using simulated reflectance data only.
The simulated data has volume fraction $\phi = 0.59$, particle radius = $119$ nm, thickness = $120$ um, and losses $l_0=0.1$, $l_1=0$. It was generated by a single run of the multiple scattering code. We perform a basic MCMC computation with 14 walkers each taking 10 steps. This is not nearly sufficient to achieve a representative sampling of parameter values or even to 'burn-in', so the output will be highly dependent on the initial parameter distributions.
We show traces of each parameter, as well as the log-probability over the walkers' trajectories. The traces are not particularly informative for such a small number of short trajectories. The main takeaway is that the walkers have not equilibrated over 10 steps (most visible in the plots of $l_0$ and $\log(p)$), and so the sparse distributions obtained do not truly represent the posterior probability. For this reason, we do not report most probable values from this calculation, since they would not be physically meaningful.
In [3]:
# Best guess for the parameters to be inferred
theta_guess = {'phi':0.59, 'radius':119, 'thickness':120, 'l0_r':0.1, 'l1_r':0}
# Range of values that contain the parameters according to the user's knowledge of the system.
theta_range = {'min_phi':0.35, 'max_phi':0.74, 'min_radius':70, 'max_radius': 120,
'min_thickness':10, 'max_thickness':300}
#convert test data into a spectrum file and load it
directory= os.path.join(os.getcwd(), 'infer_structcol', 'tests', 'test_data', 'simulated_data', 'reflection')
ifs.convert_data([450,500,550,600,650,700,750,800], 'ref.txt', 'dark.txt', directory)
spect = ifs.load_spectrum(refl_filepath = os.path.join(directory , 'converted', '0_data_file.txt'))
#define sample values
samp = ifs.Sample(spect.wavelength, particle_index=1.59, matrix_index=1)
#perform inference calculation
t0 = time.time()
walkers = ifs.run_mcmc(spect, samp, nwalkers=14, nsteps=10, theta_guess=theta_guess, theta_range=theta_range, seed=2)
print(time.time() - t0)
In [5]:
#code to visualize traces modified from notebook_week_09.ipynb
fig, (ax_l1, ax_l0, ax_vf, ax_radius, ax_thickness, ax_lnprob) = plt.subplots(6, figsize=(8,12))
ax_vf.set(ylabel='vf')
ax_radius.set(ylabel='radius')
ax_thickness.set(ylabel='thickness')
ax_l0.set(ylabel='l_0')
ax_l1.set(ylabel='l_1')
ax_lnprob.set(ylabel='ln(p)')
ax_vf.set(ylim=[0.3,0.7])
ax_l0.set(ylim=[0,0.2])
ax_l1.set(ylim=[-0.2,0.2])
ax_lnprob.set(ylim=[-100,20])
for i in range(6):
sns.tsplot(walkers.chain[i,:,0], ax=ax_vf)
sns.tsplot(walkers.chain[i,:,1], ax=ax_radius)
sns.tsplot(walkers.chain[i,:,2], ax=ax_thickness)
sns.tsplot(walkers.chain[i,:,3], ax=ax_l0)
sns.tsplot(walkers.chain[i,:,4], ax=ax_l1)
sns.tsplot(walkers.lnprobability[i,:], ax=ax_lnprob)
plt.show()