{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c9023fe-f86d-428a-8072-4264ff672025",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "from rich.console import Console\n",
    "from rich.table import Table\n",
    "\n",
    "from numpy import loadtxt\n",
    "\n",
    "from numpy import float32, bool as np_bool\n",
    "from numpy.typing import NDArray\n",
    "from typing import Any\n",
    "from collections.abc import Callable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d82f94d-b69a-4ad0-b042-08d045925545",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import mean, diff, sqrt, pi, linspace\n",
    "from numpy.random import default_rng\n",
    "\n",
    "from backscattering_analyzer.experiment import Experiment\n",
    "from backscattering_analyzer.acquisition import AcquisitionType\n",
    "\n",
    "from matplotlib.pyplot import figure, show\n",
    "\n",
    "from science_signal import Signal\n",
    "from science_signal.generator import sin as sin_gen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d91d36b8-28fb-40d8-a4e6-77b659839fbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "from labellines import labelLines  # from matplotlib-label-lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62328640-6c4f-47a1-8461-9637a50f52e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_over(signal1: Signal, signal2: Signal) -> bool:\n",
    "    \"\"\"\n",
    "    True if signal1 is over signal2\n",
    "    \"\"\"\n",
    "    return not ((signal1 - signal2).y < 0).any()\n",
    "\n",
    "\n",
    "def get_speed(signal: Signal) -> NDArray[float32]:\n",
    "    \"\"\"\n",
    "    get speed from position signal\n",
    "    \"\"\"\n",
    "    return 1 / mean(diff(signal.x)) * diff(signal.y)\n",
    "\n",
    "\n",
    "def compute_rms(signal: NDArray[float32]) -> float:\n",
    "    \"\"\"\n",
    "    compute RMS of data\n",
    "    \"\"\"\n",
    "    return sqrt(mean(signal**2))  # pyright: ignore[reportAny]\n",
    "\n",
    "\n",
    "def accelerate_timesignal(factor: float, signal: Signal) -> Signal:\n",
    "    \"\"\"\n",
    "    accelerate a signal by a factor\n",
    "    \"\"\"\n",
    "    return Signal(\n",
    "        signal.x[0] + (signal.x - signal.x[0]) * (1 / factor),\n",
    "        signal.y,\n",
    "    )\n",
    "\n",
    "\n",
    "def generate_movement(max_speed: float, start: int, end: int) -> Signal:\n",
    "    \"\"\"\n",
    "    generate a sine movement with a given maximum speed\n",
    "    \"\"\"\n",
    "    amplitude: float = sqrt(max_speed / (2 * pi))\n",
    "    frequency: float = sqrt(max_speed / (2 * pi))\n",
    "    signal = sin_gen(end - start, 1000, frequency, amplitude)\n",
    "    return Signal(\n",
    "        start + signal.x,\n",
    "        signal.y,\n",
    "    )\n",
    "\n",
    "\n",
    "def fit_value(\n",
    "    start: float,\n",
    "    end: float,\n",
    "    nb_loop: int,\n",
    "    experiment: Experiment,\n",
    "    reference: Signal,\n",
    "    generate: Callable[[float, Any], Signal],\n",
    "    condition: Callable[[Signal, Signal], np_bool | bool],\n",
    "    generate_args: list[Any],\n",
    ") -> float:\n",
    "    \"\"\"\n",
    "    find the value that is the closest to NOT respecting the condition function\n",
    "    Considering it start by not respecting the condition, and after some time, it will\n",
    "    \"\"\"\n",
    "    for _ in range(nb_loop):\n",
    "        values = linspace(start, end, 10, dtype=float32)\n",
    "        detected = False\n",
    "        for j in range(len(values)):\n",
    "            value = values[j]\n",
    "            experiment.reference_movement = generate(value, *generate_args)\n",
    "            experiment.projection_reference = experiment.compute_projection(\n",
    "                AcquisitionType.REFERENCE\n",
    "            )\n",
    "            if (\n",
    "                not condition(reference, experiment.projection_reference)\n",
    "                and not detected\n",
    "            ):\n",
    "                detected = True\n",
    "                end = value\n",
    "                start = values[j - 1]\n",
    "        if not detected:\n",
    "            raise Exception(\"not in range\")\n",
    "    return start + (end - start) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "587e2e3a-767b-4a6e-80df-188b6152bbb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "console = Console()\n",
    "o5_mat = (\n",
    "    Path(\"/home/demagny/data\")\n",
    "    / \"simulation\"\n",
    "    / \"optickle\"\n",
    "    / \"transfer_function\"\n",
    "    / \"O5.mat\"\n",
    ")\n",
    "sensitivities: dict[str, Signal] = dict()\n",
    "for name in [\"high\", \"low\"]:\n",
    "    data = loadtxt(\n",
    "        Path(\n",
    "            \"/home/demagny/data/sensitivity/O5/23932_O5{}SensASD.txt\".format(\n",
    "                name.capitalize()\n",
    "            )\n",
    "        ),\n",
    "        dtype=float32,\n",
    "    )\n",
    "    sensitivities[name] = Signal(\n",
    "        data[:, 0],\n",
    "        data[:, 1],\n",
    "    )\n",
    "sensitivities[\"high\"].y /= 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7f534cc-0b50-4fc1-ba62-8a36747eb646",
   "metadata": {},
   "outputs": [],
   "source": [
    "C_BENCH = \"SDB1\"\n",
    "C_COUPLING = \"../main_script_virgo/values-coupling.toml\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a72635ea-8395-4783-9be1-8a89e82b7e21",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_experiment = Experiment(C_BENCH, \"2024_06_07\", C_COUPLING, 0.1)\n",
    "temp_experiment = Experiment(C_BENCH, \"2024_06_07\", C_COUPLING, 0.1)\n",
    "base_experiment.factors = {\"pre\": 2e-10 * 1e6, \"true\": 2e-10}\n",
    "temp_experiment.factors = {\"pre\": 2e-10 * 1e6, \"true\": 2e-10}\n",
    "base_experiment.modelisation_file = str(o5_mat)\n",
    "temp_experiment.modelisation_file = str(o5_mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "850dd2af-2b7b-43b3-8758-6b36d9bced20",
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = default_rng().random(1) * 10\n",
    "# check if frequency is multiplied by the factor after going into accelerate_timesignal function\n",
    "assert (\n",
    "    abs(\n",
    "        factor / mean(diff(base_experiment.reference_movement.x))\n",
    "        - 1\n",
    "        / mean(\n",
    "            diff(accelerate_timesignal(factor, base_experiment.reference_movement).x)\n",
    "        )\n",
    "    )\n",
    "    < 1e-1\n",
    "), abs(\n",
    "    factor / mean(diff(base_experiment.reference_movement.x))\n",
    "    - 1\n",
    "    / mean(diff(accelerate_timesignal(factor, base_experiment.reference_movement).x))\n",
    ")\n",
    "# check if amplitude is not modified after going into accelerate_timesignal function\n",
    "assert (\n",
    "    base_experiment.reference_movement.y\n",
    "    == accelerate_timesignal(factor, base_experiment.reference_movement).y\n",
    ").all(), (\n",
    "    base_experiment.reference_movement.y\n",
    "    - accelerate_timesignal(factor, base_experiment.reference_movement).y\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dbb68d6-fba3-4488-b951-b5bf860152d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_fitted_value = fit_value(\n",
    "    1e-6,\n",
    "    1e-5,\n",
    "    5,\n",
    "    temp_experiment,\n",
    "    sensitivities[\"high\"],\n",
    "    generate_movement,\n",
    "    is_over,\n",
    "    [\n",
    "        base_experiment.reference_movement.x[0],\n",
    "        base_experiment.reference_movement.x[-1],\n",
    "    ],\n",
    ")\n",
    "low_fitted_value = fit_value(\n",
    "    1e-6,\n",
    "    1e-5,\n",
    "    5,\n",
    "    temp_experiment,\n",
    "    sensitivities[\"low\"],\n",
    "    generate_movement,\n",
    "    is_over,\n",
    "    [\n",
    "        base_experiment.reference_movement.x[0],\n",
    "        base_experiment.reference_movement.x[-1],\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "132a82c5-b69b-402c-b278-c565ccf6e504",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_movement = generate_movement(\n",
    "    high_fitted_value,\n",
    "    base_experiment.reference_movement.x[0],\n",
    "    base_experiment.reference_movement.x[-1],\n",
    ")\n",
    "temp_experiment.reference_movement = high_movement\n",
    "high_projection = temp_experiment.compute_projection(AcquisitionType.REFERENCE)\n",
    "low_movement = generate_movement(\n",
    "    low_fitted_value,\n",
    "    base_experiment.reference_movement.x[0],\n",
    "    base_experiment.reference_movement.x[-1],\n",
    ")\n",
    "temp_experiment.reference_movement = low_movement\n",
    "low_projection = temp_experiment.compute_projection(AcquisitionType.REFERENCE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2966a38e-9020-4cdd-a186-70519afb85de",
   "metadata": {},
   "outputs": [],
   "source": [
    "Figure = figure()\n",
    "_ = Figure.gca().loglog(\n",
    "    high_projection.x, high_projection.y, label=\"projection for high\"\n",
    ")\n",
    "_ = Figure.gca().loglog(\n",
    "    sensitivities[\"high\"].x,\n",
    "    sensitivities[\"high\"].y,\n",
    "    label=\"one order below high sensitivity\",\n",
    ")\n",
    "_ = Figure.gca().loglog(low_projection.x, low_projection.y, label=\"projection for low\")\n",
    "_ = Figure.gca().loglog(\n",
    "    sensitivities[\"low\"].x, sensitivities[\"low\"].y, label=\"low sensitivity\"\n",
    ")\n",
    "_ = Figure.gca().legend()\n",
    "Figure.gca().grid(True, \"both\", \"both\")\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf7aa342-1ffe-4262-80fd-e900acf590c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "table = Table(title=\"Vérification des calculs: haute sensibilité\")\n",
    "table.add_column(\"nom\")\n",
    "table.add_column(\"RMS\")\n",
    "table.add_column(\"vitesse max\")\n",
    "\n",
    "table.add_row(\n",
    "    \"calculé\",\n",
    "    \"{:.2E}\".format(high_fitted_value / sqrt(2)),\n",
    "    \"{:.2E}\".format(high_fitted_value),\n",
    ")\n",
    "table.add_row(\n",
    "    \"mesuré\",\n",
    "    \"{:.2E}\".format(\n",
    "        compute_rms(get_speed(high_movement)),\n",
    "    ),\n",
    "    \"{:.2E}\".format(\n",
    "        max(get_speed(high_movement)),\n",
    "    ),\n",
    ")\n",
    "\n",
    "console.print(table)\n",
    "\n",
    "table = Table(title=\"Vérification des calculs: basse sensibilité\")\n",
    "table.add_column(\"nom\")\n",
    "table.add_column(\"RMS\")\n",
    "table.add_column(\"vitesse max\")\n",
    "\n",
    "table.add_row(\n",
    "    \"calculé\",\n",
    "    \"{:.2E}\".format(low_fitted_value / sqrt(2)),\n",
    "    \"{:.2E}\".format(low_fitted_value),\n",
    ")\n",
    "table.add_row(\n",
    "    \"mesuré\",\n",
    "    \"{:.2E}\".format(\n",
    "        compute_rms(get_speed(low_movement)),\n",
    "    ),\n",
    "    \"{:.2E}\".format(\n",
    "        max(get_speed(low_movement)),\n",
    "    ),\n",
    ")\n",
    "\n",
    "console.print(table)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7ca8aea-80f5-4306-916d-a4cb9240bd1b",
   "metadata": {},
   "source": [
    "# Utilisation du mouvement réel du banc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c61da9e2-0339-4345-b601-395ef00b380c",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_fitted_value = fit_value(\n",
    "    1,\n",
    "    3,\n",
    "    5,\n",
    "    temp_experiment,\n",
    "    sensitivities[\"high\"],\n",
    "    accelerate_timesignal,\n",
    "    is_over,\n",
    "    [\n",
    "        base_experiment.reference_movement,\n",
    "    ],\n",
    ")\n",
    "low_fitted_value = fit_value(\n",
    "    1,\n",
    "    6,\n",
    "    5,\n",
    "    temp_experiment,\n",
    "    sensitivities[\"low\"],\n",
    "    accelerate_timesignal,\n",
    "    is_over,\n",
    "    [\n",
    "        base_experiment.reference_movement,\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13cbbcf4-4753-4cef-90a3-83477177a102",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_movement = accelerate_timesignal(\n",
    "    high_fitted_value, base_experiment.reference_movement\n",
    ")\n",
    "low_movement = accelerate_timesignal(\n",
    "    low_fitted_value, base_experiment.reference_movement\n",
    ")\n",
    "temp_experiment.reference_movement = high_movement\n",
    "high_projection = temp_experiment.compute_projection(AcquisitionType.REFERENCE)\n",
    "temp_experiment.reference_movement = low_movement\n",
    "low_projection = temp_experiment.compute_projection(AcquisitionType.REFERENCE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12b9b3ee-a882-4088-9705-5c8a4c38834c",
   "metadata": {},
   "outputs": [],
   "source": [
    "table = Table(title=\"Vérification des calculs: haute sensibilité\")\n",
    "table.add_column(\"nom\")\n",
    "table.add_column(\"RMS\")\n",
    "table.add_column(\"vitesse max\")\n",
    "\n",
    "table.add_row(\n",
    "    \"calculé\",\n",
    "    \"{:.2E}\".format(\n",
    "        compute_rms(get_speed(base_experiment.reference_movement)) * high_fitted_value\n",
    "    ),\n",
    "    \"{:.2E}\".format(\n",
    "        max(get_speed(base_experiment.reference_movement)) * high_fitted_value\n",
    "    ),\n",
    ")\n",
    "table.add_row(\n",
    "    \"mesuré\",\n",
    "    \"{:.2E}\".format(\n",
    "        compute_rms(get_speed(high_movement)),\n",
    "    ),\n",
    "    \"{:.2E}\".format(\n",
    "        max(get_speed(high_movement)),\n",
    "    ),\n",
    ")\n",
    "\n",
    "console.print(table)\n",
    "\n",
    "table = Table(title=\"Vérification des calculs: basse sensibilité\")\n",
    "table.add_column(\"nom\")\n",
    "table.add_column(\"RMS\")\n",
    "table.add_column(\"vitesse max\")\n",
    "\n",
    "table.add_row(\n",
    "    \"calculé\",\n",
    "    \"{:.2E}\".format(\n",
    "        compute_rms(get_speed(base_experiment.reference_movement)) * low_fitted_value\n",
    "    ),\n",
    "    \"{:.2E}\".format(\n",
    "        max(get_speed(base_experiment.reference_movement)) * low_fitted_value\n",
    "    ),\n",
    ")\n",
    "table.add_row(\n",
    "    \"mesuré\",\n",
    "    \"{:.2E}\".format(\n",
    "        compute_rms(get_speed(low_movement)),\n",
    "    ),\n",
    "    \"{:.2E}\".format(\n",
    "        max(get_speed(low_movement)),\n",
    "    ),\n",
    ")\n",
    "\n",
    "console.print(table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ff45ff8-8d0f-4f9a-a060-839700121e43",
   "metadata": {},
   "outputs": [],
   "source": [
    "Figure = figure(figsize=(8, 5))\n",
    "_ = Figure.gca().loglog(\n",
    "    high_projection.x,\n",
    "    high_projection.y,\n",
    "    label=\"rms = {:.2E}, max = {:.2E}\".format(\n",
    "        compute_rms(get_speed(high_movement)), max(get_speed(high_movement))\n",
    "    ),\n",
    ")\n",
    "_ = Figure.gca().loglog(\n",
    "    low_projection.x,\n",
    "    low_projection.y,\n",
    "    label=\"rms= {:.2E} max= {:.2E}\".format(\n",
    "        compute_rms(get_speed(low_movement)), max(get_speed(low_movement))\n",
    "    ),\n",
    ")\n",
    "_ = Figure.gca().loglog(\n",
    "    sensitivities[\"low\"].x, sensitivities[\"low\"].y, label=\"low sensitivity\"\n",
    ")\n",
    "_ = Figure.gca().loglog(\n",
    "    sensitivities[\"high\"].x,\n",
    "    sensitivities[\"high\"].y,\n",
    "    label=\"one order below high sensitivity\",\n",
    ")\n",
    "\n",
    "labelLines(Figure.gca().get_lines(), xvals=[9, 20, 50, 50], align=False)\n",
    "\n",
    "# _ = Figure.gca().legend(loc=\"upper right\")\n",
    "_ = Figure.gca().set_xlim(5, 100)\n",
    "_ = Figure.gca().set_ylim(1e-26, 1e-18)\n",
    "_ = Figure.gca().set_xlabel(\"Frequencies Hz\")\n",
    "_ = Figure.gca().set_ylabel(\"Sensitivity $\\\\frac { 1 } { \\\\sqrt{ Hz } }$\")\n",
    "_ = Figure.gca().set_title(\"projection of noise from backscatterd light on SRB\")\n",
    "Figure.gca().grid(True, \"both\", \"both\")\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e985e2c-ac0d-4ba7-8527-cadf16fcb22d",
   "metadata": {},
   "outputs": [],
   "source": [
    "console.print(\n",
    "    \"base rms: [repr.number]{:.2E}[/repr.number]\\nbase max speed: [repr.number]{:.2E}[/repr.number]\".format(\n",
    "        compute_rms(get_speed(base_experiment.reference_movement)),\n",
    "        max(get_speed(base_experiment.reference_movement)),\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f447e42-d7a9-4c69-8721-d6c21ba0cf15",
   "metadata": {},
   "outputs": [],
   "source": [
    "console.print(high_fitted_value, low_fitted_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ef341b0-8d63-4130-81f7-82cdad1b0052",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}