[Python] overflow encountered in matmul, what can I do?

I am currently using an open source Python package called snompy (https://github.com/TomVincentUK/snompy). While implementing the package, a function that defines the transfer matrix returns an overflow error in the line that uses numpy's matrix product (matmul). The part of the code that raises the problem is as follows:

def transfer_matrix_qs(self, q=0.0):
"""Return the transfer matrix for the sample in the quasistatic
limit.

Parameters
----------
q : float, default 0.0
In-plane electromagnetic wave momentum.
Must be broadcastable with all `eps_stack[i, ...]` and
`t_stack[i, ...]`. Either `q` or `theta_in` must be None.

Returns
-------
M_qs : complex
The quasistatic transfer matrix of the sample.

Notes
-----
This implementation of the transfer matrix method is based on the
description given in reference [1]_.

References
----------
.. [1] T. Zhan, X. Shi, Y. Dai, X. Liu, and J. Zi, “Transfer matrix
method for optics in graphene layers,” J. Phys. Condens. Matter,
vol. 25, no. 21, p. 215301, May 2013,
doi: 10.1088/0953-8984/25/21/215301.

"""
trans_factor = self.eps_stack[:-1] / self.eps_stack[1:]
trans_matrices = (
1 + np.array([[1, -1], [-1, 1]]) * trans_factor[..., np.newaxis, np.newaxis]
) / 2

# Convert stack to single transfer matrix (accounting for propagation if needed)
M_qs = trans_matrices[0]
if self.multilayer:
# Optical path length of internal layers
prop_factor = np.array([q * t for t in self.t_stack])

# Avoid overflow by clipping (there's probably a more elegant way)
prop_factor = np.clip(prop_factor, -_MAX_EXP_ARG, _MAX_EXP_ARG)
prop_matrices = np.exp(
np.array([1, -1]) * prop_factor[..., np.newaxis, np.newaxis]
) * np.eye(2)

for T, P in zip(trans_matrices[1:], prop_matrices):
M_qs = M_qs @ P @ T
else:
M_qs = M_qs * np.ones_like(q)[..., np.newaxis, np.newaxis]

return M_qs


The error then propagates to functions that depend on the transfer_matrix_qs results. The error window on my side is displayed as follows:

C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:228: RuntimeWarning: overflow encountered in matmul M_qs = M_qs @ P @ T C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:228: RuntimeWarning: invalid value encountered in matmul M_qs = M_qs @ P @ T C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:252: RuntimeWarning: overflow encountered in divide return M[..., 1, 0] / M[..., 0, 0] C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:252: RuntimeWarning: invalid value encountered in divide return M[..., 1, 0] / M[..., 0, 0] C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:723: RuntimeWarning: invalid value encountered in multiply phi = np.sum(w_lag * beta_q, axis=0) / (2 * z_Q) C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:723: RuntimeWarning: invalid value encountered in divide phi = np.sum(w_lag * beta_q, axis=0) / (2 * z_Q) C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:724: RuntimeWarning: invalid value encountered in multiply E = np.sum(w_lag * x_lag * beta_q, axis=0) / (4 * z_Q**2) C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:724: RuntimeWarning: invalid value encountered in divide E = np.sum(w_lag * x_lag * beta_q, axis=0) / (4 * z_Q**2) C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:796: RuntimeWarning: invalid value encountered in divide d_image = np.abs(phi / E) - z_Q C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\sample.py:797: RuntimeWarning: invalid value encountered in divide beta_image = phi**2 / E C:\Spyder\envs\spyder-runtime\Lib\site-packages\snompy\fdm.py:261: RuntimeWarning: invalid value encountered in divide alpha_eff = 1 + (beta_0 * f_0) / (2 * (1 - beta_1 * f_1))
​

I have prepared a test code you can copy and run yourself to reproduce the error:

import numpy as np
import snompy

# Set some experimental parameters
A_tip = 90e-9 # AFM tip tapping amplitude
r_tip = 10e-9 # AFM tip radius of curvature
L_tip = 300e-9 # Semi-major axis length of ellipsoid tip model
n = 2 # Harmonic for demodulation
theta_in = np.deg2rad(60) # Light angle of incidence
c_r = 0.9 # Experimental weighting factor
nu_vac = np.linspace(600,1000, 1000) * 1e2 # Vacuum wavenumber
method = "multi" # The FDM method to use
t_aln= 250e-9 #AlN layer thickness
t_sample=np.array([t_aln])
eps_air = 1.0 # Air permittivity
eps_Si = 11.7 # Si permitivitty in the mid-infrared
### GaN Lorentzian model parameters
einf_gan = 5.35 #epsilon infinity GaN
to_gan = 530 #TO phonon frequency in cm-1
lo_gan = 735 #LO phonon frequency in cm-1
gamma_gan = 3.5 #relaxation factor
### AlN Lorentzian model parameters
einf_aln=4.6
to_aln=649.5
lo_aln=880.3
gamma_aln= 9

#Dielectric permittivity from a single Lorentzian oscillator:
eps_gan = einf_gan*(1+ (lo_gan**2 -to_gan**2)/(to_gan**2 - (nu_vac*1e-2)**2 - (nu_vac*1e-2)*1j*gamma_gan))
eps_aln= einf_aln*(1+ (lo_aln**2 -to_aln**2)/(to_aln**2 - (nu_vac*1e-2)**2 - (nu_vac*1e-2)*1j*gamma_aln))

#Model of AlN on GaN
sample_aln = snompy.Sample(
eps_stack=(eps_air,eps_aln, eps_gan),
t_stack = t_sample,
nu_vac=nu_vac,)

# Model of Si
sample_si = snompy.bulk_sample(eps_sub=eps_Si, eps_env=eps_air, nu_vac=nu_vac)

# Measurement
alpha_eff_aln = snompy.fdm.eff_pol_n(
sample=sample_aln, A_tip=A_tip, n=n, r_tip=r_tip, L_tip=L_tip, method=method
)
r_coef_aln = sample_aln.refl_coef(theta_in=theta_in)
sigma_aln = (1 + c_r * r_coef_aln) ** 2 * alpha_eff_aln

# Silicon reference
alpha_eff_si = snompy.fdm.eff_pol_n(
sample=sample_si, A_tip=A_tip, n=n, r_tip=r_tip, L_tip=L_tip, method=method
)
r_coef_si = sample_si.refl_coef(theta_in=theta_in)
sigma_si = (1 + c_r * r_coef_si) ** 2 * alpha_eff_si

# Normalised complex scattering
eta_n = sigma_aln / sigma_si


I'm interested in eta_n's amplitude and phase, but it has a bunch of nans caused by the overflow. Reducing the AlN layer thickness (t_aln) helps in reducing the amount of elements in the array that overflow, but I want the code to work for the actual thickness.

Before asking directly to the author of the code, I was wondering if there's anything on my end that I can do to avoid the function to overflow. The function uses complex numbers, so very naïvely I tried to change the type from complex128 to complex256 but apparently that type is not supported on my machine. Any further suggestions are highly appreciated.

Continue reading...