{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "75b8eca5",
   "metadata": {},
   "source": [
    "# pyOMA – Single-Setup Operational Modal Analysis\n",
    "\n",
    "This notebook walks through a complete Operational Modal Analysis (OMA) workflow using\n",
    "a single measurement setup.  It uses the example data bundled with the repository so\n",
    "you can run it immediately after installation.\n",
    "\n",
    "**Workflow overview**\n",
    "\n",
    "| Step | Class / function | Purpose |\n",
    "|------|------------------|---------|\n",
    "| 1 | `GeometryProcessor` | Load structural geometry (nodes, lines) |\n",
    "| 2 | `PreProcessSignals` | Load measurement data and assign metadata |\n",
    "| 3 | `PreProcessSignals` | Decimate, filter, compute correlations |\n",
    "| 4 | `BRSSICovRef` (or other) | Identify modal parameters |\n",
    "| 5 | `StabilCalc` / `StabilCluster` | Compute stabilisation diagram |\n",
    "| 6 | Results & export | Inspect frequencies, damping, mode shapes |\n",
    "\n",
    "> **Tip:** Run cells top-to-bottom with **Shift+Enter**.  Each section begins with a\n",
    "> short explanation of what it does and why."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3ac620f",
   "metadata": {},
   "source": [
    "## 0  Imports and display setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a825d0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "\n",
    "# Use the interactive widget backend when available, fall back to inline.\n",
    "try:\n",
    "    import ipympl  # noqa: F401 – presence check only\n",
    "    %matplotlib widget\n",
    "except ImportError:\n",
    "    %matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# --- pyOMA imports -----------------------------------------------------------\n",
    "from pyOMA.core import (\n",
    "    GeometryProcessor,\n",
    "    PreProcessSignals,\n",
    "    SignalPlot,\n",
    "    BRSSICovRef,\n",
    "    SSIData,\n",
    "    PLSCF,\n",
    "    VarSSIRef,\n",
    "    StabilCalc,\n",
    "    StabilCluster,\n",
    "    StabilPlot,\n",
    "    ModeShapePlot,\n",
    ")\n",
    "\n",
    "print('pyOMA imported successfully')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d305ffd7",
   "metadata": {},
   "source": [
    "## 1  Paths\n",
    "\n",
    "The cell below locates the **example data** that ships with the repository\n",
    "(`tests/files/`).  When you apply this workflow to your own project, replace\n",
    "`EXAMPLE_DATA` with the folder that contains your measurement files and set\n",
    "`SETUP_DIR` to the folder for the specific setup you want to analyse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "369f791c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyOMA\n",
    "\n",
    "# Root of the repository — works whether pyOMA is installed or run in-place\n",
    "REPO_ROOT    = Path(pyOMA.__file__).parent.parent\n",
    "EXAMPLE_DATA = REPO_ROOT / 'tests' / 'files'\n",
    "SETUP_DIR    = EXAMPLE_DATA / 'measurement_1'\n",
    "\n",
    "assert EXAMPLE_DATA.exists(), f'Example data not found at {EXAMPLE_DATA}'\n",
    "print(f'Example data : {EXAMPLE_DATA}')\n",
    "print(f'Setup folder : {SETUP_DIR}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adc77bb4",
   "metadata": {},
   "source": [
    "---\n",
    "## 2  Structural geometry\n",
    "\n",
    "`GeometryProcessor` stores the 3-D coordinates of measurement nodes, the\n",
    "connectivity lines between them, and optional parent–child relationships used\n",
    "when sensors do not sit directly on structural nodes.\n",
    "\n",
    "| File | Content | Required |\n",
    "|------|---------|----------|\n",
    "| `grid.txt` | Node names and (x, y, z) coordinates | Yes |\n",
    "| `lines.txt` | Pairs of connected node names | No |\n",
    "| `parent_child_assignments.txt` | Sensor-to-node offsets | No |\n",
    "\n",
    "See **Input File Formats** in the documentation for the exact column layout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57d97be2",
   "metadata": {},
   "outputs": [],
   "source": [
    "geometry_data = GeometryProcessor.load_geometry(\n",
    "    nodes_file=EXAMPLE_DATA / 'grid.txt',\n",
    "    lines_file=EXAMPLE_DATA / 'lines.txt',\n",
    "    parent_childs_file=EXAMPLE_DATA / 'parent_child_assignments.txt',\n",
    ")\n",
    "\n",
    "print(f'{len(geometry_data.nodes)} nodes, {len(geometry_data.lines)} lines loaded')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5082a233",
   "metadata": {},
   "source": [
    "---\n",
    "## 3  Loading measurement signals\n",
    "\n",
    "### 3.1  Direct construction (recommended for notebooks)\n",
    "\n",
    "The most transparent approach is to construct `PreProcessSignals` directly\n",
    "from a NumPy array.  All parameters are set in Python — no config files needed:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66d95bab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ── Load your measurement data ───────────────────────────────────────────────\n",
    "# Replace this line with whatever loads your file format:\n",
    "signals = np.load(SETUP_DIR / 'measurement_1.npy')  # shape (n_samples, n_channels)\n",
    "\n",
    "# ── Construct PreProcessSignals ───────────────────────────────────────────────\n",
    "prep_signals = PreProcessSignals(\n",
    "    signals,\n",
    "    sampling_rate=256,            # Hz — from your data acquisition system\n",
    "    ref_channels=[3, 4],          # column indices of reference sensors\n",
    "    accel_channels=[3, 4],        # columns with accelerometer data\n",
    "    velo_channels=[0, 1, 2],      # columns with velocimeter data\n",
    "    disp_channels=[],             # columns with displacement sensor data\n",
    "    setup_name='measurement_1',\n",
    ")\n",
    "\n",
    "# ── Assign channel DOFs for mode-shape plotting ───────────────────────────────\n",
    "# Format: (channel_number, node_name, azimuth_deg, elevation_deg, label)\n",
    "prep_signals.chan_dofs = [\n",
    "    (0, '5',   28.9,  -8.7, 'vib_l'),\n",
    "    (1, '5',   81.0,  -7.3, 'vib_r'),\n",
    "    (2, '5',   55.9,  29.3, 'vib_t'),\n",
    "    (3, '24',   0.0, 180.0, 'ref_x'),\n",
    "    (4, '24', -90.0,   0.0, 'ref_y'),\n",
    "]\n",
    "\n",
    "print(\n",
    "    f'Loaded   : {prep_signals.total_time_steps} samples'\n",
    "    f'  x {prep_signals.num_analised_channels} channels'\n",
    "    f'  ({prep_signals.duration:.1f} s  @  {prep_signals.sampling_rate} Hz)'\n",
    ")\n",
    "print(f'Ref. chans : {prep_signals.ref_channels}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89c3d0be",
   "metadata": {},
   "source": [
    "> **Alternative – config file loader**\n",
    "> If you prefer to keep metadata in text files (useful for batch processing many\n",
    "> setups with identical settings), set `PreProcessSignals.load_measurement_file`\n",
    "> and call `init_from_config` instead:\n",
    "> ```python\n",
    "> PreProcessSignals.load_measurement_file = np.load\n",
    "> prep_signals = PreProcessSignals.init_from_config(\n",
    ">     conf_file=SETUP_DIR / 'setup_info.txt',\n",
    ">     meas_file=SETUP_DIR / 'measurement_1.npy',\n",
    ">     chan_dofs_file=SETUP_DIR / 'channel_dofs.txt',\n",
    "> )\n",
    "> ```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c49741e",
   "metadata": {},
   "source": [
    "### 3.2  Inspect raw signals and spectra\n",
    "\n",
    "`plot_signals` shows time histories and auto-spectral densities for each channel.\n",
    "Use this to check for offsets, spikes, and to identify the frequency band that\n",
    "contains structural modes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59ccb007",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_signals = SignalPlot(prep_signals)\n",
    "plot_signals.plot_signals(per_channel_axes=True, n_lines=512)\n",
    "plt.suptitle('Raw signals', fontsize=12)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "geo-viz-header",
   "metadata": {},
   "source": [
    "### 3.3  Geometry with sensor positions\n",
    "\n",
    "`ModeShapePlot` can visualise the structural geometry together with the channel-DOF arrows showing each sensor's measurement direction before any identification is run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "geo-viz-code",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nodes (grey dots), structural lines and channel-DOF direction arrows\n",
    "geo_plot = ModeShapePlot(geometry_data=geometry_data, prep_signals=prep_signals)\n",
    "geo_plot.reset_view()\n",
    "geo_plot.draw_nodes()\n",
    "geo_plot.draw_lines()\n",
    "geo_plot.draw_chan_dofs()\n",
    "geo_plot.refresh_parent_childs(False)\n",
    "geo_plot.subplot.view_init(elev=20, azim=110)\n",
    "display(geo_plot.fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "corr-header",
   "metadata": {},
   "source": [
    "### 3.4  Correlation functions and power spectral densities\n",
    "\n",
    "After calling `prep_signals.correlation()` you can inspect the estimated correlation functions and their PSDs — useful to verify ``m_lags`` is large enough and the frequency content looks reasonable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "corr-code",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run after pre-processing (Step 4) to inspect correlation functions\n",
    "sig_plot2 = SignalPlot(prep_signals)\n",
    "sig_plot2.plot_signals(per_channel_axes=True, timescale='lags', psd_scale='db')\n",
    "plt.suptitle('Correlation functions and PSDs', fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fa0a3e0",
   "metadata": {},
   "source": [
    "---\n",
    "## 4  Signal pre-processing\n",
    "\n",
    "### 4.1  Decimation\n",
    "\n",
    "Reduces the sampling rate by an integer factor.  An anti-aliasing filter is\n",
    "applied automatically before down-sampling.\n",
    "\n",
    "The example data was recorded at 256 Hz.  Decimating twice by 3 gives\n",
    "≈ 28.4 Hz — sufficient for modes below ~12 Hz.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3822952",
   "metadata": {},
   "outputs": [],
   "source": [
    "prep_signals.decimate_signals(3)\n",
    "prep_signals.decimate_signals(3)\n",
    "\n",
    "print(f'After decimation: {prep_signals.sampling_rate:.2f} Hz  '\n",
    "      f'({prep_signals.total_time_steps} samples  '\n",
    "      f'{prep_signals.duration:.1f} s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e49cb18e",
   "metadata": {},
   "source": [
    "### 4.2  Correlation functions\n",
    "\n",
    "Covariance-driven SSI and PLSCF need the cross-correlation matrix between all\n",
    "channels and the reference channels.  pyOMA provides two estimators:\n",
    "\n",
    "| Method | Function | Notes |\n",
    "|--------|----------|-------|\n",
    "| Welch | `corr_welch` | Fast; slight damping bias from windowing |\n",
    "| Blackman–Tukey | `corr_blackman_tukey` | Better resolution; recommended for SSI-cov |\n",
    "\n",
    "`m_lags` must exceed `num_block_columns + num_block_rows` (set in Step 5).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29572fc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "prep_signals.corr_blackman_tukey(m_lags=200)\n",
    "\n",
    "print(f'Correlation matrix shape: {prep_signals.corr_matrix.shape}'\n",
    "      f'  (channels x refs x lags)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6351dfab",
   "metadata": {},
   "source": [
    "---\n",
    "## 5  System identification\n",
    "\n",
    "pyOMA implements several classical OMA algorithms.  Pick the class and\n",
    "call the two core methods — no config file required.\n",
    "\n",
    "| Class | Algorithm | Call sequence |\n",
    "|-------|-----------|---------------|\n",
    "| `BRSSICovRef` | Covariance-driven SSI | `build_toeplitz_cov` → `compute_modal_params` |\n",
    "| `SSIData` | Data-driven SSI | `build_block_hankel` → `compute_modal_params` |\n",
    "| `VarSSIRef` | SSI + uncertainty | `build_subspace_mat` → `compute_state_matrices` → `prepare_sensitivities` → `compute_modal_params` |\n",
    "| `PLSCF` | Poly-ref. LS Complex Freq. | `build_half_spectra` → `compute_modal_params` |\n",
    "\n",
    "**Key parameters**\n",
    "\n",
    "- `num_block_columns` — number of block columns in the Toeplitz / Hankel matrix.\n",
    "  Must satisfy `num_block_columns + num_block_rows < m_lags` (here < 200).\n",
    "  Typical range: 50–150.\n",
    "- `max_model_order` — modal parameters are estimated at every order 1 … max.\n",
    "  Set to roughly 2× the number of expected modes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7149734",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ── Choose a method by uncommenting one block ─────────────────────────────────\n",
    "\n",
    "# --- BRSSICovRef (Covariance-driven SSI — most common) -----------------------\n",
    "modal_data = BRSSICovRef(prep_signals)\n",
    "modal_data.build_toeplitz_cov(num_block_columns=100)\n",
    "modal_data.compute_modal_params(max_model_order=40)\n",
    "\n",
    "# --- SSIData (Data-driven SSI) -----------------------------------------------\n",
    "# modal_data = SSIData(prep_signals)\n",
    "# modal_data.build_block_hankel(num_block_rows=100)\n",
    "# modal_data.compute_modal_params(max_model_order=40)\n",
    "\n",
    "# --- PLSCF (frequency-domain) -------------------------------------------------\n",
    "# modal_data = PLSCF(prep_signals)\n",
    "# modal_data.build_half_spectra(\n",
    "#     nperseg=200,           # same as m_lags above\n",
    "#     begin_frequency=0.0,   # Hz\n",
    "#     end_frequency=12.0,    # Hz\n",
    "# )\n",
    "# modal_data.compute_modal_params(max_model_order=40)\n",
    "# ─────────────────────────────────────────────────────────────────────────────\n",
    "\n",
    "print(f'Method   : {type(modal_data).__name__}')\n",
    "print(f'Orders   : 1 - {modal_data.max_model_order}')\n",
    "print(f'Channels : {modal_data.num_analised_channels}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df5c37d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "stabil_calc = StabilCluster(modal_data)\n",
    "\n",
    "stabil_calc.calculate_stabilization_masks(\n",
    "    order_range=(0, 1, modal_data.max_model_order),\n",
    "    d_range=(0, 0.10),   # damping 0–10 %\n",
    "    df_max=0.01,          # max relative frequency change between orders\n",
    "    dd_max=0.05,          # max relative damping change between orders\n",
    "    dmac_max=0.05,        # max MAC difference between orders\n",
    ")\n",
    "\n",
    "n_stable = stabil_calc.masks['mask_stable'].sum()\n",
    "print(f'{n_stable} stable poles found across all model orders')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29a91e9b",
   "metadata": {},
   "source": [
    "### 6.1  Static stabilisation diagram\n",
    "\n",
    "Marker colours indicate stability status:\n",
    "\n",
    "- **Light grey dots** — all poles\n",
    "- **Coloured markers** — poles stable in individual criteria\n",
    "- **Black stars / filled markers** — stable in *all* active criteria simultaneously"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed845d05",
   "metadata": {},
   "outputs": [],
   "source": [
    "stabil_plot = StabilPlot(stabil_calc)\n",
    "stabil_plot.plot_stabil('plot_stable')\n",
    "plt.tight_layout\n",
    "display(stabil_plot.fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9df866ca",
   "metadata": {},
   "source": [
    "### 6.2  Interactive stabilisation GUI (requires `pyOMA[jupyter]`)\n",
    "\n",
    "Install the optional Jupyter dependencies with:\n",
    "```\n",
    "pip install \"pyOMA[jupyter]\"\n",
    "```\n",
    "then run the cell below.  The widget lets you **click on poles** to select them\n",
    "and see the corresponding mode shape animated in 3-D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7aaa729b",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    from pyOMA.GUI.JupyterGUI import StabilGUIWeb\n",
    "    stabil_widget, snap_cursor = StabilGUIWeb(stabil_plot)\n",
    "    display(stabil_widget)\n",
    "except ImportError:\n",
    "    print('ipywidgets / ipympl not installed.  '\n",
    "          'Run: pip install \"pyOMA[jupyter]\"')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5c0b365",
   "metadata": {},
   "source": [
    "---\n",
    "## 7  Automated mode selection and result export\n",
    "\n",
    "`automatic_clearing` removes poles that are clearly spurious by applying soft\n",
    "statistical criteria.  For a fully headless workflow you can follow it with\n",
    "`automatic_classification` and `automatic_selection` (see the API reference).\n",
    "\n",
    "`export_results` writes a tab-separated summary table with frequencies, damping\n",
    "ratios, and quality metrics for all selected modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c05b8b36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Automated clearing\n",
    "try:\n",
    "    stabil_calc.automatic_clearing()\n",
    "    print(f'Poles remaining after clearing: '\n",
    "          f'{stabil_calc.masks[\"mask_stable\"].sum()}')\n",
    "except Exception as exc:\n",
    "    print(f'Automatic clearing skipped ({exc})')\n",
    "\n",
    "# Export to a text file\n",
    "import tempfile, os\n",
    "with tempfile.NamedTemporaryFile(suffix='.txt', delete=False, mode='w') as tf:\n",
    "    export_path = tf.name\n",
    "\n",
    "stabil_calc.export_results(export_path)\n",
    "print(f'\\nResults written to {export_path}')\n",
    "\n",
    "# Print the first lines\n",
    "with open(export_path) as f:\n",
    "    for line in f.readlines()[:20]:\n",
    "        print(line, end='')\n",
    "os.unlink(export_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c398e8",
   "metadata": {},
   "source": [
    "---\n",
    "## 8  Mode shape visualisation\n",
    "\n",
    "`ModeShapePlot` animates selected mode shapes on the 3-D geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0da3a40b",
   "metadata": {},
   "outputs": [],
   "source": [
    "msh_plot = ModeShapePlot(\n",
    "    geometry_data=geometry_data,\n",
    "    stabil_calc=stabil_calc,\n",
    "    modal_data=modal_data,\n",
    "    prep_signals=prep_signals,\n",
    "    amplitude=20,\n",
    ")\n",
    "\n",
    "# Interactive Jupyter version (requires pyOMA[jupyter])\n",
    "try:\n",
    "    from pyOMA.GUI.JupyterGUI import PlotMSHWeb\n",
    "    display(PlotMSHWeb(msh_plot))\n",
    "except ImportError:\n",
    "    print('Jupyter widget not available.  '\n",
    "          'For the desktop GUI: from pyOMA.GUI.PlotMSHGUI import start_msh_gui; '\n",
    "          'start_msh_gui(msh_plot)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ce5ceb7",
   "metadata": {},
   "source": [
    "---\n",
    "## 9  Saving and resuming work\n",
    "\n",
    "Every major object serialises to a compressed NumPy archive.  Checkpoint\n",
    "between stages to avoid recomputing expensive steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46d4ba7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tempfile\n",
    "\n",
    "with tempfile.TemporaryDirectory() as tmpdir:\n",
    "    tmp = Path(tmpdir)\n",
    "\n",
    "    # --- save ----------------------------------------------------------------\n",
    "    prep_signals.save_state(tmp / 'prep_signals.npz')\n",
    "    modal_data.save_state(tmp / 'modal_data.npz')\n",
    "    stabil_calc.save_state(tmp / 'stabil_data.npz')\n",
    "\n",
    "    # --- reload --------------------------------------------------------------\n",
    "    prep2   = PreProcessSignals.load_state(tmp / 'prep_signals.npz')\n",
    "    modal2  = BRSSICovRef.load_state(tmp / 'modal_data.npz', prep2)\n",
    "    stabil2 = StabilCalc.load_state(tmp / 'stabil_data.npz', modal2)\n",
    "\n",
    "    import numpy.testing as npt\n",
    "    npt.assert_array_equal(prep2.signals, prep_signals.signals)\n",
    "    npt.assert_array_equal(modal2.modal_frequencies,\n",
    "                           modal_data.modal_frequencies)\n",
    "    print('Save / load round-trip: OK')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28962c98",
   "metadata": {},
   "source": [
    "---\n",
    "## Next steps\n",
    "\n",
    "- **Multi-setup analysis** — see `PogerSSICovRef` and `MergePoSER` for the\n",
    "  PoGer / PoSER methodology when all channels cannot be measured simultaneously.\n",
    "- **Uncertainty quantification** — `VarSSIRef` adds frequency and damping\n",
    "  uncertainty estimates to the identified modes.\n",
    "- **Desktop GUI** — `from pyOMA.GUI.StabilGUI import start_stabil_gui` after\n",
    "  `pip install \"pyOMA[gui]\"`.\n",
    "- **API reference** — https://py-oma.readthedocs.io"
   ]
  }
 ],
 "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.14.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}