ramansep

ramansep is a Python package that tells strain apart from charge in Raman maps of 2D materials (sheets only one or a few atoms thick, such as graphene or monolayer MoS2).

A Raman peak moves when the material is stretched, and it also moves when extra electrons or holes are added. One peak alone cannot say which of the two happened. Two peaks that respond differently can: if you know how much each peak moves per unit of strain and per unit of charge, the two measured shifts give two equations with two unknowns. ramansep solves those equations at every pixel of a map and gives a strain map and a charge-density map, each with an error bar. It also covers the steps around that:

How do I read the text files my instrument exports, and get the peak shifts out of the measured spectra? Which pixels failed or carry no peak at all?
What if two peaks overlap?
Which pair of peaks separates the two causes best?
With three or more peaks, does "strain plus charge" actually explain what I measured, pixel by pixel?
Can I trade some resolution for smaller error bars if the fields are smooth, and fill the pixels whose fits failed?
How do I measure the response numbers on my own instrument, even when the starting strain and charge of my sample are unknown, and how should I plan that measurement?
How much of my final error bar comes from that calibration, with two peaks or more?
What if laser heating also shifts the peaks?

The results are checked by automated tests against independent calculations and exact limits (see How the results are checked). When an input cannot give a trustworthy answer, the package stops with an error that says why, rather than returning a number that looks fine.

A short guide to the words used here
Install, requirements and units
Examples (each with the output it prints)
What is in the package
Cited coefficient sets
When it refuses, and why
How the results are checked
Corrections in earlier versions
Limits, and what is deliberately not included
Where it comes from
Citing, support and license

A short guide to the words used here

Raman spectrum -- the light scattered by a material, plotted against its energy loss in cm^-1 (wavenumbers). Each lattice vibration (phonon, or mode) gives a peak.
Mode names -- labels such as A'1, E', 2LA(M) (MoS2) or G, 2D (graphene). They only name which peak is meant.
Shift (dw) -- how far a peak has moved from its position in the undisturbed ("pristine") material, in cm^-1. Negative means the peak moved to lower wavenumber (the mode "softened").
Strain -- stretching of the sheet, in percent. Biaxial means stretched equally in both directions; uniaxial means along one direction. Positive is tension.
Carrier density (or doping, charge) -- extra electrons or holes (missing electrons) per area, in the MoS2 sets in units of 1e13 cm^-2.
Gate, gated -- an electrode that adds a known amount of charge to the sheet when a voltage is applied; a strain stage stretches the sample by a known amount. Both are used to measure how far a peak moves per unit of charge or strain.
Lever arm -- how much one peak shifts per unit of one cause, for example -5.1 cm^-1 per percent strain. The lever arms of all the peaks form the lever-arm matrix K (one row per peak, one column per cause). The package assumes the shifts are linear in the causes: shifts = K @ (strain, density).
Condition number -- a single number that says how close the matrix is to being unsolvable. Large values mean the peaks respond too similarly, so small shift errors become large strain or density errors.
Sigma / error bar -- one standard deviation of a value. Correlation between two estimates near +1 or -1 means they are hard to tell apart. The covariance matrix collects the squared error bars (on its diagonal) and how each pair of estimates varies together (off the diagonal).
Weighted least squares -- the standard way to solve more equations than unknowns: find the values that make the weighted sum of squared mismatches smallest, each measurement weighted by 1/sigma^2 so that more precise peaks count more.
Chi-square and p-value -- with more peaks than causes, the equations cannot all be met exactly when something else is going on. Chi-square measures the mismatch, weighted by the error bars; the p-value is the probability of a mismatch at least that large if strain and charge alone were at work. A tiny p-value flags a pixel. The degrees of freedom (dof) of this check are the number of peaks minus the number of causes solved for.
Hyperspectral map (cube) -- one spectrum per pixel, stored as an array of shape (rows H, columns W, spectral points L).
Lorentzian / Voigt -- peak shapes. A Lorentzian is the natural shape of a phonon peak; a Voigt is a Lorentzian blurred by a Gaussian, for peaks broadened by the instrument.
FWHM -- full width at half maximum, the width of a peak measured halfway up its height.
Baseline -- the background under a peak. A sloped background (for example from fluorescence) pulls a fitted peak centre sideways unless the slope is fitted too.
Stokes / anti-Stokes -- Raman peaks on the energy-loss side and the energy-gain side. The ratio of their heights depends on temperature, so it works as a thermometer.
DFT (density functional theory) -- a standard quantum-mechanical computer calculation of a material. "Frozen-phonon" DFT computes a vibration's frequency by displacing the atoms along that vibration. The source paper used it to compute the MoS2 strain lever arms.
Grueneisen parameter (gamma) -- a number that says how strongly a vibration's frequency w changes when the sheet is stretched; for biaxial strain the package uses gamma = -(1/(2w)) dw/d(strain), with strain as a fraction (not percent).

Install, requirements and units

pip install ramansep

For development: clone the repository and pip install -e .[test].

It needs Python 3.9 or newer, NumPy 1.22 or newer and SciPy 1.8 or newer. matplotlib is optional (pip install ramansep[plot]) and is used only by the demo script examples/synthetic_map.py.

Units and signs:

Shifts are in cm^-1, relative to the pristine reference frequency of each mode (you supply those references, e.g. ref1, ref2).
The inversion returns strain and density in whatever units your lever arms use. The shipped MoS2 sets use percent biaxial strain (positive = tension) and electron density in 1e13 cm^-2. The graphene set uses percent uniaxial strain and, on purpose, a doping coordinate in cm^-1 (see Cited coefficient sets).
Every ModeCoefficients must carry a reference string that says where its numbers came from; it is a required field.

Examples

Each example below runs as written, and the output shown is what it printed with ramansep 0.13.0. Unless a value is said to come from a cited source, it is an illustrative value chosen for the example. The repository also has a longer demo, examples/synthetic_map.py, that inverts a noisy synthetic 128 x 128 map and plots it with matplotlib.

1. Strain and charge at two spots of a MoS2 ribbon

The shifts here are the ones quoted by the source paper (Mahim and Rahman, under review) from tip-enhanced Raman maps (maps taken with a sharp metal tip that gives nanometre resolution; Krayev et al., Appl. Phys. Lett. 128, 203102 (2026)); the 0.1 cm^-1 shift error is illustrative.

import numpy as np
import ramansep as rs

model = rs.SeparationModel(rs.mos2_a1_2la())    # monolayer MoS2, A'1 + 2LA(M)

# Peak shifts (cm^-1) quoted by the source paper at two spots:
# the ribbon edge (A'1 -0.5, 2LA(M) 0.0) and an interior spot (-0.6, -2.8).
dw_a1 = np.array([-0.5, -0.6])
dw_2la = np.array([0.0, -2.8])
res = model.invert(dw_a1, dw_2la, sigma1=0.1, sigma2=0.1)   # 0.1 cm^-1: illustrative

for name, i in [("edge", 0), ("interior", 1)]:
    print(f"{name:8s} strain {res.strain[i]:+.3f} +/- {res.strain_sigma[i]:.3f} %   "
          f"electrons {res.density[i] * 1e13:+.2e} +/- {res.density_sigma[i] * 1e13:.1e} cm^-2")
print(f"condition number of the matrix: {res.condition_number:.1f}")

edge     strain +0.000 +/- 0.005 %   electrons +2.27e+12 +/- 4.7e+11 cm^-2
interior strain +0.134 +/- 0.005 %   electrons -3.78e+11 +/- 4.7e+11 cm^-2
condition number of the matrix: 10.1

The edge carries about 2.3e12 electrons per cm^2 and no strain; the interior spot is under 0.134 % tension, and its density is within one error bar of zero. These are the paper's published separations; the test suite checks them (see below). res.correlation gives the per-pixel correlation of the two estimates.

2. Choosing a mode pair

import numpy as np
import ramansep as rs

sig = 0.1                                        # illustrative shift error, cm^-1
paper = rs.SeparationModel(rs.mos2_a1_2la()).invert(0.0, 0.0, sig, sig)
optical = rs.SeparationModel(rs.mos2_eprime_a1()).invert(0.0, 0.0, sig, sig)
print(f"A'1 + 2LA(M): strain error {float(paper.strain_sigma):.4f} %")
print(f"E'  + A'1   : strain error {float(optical.strain_sigma):.4f} %")
print(f"ratio: {float(optical.strain_sigma / paper.strain_sigma):.1f}")

A'1 + 2LA(M): strain error 0.0048 %
E'  + A'1   : strain error 0.0257 %
ratio: 5.4

Both optical modes E' and A'1 move only a little with strain, so with the same shift errors that pair gives a strain error about five times larger than the A'1 + 2LA(M) pair. The error bars depend only on the lever arms and the shift errors, so this comparison can be made before measuring anything.

3. From a map file to strain and charge maps

import os, tempfile
import numpy as np
import ramansep as rs
from ramansep import lorentzian

coeffs = rs.mos2_a1_2la()
model = rs.SeparationModel(coeffs)

