Update on refactoring

This commit is contained in:
linarphy 2024-06-12 15:17:16 +02:00
parent 2dea7090e3
commit 2283a9a778
No known key found for this signature in database
GPG key ID: E61920135EFF2295
5 changed files with 192 additions and 134 deletions

View file

@ -1,17 +1,20 @@
from pathlib import Path from pathlib import Path
from tomllib import load from tomllib import load
from numpy.lib.npyio import loadtxt from numpy.lib.npyio import loadtxt
from science_signal import to_signal from science_signal import to_signal # type: ignore[reportMissingTypeStubs]
from science_signal.signal import Signal from science_signal.signal import Signal # type: ignore[reportMissingTypeStubs]
from scipy.io.matlab import loadmat from scipy.io.matlab import loadmat # type: ignore[reportMissingTypeStubs]
from backscattering_analyzer import logger from backscattering_analyzer import logger
from backscattering_analyzer.acquisition import Acquisition, AcquisitionType from backscattering_analyzer.acquisition import (
Acquisition,
AcquisitionType,
)
from backscattering_analyzer.movements import Movements from backscattering_analyzer.movements import Movements
class Data: class Data:
def __init__(self, folder: Path, bench: str, date: str): def __init__(self, folder: Path, bench: str, date: str):
logger.info( logger.debug(
_("loading data for {bench} in {date}").format( _("loading data for {bench} in {date}").format(
bench=bench, date=date bench=bench, date=date
) )
@ -27,7 +30,7 @@ class Data:
try: try:
with open(files["metadata"], "rb") as file: with open(files["metadata"], "rb") as file:
metadata = load(file) metadata = load(file)
logger.info( logger.debug(
_("successfully load {metadata}").format( _("successfully load {metadata}").format(
metadata=files["metadata"], metadata=files["metadata"],
) )
@ -136,53 +139,52 @@ class Data:
"reference" "reference"
] / "{channel}.csv".format(channel=channel_names["sensibility"]) ] / "{channel}.csv".format(channel=channel_names["sensibility"])
self.reference = Acquisition(type = AcquisitionType(0)) time: float
duration: float
try: try:
self.reference.movements = Movements( movements = Movements(
bench=to_signal(loadtxt(files["ref_bench"]).T), bench=to_signal(loadtxt(files["ref_bench"]).T),
mirror=to_signal(loadtxt(files["ref_mirror"]).T), mirror=to_signal(loadtxt(files["ref_mirror"]).T),
) )
self.reference.sensibility = to_signal( sensibility = to_signal(loadtxt(files["ref_sensibility"]).T)
loadtxt(files["ref_sensibility"]).T time = sensibility.x[0]
duration = sensibility.x[-1] - sensibility.x[0]
logger.debug(
_(
"successfully loaded reference data files: "
+ "\n {bench}"
+ "\n {mirror}"
+ "\n {sensibility}"
).format(
bench=files["ref_bench"],
mirror=files["ref_mirror"],
sensibility=files["ref_sensibility"],
)
) )
self.reference.time = self.reference.sensibility[0][0]
self.reference.duration = (
self.reference.sensibility[0][-1]
- self.reference.sensibility[0][0]
)
logger.debug(_("successfully loaded reference data"))
try: try:
self.reference.time = metadata_reference["time"] time = metadata_reference["time"]
self.reference.duration = metadata_reference["duration"] duration = metadata_reference["duration"]
if ( if sensibility.x[0] != time:
self.reference.sensibility[0][0]
!= self.reference.time
):
logger.warning( logger.warning(
_( _(
"time between metadata and data does not match: " "time between metadata and data does not match: "
+ "{data} != {metadata}" + "{data} != {metadata}"
).format( ).format(
data=self.reference.sensibility[0][0], data=sensibility.x[0],
metadata=self.reference.time, metadata=time,
) )
) )
if ( if abs(
self.reference.sensibility[0][-1] sensibility.x[-1] - sensibility.x[0] - duration
- self.reference.sensibility[0][0] ) > 1e-2:
!= self.reference.duration
):
logger.warning( logger.warning(
_( _(
"duration between metadata and data does not " "duration between metadata and data does not "
+ "match: {data} != {metadata}" + "match: {data} != {metadata}"
).format( ).format(
data=( data=(sensibility.x[-1] - sensibility.x[0]),
self.reference.sensibility[0][-1] metadata=duration,
- self.reference.sensibility[0][0]
),
metadata=self.reference.duration,
) )
) )
except KeyError: except KeyError:
@ -202,9 +204,15 @@ class Data:
) )
raise e raise e
logger.debug(_("reference data loaded and checked")) self.reference = Acquisition(
time=time,
duration=duration,
movements=movements,
sensibility=sensibility,
type=AcquisitionType.REFERENCE,
)
self.excited = Acquisition(type = AcquisitionType(1)) logger.debug(_("reference data loaded and checked"))
if excited: if excited:
files["exc_bench"] = folders[ files["exc_bench"] = folders[
@ -219,50 +227,52 @@ class Data:
channel=channel_names["sensibility"] channel=channel_names["sensibility"]
) )
try: try:
self.excited.movements = Movements( movements = Movements(
bench=to_signal(loadtxt(files["exc_bench"]).T), bench=to_signal(loadtxt(files["exc_bench"]).T),
mirror=to_signal(loadtxt(files["exc_mirror"]).T), mirror=to_signal(loadtxt(files["exc_mirror"]).T),
) )
self.excited_sensibility = to_signal( sensibility = to_signal(
loadtxt(files["exc_sensibility"]).T loadtxt(files["exc_sensibility"]).T
) )
self.excited.time = self.excited.sensibility[0][0] time = sensibility.x[0]
self.excited.duration = ( duration = sensibility.x[-1] - sensibility.x[0]
self.excited.sensibility[0][-1] logger.debug(
- self.excited.sensibility[0][0] _(
"successfully loaded excited data files: "
+ "\n {bench}"
+ "\n {mirror}"
+ "\n {sensibility}"
).format(
bench=files["exc_bench"],
mirror=files["exc_mirror"],
sensibility=files["exc_sensibility"],
)
) )
logger.debug(_("successfully loadd excited data"))
try: try:
self.excited.time = metadata_excited["time"] time = metadata_excited["time"]
self.excited.duration = metadata_excited["duration"] duration = metadata_excited["duration"]
if ( if sensibility.x[0] != time:
self.excited.sensibility[0][0]
!= self.excited.time
):
logger.warning( logger.warning(
_( _(
"time between metadata and data does not " "time between metadata and data does not "
+ "match: {data} != {metadata}" + "match: {data} != {metadata}"
).format( ).format(
data=self.excited.sensibility[0][0], data=sensibility.x[0],
metadata=self.excited.time, metadata=time,
) )
) )
if ( if abs(
self.excited.sensibility[0][-1] sensibility.x[-1] - sensibility.x[0] - duration
- self.excited.sensibility[0][0] ) > 1e-2:
!= self.excited.duration
):
logger.warning( logger.warning(
_( _(
"duration between metadata and data does " "duration between metadata and data does "
+ "not match: {data} != {metadata}" + "not match: {data} != {metadata}"
).format( ).format(
data=( data=(
self.excited.sensibility[0][-1] sensibility.x[-1] - sensibility.x[0]
- self.excited.sensibility[0][0]
), ),
metadata=self.excited.duration, metadata=duration,
) )
) )
except KeyError: except KeyError:
@ -281,6 +291,13 @@ class Data:
) )
) )
) )
self.excited = Acquisition(
time=time,
duration=duration,
movements=movements,
sensibility=sensibility,
type=AcquisitionType.EXCITED,
)
logger.debug(_("excited data loaded and checked")) logger.debug(_("excited data loaded and checked"))
@ -313,7 +330,7 @@ def load_coupling(
""" """
Load and correct coupling date from modelisation Load and correct coupling date from modelisation
""" """
logger.info( logger.debug(
_("loading coupling for {bench} in {date}").format( _("loading coupling for {bench} in {date}").format(
bench=bench, date=date bench=bench, date=date
) )
@ -321,7 +338,7 @@ def load_coupling(
modelisation_file = folder / file modelisation_file = folder / file
try: try:
modelisation_data = loadmat(modelisation_file) modelisation_data = loadmat(modelisation_file)
logger.info( logger.debug(
_("file successfully loaded:\n{file}").format( _("file successfully loaded:\n{file}").format(
file=modelisation_file file=modelisation_file
) )
@ -378,7 +395,7 @@ def load_coupling(
/ model_powers["laser"] / model_powers["laser"]
) # radiation ) # radiation
elif bench in ["SIB2", "SPRB"]: elif bench in ["SIB2", "SPRB"]:
logger.info( logger.warning(
_( _(
"cannot fix projection power for SIB2 and SPRB coupling" "cannot fix projection power for SIB2 and SPRB coupling"
) )

View file

@ -16,44 +16,38 @@ def show_projection(
show, show,
legend, legend,
close, close,
vlines,
title, title,
xlabel, xlabel,
ylabel, ylabel,
) )
excited = experiment.signals["excited"].psd().sqrt() excited = experiment.data.excited.sensibility.psd().sqrt()
reference = experiment.signals["reference"].psd().sqrt() reference = experiment.data.reference.sensibility.psd().sqrt()
loglog( loglog(
experiment.projection.x, experiment.projection_excited.x,
experiment.projection.y, experiment.projection_excited.y,
label="projection", label=_("projection with excitation"),
) # type: ignore[ReportUnusedCallResult] ) # type: ignore[ReportUnusedCallResult]
loglog(excited.x, excited.y, label=_("excited bench")) # type: ignore[ReportUnusedCallResult]
loglog(reference.x, reference.y, label=_("reference bench")) # type: ignore[ReportUnusedCallResult]
loglog( loglog(
(experiment.projection + reference).x, experiment.projection_reference.x,
(experiment.projection + reference).y, experiment.projection_reference.y,
label=_("sum reference + excited"), label=_("projection without excitation"),
) # type: ignore[ReportUnunsedCallResult]
loglog(
excited.x,
excited.y,
label=_("detected sensibility with excitation"),
) # type: ignore[ReportUnusedCallResult]
loglog(
reference.x,
reference.y,
label=_("detected sensibility without excitation"),
) # type: ignore[ReportUnusedCallResult]
loglog(
(experiment.projection_excited + reference).x,
(experiment.projection_excited + reference).y,
label=_("sum: no excitation and projection with excitation"),
) # type: ignore[ReportUnusedCallResult] ) # type: ignore[ReportUnusedCallResult]
if start_frequency is not None:
vlines(
[
start_frequency,
],
min(experiment.projection.y),
max(experiment.projection.y),
color="k",
) # type: ignore[ReportUnusedCallResult]
if end_frequency is not None:
vlines(
[
end_frequency,
],
min(experiment.projection.y),
max(experiment.projection.y),
color="k",
) # type: ignore[ReportUnusedCallResult]
legend() # type: ignore[ReportUnusedCallResult] legend() # type: ignore[ReportUnusedCallResult]
title(experiment.bench) #  type: ignore[ReportUnusedCallResult] title(experiment.bench) #  type: ignore[ReportUnusedCallResult]
xlabel(_("frequency (Hz)")) # type: ignore[ReportUnusedCallResult] xlabel(_("frequency (Hz)")) # type: ignore[ReportUnusedCallResult]

View file

@ -7,6 +7,7 @@ from numpy.core.multiarray import array
from science_signal import interpolate from science_signal import interpolate
from science_signal.signal import Signal from science_signal.signal import Signal
from backscattering_analyzer import logger from backscattering_analyzer import logger
from backscattering_analyzer.acquisition import AcquisitionType
from backscattering_analyzer.data import ( from backscattering_analyzer.data import (
Data, Data,
load_coupling, load_coupling,
@ -40,7 +41,7 @@ class Experiment:
try: try:
with open(values_file, "rb") as file: with open(values_file, "rb") as file:
values = load(file) values = load(file)
logger.info( logger.debug(
_("successfully loaded values from {file}").format( _("successfully loaded values from {file}").format(
file=values_file file=values_file
) )
@ -254,25 +255,74 @@ class Experiment:
) )
self.factors = self.fit_factors() self.factors = self.fit_factors()
self.projection = self.compute_projection() self.projection_excited = self.compute_projection(
AcquisitionType.EXCITED,
)
self.projection_reference = self.compute_projection(
AcquisitionType.REFERENCE,
)
logger.debug(_("end of experiment initialisation")) logger.debug(_("end of experiment initialisation"))
def get_factors( def get_factors(
self, self,
phase: float,
*signals: Signal,
start: float | None = None, start: float | None = None,
end: float | None = None, end: float | None = None,
) -> tuple[Signal, Signal, Signal, Signal, Signal, Signal, float]: type: AcquisitionType = AcquisitionType.EXCITED,
) -> tuple[Signal, ...]:
""" """
Get factor to compute the projection for a given scatter factor Get factor to compute the projection for a given scatter factor
""" """
phase = 4 * pi / self.wavelength
factor_n = (self.excited_movement * phase).sin().psd().sqrt()
coupling_n = self.coupling[0] coupling_n = self.coupling[0]
factor_d = (
(self.excited_movement["true"] * phase).cos().psd().sqrt()
)
coupling_d = self.coupling[1] coupling_d = self.coupling[1]
if type is AcquisitionType.EXCITED:
if self.excited_movement is not None:
factor_n = (
(self.excited_movement * phase).sin().psd().sqrt()
)
factor_d = (
(self.excited_movement * phase).cos().psd().sqrt()
)
else:
logger.critical(
_(
"no excited signal given, cannot get factors "
+ "of the excited signal"
)
)
raise Exception(
_(
"no excited signal given, cannot get facotrs "
+ "of the excited signal"
)
)
elif type is AcquisitionType.REFERENCE:
if self.reference_movement is not None:
factor_n = (
(self.reference_movement * phase).sin().psd().sqrt()
)
factor_d = (
(self.reference_movement * phase).cos().psd().sqrt()
)
else:
logger.critical(
_(
"no reference signal given, cannot get factors "
+ "of the reference signal"
)
)
raise Exception(
_(
"no reference signal given, cannot get factors "
+ "of the reference signal"
)
)
else:
logger.critical(_("unknown acquisition type"))
raise Exception(_("unknown acquisition type"))
if start is None: if start is None:
start = coupling_n.x[0] start = coupling_n.x[0]
if end is None: if end is None:
@ -280,30 +330,12 @@ class Experiment:
coupling_n = coupling_n.cut(start, end) coupling_n = coupling_n.cut(start, end)
( return interpolate(
factor_n, factor_n,
coupling_n, coupling_n,
factor_d, factor_d,
coupling_d, coupling_d,
excited, *signals,
reference,
) = interpolate(
factor_n,
coupling_n,
factor_d,
coupling_d,
self.signals["excited"].psd().sqrt(),
self.signals["reference"].psd().sqrt(),
)
return (
factor_n,
coupling_n,
factor_d,
coupling_d,
excited,
reference,
phase,
) )
def fit_factors( def fit_factors(
@ -318,6 +350,7 @@ class Experiment:
Find the best factor (first order only) to get the projection Find the best factor (first order only) to get the projection
that best match the excited signal that best match the excited signal
""" """
phase = 4 * pi / self.wavelength
( (
factor_n, factor_n,
coupling_n, coupling_n,
@ -325,8 +358,14 @@ class Experiment:
coupling_d, coupling_d,
excited, excited,
reference, reference,
) = self.get_factors(
phase, phase,
) = self.get_factors(start=start_frequency, end=end_frequency) self.data.excited.sensibility.psd().sqrt(),
self.data.reference.sensibility.psd().sqrt(),
start=start_frequency,
end=end_frequency,
type=AcquisitionType.EXCITED,
)
for index in range(precision): for index in range(precision):
logger.debug(_("search for a local minimum")) logger.debug(_("search for a local minimum"))
@ -358,17 +397,17 @@ class Experiment:
if argmin(y) == 0: if argmin(y) == 0:
logger.warning(_("smaller than current range allows")) logger.warning(_("smaller than current range allows"))
scatter_max: float = x[1] scatter_max = x[1]
elif argmin(y) == len(x) - 1: elif argmin(y) == len(x) - 1:
logger.warning(_("bigger than current range allows")) logger.warning(_("bigger than current range allows"))
scatter_min: float = x[-2] scatter_min = x[-2]
else: else:
scatter_min: float = x[argmin(y) - 1] scatter_min = x[argmin(y) - 1]
scatter_max: float = x[argmin(y) + 1] scatter_max = x[argmin(y) + 1]
logger.debug(_("local minimum found")) logger.debug(_("local minimum found"))
pre_scatter_factor: float = scatter_min / 2 + scatter_max / 2 pre_scatter_factor = scatter_min / 2 + scatter_max / 2
logger.info( logger.info(
_("found a scattering factor of {factor}").format( _("found a scattering factor of {factor}").format(
@ -381,19 +420,19 @@ class Experiment:
"true": pre_scatter_factor * 1e-6, "true": pre_scatter_factor * 1e-6,
} }
def compute_projection(self) -> Signal: def compute_projection(
self, type: AcquisitionType = AcquisitionType.EXCITED
) -> Signal:
""" """
Compute the projection of the noise from current values Compute the projection of the noise from current values
""" """
phase = 4 * pi / self.wavelength
( (
factor_n, factor_n,
coupling_n, coupling_n,
factor_d, factor_d,
coupling_d, coupling_d,
_, ) = self.get_factors(phase, type=type)
_,
phase,
) = self.get_factors()
return Signal( return Signal(
factor_n.x, factor_n.x,

View file

@ -8,5 +8,9 @@ class Movements:
Container for the movement of the bench and the mirror Container for the movement of the bench and the mirror
""" """
bench: str bench: Signal
mirror: str mirror: Signal
def __init__(self, bench: Signal, mirror: Signal):
self.bench = bench * 1e-6
self.mirror = mirror * 1e-6

View file

@ -14,7 +14,7 @@ def process_movement(
""" """
Clean and compute relative movement between the bench and the mirror Clean and compute relative movement between the bench and the mirror
""" """
logger.info(_("computing relative movement")) logger.debug(_("computing relative movement"))
return ( return (
( (
@ -28,7 +28,7 @@ def process_movement(
def compute_light( def compute_light(
pre_scatter_factor: float, pre_scatter_factor: float | None,
factor_n: Signal, factor_n: Signal,
coupling_n: Signal, coupling_n: Signal,
factor_d: Signal, factor_d: Signal,
@ -38,6 +38,10 @@ def compute_light(
""" """
Optimized computing of light with pre-computed factor Optimized computing of light with pre-computed factor
""" """
if pre_scatter_factor is None:
raise Exception(_(
"Cannot compute projection without backscattering factor"
))
return ( return (
sqrt(pre_scatter_factor) sqrt(pre_scatter_factor)
/ phase / phase