Showcasing of AMBER utilizing synthetic neutron scattering data

import numpy as np
from AMBER.background import background
import matplotlib.pyplot as plt

Generate signal data

The data for this tutorial is generated using the expected neutron scattering signal in an inelastic exmperiment measuring MnF\(_2\). The dispersion relation, i.e. the energy at a give (H,K,L) position is given by the analytical formula Yamany et al. 2010:

\[\omega(H,K,L) = \sqrt{\left(2S z_2 J_2+D+2 z_1 S J_1 \sin{L \pi}^2\right)^2-\left(2 S z_2 J_2 \cos(H \pi)\cos(K \pi) \cos(L \pi)\right)^2}\]

In the above \(J_1\), \(J_2\), and \(D\) are the magnetic coupling strengths and single-ion anisotropy which determine the amplitude of the dispersion. The parameters \(S = 5/2\), \(z_1 = 2\), and \(z_2 = 8\) corresponds the the size of the magnetic spin, and the number of nearest and next-nearest neighbours in the MnF\(_2\) lattice.

From the above equation, the imporant fact is that it depends on the 3D vector \(Q = (H,K,L)\). For simplicity, we assume that \(K\) is zero and thus reduce the data to 3D, i.e. \((H,L,\omega)\)

Due to how neutron scattering experiments function, the above dispersion will not be infinitely thin but rather extended. Here, this is replicated by a simple gaussian smearing along \(\omega\)

# Define parameters, dispersion, and Smearing
S = 5/2
z1 = 2
z2 = 8

def SpinWave(Q,J1,J2,D):
    return np.sqrt(np.power(2*S*z2*J2+D+2*z1*S*J1*np.sin(Q[2]*np.pi)**2,2.0)-

            np.power(2*S*z2*J2*np.cos(Q[0]*np.pi)*np.cos(Q[1]*np.pi)*np.cos(Q[2]*np.pi),2.0))

def Intensity(H,K,L,E):
    sigmaE = 0.25
    omega = SpinWave(np.array([H,K,L]),  J1=0.0354,J2=0.1499,D=0.131 )
    I = np.exp(-np.power(omega-E,2.0)/(2*sigmaE**2))

    return I

Define 3D grid domain

The input to AMBER is a 3D cube of intensities which is define below

# Data will ge simulated along H and L but with K = 0

h = np.linspace(-0.1,2.1,101)
k = 0
l = np.linspace(-0.1,2.1,101)

# Choose a sufficient energy range
e = np.linspace(0.5,8,61)

# Generate the grid upon which the dispersion is calculated

H,L,E = np.meshgrid(h,l,e)
K = np.zeros_like(H)

# Intensities are calculated with the smearing and scaled
I = 30*Intensity(H,K,L,E)

# Poisson noice is added to the intensity
I = np.random.poisson(I).astype(float)

Introduction of NaN-values

Triple axis instruments cannot measure all lengths of \(Q\). this we can mimic by exchaning the intensity at these points by NaNs

QLength = np.linalg.norm([H,L],axis=0)

I[QLength<0.35] = np.nan
I[QLength>2.5] = np.nan

Add background

The main feature of AMBER is the background segmentation, which requires the data to have a background. As describe in the article, background is assumed to:

• Rotation independence of the background.

• Smooth change of background along energy and \(|\vec{Q}|\).

• The signal is sparse but continuous in energy and \(|\vec{Q}|\).

In the following, a background is generated with a higher amplitude for low and high \(|\vec{Q}|\) mimicing instrumental background artefacts from triple axis neutron experiments - these corresponds to the direct beam contribution and to increased background at larger scattering angles.

# Background definition
def background_simulation(q,amplitude,gamma,mu,amplitude2,gamma2):
    return amplitude*((gamma/(q**2+gamma**2))+amplitude2*np.exp(-np.power(q-mu,2.0)/(2*gamma2**2)))


Q = np.linalg.norm([H,K,L],axis=0)

# Choose suitable valies
gamma = 0.5
mu = 3.0
amplitude = 20
amplitude2 = 1.0
gamma2 = 0.5

# Generate an example of the background for visual inspection
q = np.linspace(0.25,Q.max(),201)
bg_test = background_simulation(q,amplitude,gamma,mu,amplitude2,gamma2)

# display background amplitude
fig,ax = plt.subplots()
ax.plot(q,bg_test)
ax.set_xlabel('|Q|')
ax.set_ylabel('Background intensity')