# An illustrative 4 x 5 map: known strain (%) and electron density (1e13 cm^-2).
rng = np.random.default_rng(3)
strain = rng.uniform(-0.1, 0.1, (4, 5))
density = rng.uniform(0.0, 0.2, (4, 5))
dw1, dw2 = model.forward(strain, density)

# Make a spectrum for every pixel: two Lorentzian peaks plus noise.
x = np.linspace(380.0, 480.0, 401)
cube = np.empty((4, 5, x.size))
for i in range(4):
    for j in range(5):
        cube[i, j] = (lorentzian(x, 404.7 + dw1[i, j], 3.0, 800.0, 50.0)
                      + lorentzian(x, 452.0 + dw2[i, j], 4.0, 500.0, 0.0)
                      + rng.normal(0.0, 3.0, x.size))
cube[3, 4, 100] = np.nan                         # one broken pixel

path = os.path.join(tempfile.mkdtemp(), "map.csv")
rs.save_map_csv(path, x, cube)                   # the documented file format

x2, cube2 = rs.load_map_csv(path)
fits = rs.fit_map(x2, cube2, window1=(395, 415), window2=(440, 465),
                  ref1=404.7, ref2=452.0)
print("masked pixels:", fits.n_masked, "at", np.argwhere(~fits.ok).tolist())

inv = model.invert(fits.dw1, fits.dw2, fits.sigma1, fits.sigma2)
ok = fits.ok
z_s = np.abs(inv.strain - strain)[ok] / inv.strain_sigma[ok]
z_d = np.abs(inv.density - density)[ok] / inv.density_sigma[ok]
print(f"typical error bars: strain {np.median(inv.strain_sigma[ok]):.4f} %, "
      f"density {np.median(inv.density_sigma[ok]) * 1e13:.1e} cm^-2")
print(f"pixels whose error is within 3 error bars: strain {(z_s < 3).sum()} of {ok.sum()}, "
      f"density {(z_d < 3).sum()} of {ok.sum()}")
print("strain at the broken pixel:", inv.strain[3, 4])

masked pixels: 1 at [[3, 4]]
typical error bars: strain 0.0002 %, density 1.3e+10 cm^-2
pixels whose error is within 3 error bars: strain 19 of 19, density 19 of 19
strain at the broken pixel: nan

The map file is plain CSV with the header row,col,wavenumber_cm1,counts (0-based pixel indices, one shared wavenumber axis). A single spectrum uses wavenumber_cm1,counts (save_spectrum_csv, load_spectrum_csv). To read the plain text tables that instruments export, see example 16. fit_map fits both peaks at every pixel. A pixel is masked (NaN in all four output maps, False in fits.ok) when its counts are not finite, a fit stops with an error or does not converge, a fitted centre ends up outside its window, or a fitted peak height is zero or negative (a dip is not a Raman peak; this check is new in 0.12.0). fits.reason records, per pixel, why it was masked (rs.MASK_REASONS spells out the codes). Maps with pixels that carry no peak at all need one more check, min_snr (example 12). Mistakes in the settings themselves (window limits, baseline name, a wavenumber axis containing NaN) are refused with an error before any pixel is fitted. The per-pixel inversion keeps masked pixels as NaN and does not touch their neighbours. For spectra on a sloped background, pass baseline="linear"; for peaks close enough to overlap, pass joint=True (example 13).

4. Three peaks, and a check that the model holds

import numpy as np
import ramansep as rs

# Lever arms (cm^-1 per unit strain, per unit density) of three modes.
# Rows 1-2 are the MoS2 A'1 and 2LA(M) values; row 3 is an invented,
# illustrative third mode.
K = np.array([[-5.1, -2.2],
              [-20.9, 0.0],
              [-8.0, -1.1]])
sig = [0.15, 0.10, 0.12]                         # illustrative shift errors, cm^-1
mm = rs.MultiModeModel(K, mode_names=["A'1", "2LA(M)", "X"])

rng = np.random.default_rng(3)
strain = rng.normal(0, 0.2, (30, 30))
density = rng.normal(0, 1.0, (30, 30))
shifts = [s + rng.normal(0, e, s.shape) for s, e in zip(mm.forward(strain, density), sig)]
shifts[2][10:15, :] += 2.0                       # rows 10-14: a cause the model does not know

res = mm.invert(shifts, sigmas=sig)
print("degrees of freedom of the check:", res.dof)
print(f"median p-value, rows 10-14 : {np.median(res.p_value[10:15]):.1e}")
print(f"median p-value, other rows : {np.median(res.p_value[:10]):.2f}")

for r in rs.compare_mode_sets(K, sig, mode_names=["A'1", "2LA(M)", "X"]):
    print(r["names"], f"var_strain {r['var_strain']:.2e}  var_density {r['var_density']:.2e}")
print(f"all three     var_strain {res.strain_sigma[0, 0]**2:.2e}  var_density {res.density_sigma[0, 0]**2:.2e}")

degrees of freedom of the check: 1
median p-value, rows 10-14 : 1.0e-43
median p-value, other rows : 0.50
["A'1", '2LA(M)'] var_strain 2.29e-05  var_density 4.77e-03
['2LA(M)', 'X'] var_strain 2.29e-05  var_density 1.31e-02
["A'1", 'X'] var_strain 6.74e-04  var_density 1.26e-02
all three     var_strain 2.21e-05  var_density 3.65e-03

With more peaks than causes, MultiModeModel uses weighted least squares (each peak weighted by 1/sigma^2). The extra peak shrinks the error bars and makes the chi-square check possible: the rows with the hidden extra shift stand out with p-values near zero, while elsewhere the p-values scatter around 0.5 as they should. With exactly two peaks there is no check (dof is 0, chi2_map and p_value are None). compare_mode_sets ranks every subset of the modes (pairs by default) that can separate the two causes by the sum of the two variances (squared error bars), smallest first.

5. Smoothing a whole map

import numpy as np
import ramansep as rs

K = np.array([[-2.3, 1.1], [-0.9, -1.7], [0.4, 2.2]])   # illustrative, non-physical
sig = np.array([0.03, 0.05, 0.02])                      # cm^-1

# A smooth illustrative truth on a 12 x 12 map, measured with noise.
yy, xx = np.mgrid[0:12, 0:12]
strain = 0.002 * xx                                     # gentle gradient
density = 0.01 + 0.001 * yy
rng = np.random.default_rng(7)
shifts = np.einsum("mj,jhw->mhw", K, np.stack([strain, density])) \
    + rng.normal(size=(3, 12, 12)) * sig[:, None, None]

for lam in (0.0, 1000.0, 10000.0):
    r = rs.bayesian_map_inversion(K, shifts, sig, lam_strain=lam, posterior_sigma=True)
    err = np.sqrt(np.mean((r.strain - strain) ** 2))
    print(f"lam = {lam:6.0f}: mean strain sigma {r.strain_sigma.mean():.4f}, "
          f"RMS strain error {err:.4f}")

shifts[0, 2, 3] = np.nan                                # a masked pixel
try:
    rs.bayesian_map_inversion(K, shifts, sig, lam_strain=100.0)
except ValueError as err:
    print("refused:", str(err)[:60], "...")

lam =      0: mean strain sigma 0.0123, RMS strain error 0.0115
lam =   1000: mean strain sigma 0.0101, RMS strain error 0.0078
lam =  10000: mean strain sigma 0.0057, RMS strain error 0.0039
refused: shifts contain non-finite values; the smoothness prior coupl ...

bayesian_map_inversion solves the whole map at once, adding a penalty lam times the squared differences between neighbouring pixels (a "smoothness prior"; lam_strain and lam_density can differ). It is one sparse linear system, solved exactly. At lam = 0 it gives the same answer as the pixel-by-pixel inversion, and a larger lam never increases the error bars. Smoothing helps only when the real fields are smooth, as they are here by construction; you choose lam, and the package does not estimate it for you. Masked (NaN) pixels are refused by default because the prior links every pixel to its neighbours; fill_masked=True fills them from their neighbours instead, and sigmas may be given per pixel (example 14). Exact per-pixel error bars (posterior_sigma=True) need a dense matrix inverse and are refused above max_dense = 4096 unknowns (two per pixel) unless you raise that limit.

6. Measuring the lever arms on your own instrument

import numpy as np
import ramansep as rs

# Reference states you prepared: strain (%) from a strain stage,
# electron density (1e13 cm^-2) from a gate. Illustrative values.
strain = np.array([0.0, 0.4, 0.8, 0.0, 0.3, 0.6, 0.2, 0.0])
density = np.array([0.0, 0.0, 0.0, 0.9, 0.6, 0.3, 1.2, 0.5])

# Pretend the true lever arms are the MoS2 A'1 / 2LA(M) values and
# measure the shifts with 0.1 cm^-1 noise.
K_true = np.array([[-5.1, -2.2], [-20.9, 0.0]])
rng = np.random.default_rng(1)
shifts = np.column_stack([strain, density]) @ K_true.T + rng.normal(0, 0.1, (8, 2))

cal = rs.calibrate_lever_arms(strain, density, shifts, sigmas=0.1,
                              mode_names=["A'1", "2LA(M)"])
np.set_printoptions(precision=2, suppress=True)
print("fitted lever arms:\n", cal.K)
print("their error bars:\n", cal.K_sigma)
print("chi2 / dof per mode:", cal.chi2 / cal.chi2_dof)

