import logging
from multiprocessing.sharedctypes import Value
log = logging.getLogger(__name__)
import numpy as np
from .prsdata import PrsData, TRANSMISSION, REFLECTION, ABSORPTION
from . import functions as f
from ..version import version as prsctrl_current_version
from devctrl.data.spectrum import transforms as st
from devctrl.data import transforms as tf, NDArray, Columns
[docs]
def remove_invalid_points_from_spectrum(spectrum_data: np.ndarray, columns: list[str], I_low_threshold=0.010, I_high_threshold=10.0):
"""
Remove some points that are very likely not physical:
- I (DC) is lower than 10 mV (default)
- I (DC) is larger than 10 V (default)
"""
invalid_idx = []
if "I" in columns:
I_idx = columns.index("I")
for i in range(spectrum_data.shape[0]):
Ival = spectrum_data[i,I_idx]
if Ival < I_low_threshold:
invalid_idx.append(i)
elif Ival > I_high_threshold:
invalid_idx.append(i)
remove_columns = ["I", "dI", "dI-X", "dI-Y", "dI_I", "dI-X_I", "dI-Y_I", "theta-corr"]
remove_columns += [f"s{col}" for col in remove_columns]
remove_columns_idx = [columns.index(col) for col in remove_columns if col in columns]
log.debug(f"Removing invalid data points {invalid_idx}")
for i in invalid_idx:
for c in remove_columns_idx:
spectrum_data[i,c] = np.nan
return spectrum_data, columns
[docs]
def remove_noisy_points_from_spectrum(spectrum_data: np.ndarray, columns: list[str], dI_I_min_amplitude=1e-6, I_noise_multiplier=5, dI_I_noise_multiplier=2):
"""
Remove points that are dominated by noise.
- # not used anymore: theta-corr if dI_I, dI-X_I or dI-Y_I are close to 0
- theta-corr if abs(dI_I) <= sdI_I * noise multiplier
- I (DC) <= sI * noise multiplier
NOTE: Make sure data for I and sI is included!
"""
# THETA-CORR
invalid_idx = []
if "theta-corr" in columns:
# Approach 1: dI_I Amplitude based
# dI_I_cols = ["dI_I", "dI-X_I", "dI-Y_I"]
# try:
# dI_I_col = next(col for col in dI_I_cols)
# dI_I_idx = columns.index(dI_I_col)
# for i in range(spectrum_data.shape[0]):
# dI_I = spectrum_data[i, dI_I_idx]
# if dI_I < dI_I_min_amplitude:
# invalid_idx.append(i)
# except StopIteration:
# log.warning(f"Can not remove invalid 'theta-corr' points from spectrum because none of the following columns are present: {dI_I_cols}")
# Approach 2: Noise based: dI_I near 0 will have noise > value, these values will produce noisy phase
dI_I_col = None
dI_I_idx = 0
sdI_I_idx = 0
dI_I_cols = ["dI_I", "dI-X_I", "dI-Y_I"]
for col in dI_I_cols:
if col in columns and "s"+col in columns:
dI_I_col = col
dI_I_idx = columns.index(col)
sdI_I_idx = columns.index("s"+col)
break
if dI_I_col is not None:
for i in range(spectrum_data.shape[0]):
dI_I = spectrum_data[i, dI_I_idx]
sdI_I = spectrum_data[i, sdI_I_idx]
if np.abs(dI_I) < sdI_I * dI_I_noise_multiplier:
invalid_idx.append(i)
else:
log.warning(f"Can not remove invalid 'theta-corr' points from spectrum because none of the following columns are present: {dI_I_cols}")
remove_columns = ["theta-corr", "stheta-corr"]
remove_columns_idx = [columns.index(col) for col in remove_columns if col in columns]
for i in invalid_idx:
for c in remove_columns_idx:
spectrum_data[i,c] = np.nan
# I
invalid_idx = []
if "I" in columns and "sI" in columns:
I_idx = columns.index("I")
sI_idx = columns.index("sI")
for i in range(spectrum_data.shape[0]):
Ival = spectrum_data[i,I_idx]
sIval = spectrum_data[i,sI_idx]
if Ival <= sIval * I_noise_multiplier:
invalid_idx.append(i)
else:
log.warning(f"No I or sI columns in spectrum data.")
remove_columns = ["I", "dI", "dI-X", "dI-Y", "dI_I", "dI-X_I", "dI-Y_I", "theta-corr"]
remove_columns += [f"s{col}" for col in remove_columns]
remove_columns_idx = [columns.index(col) for col in remove_columns if col in columns]
for i in invalid_idx:
for c in remove_columns_idx:
spectrum_data[i,c] = np.nan
return spectrum_data, columns
[docs]
def remove_spectrum_outliers(data: np.ndarray, columns: list[str], check_cols: list[str]|None=None, threshold: float=0.6, min_points=20):
"""
Set outlier points to nan.
An outlier point is a point that deviates more that full_range * threshold from the left AND right neighbor.
full_range is difference between the maximum and minimum value
"""
if not check_cols:
check_cols = columns
outliers = set()
for col in check_cols:
try:
col_idx = columns.index(col)
except ValueError as e:
log.error(f"Can not remove outliers for column '{col}' - no such column (columns={columns})")
outliers |= set(find_outliers(data[:,col_idx]))
for idx in outliers:
for col in columns:
if col in ["E", "wl"]: continue
col_idx = columns.index(col)
data[idx, col_idx] = np.nan
return data, columns
if data.shape[0] < min_points: return data
full_range = np.max(data) - np.min(data)
def check_left(data, i):
return np.abs(data[i] - data[i-1]) > full_range * threshold
def check_right(data, i):
return np.abs(data[i+1] - data[i]) > full_range * threshold
if check_right(data, 0):
data[0] = np.nan
for i in range(1, data.shape[0]-1):
if not np.isnan(data[i-1]) and check_left(data, i) and check_right(data, i):
print(f"Removing point {i}: {data[i]}")
data[i] = np.nan
if check_left(data, data.shape[0]-1):
data[-1] = np.nan
return data
sane_transforms = [remove_spectrum_outliers, remove_invalid_points_from_spectrum]
[docs]
def find_outliers(data: np.ndarray, threshold: float=0.6, min_points=20) -> list[int]:
"""
Set outlier points to nan.
An outlier point is a point that deviates more that full_range * threshold from the left AND right neighbor.
full_range is difference between the maximum and minimum value
"""
outlier_idxs = []
if data.shape[0] < min_points: return outlier_idxs
full_range = np.max(data) - np.min(data)
def check_left(data, i):
return np.abs(data[i] - data[i-1]) > full_range * threshold
def check_right(data, i):
return np.abs(data[i+1] - data[i]) > full_range * threshold
if check_right(data, 0):
outlier_idxs.append(0)
for i in range(1, data.shape[0]-1):
if not np.isnan(data[i-1]) and check_left(data, i) and check_right(data, i):
# log.debug(f"Found outlier point {i}: {data[i]}")
outlier_idxs.append(i)
if check_left(data, data.shape[0]-1):
outlier_idxs.append(data.shape[0]-1)
return outlier_idxs
[docs]
def remove_outliers(data: np.ndarray, threshold: float=0.6, min_points=20):
"""
Set outlier points to nan.
An outlier point is a point that deviates more that full_range * threshold from the left AND right neighbor.
full_range is difference between the maximum and minimum value
"""
if data.shape[0] < min_points: return data
full_range = np.max(data) - np.min(data)
def check_left(data, i):
return np.abs(data[i] - data[i-1]) > full_range * threshold
def check_right(data, i):
return np.abs(data[i+1] - data[i]) > full_range * threshold
if check_right(data, 0):
data[0] = np.nan
for i in range(1, data.shape[0]-1):
if not np.isnan(data[i-1]) and check_left(data, i) and check_right(data, i):
# print(f"Removing point {i}: {data[i]}")
data[i] = np.nan
if check_left(data, data.shape[0]-1):
data[-1] = np.nan
return data, columns
[docs]
def get_scale(data: PrsData):
"""
Calculate how the data should be scaled, so that all measurements can be compared against each other.
Data recorded with an pre-amplifier gain of 10^6 V/A corresponds to a scaling of 1
10^6 V/A -> scale = 1
10^7 V/A -> scale = 0.1
"""
scale = 1.
gain = -1
try:
gain = int(data.metadata["pre-amp"]["gain"])
except:
pass
if gain == -1:
log.warning(f"Failed to determine data scale. (data.mode={data.mode}, data.name={data.name})")
else:
scale = 10 ** (6-gain)
# elif gain == 10:
# scale = 0.0001
# elif gain == 9:
# scale = 0.001
# elif gain == 8:
# scale = 0.01
# elif gain == 7:
# scale = 0.1
# elif gain == 6:
# scale = 1.0
# elif gain == 5:
# scale = 10.0
# elif gain == 4:
# scale = 100.0
# elif gain == 3:
# scale = 1000.0
# else:
# raise ValueError(f"Unexpected gain: {gain} (data.mode={data.mode}, data.name={data.name})")
return scale
[docs]
def calc_absorption_data(ref_data: PrsData, tra_data: PrsData, norm_data: PrsData|None=None, shift_to_sane_values=True, **prsdata_kw) -> PrsData:
"""
Calculate the absorption data from reflection, transmission and reference data (normalization data).
The reference data is interpolated and sampled at the wavelengths that are shared by ref_data and tra_data.
:param shift_to_sane_values: If True, shift the DC intensity (A) values so that no value is negative (since that is unphysical and likely due to losses in reference spectrum).
:param prsdata_kw: Keywords for the PrsData constructor. Pass file_mode="w" if you want to save the absorption data later.
"""
md = {
"data_info": "Absorption data calculated from reflection and transmission data",
"data_ref_dir": ref_data.dirname,
"data_ref_name": ref_data.name,
"data_tra_dir": tra_data.dirname,
"data_tra_name": tra_data.name,
"name": ref_data.name,
}
if norm_data:
md |= {
"data_norm_dir": norm_data.dirname,
"data_norm_name": norm_data.name,
}
else:
md |= {
"data_norm": "No normalization data with 'A' values given -> the absorption data is in arbitrary units!"
}
if "prsctrl_version" in ref_data.metadata:
md["prsctrl_version"] = prsctrl_current_version
abs_data = PrsData(data={}, metadata=md, exp_mode=ABSORPTION, **(dict(file_mode="") | prsdata_kw))
wls = np.intersect1d(ref_data.wavelengths, tra_data.wavelengths).tolist()
# get the scaling factor of the data (from pre-amp settings)
ref_scale = get_scale(ref_data)
tra_scale = get_scale(tra_data)
norm_scale = get_scale(norm_data) if norm_data else 1.0
log.info(f"{ref_data.name}: Scaling for absorption calculation (r,t,n): {ref_scale}, {tra_scale}, {norm_scale}")
warnings = []
if norm_data:
norm_cols = ["wl", "I", "sI"]
sdata_norm, _ = norm_data.get_spectrum_data(columns=norm_cols)
sdata_norm, _ = st.interpolate_to_wavelengths(sdata_norm, norm_cols, wls)
else:
sdata_norm = np.array([[wl, 10, 10] for wl in wls])
I0 = sdata_norm[:,1] * norm_scale
sI0 = sdata_norm[:,2] * norm_scale
rt_cols = ["wl", "I", "sI", "dI", "sdI", "dI-X", "sdI-X", "dI-Y", "sdI-Y"]
sdata_ref, rt_cols = ref_data.get_spectrum_data(wavelengths=wls, columns=rt_cols)
sdata_tra, rt_cols = tra_data.get_spectrum_data(wavelengths=wls, columns=rt_cols)
if sdata_ref.shape != sdata_tra.shape:
raise ValueError(f"Ref and Tra data have different shapes: ref={sdata_ref.shape}, tra={sdata_tra.shape}.")
R = sdata_ref[:,rt_cols.index("I")] * ref_scale
sR = sdata_ref[:,rt_cols.index("sI")] * ref_scale
T = sdata_tra[:,rt_cols.index("I")] * tra_scale
sT = sdata_tra[:,rt_cols.index("sI")] * tra_scale
A, sA = f.A(I0, sI0, R, sR, T, sT)
if shift_to_sane_values:
_min = np.min(A)
if _min < 0:
extra_offset = 1e-6 # avoid having A=0
A += np.abs(_min) + extra_offset
sA += np.abs(_min) + extra_offset
log.warning(f"Shifting A values by {np.abs(_min)} V")
if not "messages" in abs_data.metadata:
abs_data.metadata["messages"] = []
abs_data.metadata["messages"].append(f"Shifted A values by {np.abs(_min)} V + {extra_offset:.2e} V to avoid negative/zero values")
dAs = {}
sdAs = {}
calc_dA_keys = ["dI", "dI-X", "dI-Y"]
for k in calc_dA_keys:
dR = sdata_ref[:,rt_cols.index(k)] * ref_scale
sdR = sdata_ref[:,rt_cols.index("s"+k)] * ref_scale
dT = sdata_tra[:,rt_cols.index(k)] * tra_scale
sdT = sdata_tra[:,rt_cols.index("s"+k)] * tra_scale
dA, sdA = f.dA(dR, sdR, dT, sdT)
dAs[k] = dA
sdAs[k] = sdA
for i, wl in enumerate(wls):
abs_data[wl] = {}
A_val = A[i]
if A_val <= 0:
warnings.append((wl, "A", A_val))
abs_data[wl]["I"] = A_val
# abs_data[wl]["lock-in-aux"] = A_val
abs_data[wl]["sI"] = sA[i]
for k in calc_dA_keys:
abs_data[wl][k] = dAs[k][i]
abs_data[wl]["s"+k] = sdAs[k][i]
abs_data.get_for_wl(wl, "dI-X_I") # force calculation of values
if warnings:
log.warning(f"calc_absorption_data: Invalid 'A' values at wavelengths (nm): {[w[2] for w in warnings]}")
return abs_data
[docs]
def calc_modulus_spectrum(data: PrsData, which_real="dI-X_I", extra_cols=["wl", "E"], **get_spectrum_kw):
"""
Calculate the modulus spectrum.
The modulus spectrum is the modulus of the real spectrum and an imaginary spectrum, which is obtained from the real spectrum via Kramers-Kronig transformation.
See T. J. C. HOSEA in phys. stat. sol. (b) 182, K43 (1994)
The function returns spectrum data with two additional columns: `<which_real>-imag`, the imaginary spectrum from KK, and `L`, the modulus spectrum.
:param which_real: Which column to use as real spectrum. Use either dI_I, dI-X_I (default) or dI-Y_I
:return: spectrum data, columns
"""
try:
import pykk
except ImportError as e:
raise ImportError("Module pykk is required for Kramers-Kronig transformation during modulus spectrum calculation. Please install pykk.") from e
cols = list(set(extra_cols + ["wl", "E", which_real, "s" + which_real]))
sdata, cols = data.get_spectrum_data(columns=cols, **get_spectrum_kw)
sdata, cols = tf.remove_nans(sdata, cols)
sdata, cols = st.make_const_energy_interval_spectrum(sdata, cols)
print(sdata)
Es = sdata[:,cols.index("E")]
dR_Rs = sdata[:,cols.index(which_real)]
try:
kk_data = np.array(pykk.real2imag(Es, dR_Rs))
except ValueError as e:
raise ValueError(f"ValueError occured in pykk.real2imag. Es={Es}") from e
print(sdata.shape, kk_data.shape)
L = np.sqrt(dR_Rs**2 + kk_data**2)
sdata = np.vstack([sdata.T, kk_data, L]).T
cols += ["imag", "L"]
return sdata, cols
