Skip to content

Commit

Permalink
Adjust bounds calculation for 4-peak-fit
Browse files Browse the repository at this point in the history
  • Loading branch information
@raimund-schluessler
raimund-schluessler committed Sep 28, 2022
1 parent a168083 commit 41460dd
Show file tree
Hide file tree
Showing 2 changed files with 112 additions and 11 deletions.
48 changes: 37 additions & 11 deletions bmlab/controllers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -703,24 +703,49 @@ def create_bounds(self, brillouin_regions, rayleigh_peaks):
if bounds is None:
return None

# We need a Rayleigh peak position for
# every brillouin_region
if rayleigh_peaks.shape[0] != (len(brillouin_regions)):
# We need two Rayleigh peak positions to determine the FSR
# and decide whether a peak is Stokes or Anti-Stokes
if rayleigh_peaks.shape[0] != 2:
return None

w0_bounds = []
# We have to create a separate bound for every region
for region_idx, region in enumerate(brillouin_regions):
local_time = []
for f_rayleigh in rayleigh_peaks[region_idx, :]:
# In case this is an Anti-Stokes peak, we find the peak
# with the higher frequency on the left hand side and
# have to flip the bounds
is_anti_stokes = np.mean(region) <\
f_rayleigh
for rayleigh_idx in range(rayleigh_peaks.shape[1]):
anti_stokes_limit = np.nanmean(rayleigh_peaks[:, rayleigh_idx])
# We need to figure out, whether a peak is a
# Stokes or Anti-Stokes peak in order to decide
# to which Rayleigh peak we need to relate to.
tmp = region >= anti_stokes_limit
is_pure_region = (tmp == tmp[0]).all()

is_anti_stokes_region = np.mean(region) >\
anti_stokes_limit

local_bound = []
for bound in bounds:
parsed_bound = []
for limit in bound:
try:
parsed_bound.append(float(limit))
except ValueError:
parsed_bound.append(np.NaN)
# We don't treat Inf as a value
with warnings.catch_warnings():
warnings.filterwarnings(
action='ignore',
message='Mean of empty slice'
)
is_anti_stokes_peak =\
np.nanmean(
np.array(parsed_bound)[
np.isfinite(parsed_bound)]
) < 0

is_anti_stokes = is_anti_stokes_region if\
is_pure_region else is_anti_stokes_peak

local_limit = []
for limit in bound:
if limit.lower() == 'min':
Expand All @@ -739,10 +764,11 @@ def create_bounds(self, brillouin_regions, rayleigh_peaks):
# a value in pixel depending on the time
try:
val = ((-1) ** is_anti_stokes)\
* 1e9 * float(limit) + f_rayleigh
* 1e9 * abs(float(limit))\
+ rayleigh_peaks[
int(is_anti_stokes), rayleigh_idx]
except BaseException:
val = np.Inf
# print(val)
local_limit.append(val)
# Check that the bounds are sorted ascendingly
# (for anti-stokes, they might not).
Expand Down
75 changes: 75 additions & 0 deletions tests/test_evaluation_controller.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -635,6 +635,8 @@ def test_create_bounds(mocker):
fit_bounds = EvaluationController().create_bounds(None, None)
assert fit_bounds is None

# Test that we get correct bounds for regions
# that only contain one type of peaks (either Stokes or Anti-Stokes)
evm.bounds = [['min', '5'], ['5.5', 'Inf'], ['-inf', 'max']]

brillouin_regions = [(3.3e9, 7.0e9), (8.0e9, 12.1e9)]
Expand All @@ -652,3 +654,76 @@ def test_create_bounds(mocker):
[[10.3e9, 12.1e9], [-np.inf, 9.8e9], [8.0e9, np.inf]]
]
], fit_bounds, atol=1e6)

# Test that we get the correct bounds for regions
# that contain both Stokes and Anti-Stokes peaks
evm.bounds = [['min', '5'], ['5.5', 'Inf'], ['-inf', 'max'],
['min', '-5'], ['-5.5', 'Inf'], ['-inf', 'max']]