model = rs.SeparationModel(cal.coefficients(reference="stage + gate calibration, illustrative"))
print(model.coefficients.reference)

try:                                             # a strain-only sweep
    rs.calibrate_lever_arms(strain[:3], density[:3], shifts[:3], sigmas=0.1)
except ValueError as err:
    print("refused:", err)

fitted lever arms:
 [[ -5.    -2.27]
 [-20.9    0.03]]
their error bars:
 [[0.09 0.06]
 [0.09 0.06]]
chi2 / dof per mode: [0.09 0.56]
stage + gate calibration, illustrative
refused: reference states are collinear in the (strain, density) plane, so the strain and density lever arms are not separately identifiable; add a reference point off that line (e.g. a gated point to a strain-only sweep)

calibrate_lever_arms fits each mode's two lever arms by weighted least squares from reference states of known strain and density. With sigmas given, the error bars are the exact result for that noise level, and chi2 / dof near 1 says the straight-line model and your stated errors agree. Without sigmas, the error bars are estimated from the scatter of the fit, which needs at least 3 reference points. A sweep that changes strain only (or strain and charge in fixed proportion) cannot tell the two lever arms apart, so it is refused. This fit assumes the shifts are measured from the true pristine frequencies and that strain and density are known in absolute terms. If they are known only relative to the starting state of the sample, use fit_offset=True (example 17). cal.K of a calibration with three or more modes goes straight into MultiModeModel.

7. Planning the calibration before measuring it

import numpy as np
from ramansep import plan_calibration, design_references, repeats_for_sigma

# Four planned reference states, 0.15 cm^-1 expected shift error (illustrative).
plan = plan_calibration(strain=[0.0, 0.3, 0.6, 0.0],
                        density=[0.0, 0.0, 0.1, 0.8], sigmas=0.15)
print("identifiable:", plan["identifiable"])
print("expected error bars (strain arm, density arm):", plan["K_sigma"].round(3))

# States the stage and gate can reach; which 4 are worth measuring?
eps = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.0, 0.0, 0.25])
rho = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.8, 0.8])
pick = design_references(eps, rho, n_pick=4, sigmas=0.15)
print("measure candidates:", pick["indices"])

r, plan_r = repeats_for_sigma(0.1, eps[pick["indices"]], rho[pick["indices"]], sigmas=0.15)
print(f"repeat that set {r} times -> worst error bar {plan_r['K_sigma'].max():.3f}")

line = plan_calibration([0.1, 0.2, 0.4], [0.0, 0.0, 0.0], sigmas=0.15)
print("strain-only sweep identifiable:", line["identifiable"])

identifiable: True
expected error bars (strain arm, density arm): [0.225 0.187]
measure candidates: [8, 5, 7, 4]
repeat that set 6 times -> worst error bar 0.092
strain-only sweep identifiable: False

Because the calibration is a linear fit, the error bars it will report depend only on which reference states you prepare and on the shift errors, not on the measured shifts. plan_calibration computes them in advance. design_references picks states one at a time, each time taking the one that shrinks the joint uncertainty most (the "D-optimal" rule; F. Pukelsheim, Optimal Design of Experiments, SIAM (2006)); this is a good heuristic, not a guarantee of the best possible subset. repeats_for_sigma uses the fact that repeating the whole set r times divides the covariance by r. design_references and repeats_for_sigma refuse a candidate set that lies on one line.

8. How much of the error bar comes from the calibration

import numpy as np
from ramansep import calibrate_lever_arms, separation_with_calibration

# Five reference states (illustrative) and noise-free shifts from the
# MoS2 A'1 / 2LA(M) lever arms; the calibration is told they carry 0.2 cm^-1 errors.
K_true = np.array([[-5.1, -2.2], [-20.9, 0.0]])
eps = np.array([0.0, 0.4, 0.8, 0.0, 0.3])
rho = np.array([0.0, 0.0, 0.2, 0.9, 0.6])
cal = calibrate_lever_arms(eps, rho, np.column_stack([eps, rho]) @ K_true.T, sigmas=0.2)

out = separation_with_calibration(cal, dw1=-2.0, dw2=-8.0, sigma1=0.1, sigma2=0.1)
print(f"strain  {float(out['strain']):.4f} %")
print(f"  error from the spectra     {float(out['strain_sigma_shifts']):.4f}")
print(f"  error from the calibration {float(out['strain_sigma_calibration']):.4f}")
print(f"  total                      {float(out['strain_sigma']):.4f}")
print(f"density {float(out['density']):.4f} x 1e13 cm^-2")
print(f"  error from the spectra     {float(out['density_sigma_shifts']):.4f}")
print(f"  error from the calibration {float(out['density_sigma_calibration']):.4f}")
print(f"  total                      {float(out['density_sigma']):.4f}")

strain  0.3828 %
  error from the spectra     0.0048
  error from the calibration 0.0040
  total                      0.0063
density 0.0217 x 1e13 cm^-2
  error from the spectra     0.0468
  error from the calibration 0.0396
  total                      0.0613

Other inversions treat the lever arms as exact. separation_with_calibration also carries the calibration's own uncertainty into the result, and reports the two parts separately; the total is their quadrature sum (square root of the sum of squares). It uses the first-order (linearised) error propagation, based on the identity dx/dK_mj = -K^-1 E_mj x for the two-mode solution x = K^-1 y (E_mj is the matrix with a single 1 at row m, column j). It takes two-mode calibrations only; for three or more modes use multimode_with_calibration (example 15). The calibration part comes from one set of lever arms shared by every pixel, so it is an error common to the whole map: unlike the spectral noise, it does not shrink when you average many pixels.

9. Temperature from the anti-Stokes / Stokes ratio

from ramansep import (HC_OVER_KB_CM_K, anti_stokes_ratio,
                      calibrate_anti_stokes, temperature_from_anti_stokes)

omega = 385.0                  # the mode's Raman shift, cm^-1
print(f"hc/kB = {HC_OVER_KB_CM_K:.7f} cm K")

# Calibrate once, on a spot you know is at 295 K (illustrative ratio).
C = calibrate_anti_stokes(ratio_measured=0.13, omega_cm1=omega, temperature_K=295.0)
print(f"calibration constant C = {C:.4f}")

# Later, under stronger laser power, the ratio has risen.
out = temperature_from_anti_stokes(0.17, omega, C, ratio_sigma=0.005)
print(f"T = {out['temperature']:.1f} +/- {out['temperature_sigma']:.1f} K")
print(f"with C wrongly set to 1: T = {temperature_from_anti_stokes(0.17, omega, 1.0):.1f} K")

try:
    temperature_from_anti_stokes(1.2 * C, omega, C)
except ValueError as err:
    print("refused:", str(err)[:60], "...")

hc/kB = 1.4387769 cm K
calibration constant C = 0.8500
T = 344.2 +/- 6.3 K
with C wrongly set to 1: T = 312.6 K
refused: ratio >= calibration constant: the Bose factor is strictly < ...

The model is I_AS / I_S = C exp(-hc omega / (kB T)), where omega is the mode's Raman shift. The constant C absorbs the instrument's response and the scattering factors; it is close to, but not exactly, 1, and assuming 1 gives a wrong temperature (as the last printed line shows). So C is measured once at a known temperature. The constant hc/kB is computed from the exact SI values of h, c and kB, not typed in by hand. A ratio at or above C has no finite positive temperature and is refused; it points to a wrong calibration or background subtraction. anti_stokes_ratio is the forward model.

10. Strain, charge and temperature together

import numpy as np
from ramansep import ThreeCauseModel

# Four modes; columns: shift per % strain, per 1e13 cm^-2, per kelvin.
# Invented, illustrative values -- no temperature coefficients ship with the package.
K = np.array([[-2.1, -0.8, -0.011],
              [-1.0, -2.3, -0.016],
              [-3.2, -0.4, -0.007],
              [-0.6, -1.1, -0.021]])
model = ThreeCauseModel(K)
shifts = model.forward(strain=0.2, density=0.5, temperature=40.0)   # 40 K of heating
res = model.invert(shifts, sigmas=[0.05, 0.05, 0.05, 0.05])
print(f"strain {float(res.strain):.3f} +/- {float(res.strain_sigma):.3f} %")
print(f"density {float(res.density):.3f} +/- {float(res.density_sigma):.3f} x 1e13 cm^-2")
print(f"temperature {float(res.temperature):.1f} +/- {float(res.temperature_sigma):.1f} K")

strain 0.200 +/- 0.016 %
density 0.500 +/- 0.042 x 1e13 cm^-2
temperature 40.0 +/- 4.2 K

Heating shifts peaks too, and a two-cause analysis would count that shift as strain or charge. ThreeCauseModel adds a third column of lever arms, the temperature coefficients (cm^-1 per kelvin) that you have calibrated for your material, and needs at least three modes. Because shifts are measured from the reference frequencies, the "temperature" it returns is the change from the temperature at which those references hold. With four or more modes it also returns the chi-square check (chi2_map, p_value) and the full 3 x 3 covariance per pixel. A pixel with a NaN shift or error bar (a pixel masked by fit_map) comes out as NaN without affecting the others. Measuring the temperature with example 9 and comparing it with this estimate is a direct cross-check.