## Add background to data
bg_tmp = background_simulation(Q,amplitude,gamma,mu,amplitude2,gamma2)

I+=np.random.poisson(bg_tmp)

Plot a constant energy plot

The now background affected data is plotted to show a “before” picture. This is done by choosing a specific energy (\(\omega \sim 4.5\) meV) and plotting the intensity as a color map

fig,ax = plt.subplots()

# Find the closest energy slice
energy = 4.5
EIdx = np.argmin(np.abs(E-energy))
ax.pcolormesh(H[:, :, EIdx], L[:, :, EIdx], I[:, :, EIdx], vmin=0, vmax = 50)

ax.axis('equal')
ax.set_title('E = {:.2f} meV'.format(e[EIdx]))
ax.set_xlabel('H [r.l.u.]')
ax.set_ylabel('L [r.l.u.]')

Text(0, 0.5, 'L [r.l.u.]')

Plot a constant H map

fig,ax = plt.subplots()

h_value = 0.25
HIdx = np.argmin(np.abs(h-h_value))
ax.pcolormesh(L[:, HIdx, :],E[:, HIdx, :],I[:, HIdx, :], vmin=0, vmax=50)

ax.set_xlabel('L [r.l.u.]')
ax.set_ylabel('E [meV]')

Text(0, 0.5, 'E [meV]')

Run denoising algorithm

# Initialize the AMBER background object

AMBER = background(dtype=np.float32)

# Set the grid sizes
AMBER.set_gridcell_size(dqx = 0.022, dqy = 0.022, dE = 0.125)

# Alternatively these can be set from the h, l, and e arrays like
# AMBER.set_volume_from_limits([h[0],l[0],e[0]],[h[-1],l[-1],e[-1]],)

# Input the data
AMBER.set_binned_data(h, l, e, I)

bins = int((q.max()-q.min())/0.022)

# define maximum radius and number of bins
AMBER.set_radial_bins(q.max(),n_bins=bins)

Set algorithm parameters (\(\lambda, \beta, \mu\))

lambda and mu will be determined as described in the paper and we select beta using cross validation

beta_ is selected using cross validation, i.e. mask out the top q quantile of intensity. The beta value for the lowest Root-Mean-Square-Error is then chosen

lambda_tmp = AMBER.MAD_lambda()
mu_tmp = AMBER.mu_estimator()


beta_range_tmp = np.array([0.1,1.0,10.0,100.0,200.0,300.0,400.0,500.0])
rmse = AMBER.cross_validation(q=0.3,beta_range= beta_range_tmp, lambda_=lambda_tmp, mu_=mu_tmp,n_epochs=15,verbose=False)


beta_tmp = beta_range_tmp[np.argmin(rmse)]

Test - ( 0.1 )
RMSE - ( 4.4478 0.1 55.79589160155072 ) :  3.4932678
Test - ( 1.0 )
RMSE - ( 4.4478 1.0 55.79589160155072 ) :  3.4923534
Test - ( 10.0 )
RMSE - ( 4.4478 10.0 55.79589160155072 ) :  3.4850137
Test - ( 100.0 )
RMSE - ( 4.4478 100.0 55.79589160155072 ) :  3.459206
Test - ( 200.0 )
RMSE - ( 4.4478 200.0 55.79589160155072 ) :  3.452546
Test - ( 300.0 )
RMSE - ( 4.4478 300.0 55.79589160155072 ) :  3.4521124
Test - ( 400.0 )
RMSE - ( 4.4478 400.0 55.79589160155072 ) :  3.454414
Test - ( 500.0 )
RMSE - ( 4.4478 500.0 55.79589160155072 ) :  3.4581234

Run the denoising algorithm using the parameters obtained using the heuristic

# set number of epochs
n_epochs = 20

AMBER.denoising(AMBER.Ygrid,lambda_tmp,beta_tmp,mu_tmp,n_epochs,verbose=True)

 Iteration  1
 Loss function:  24658333.70478791
 Iteration  2
 Loss function:  18359579.755618207
 Iteration  3
 Loss function:  15629052.90024966
 Iteration  4
 Loss function:  14376783.00374195
 Iteration  5
 Loss function:  13776632.868500363
 Iteration  6
 Loss function:  13504604.810678879
 Iteration  7
 Loss function:  13392707.644277826
 Iteration  8
 Loss function:  13349203.821932396
 Iteration  9
 Loss function:  13332601.151676292
 Iteration  10
 Loss function:  13326251.723918369
 Iteration  11
 Loss function:  13323807.99708973
 Iteration  12
 Loss function:  13322858.198862182
 Iteration  13
 Loss function:  13322487.582507188
 Iteration  14
 Loss function:  13322342.369417826
 Iteration  15
 Loss function:  13322285.33291491
 Iteration  16
 Loss function:  13322262.891915664
 Iteration  17
 Loss function:  13322254.058312038
 Iteration  18
 Loss function:  13322250.579985036
 Iteration  19
 Loss function:  13322249.209292663
 Iteration  20
 Loss function:  13322248.66973307