brillouin_regions = [(2.3e9, 12.1e9), (3.3e9, 12.1e9)]
rayleigh_peaks = 1e9 * np.array([[0.0, 0.1, -0.2], [14.8, 15.0, 15.3]])
fit_bounds = evc.create_bounds(brillouin_regions, rayleigh_peaks)

np.testing.assert_allclose([
[
[[2.3e9, 5.0e9], [5.5e9, np.inf], [-np.inf, 12.1e9],
[9.8e9, 12.1e9], [-np.inf, 9.3e9], [-np.inf, 12.1e9]],
[[2.3e9, 5.1e9], [5.6e9, np.inf], [-np.inf, 12.1e9],
[10.0e9, 12.1e9], [-np.inf, 9.5e9], [-np.inf, 12.1e9]],
[[2.3e9, 4.8e9], [5.3e9, np.inf], [-np.inf, 12.1e9],
[10.3e9, 12.1e9], [-np.inf, 9.8e9], [-np.inf, 12.1e9]]
],
[
[[3.3e9, 5.0e9], [5.5e9, np.inf], [-np.inf, 12.1e9],
[9.8e9, 12.1e9], [-np.inf, 9.3e9], [-np.inf, 12.1e9]],
[[3.3e9, 5.1e9], [5.6e9, np.inf], [-np.inf, 12.1e9],
[10.0e9, 12.1e9], [-np.inf, 9.5e9], [-np.inf, 12.1e9]],
[[3.3e9, 4.8e9], [5.3e9, np.inf], [-np.inf, 12.1e9],
[10.3e9, 12.1e9], [-np.inf, 9.8e9], [-np.inf, 12.1e9]]
]
], fit_bounds, atol=1e6)

# Test real values for a two peak fit
evm.bounds = [['3.5', '5'], ['5.0', '7.5']]

brillouin_regions = [(3.3e9, 7.5e9), (8.0e9, 12.1e9)]
rayleigh_peaks = 1e9 * np.array([[0.0, 0.1, -0.2], [14.8, 15.0, 15.3]])
fit_bounds = evc.create_bounds(brillouin_regions, rayleigh_peaks)

np.testing.assert_allclose([
[
[[3.5e9, 5.0e9], [5.0e9, 7.5e9]],
[[3.6e9, 5.1e9], [5.1e9, 7.6e9]],
[[3.3e9, 4.8e9], [4.8e9, 7.3e9]]
], [
[[9.8e9, 11.3e9], [7.3e9, 9.8e9]],
[[10.0e9, 11.5e9], [7.5e9, 10.0e9]],
[[10.3e9, 11.8e9], [7.8e9, 10.3e9]]
]
], fit_bounds, atol=1e6)

# Test real values for a four peak fit over the whole spectrum
evm.bounds = [['3.5', '5'], ['5.0', '7.5'],
['-7.5', '-5.0'], ['-5', '-3.5']]

brillouin_regions = [(2.3e9, 12.1e9), (3.3e9, 12.1e9)]
rayleigh_peaks = 1e9 * np.array([[0.0, 0.1, -0.2], [14.8, 15.0, 15.3]])
fit_bounds = evc.create_bounds(brillouin_regions, rayleigh_peaks)

np.testing.assert_allclose([
[
[[3.5e9, 5.0e9], [5.0e9, 7.5e9],
[7.3e9, 9.8e9], [9.8e9, 11.3e9]],
[[3.6e9, 5.1e9], [5.1e9, 7.6e9],
[7.5e9, 10.0e9], [10.0e9, 11.5e9]],
[[3.3e9, 4.8e9], [4.8e9, 7.3e9],
[7.8e9, 10.3e9], [10.3e9, 11.8e9]]
], [
[[3.5e9, 5.0e9], [5.0e9, 7.5e9],
[7.3e9, 9.8e9], [9.8e9, 11.3e9]],
[[3.6e9, 5.1e9], [5.1e9, 7.6e9],
[7.5e9, 10.0e9], [10.0e9, 11.5e9]],
[[3.3e9, 4.8e9], [4.8e9, 7.3e9],
[7.8e9, 10.3e9], [10.3e9, 11.8e9]]
]
], fit_bounds, atol=1e6)

0 comments on commit 41460dd

Please sign in to comment.