11. Graphene: the G + 2D decomposition

import warnings
import ramansep as rs

coeffs = rs.graphene_g_2d_lee2012()
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    rs.SeparationModel(coeffs)                   # default condition_warn = 30
print("warning:", str(caught[0].message)[:62], "...")

model = rs.SeparationModel(coeffs, condition_warn=100)
res = model.invert(dw1=0.3, dw2=-6.84)          # G and 2D shifts, cm^-1 (illustrative)
print(f"strain {float(res.strain):.3f} %, hole-doping G shift {float(res.density):.2f} cm^-1")

warning: coefficient matrix condition number is 91.5; the mode pair is  ...
strain 0.200 %, hole-doping G shift 5.00 cm^-1

The G peak has barely moved here, but only because strain and doping pushed it in opposite directions; the 2D peak separates the two. The condition number of this set (91.5) is high mainly because its two columns are in different units (percent and cm^-1), so the docstring advises condition_warn=100. The noise amplification is moderate: a 1 cm^-1 error on both shifts costs about 0.035 % of strain and 1.6 cm^-1 of doping coordinate (asserted in the tests). The "density" here is a G-band shift, not a carrier density; see below.

12. Pixels with no peak at all

import warnings
import numpy as np
import ramansep as rs
from ramansep import lorentzian

# An illustrative 2 x 20 map: row 0 is off the flake (background and
# noise only, no Raman peak), row 1 carries both MoS2 peaks.
x = np.linspace(380.0, 480.0, 401)
rng = np.random.default_rng(6)
cube = np.empty((2, 20, x.size))
cube[0] = 50.0 + rng.normal(0.0, 3.0, (20, x.size))
cube[1] = (lorentzian(x, 404.7, 3.0, 40.0, 50.0) + lorentzian(x, 452.0, 4.0, 25.0, 0.0)
           + rng.normal(0.0, 3.0, (20, x.size)))

with warnings.catch_warnings():
    warnings.simplefilter("ignore", RuntimeWarning)   # overflow in runaway noise fits
    plain = rs.fit_map(x, cube, (395, 415), (440, 465), 404.7, 452.0)
    checked = rs.fit_map(x, cube, (395, 415), (440, 465), 404.7, 452.0, min_snr=3.0)

for name, f in (("without min_snr", plain), ("with min_snr=3 ", checked)):
    print(f"{name}: kept off-flake {f.ok[0].sum()} of 20, on-flake {f.ok[1].sum()} of 20")
for j in np.flatnonzero(plain.ok[0]):
    print(f"  off-flake pixel {j:2d} kept without min_snr: "
          f"dw1 {plain.dw1[0, j]:+.2f} +/- {plain.sigma1[0, j]:.2f} cm^-1")
for code in np.unique(checked.reason[0]):
    print(f"  off-flake reason {code}: {(checked.reason[0] == code).sum():2d} x {rs.MASK_REASONS[code]}")

without min_snr: kept off-flake 3 of 20, on-flake 20 of 20
with min_snr=3 : kept off-flake 0 of 20, on-flake 20 of 20
  off-flake pixel 11 kept without min_snr: dw1 +9.92 +/- 0.05 cm^-1
  off-flake pixel 15 kept without min_snr: dw1 -7.72 +/- 0.11 cm^-1
  off-flake pixel 17 kept without min_snr: dw1 -0.32 +/- 1.14 cm^-1
  off-flake reason 3: 13 x fit did not converge
  off-flake reason 4:  1 x fitted centre outside its window
  off-flake reason 5:  3 x fitted peak height not positive (a dip or no peak)
  off-flake reason 6:  3 x peak height below min_snr times its error bar

A fit to a spectrum that holds no peak (off the flake, over a hole) sometimes converges on a noise spike and returns a shift that looks like a measurement: pixel 11 above reports +9.92 +/- 0.05 cm^-1 from pure noise. min_snr=3 masks a pixel unless each fitted peak height is at least 3 times its own error bar (a signal-to-noise test), and removes all such pixels here while keeping every real one; the right threshold for your data is your choice. fwhm_range=(lo, hi) adds a check on the fitted widths. Both are off by default. Before 0.12.0, peaks with a negative fitted height were kept too: on this map 6 off-flake pixels were kept instead of 3.

13. Two peaks that overlap

import numpy as np
import ramansep as rs
from ramansep import lorentzian

# Two broad, overlapping peaks at known positions 385 and 404 cm^-1
# (illustrative widths), no noise: both shifts should come out as zero.
x = np.linspace(360.0, 430.0, 561)
y = lorentzian(x, 385.0, 5.0, 600.0, 20.0) + lorentzian(x, 404.0, 6.0, 1000.0, 0.0)

for joint in (False, True):
    dw1, dw2, s1, s2, f1, f2 = rs.fit_two_modes(x, y, (375, 394), (395, 415),
                                                ref1=385.0, ref2=404.0, joint=joint)
    print(f"joint={joint!s:5}: dw1 {dw1:+.4f}  dw2 {dw2:+.4f} cm^-1")
print(f"correlation of the two fitted centres: {f1.center_correlation:.3f}")

joint=False: dw1 +0.0795  dw2 -0.0279 cm^-1
joint=True : dw1 +0.0000  dw2 +0.0000 cm^-1
correlation of the two fitted centres: 0.014

Fitting each peak in its own window leaves the tail of the other peak in that window, and the tail pulls the fitted centre: here by 0.08 cm^-1, comparable to the 0.1 cm^-1 shift error used in the examples above, and as a systematic offset rather than random scatter. joint=True fits both peaks together, as two Lorentzians on one shared baseline, over the whole range from the lower edge of the first window to the upper edge of the second; that range must not contain a third peak. The two fitted centres are then slightly correlated (center_correlation); the inversion treats the two shifts as independent, which is a good approximation when this number is small. fit_map(..., joint=True) does the same at every pixel.

14. Smoothing a map that has gaps

import numpy as np
import ramansep as rs

K = np.array([[-2.3, 1.1], [-0.9, -1.7], [0.4, 2.2]])   # illustrative, non-physical
rng = np.random.default_rng(7)
yy, xx = np.mgrid[0:12, 0:12]
strain = 0.002 * xx
density = 0.01 + 0.001 * yy

# One error bar per mode AND pixel (as fit_map returns them), and a
# 3 x 3 patch of pixels whose fits failed.
sig = np.array([0.03, 0.05, 0.02])[:, None, None] * rng.uniform(0.5, 2.0, (3, 12, 12))
shifts = np.einsum("mj,jhw->mhw", K, np.stack([strain, density])) \
    + rng.normal(size=(3, 12, 12)) * sig
shifts[:, 4:7, 4:7] = np.nan

r = rs.bayesian_map_inversion(K, shifts, sig, lam_strain=1000.0,
                              fill_masked=True, posterior_sigma=True)
hole = np.zeros((12, 12), bool)
hole[4:7, 4:7] = True
print(f"RMS strain error, measured pixels: {np.sqrt(np.mean((r.strain - strain)[~hole] ** 2)):.4f}")
print(f"RMS strain error, filled pixels  : {np.sqrt(np.mean((r.strain - strain)[hole] ** 2)):.4f}")
print(f"strain sigma: hole centre {r.strain_sigma[5, 5]:.4f}, "
      f"measured pixels (median) {np.median(r.strain_sigma[~hole]):.4f}")
print("filled centre = mean of its 4 neighbours:",
      np.isclose(r.strain[5, 5], (r.strain[4, 5] + r.strain[6, 5] + r.strain[5, 4] + r.strain[5, 6]) / 4))

RMS strain error, measured pixels: 0.0077
RMS strain error, filled pixels  : 0.0064
strain sigma: hole centre 0.0197, measured pixels (median) 0.0116
filled centre = mean of its 4 neighbours: True

bayesian_map_inversion accepts (since 0.12.0) one error bar per mode and per pixel (an array of the same shape as shifts), such as the sigma1, sigma2 maps of fit_map. With fill_masked=True, a NaN shift (or NaN error bar) counts as a missing measurement: it gets zero weight, and the smoothness prior fills the gap from the neighbouring pixels. A pixel with no data at all ends up exactly at the average of its neighbours, and its error bar is larger than that of a measured pixel. The filled values are interpolations, not measurements; they were accurate here because the true fields are smooth by construction. Filling needs both smoothness weights above zero.

15. Calibration uncertainty with three or more modes

import numpy as np
import ramansep as rs

# Three modes calibrated together on six reference states (illustrative);
# the shifts are generated from K_true with 0.1 cm^-1 noise.
K_true = np.array([[-5.1, -2.2], [-20.9, 0.0], [-8.0, -1.1]])
eps = np.array([0.0, 0.4, 0.8, 0.0, 0.3, 0.6])
rho = np.array([0.0, 0.0, 0.2, 0.9, 0.6, 0.3])
rng = np.random.default_rng(4)
ref_shifts = np.column_stack([eps, rho]) @ K_true.T + rng.normal(0, 0.1, (6, 3))
cal = rs.calibrate_lever_arms(eps, rho, ref_shifts, sigmas=0.1)

