344 lines
10 KiB
Python
Executable file
344 lines
10 KiB
Python
Executable file
#!/usr/bin/env python3
|
|
# pyright: reportUnusedCallResult=false
|
|
|
|
from pathlib import Path
|
|
from numpy import loadtxt
|
|
from matplotlib.pyplot import figure, show
|
|
|
|
from backscattering_analyzer.acquisition import AcquisitionType
|
|
from backscattering_analyzer.experiment import Experiment
|
|
|
|
from science_signal import Signal
|
|
|
|
data_folder = Path("/home/demagny/data")
|
|
o5_mat = data_folder / "simulation" / "optickle" / "transfer_function" / "O5.mat"
|
|
|
|
files: dict[str, Path] = {
|
|
"high": data_folder / "sensitivity" / "O5" / "23932_O5HighSensASD.txt",
|
|
"low": data_folder / "sensitivity" / "O5" / "23932_O5LowSensASD.txt",
|
|
"design": data_folder / "sensitivity" / "O5" / "23932_O5DesignASD.txt",
|
|
"planB": data_folder / "sensitivity" / "O5" / "23932_O5PlanBASD.txt",
|
|
}
|
|
sensitivities: dict[str, Signal] = {}
|
|
for name, file in files.items():
|
|
data = loadtxt(file)
|
|
sensitivities[name] = Signal(data[:, 0], data[:, 1])
|
|
data = loadtxt(files["high"])
|
|
sensitivities["one order below high"] = Signal(data[:, 0], data[:, 1] / 10)
|
|
|
|
benches = ["SDB1", "SDB2", "SNEB", "SWEB"]
|
|
# benches = ["SDB1"]
|
|
|
|
experiments: dict[str, list[Experiment]] = {
|
|
"good weather": [
|
|
Experiment(bench, "2024_06_07", "values-coupling.toml", 0.1)
|
|
for bench in benches
|
|
],
|
|
"bad weather": [
|
|
Experiment(bench, "2024_08_15", "values-coupling.toml", 0.1)
|
|
for bench in benches
|
|
],
|
|
}
|
|
|
|
# fscs = {
|
|
# "best": {
|
|
# "SDB1": 5e-15,
|
|
# "SDB2": 5e-16,
|
|
# "SNEB": 7e-15,
|
|
# "SWEB": 5e-15,
|
|
# },
|
|
# "'realistic'": {
|
|
# "SDB1": 2e-10,
|
|
# "SDB2": 2e-12,
|
|
# "SNEB": 2e-10,
|
|
# "SWEB": 2e-10,
|
|
# },
|
|
# "current": {
|
|
# "SDB1": 2e-9,
|
|
# "SDB2": 2e-12,
|
|
# "SNEB": 5e-9,
|
|
# "SWEB": 3e-9,
|
|
# },
|
|
# "wow": {
|
|
# "SDB1": ratio["SDB1"] * best_base_fsc,
|
|
# "SDB2": ratio["SDB2"] * best_base_fsc,
|
|
# "SNEB": ratio["SNEB"] * best_base_fsc,
|
|
# "SWEB": ratio["SWEB"] * best_base_fsc,
|
|
# },
|
|
# }
|
|
|
|
# ratio = {
|
|
# "SDB1": 1,
|
|
# "SDB2": 1e-2,
|
|
# "SNEB": 1,
|
|
# "SWEB": 1,
|
|
# }
|
|
# best_base_fsc = 5e-13
|
|
|
|
ratio = {
|
|
"SDB1": 2e-10,
|
|
"SDB2": 2e-12,
|
|
"SNEB": 2e-10,
|
|
"SWEB": 2e-10,
|
|
}
|
|
best_base_fsc = 1
|
|
|
|
fscs: list[dict[str, float | None]] = []
|
|
for bench in benches:
|
|
fscs.append(
|
|
{
|
|
"pre": ratio[bench] * best_base_fsc * 1e6,
|
|
"true": ratio[bench] * best_base_fsc,
|
|
}
|
|
)
|
|
|
|
Figure = figure()
|
|
results: dict[str, Signal] = dict()
|
|
for name, experiment in experiments.items():
|
|
projection = None
|
|
for index in range(len(benches)):
|
|
experiment[index].factors = fscs[index]
|
|
experiment[index].projection_reference = experiment[index].compute_projection(
|
|
AcquisitionType.REFERENCE
|
|
)
|
|
"""
|
|
if name == "bad weather":
|
|
Figure.gca().loglog(
|
|
experiment[index].projection_reference.x,
|
|
experiment[index].projection_reference.y,
|
|
linestyle="--",
|
|
label="projection of {} for a {}".format(benches[index], name),
|
|
)
|
|
"""
|
|
if projection is None:
|
|
projection = experiment[index].projection_reference
|
|
sensitivity = experiment[index].data.reference.sensibility.psd().sqrt()
|
|
Figure.gca().loglog(
|
|
sensitivity.x,
|
|
sensitivity.y,
|
|
label="{} sensitivity".format(name),
|
|
)
|
|
else:
|
|
projection += experiment[index].projection_reference
|
|
if projection is None:
|
|
raise Exception("wut")
|
|
results[name] = projection
|
|
if name == "bad weather":
|
|
Figure.gca().loglog(
|
|
projection.x,
|
|
projection.y,
|
|
linestyle="--",
|
|
label="projection for a {}".format(name),
|
|
)
|
|
else:
|
|
Figure.gca().loglog(
|
|
projection.x,
|
|
projection.y,
|
|
linestyle="--",
|
|
label="projection for a {}".format(name),
|
|
)
|
|
for name, sensitivity in sensitivities.items():
|
|
Figure.gca().loglog(
|
|
sensitivity.x, sensitivity.y, label="{} sensitivity".format(name)
|
|
)
|
|
Figure.gca().legend()
|
|
Figure.gca().grid(True, "both", "both")
|
|
Figure.gca().set_xlim(10, 40)
|
|
Figure.gca().set_ylim(1e-25, 1e-18)
|
|
Figure.gca().set_xlabel("Frequencies (Hz)")
|
|
Figure.gca().set_title(
|
|
'Sensitivities and projection curves to respect "low sensitivity" curve goal in bad condition'
|
|
)
|
|
show()
|
|
|
|
"""
|
|
from pathlib import Path
|
|
from numpy import argmin, linspace
|
|
from rich.pretty import pprint
|
|
|
|
from backscattering_analyzer.acquisition import Acquisition
|
|
from backscattering_analyzer.projection import Projection
|
|
from matplotlib.pyplot import subplot, show
|
|
|
|
from backscattering_analyzer import Experiment
|
|
from backscattering_analyzer.sensitivity import Sensitivity
|
|
from science_signal import interpolate, Signal
|
|
|
|
data_folder = Path("/home/demagny/data")
|
|
o5_mat = data_folder / "simulation" / "optickle" / "transfer_function" / "O5.mat"
|
|
|
|
experiments = {
|
|
"bonne météo": Experiment(
|
|
"2024_06_07",
|
|
data_folder,
|
|
o5_mat,
|
|
),
|
|
"mauvaise météo": Experiment(
|
|
"2024_08_15",
|
|
data_folder,
|
|
o5_mat,
|
|
),
|
|
# "moins mauvaise météo": Experiment(
|
|
# "2024_06_09",
|
|
# data_folder,
|
|
# o5_mat,
|
|
# ),
|
|
}
|
|
|
|
acquisitions: dict[str, Acquisition] = dict()
|
|
|
|
for name, experiment in experiments.items():
|
|
acquisitions[name] = experiment.get_measure("reference").get_acquisition(
|
|
"reference"
|
|
)
|
|
|
|
sensitivities: dict[str, Sensitivity] = {
|
|
"high": Sensitivity(
|
|
acquisitions["bonne météo"],
|
|
from_txt=Path("/home/demagny/data/sensitivity/O5/23932_O5HighSensASD.txt"),
|
|
),
|
|
"one order below high": Sensitivity(
|
|
acquisitions["bonne météo"],
|
|
from_txt=Path("/home/demagny/data/sensitivity/O5/23932_O5HighSensASD.txt"),
|
|
),
|
|
"low": Sensitivity(
|
|
acquisitions["bonne météo"],
|
|
from_txt=Path("/home/demagny/data/sensitivity/O5/23932_O5LowSensASD.txt"),
|
|
),
|
|
"design": Sensitivity(
|
|
acquisitions["bonne météo"],
|
|
from_txt=Path("/home/demagny/data/sensitivity/O5/23932_O5DesignASD.txt"),
|
|
),
|
|
"plan B": Sensitivity(
|
|
acquisitions["bonne météo"],
|
|
from_txt=Path("/home/demagny/data/sensitivity/O5/23932_O5PlanBASD.txt"),
|
|
),
|
|
}
|
|
if sensitivities["one order below high"].data is None:
|
|
raise Exception("cannot load data for high sensitivity")
|
|
sensitivities["one order below high"].data.y = (
|
|
sensitivities["one order below high"].data.y / 10
|
|
)
|
|
|
|
ratio = {
|
|
"SDB1": 1,
|
|
"SDB2": 1e-2,
|
|
"SNEB": 1,
|
|
"SWEB": 1e-1,
|
|
}
|
|
best_base_fsc = 1.57e-10
|
|
|
|
pprint("début du fit")
|
|
|
|
max_value, min_value = 1e-8, 1e-13
|
|
nb_loop = 5
|
|
nb_range = 10
|
|
for i in range(nb_loop):
|
|
pprint("zoom")
|
|
base_fscs = linspace(max_value, min_value, nb_range)
|
|
differences: list[float] = []
|
|
for base_fsc in base_fscs:
|
|
fscs: dict[str, float] = dict()
|
|
for (
|
|
name,
|
|
value,
|
|
) in ratio.items(): # compute "used" backscatterfactor for every benches
|
|
fscs[name] = base_fsc * value
|
|
acquisitions["mauvaise météo"].projection = Projection(
|
|
acquisitions["mauvaise météo"],
|
|
fscs,
|
|
)
|
|
acquisitions["mauvaise météo"].projection.compute()
|
|
pprint("projection calculée")
|
|
if acquisitions["mauvaise météo"].projection.data is None:
|
|
raise Exception("la projection n'a pas été calculé")
|
|
if sensitivities["low"].data is None:
|
|
raise Exception("la sensibilité n'a pas été chargé")
|
|
low_sensitivity, bad_weather_projection = interpolate(
|
|
sensitivities["low"].data, acquisitions["mauvaise météo"].projection.data
|
|
)
|
|
difference = abs(
|
|
sensitivities["low"].data
|
|
- acquisitions["mauvaise météo"].projection.data.cut(20, 40)
|
|
)
|
|
differences.append(sum(difference.y))
|
|
index_min = argmin(differences)
|
|
if index_min == 0:
|
|
pprint(
|
|
"Attention, il semblerait que le fsc demandé soit plus haut que ce qui est permit (bizarre)"
|
|
)
|
|
min_value = base_fscs[0]
|
|
max_value = base_fscs[1]
|
|
elif index_min == nb_range - 1:
|
|
pprint(
|
|
"Attention, il semblerait que le fsc demandé soit plus bas que ce qui est permis (Aïe Aïe Aïe)"
|
|
)
|
|
min_value = base_fscs[-2]
|
|
max_value = base_fscs[-1]
|
|
else:
|
|
min_value = base_fscs[index_min - 1]
|
|
max_value = base_fscs[index_min + 1]
|
|
|
|
print(min_value, max_value)
|
|
|
|
fscs = {
|
|
# "best": {
|
|
# "SDB1": 5e-15,
|
|
# "SDB2": 5e-16,
|
|
# "SNEB": 7e-15,
|
|
# "SWEB": 5e-15,
|
|
# },
|
|
# "'realistic'": {
|
|
# "SDB1": 1e-10,
|
|
# "SDB2": 1e-12,
|
|
# "SNEB": 1e-10,
|
|
# "SWEB": 1e-11,
|
|
# },
|
|
"current": {
|
|
"SDB1": 2e-9,
|
|
"SDB2": 2e-12,
|
|
"SNEB": 5e-9,
|
|
"SWEB": 3e-9,
|
|
},
|
|
# "wow": {
|
|
# "SDB1": ratio["SDB1"] * best_base_fsc,
|
|
# "SDB2": ratio["SDB2"] * best_base_fsc,
|
|
# "SNEB": ratio["SNEB"] * best_base_fsc,
|
|
# "SWEB": ratio["SWEB"] * best_base_fsc,
|
|
# },
|
|
}
|
|
|
|
results: dict[str, dict[str, Signal]] = dict()
|
|
|
|
i = 1
|
|
for name_fsc, fsc in fscs.items():
|
|
axes = subplot(len(fscs), 1, i)
|
|
results[name_fsc] = dict()
|
|
for name, acquisition in acquisitions.items():
|
|
acquisition.projection = Projection(
|
|
acquisition,
|
|
fsc,
|
|
)
|
|
acquisition.show(axes, True)
|
|
# acquisition.projection.compute()
|
|
# acquisition.sensitivities[0].show(axes, True)
|
|
if acquisition.projection.data is None:
|
|
raise Exception(
|
|
"could not compute {} for the {} fsc".format(name, name_fsc)
|
|
)
|
|
results[name_fsc][name] = acquisition.projection.data
|
|
result = results[name_fsc][name]
|
|
axes.loglog(result.x, result.y, label=name)
|
|
for name, sensitivity in sensitivities.items():
|
|
if sensitivity.data is None:
|
|
raise Exception("could not retrieve sensitivity {}".format(name))
|
|
axes.loglog(sensitivity.data.x, sensitivity.data.y, label=name, linestyle="--")
|
|
axes.set_title("sensitivities for {} fsc".format("current"))
|
|
axes.set_xlabel("frequencies (Hz)")
|
|
axes.set_xlim(10, 60)
|
|
axes.set_ylim(1e-24, 1e-18)
|
|
axes.grid(True, "both", "both")
|
|
axes.legend()
|
|
i += 1
|
|
show()
|
|
"""
|