Compute substracted signal

# The subtracted signal is given by

Y_sub = AMBER.Ygrid - AMBER.b_grid
Y_back = AMBER.b_grid

# reshape the data to fit
Y_sub = Y_sub.reshape(AMBER.E_size,AMBER.Qx_size,AMBER.Qy_size).T
Y_back = Y_back.reshape(AMBER.E_size,AMBER.Qx_size,AMBER.Qy_size).T

# reshape the observation
Y_obs =  AMBER.Ygrid.reshape(AMBER.E_size,AMBER.Qx_size,AMBER.Qy_size).T

Display substracted signal

# reshape the observation
Y_obs =  AMBER.Ygrid.reshape(AMBER.E_size,AMBER.Qx_size,AMBER.Qy_size).T


fig0 = plt.figure(figsize=(20, 9))

# Plot 1: Observations Y
ax0 = fig0.add_subplot(1, 2, 1)

energy = 4.5
EIdx = np.argmin(np.abs(E-energy))
ax0.pcolormesh(H[:,:,EIdx],L[:,:,EIdx],Y_obs[:,:,EIdx],vmin=-2,vmax=50)

ax0.axis('equal')
ax0.set_xlabel('H [r.l.u.]')
ax0.set_ylabel('L [r.l.u.]')
ax0.set_title('Observations')

# Plot 2: Subtracted Y
ax1 = fig0.add_subplot(1, 2, 2)

energy = 4.5
EIdx = np.argmin(np.abs(E-energy))
p = ax1.pcolormesh(H[:,:,EIdx],L[:,:,EIdx],Y_sub[:,:,EIdx],vmin=-2,vmax=50)

ax1.axis('equal')
ax1.set_xlabel('H [r.l.u.]')
ax1.set_ylabel('L [r.l.u.]')
ax1.set_title('Subtracted signal')

fig0.colorbar(p)
plt.tight_layout()
plt.show()

fig0,axx = plt.subplots(figsize=(20, 9),ncols=2,nrows=1)

# Plot 1: Observations Y

ax0 = axx[0]

h_value = 0.25
HIdx = np.argmin(np.abs(h-h_value))
ax0.pcolormesh(L[:,HIdx,:],E[:,HIdx,:],Y_obs[:,HIdx,:],vmin=-2,vmax=50)

ax0.set_xlabel('H [r.l.u.]')
ax0.set_ylabel('L [r.l.u.]')
ax0.set_title('Observations')

# Plot 2: Subtracted Y
ax1 = axx[1]


HIdx = np.argmin(np.abs(h-h_value))
p = ax1.pcolormesh(L[:,HIdx,:],E[:,HIdx,:],Y_sub[:,HIdx,:],vmin=-2,vmax=50)

ax1.set_xlabel('L [r.l.u.]')
ax1.set_ylabel('E [meV]')
ax1.set_title('Subtracted signal')


fig0.colorbar(p,ax =axx, shrink=0.2, location='bottom')

#plt.tight_layout()
plt.show()

Compare the signal with the subtracted signal in a 1D cut

fig0 = plt.figure()

# Plot 1: Observations Y
ax0 = fig0.add_subplot(1, 1, 1)

h_value = 0.25
e_value = 5.5 # meV
EIdx = np.argmin(np.abs(e-e_value))
HIdx = np.argmin(np.abs(h-h_value))
ax0.scatter(L[:,HIdx,EIdx],Y_obs[:,HIdx,EIdx],label='Observations')
ax0.scatter(L[:,HIdx,EIdx],Y_sub[:,HIdx,EIdx],label='Subtracted Signal')
ax0.scatter(L[:,HIdx,EIdx],Y_back[:,HIdx,EIdx],label='Background')



ax0.set_xlabel('L [r.l.u.]')
ax0.set_ylabel('I [Arb. unit]')
ax0.legend()


plt.tight_layout()
plt.show()