out = rs.multimode_with_calibration(cal, shifts=[-2.0, -8.0, -3.4],
                                    sigmas=[0.15, 0.10, 0.12])
print(f"strain  {float(out['strain']):.4f} %")
print(f"  error from the spectra     {float(out['strain_sigma_shifts']):.4f}")
print(f"  error from the calibration {float(out['strain_sigma_calibration']):.4f}")
print(f"  total                      {float(out['strain_sigma']):.4f}")
print(f"density {float(out['density']):.4f} x 1e13 cm^-2 +/- {float(out['density_sigma']):.4f}")
print(f"p-value of the model check: {float(out['p_value']):.2f}")

strain  0.3821 %
  error from the spectra     0.0046
  error from the calibration 0.0016
  total                      0.0049
density 0.1036 x 1e13 cm^-2 +/- 0.0642
p-value of the model check: 0.06

multimode_with_calibration is the counterpart of example 8 for any number of modes: the weighted least-squares separation of example 4, with the calibration's own uncertainty carried into the error bars, to first order. With more modes than causes the derivative of the answer with respect to a lever arm also depends on the leftover mismatch (residual) of the fit; the docstring gives the formula, and the tests check it against a direct numerical derivative. The chi-square check (chi2_map, p_value) is the one of MultiModeModel, which treats the lever arms as exact.

16. Reading a map exported by the instrument

import os, tempfile
import numpy as np
import ramansep as rs
from ramansep import lorentzian

# Write an illustrative 3 x 4 map the way mapping software often exports
# it: tab-separated "x y wavenumber counts", stage positions in
# micrometers, a serpentine scan, and the wavenumber running high to low.
rng = np.random.default_rng(2)
x_um = np.array([0.0, 0.5, 1.0, 1.5])
y_um = np.array([10.0, 10.5, 11.0])
wn = np.linspace(480.0, 380.0, 401)
path = os.path.join(tempfile.mkdtemp(), "map_export.txt")
with open(path, "w") as fh:
    fh.write("#X\t#Y\t#Wave\t#Intensity\n")
    for i, y in enumerate(y_um):
        for x in (x_um if i % 2 == 0 else x_um[::-1]):
            counts = (lorentzian(wn, 404.7 + 0.4 * x, 3.0, 800.0, 50.0)
                      + lorentzian(wn, 452.0, 4.0, 500.0, 0.0)
                      + rng.normal(0.0, 3.0, wn.size))
            for w, c in zip(wn, counts):
                fh.write(f"{x}\t{y}\t{w:.2f}\t{c:.1f}\n")

wavenumber, cube, xs, ys = rs.load_map_table(path)
print("cube:", cube.shape, " axis:", wavenumber[0], "to", wavenumber[-1], "cm^-1")
print("x (um):", xs, " y (um):", ys)
fits = rs.fit_map(wavenumber, cube, (395, 415), (440, 465), ref1=404.7, ref2=452.0)
print("A'1 shift at y = 10 um (cm^-1):", " ".join(f"{v:+.2f}" for v in fits.dw1[0]))

with open(path) as fh:                           # a scan that stopped early
    lines = fh.readlines()
with open(path, "w") as fh:
    fh.writelines(lines[:-401])
try:
    rs.load_map_table(path)
except ValueError as err:
    print("refused:", str(err).split(";")[0])

cube: (3, 4, 401)  axis: 380.0 to 480.0 cm^-1
x (um): [0.  0.5 1.  1.5]  y (um): [10.  10.5 11. ]
A'1 shift at y = 10 um (cm^-1): -0.00 +0.20 +0.40 +0.60
refused: the stage positions do not form a complete rectangular grid: 1 of the 3 x 4 combinations of distinct y and x values are missing (first: x = 1.5, y = 11.0)

Mapping software rarely writes the CSV layout of example 3. A common plain-text export has four columns, the stage position x, the stage position y, the wavenumber and the counts, with one line per pixel and spectral point. load_map_table reads such a file in any line order (raster or serpentine scans) and with the wavenumber running in either direction. It rebuilds the pixel grid from the stage positions: cube[i, j] is the spectrum at x[j], y[i], and both position arrays increase. Columns may be separated by tabs, spaces or commas. Lines starting with # and one text header line are skipped. If your export has its columns in another order, say which with usecols. A scan that is not a complete rectangle (here, one stopped before the last pixel), a pixel with its own wavenumber axis, or a repeated point is refused rather than padded or interpolated. Stage positions that differ only in the last digits can be merged with decimals=. load_spectrum_table reads a single two-column spectrum the same way. Check that your export matches this layout before relying on it; the package does not read vendor binary formats.

17. Calibrating when the starting state is unknown

import numpy as np
import ramansep as rs

# A stage + gate calibration (illustrative). Strain (%) and electron
# density (1e13 cm^-2) are known only relative to the starting state, and
# the peak positions are absolute frequencies (cm^-1).
eps = np.array([0.0, 0.2, 0.4, 0.6, 0.0, 0.3, 0.1])
rho = np.array([0.0, 0.0, 0.0, 0.1, 0.8, 0.5, 1.1])
K_true = np.array([[-5.1, -2.2], [-20.9, 0.0]])
start = np.array([403.6, 451.1])        # the starting state: unknown in practice
rng = np.random.default_rng(2)
freq = start + np.column_stack([eps, rho]) @ K_true.T + rng.normal(0, 0.1, (7, 2))

np.set_printoptions(precision=2, suppress=True)
old = rs.calibrate_lever_arms(eps, rho, freq - [403.0, 452.0], sigmas=0.1)
print("through zero, shifts from 403.0 and 452.0 cm^-1:\n", old.K)
print("  chi2 / dof:", old.chi2 / old.chi2_dof)
new = rs.calibrate_lever_arms(eps, rho, freq, sigmas=0.1, fit_offset=True)
print("with a fitted offset per mode:\n", new.K)
print("  error bars:\n", new.K_sigma)
print("  chi2 / dof:", new.chi2 / new.chi2_dof)
print("  starting frequencies:", new.offset, "+/-", new.offset_sigma)
print("planned error bars:", rs.plan_calibration(eps, rho, 0.1, fit_offset=True)["K_sigma"])

through zero, shifts from 403.0 and 452.0 cm^-1:
 [[ -3.8   -1.67]
 [-22.62  -0.87]]
  chi2 / dof: [13.5  34.69]
with a fitted offset per mode:
 [[ -5.05  -2.23]
 [-20.57   0.03]]
  error bars:
 [[0.2 0.1]
 [0.2 0.1]]
  chi2 / dof: [0.95 1.35]
  starting frequencies: [403.63 450.98] +/- [0.08 0.08]
planned error bars: [0.2 0.1]

In a real calibration, the stage reports strain relative to where it started and the gate reports charge relative to the starting voltage. The starting strain and doping of the flake, and so the pristine frequencies, are unknown. The fit of example 6 goes through zero and needs them. Here the shifts were measured from guessed references (403.0 and 452.0 cm^-1, both off), and the lever arms come out wrong by up to 1.7 cm^-1 per percent. The fit itself says so: chi2 / dof is far above 1. With fit_offset=True each mode gets one more free parameter, a constant offset (the shift at the calibration's own zero). The lever arms then depend only on how the shifts change between the reference states, and they come out within two error bars of the true values, with chi2 / dof near 1. The offset is reported with its error bar, but it is not the pristine frequency, and the inversions do not use it. The error bars of the lever arms already include the uncertainty of the offset, so separation_with_calibration and multimode_with_calibration take the result as it is. The price of the extra parameter: the reference states must not lie on any one straight line (not only a line through zero), and the fit needs one more reference point. plan_calibration and repeats_for_sigma take fit_offset=True too, so the error bars can be planned in advance.

What is in the package

Every name below is exported from ramansep (its __all__). Each docstring (help(ramansep.fit_map), for example) gives inputs, units and conventions.

Two-mode inversion

SeparationModel(coefficients, condition_warn=30.0) -- inverts two shift maps (invert(dw1, dw2, sigma1, sigma2)) into a SeparationResult with strain, density, strain_sigma, density_sigma, correlation and condition_number; forward predicts shifts. Error bars are computed only when both sigma1 and sigma2 are given.
ModeCoefficients -- the four lever arms of two modes plus the required reference string.

Coefficient sets

mos2_a1_2la(), mos2_eprime_a1(), graphene_g_2d_lee2012() -- the cited sets described in the next section.
synthetic_demo() -- arbitrary, non-physical values for tests and demos only.

More modes, and whole maps

MultiModeModel(K, mode_names), MultiModeResult -- weighted least squares for any number (two or more) of modes, with chi2_map, p_value and dof.
compare_mode_sets(K, sigmas, mode_names, subset_size=2) -- ranks mode subsets by delivered variance.
bayesian_map_inversion(K, shifts, sigmas, lam_strain, lam_density, posterior_sigma, max_dense, fill_masked), BayesianMapResult -- joint inversion of a whole map with a smoothness prior (a Gaussian Markov random field, Rue and Held 2005: a standard statistical way of saying that each pixel is probably close to its four neighbours).

From spectra to shifts

lorentzian, fit_lorentzian, PeakFit -- Lorentzian peak shape and a least-squares fit by the Levenberg-Marquardt method (a standard step-by-step curve-fitting method), using exact formulas for the derivatives; baseline="constant" or "linear". Parameter errors are the standard linearised estimates, which are reliable when the residuals are dominated by uncorrelated noise.
voigt, fit_voigt, VoigtFit -- the Voigt shape, computed exactly with the Faddeeva function (a standard special function, scipy.special.wofz), normalised to peak height.
fit_two_modes(x, y, window1, window2, ref1, ref2, baseline, joint=False) -- fits both peaks, each in its own non-overlapping window, or with joint=True both together (example 13), and returns (dw1, dw2, sigma1, sigma2, fit1, fit2).
fit_map(..., joint=False, min_snr=None, fwhm_range=None), MapFitResult, MASK_REASONS -- the same at every pixel of a cube, with masking of failed pixels and a per-pixel reason code (example 12).
load_spectrum_csv, save_spectrum_csv, load_map_csv, save_map_csv -- the documented CSV formats of example 3.
load_spectrum_table(path, usecols, delimiter), load_map_table(path, usecols, delimiter, decimals): plain text tables as instruments export them, with stage positions for maps (example 16, new in 0.13.0).

Calibration

calibrate_lever_arms(strain, density, shifts, sigmas, mode_names, fit_offset=False), CalibrationResult (K, K_sigma, cov, chi2, chi2_dof, n_points, condition_number, mode_names, coefficients(reference) for two modes, and, with fit_offset=True, offset, offset_sigma and full_cov; example 17, new in 0.13.0).
plan_calibration, design_references, repeats_for_sigma -- planning, example 7 (plan_calibration and repeats_for_sigma also take fit_offset=True).
separation_with_calibration -- example 8; multimode_with_calibration -- the same for three or more modes, example 15.

Temperature

ThreeCauseModel, ThreeCauseResult -- example 10.
anti_stokes_ratio, calibrate_anti_stokes, temperature_from_anti_stokes, HC_OVER_KB_CM_K (hc/kB in cm K) -- example 9.

Cited coefficient sets

Three sets ship with full provenance in their reference field and docstring. Each docstring states which number comes from which source and the conditions it applies under. Check that your sample matches before use.

mos2_a1_2la(): the A'1 + 2LA(M) pair for monolayer MoS2. Strain arms -5.1 and -20.9 cm^-1 per percent biaxial strain, computed by frozen-phonon DFT in the source paper (Mahim and Rahman, under review; data doi:10.5281/zenodo.21778170). A'1 doping response -2.2 cm^-1 per 1e13 cm^-2 of electrons, measured with a gate by Chakraborty et al., Phys. Rev. B 85, 161403(R) (2012). The 2LA(M) doping response has never been measured and is set to zero; the source paper bounds the consequence at 16% of a recovered edge charge. Stated conditions: 1H monolayer (the common semiconducting crystal form), 532 nm laser, on an SiO2 substrate, 300 K. The tests reproduce the source paper's published separations from this set (example 1).
mos2_eprime_a1(): the historically used all-optical E' + A'1 pair (Michail et al., Appl. Phys. Lett. 108, 173102 (2016)). Strain arms -4.31 (E') and -2.50 (A'1) cm^-1 per percent, from the measured biaxial Grueneisen parameters 0.56 and 0.31 of Michail et al., ACS Appl. Mater. Interfaces 16, 49602 (2024), converted with the mode frequencies 385.0 and 403.0 cm^-1. These are rounded monolayer values (see Mignuzzi et al., Phys. Rev. B 91, 195411 (2015)); Michail et al. 2016 place their strain-doping origin at 384.6 and 402.7 cm^-1, which gives the same arms to the digits shown. Doping responses -0.33 (E') and -2.2 (A'1) cm^-1 per 1e13 cm^-2, measured by Chakraborty et al. (2012). Stated conditions: 1H monolayer, on SiO2, 300 K. Both modes respond weakly to strain, so it carries about five times the strain uncertainty (example 2) -- included for comparison.
graphene_g_2d_lee2012(): the G + 2D decomposition of Lee et al., Nat. Commun. 3, 1024 (2012). Its "doping" output is deliberately a G-band shift, not a density: graphene's doping response is nonlinear, so no universal linear rate exists to ship, and the docstring points to the gated calibrations needed for conversion. Stated for hole doping, uniaxial strain of random orientation.

Tungsten-based materials (2024-2026). For W-based monolayers the literature now provides verified strain and temperature responses -- monolayer WSe2 biaxial strain rates (Michail et al., ACS Appl. Mater. Interfaces 16, 49602 (2024)), monolayer WS2 strain rates (Roy, Yang and Gao, Sci. Rep. 14, 3860 (2024)) and temperature coefficients (Huang et al., Sci. Rep. 6, 32236 (2016)) -- while gate-calibrated doping arms remain without a clean linear coefficient (the measured behavior is an electron-only, threshold-like softening of the out-of-plane modes, the vibrations in which atoms move perpendicular to the sheet; Sohier et al., PRX 9, 031019 (2019)). That is why no W-material set ships as constants here: the strain column is citable, the doping column is not yet, and this package does not ship half a lever-arm matrix. Calibrate your own with calibrate_lever_arms -- the doping arm from your own gated reference points -- and the provenance travels with the analysis.

For monolayer MoS2 temperature coefficients, the thermal module docstring points to S. Sahoo et al., J. Phys. Chem. C 117, 9042 (2013); no temperature coefficient is shipped.

When it refuses, and why

ramansep raises ValueError instead of guessing when:

a lever-arm matrix contains NaN or infinity; a two-mode lever-arm matrix is not 2 x 2 or is exactly singular (the two modes respond identically); a merely poorly conditioned one (condition number above condition_warn, default 30) gives a warning, not an error;
a multimode or map lever-arm matrix has rank below 2 (below 3 for ThreeCauseModel, which also needs at least three modes), or the weighted problem is singular at some pixel;
the shift maps do not match in shape or number;
a shift error (sigmas) is zero or negative in MultiModeModel, ThreeCauseModel, bayesian_map_inversion, the calibration and planning tools, is zero, negative, NaN or infinite in compare_mode_sets (which also refuses a K that is not (m, 2)), or negative in SeparationModel (zero is allowed there; NaN marks a masked pixel and is passed through as NaN);
bayesian_map_inversion gets NaN or infinite shifts or error bars without fill_masked=True, missing values with a smoothness weight of zero, so little data that the unmasked shifts cannot fix both fields, per-pixel error bars of the wrong shape, a negative, NaN or infinite smoothness weight (NaN and infinity are refused since 0.13.0), or asks for exact error bars above max_dense unknowns;
MultiModeModel or ThreeCauseModel gets a number of error bars other than one per mode, or shift maps of different shapes (since 0.13.0 the message says which);
a peak fit gets too few points (5 for a Lorentzian with a constant baseline, 6 with a linear baseline or for a Voigt), non-finite data, an unknown baseline name, or a non-positive starting width;
the two fit windows are reversed, overlap, or hold too few points, the baseline name is unknown, or a reference frequency ref1, ref2 is NaN or infinite (fit_two_modes, and fit_map before fitting any pixel; the reference check is new in 0.13.0); fit_map also refuses a wavenumber axis that contains NaN or infinity, a min_snr that is not positive, or a fwhm_range that is not (lo, hi) with 0 <= lo < hi;
a CSV file has the wrong header or field count, a wavenumber axis that is not strictly increasing or contains NaN or infinity, a pixel grid with holes, or a pixel with its own wavenumber axis;
a text table (load_spectrum_table, load_map_table) has lines without the requested numeric columns, a NaN or infinite wavenumber or stage position, a wavenumber axis that repeats a value or is not monotonic, stage positions that do not form a complete rectangle, or a pixel with its own wavenumber axis or number of points;
calibration reference states lie on one line through the origin of the (strain, density) plane (the two lever arms cannot be told apart), there are fewer than 2 of them, or fewer than 3 when no sigmas are given; with fit_offset=True, the states lie on any one straight line, or there are fewer than 4 of them without sigmas; sigmas whose shape fits neither the points nor the modes;
CalibrationResult.coefficients is asked to package other than two modes, or separation_with_calibration gets such a calibration (multimode_with_calibration takes any number);
the anti-Stokes ratio is at or above the calibration constant, or a temperature, frequency, ratio or constant is not positive.

How the results are checked

114 automated tests run on every push and pull request, on Python 3.9, 3.10, 3.11, 3.12, 3.13 and 3.14, and once more on Python 3.9 with the oldest NumPy (1.22.0) and SciPy (1.8.0) the package allows. The numerical checks compare the package with an independent calculation, a closed form, a published number, or a simulation with fixed random seeds; the others check that refusals fire. The main checks, with the tolerance each test really uses. "To 1e-12" means the difference may be at most 1e-12; "plus a relative 1e-5" means that 1e-5 times the size of the expected value is also allowed. Where a test does not set both parts, NumPy's own defaults apply, and these are listed too:

Two-mode inversion and the cited sets

Forward-then-invert round trips recover strain and density to 1e-12 plus a relative 1e-7 (synthetic set, both MoS2 sets) and to 1e-12 plus a relative 1e-5 (graphene set).
Reported error bars match the scatter of 200 000 seeded noisy inversions to within 2 %.
MoS2 A'1 + 2LA(M), paper edge spot: zero strain to 1e-12, density 0.5/2.2 x 1e13 cm^-2 to a relative 1e-12, within 0.5e11 of the published 2.3e12 cm^-2. Interior spot: strain 0.134 % within 0.0005, density below 7.6e11 cm^-2 in size.
The E' + A'1 strain arms equal -2 gamma omega / 100 from the cited Grueneisen parameters within 0.005; its strain error bar is more than 4 times that of the A'1 + 2LA(M) pair.
Graphene: pure-strain and pure-doping trajectories decompose to zero doping and zero strain to 1e-12; the noise amplification equals its closed form to a relative 1e-12 plus 1e-8.

Peak fitting and maps

Analytic derivatives match finite differences (Lorentzian: 1e-7 plus a relative 1e-5; Voigt: 1e-6).
A noiseless Lorentzian is recovered from automatic starting values (centre and width to 1e-8, height and offset to 1e-6); a noiseless Voigt to 1e-8.
On 40 seeded noisy spectra, (centre error / reported sigma) has a standard deviation between 0.5 and 2.0 and a mean below 0.6 in size.
The Voigt with zero Lorentzian width equals the Gaussian to 1e-14; as the Gaussian width shrinks tenfold (1e-3 to 1e-4), the gap to the Lorentzian shrinks by a factor of more than 9 (and is below 1e-5 at width 1e-3).
A sloped background biases a constant-baseline centre by more than 3 of its sigmas; the linear baseline recovers the centre to 1e-6.
fit_map gives exactly (bit for bit) the fit_two_modes result at the pixel checked; on a noiseless synthetic cube the maps invert back to the truth (strain to 1e-4, density to 1e-3); a dead pixel is masked and the other pixels' shifts are bit-identical to the clean run; with noise, every strain error is below 5 error bars + 1e-3.
The CSV save/load round trips are exact (bit-identical arrays).
Configuration errors in fit_map (unknown baseline name, reversed or overlapping windows, too few points in a window, a NaN in the wavenumber axis) are refused, and valid settings still give the fit_two_modes result bit for bit (new in 0.11.1).

Multimode and map inversion

With two modes, MultiModeModel equals SeparationModel (maps to 1e-12 plus a relative 1e-5, error bars to a relative 1e-12 plus 1e-8).
On 1600 seeded pixels the mean chi-square is within 0.1 of its degrees of freedom and the mean p-value within 0.05 of 0.5; a hidden extra shift gives median p-values below 1e-6, against above 0.1 elsewhere.
bayesian_map_inversion at lam = 0 equals the per-pixel result (maps to 1e-12, error bars to 1e-13); a constant noiseless truth is recovered to 1e-12 at every lam tried (0, 1, 250); at lam = 1e9 it matches an independently computed pooled estimate to 1e-5; error bars never grow as lam increases (checked at lam = 0, 0.5, 2 and 10). A rank-deficient matrix is refused at lam = 0 and 1 (new in 0.11.1).

Calibration and planning

Noise-free calibration recovers the lever arms to 1e-12, and agrees with an independent numpy.linalg.lstsq solve to 1e-10.
Over 400 seeded calibrations the scatter of the fitted arms matches the reported error bars within 15 %.
The planned covariance equals the one the calibration reports to 1e-15; repeating a design 7 times divides it by 7 to a relative 1e-12; on a test set of candidates where only one is off the line, the greedy design includes it, and it is never worse than any of 30 random subsets of the same size.
separation_with_calibration: with zero calibration covariance it equals the plain inversion to 1e-12 plus a relative 1e-5; its calibration part matches 400 seeded lever-arm draws within 20 %; the derivative identity it uses matches finite differences to a relative 1e-5 plus 1e-10.

Temperature

ThreeCauseModel recovers noise-free strain and density to 1e-9 and temperature to 1e-8, agrees with an independent per-pixel least-squares solve to the same tolerances, and with a known temperature subtracted matches the two-cause solver to 1e-10.
hc/kB agrees with SciPy's CODATA table to a relative 1e-9.
The anti-Stokes thermometer round trip holds to a relative 1e-9, the calibration constant to 1e-12, and the error propagation matches a finite difference to a relative 1e-5.
With a NaN shift at one pixel and a NaN error bar at another, ThreeCauseModel returns NaN at exactly those two pixels and, at every other pixel, bit-identical results to the run without them (new in 0.12.0).

New in 0.12.0 (tests/test_v012.py, 19 tests)

Joint two-peak fit: its exact derivatives match central finite differences to 1e-6 of each derivative's largest value (plus 1e-9), with both baselines. On a noiseless pair 19 cm^-1 apart (widths 5 and 6 cm^-1) the separate-window fits are off by more than 0.05 and 0.02 cm^-1; the joint fit recovers both centres and widths to 1e-8 cm^-1 and heights and offset to 1e-6, also on a sloped background (slope to 1e-8) and with the windows given in reverse order. On 600 seeded noisy spectra of a pair 12 cm^-1 apart, the scatter of each fitted centre is within 12 % of the reported error bar, the mean error is below 4 of its standard errors, and the measured correlation of the two centres is within 0.12 of the reported one (which is above 0.05). fit_map(..., joint=True) equals fit_two_modes(..., joint=True) bit for bit at every pixel checked.
fit_map checks: on 30 seeded peak-free spectra, without min_snr at least one pixel is kept and at least one is masked for a negative height; with min_snr=3 none is kept. All 10 seeded spectra with real (weak) peaks are kept, with shifts and error bars bit-identical to the run without the checks. Constructed spectra produce the reason codes 0, 1, 5 and 7.
bayesian_map_inversion: with per-pixel error bars and lam = 0 it equals MultiModeModel given the same error bars (maps to 1e-12, error bars to 1e-13); per-pixel error bars that are constant per mode equal the scalar path (maps to 1e-12, error bars to 1e-13). With a pixel without data, a missing mode at another pixel and a missing error bar at a third, fill_masked=True agrees with an independent dense least-squares solve of the same objective (numpy.linalg.lstsq) to 1e-10. A filled pixel equals the average of its neighbours to 1e-12 (interior and corner pixel), with a larger error bar than its measured neighbours; a constant truth is recovered through a 2 x 3 hole to 1e-12 at lam = 0.5 and 250.
multimode_with_calibration: with a two-mode calibration it equals separation_with_calibration (all values and error bars, relative 1e-10 plus 1e-13); its derivative formula matches central finite differences of the MultiModeModel re-solve to a relative 1e-6 plus 1e-10 at a point with a non-zero residual; its calibration part matches 400 seeded lever-arm draws within 20 %; with zero calibration covariance it gives the MultiModeModel error bar to a relative 1e-12; masked pixels stay NaN.
New refusals fire: non-finite lever arms (SeparationModel, MultiModeModel, compare_mode_sets), a negative shift error in SeparationModel (NaN still passes as a masked pixel), a zero, negative, NaN or infinite sigma and a K that is not (m, 2) in compare_mode_sets.

New in 0.13.0 (tests/test_v013.py, 18 tests)

Each fix below has a test that fails on 0.12.0 and passes now.
Text tables: files are written line by line from arrays built independently in the test (a serpentine scan with the wavenumber running high to low; a comma-separated file with a text header and the columns in another order). load_map_table returns those arrays exactly (bit for bit), and fit_map on the loaded cube equals fit_map on the original cube bit for bit. Each refusal fires. load_spectrum_table reverses a descending export exactly.
calibrate_lever_arms(..., fit_offset=True): noise-free arms and offsets are recovered to 1e-12. With per-point error bars, the arms and offsets agree with an independent numpy.linalg.lstsq solve to 1e-10, and the covariance agrees to a relative 1e-10 with the closed form of the Frisch-Waugh-Lovell theorem (arms from the weighted-centered design; offset variance 1/sum(w) plus the propagated mean state). Without error bars the covariance equals the residual-variance formula with n - 3 degrees of freedom (relative 1e-9). The through-origin fit of offset data is off by exactly the closed-form projection of the offsets (to 1e-12), more than 1 cm^-1 per percent here; the offset fit removes it. Over 400 seeded calibrations the scatter of the arms and offsets matches the reported error bars within 15 %, and the mean chi-square is within 0.4 of its 4 degrees of freedom. The calibration part of separation_with_calibration matches 400 seeded offset calibrations within 20 %. plan_calibration(..., fit_offset=True) equals the reported covariance to 1e-15, and repeating the design r times divides it by r (relative 1e-12).

Corrections in earlier versions

0.13.0 fixed three silent failures (before and after values from running 0.12.0 and 0.13.0 on the same input).

load_spectrum_csv returned a wavenumber axis containing NaN or infinity, although it promises to refuse an axis that is not strictly increasing: a file with the axis 400, nan, 402, 403 loaded as [400, nan, 402, 403]. It is now refused. (load_map_csv refused such an axis already, but with the message that pixel (0, 0) "carries a different wavenumber axis than pixel (0, 0)"; the message now says what is wrong.)
fit_map with a NaN reference frequency reported every pixel as fitted (ok True at all pixels, n_masked 0) while every shift was NaN, and fit_two_modes returned a NaN shift. Both now refuse a reference that is not a finite number.
bayesian_map_inversion with a NaN or infinite smoothness weight returned all-NaN maps with only SciPy's "Matrix is exactly singular" warning. It now refuses such a weight and suggests a large finite one for the fully pooled limit.

Documentation corrections in 0.13.0. The docstring of MultiModeResult said that condition_number is that of the weighted design matrix when the error bars are scalars; the code has always returned the condition number of the unweighted lever-arm matrix (the value is unchanged; the docstring now says so). The docstring of fit_voigt said its default start splits the width evenly between the Gaussian and Lorentzian parts; it starts with a Gaussian FWHM equal to the estimated width and a Lorentzian FWHM of half of it. The fitting module docstring said the tests recover amplitude and baseline to 1e-8; they check 1e-6. SeparationResult.density was documented as cm^-2 by default; it is in the unit of the lever arms (1e13 cm^-2 for the MoS2 sets).

0.12.0 fixed three silent or misleading behaviours.

fit_map kept pixels whose fitted "peak" had a negative height (a dip fitted to noise). Such pixels are now masked (reason code 5). Good pixels give the same numbers as before; on the map of example 12, the off-flake pixels kept went from 6 of 20 to 3 of 20 (and to 0 with min_snr=3).
ThreeCauseModel.invert refused a whole map, with a misleading "singular" message, when one pixel carried a NaN error bar, which is how fit_map marks masked pixels. Such pixels now come out as NaN, as in MultiModeModel.
SeparationModel.invert and compare_mode_sets accepted negative shift errors and squared the sign away (for example sigma1=-0.1, sigma2=0.1 with mos2_a1_2la() gave a strain error bar of 0.0048 %). Both now refuse them. compare_mode_sets also returned NaN variances, ranked first, for a zero or NaN sigma, and silently left out every subset with an infinite lever arm; it now refuses those inputs too.

0.11.1 fixed two silent failures.

fit_map treated a mistake in its own settings (for example baseline="Linear", or overlapping windows) like a bad pixel: it masked every pixel and returned an all-NaN result without an error. It now refuses such settings before fitting, as fit_two_modes does. Valid settings give the same results as before. fit_map now also refuses a wavenumber axis that contains NaN or infinity; before, spectral points with a NaN wavenumber were silently left out of the fits.
bayesian_map_inversion did not check that the lever-arm matrix can separate the two causes. With a rank-deficient matrix it returned NaN maps (at lam = 0) or arbitrary finite maps (at lam > 0). It now refuses, as MultiModeModel does.

0.8.0: bayesian_map_inversion now refuses NaN shifts; before, it solved around them silently.

Documentation corrections in 0.11.1. The README of 0.11.0 said the tests ran on Python 3.9-3.14, but the CI matrix skipped 3.10 (added now, with an oldest-dependencies job). It called every test "pinned to an exact result" at "machine precision"; many use stated tolerances, now listed above. It said the tests reproduce published separations for all three cited sets; they do so for mos2_a1_2la only. The calibration check described as "an independent QR solve" uses numpy.linalg.lstsq. It cited the E' + A'1 pair itself to Michail et al. (2024); that paper supplies the pair's Grueneisen parameters, and the pair's historical use is Michail et al. (2016), as the docstring says.

The full history is in CHANGELOG.md.

Limits, and what is deliberately not included

The model is linear: each shift is a fixed lever arm times each cause. Where a material responds nonlinearly (graphene doping, the W-based doping behaviour above), no linear density coefficient is shipped.
No constants beyond the three cited sets. Response coefficients depend on material, mode pair, laser wavelength and substrate; a measurement tool that ships unverified constants propagates wrong results. For any other system you provide the numbers -- from the literature or from calibrate_lever_arms -- and the mandatory reference field makes their origin travel with the analysis.
The smoothing weights of the Bayesian inversion are user-chosen, because estimating them would need noise assumptions this package does not make. Pixels it fills (fill_masked=True) are interpolations from their neighbours, not measurements.
Peak-fit error bars are linearised estimates; they are reliable when the residuals are dominated by uncorrelated noise.
By default fit_two_modes fits each peak in its own window, and the tail of the other peak can shift a fitted centre (by 0.08 cm^-1 in example 13; the round-trip test allows 5e-3 cm^-1). joint=True removes this, provided no third peak lies anywhere from the lower edge of the first window to the upper edge of the second; it fits Lorentzians only (no joint Voigt fit). The inversion then ignores the correlation between the two fitted centres (reported as center_correlation).
min_snr and fwhm_range in fit_map are off by default and their values are your choice; the package does not decide what counts as "no peak" for your data.
load_map_table groups stage positions by exact value (or after rounding with decimals) and needs one shared wavenumber axis; it does not resample spectra onto a common axis, and it does not read vendor binary formats or decimal commas. design_references plans for the through-origin calibration only; check a set chosen for fit_offset=True with plan_calibration(..., fit_offset=True).
separation_with_calibration and multimode_with_calibration are first-order in the calibration errors. The calibration part is common to all pixels, so it does not average down over a map. The chi-square check of multimode_with_calibration ignores the calibration uncertainty.

Where it comes from

T. M. Mahim and M. M. Rahman, "Two Raman phonons quantify the fixed edge charge left by patterning monolayer transition metal dichalcogenides" (under review). Code for the paper: https://github.com/Tanvir-Mahmud-Mahim/Width-scaling-in-monolayer-semiconductor-nanoribbon-transistors

This package is the general-purpose, material-agnostic tool; the paper repository reproduces the specific published study.

Citing, support and license

If ramansep helps your work, please cite it with the concept DOI 10.5281/zenodo.22014913, which always resolves to the latest version; every release is archived on Zenodo. CITATION.cff has the details.

Written and maintained by Tanvir Mahmud Mahim (Department of Electrical and Electronic Engineering, BRAC University), who reviews every change and takes the final decision on scope and releases. Design questions are discussed in the open in issues and pull requests, and the standing rule of CONTRIBUTING.md binds the maintainer exactly as it binds contributors: a change that touches physics arrives with a test, and a constant arrives with its source.

Support runs through the issue tracker. Usage questions are welcome alongside bug reports; a docstring that left a unit or a sign convention unclear is treated as a documentation bug, not user error. While the version is below 1.0 the API may still move between minor versions; such changes are called out in the release notes.

Licensed under Apache-2.0.

Name		Name	Last commit message	Last commit date
.github		.github
examples		examples
src/ramansep		src/ramansep
tests		tests
.gitignore		.gitignore
CHANGELOG.md		CHANGELOG.md
CITATION.cff		CITATION.cff
CONTRIBUTING.md		CONTRIBUTING.md
LICENSE		LICENSE
README.md		README.md
pyproject.toml		pyproject.toml

Navigation Menu

Repository navigation

About

Repository files navigation

ramansep

Contents

A short guide to the words used here

Install, requirements and units

Examples

1. Strain and charge at two spots of a MoS2 ribbon

2. Choosing a mode pair

3. From a map file to strain and charge maps

4. Three peaks, and a check that the model holds

5. Smoothing a whole map

6. Measuring the lever arms on your own instrument

7. Planning the calibration before measuring it

8. How much of the error bar comes from the calibration

9. Temperature from the anti-Stokes / Stokes ratio

10. Strain, charge and temperature together

11. Graphene: the G + 2D decomposition

12. Pixels with no peak at all

13. Two peaks that overlap

14. Smoothing a map that has gaps

15. Calibration uncertainty with three or more modes

16. Reading a map exported by the instrument

17. Calibrating when the starting state is unknown

What is in the package

Cited coefficient sets

When it refuses, and why

How the results are checked

Corrections in earlier versions

Limits, and what is deliberately not included

Where it comes from

Citing, support and license

About

Releases

Sponsor this project

Packages

Contributors

Languages

Latest commit

History

Folders and files

Repository files navigation

ramansep

Contents

A short guide to the words used here

Install, requirements and units

Examples

1. Strain and charge at two spots of a MoS2 ribbon

2. Choosing a mode pair

3. From a map file to strain and charge maps

4. Three peaks, and a check that the model holds

5. Smoothing a whole map

6. Measuring the lever arms on your own instrument

7. Planning the calibration before measuring it

8. How much of the error bar comes from the calibration

9. Temperature from the anti-Stokes / Stokes ratio

10. Strain, charge and temperature together

11. Graphene: the G + 2D decomposition

12. Pixels with no peak at all

13. Two peaks that overlap

14. Smoothing a map that has gaps

15. Calibration uncertainty with three or more modes

16. Reading a map exported by the instrument

17. Calibrating when the starting state is unknown

What is in the package

Cited coefficient sets

When it refuses, and why

How the results are checked

Corrections in earlier versions

Limits, and what is deliberately not included

Where it comes from

Citing, support and license

About

Topics

Resources

Contributing

Stars

Watchers

Forks

Releases

Sponsor this project

Packages

Contributors

Languages