{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Compare dv/v from MWCS, Stretching, and WCT\n\nPlots the network-level dv/v time series produced by the three available\nmethods side by side \u2014 MWCS dt/t, Stretching, and Wavelet Cross-Transform\n(WCT) \u2014 so their agreement and systematic differences can be assessed at a\nglance.\n\nThis example is designed to run after the full MSNoise test workflow (all\ncompute steps completed).  Each ``(components, mov_stack)`` combination for\nwhich at least one method has results produces its own figure.  Methods that\nwere not computed are silently skipped.\n\nEach figure has two panels:\n\n* **Top** \u2014 dv/v time series for each available method, with a \u00b11\u03c3 shaded\n  band.\n* **Bottom** \u2014 residuals relative to the first available method (e.g.\n  MWCS \u2212 Stretching, MWCS \u2212 WCT) to highlight systematic offsets.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\n\nif \"SPHINX_DOC_BUILD\" in os.environ:\n    if \"MSNOISE_DOC\" in os.environ:\n        os.chdir(os.environ[\"MSNOISE_DOC\"])\n\nimport matplotlib\nmatplotlib.use(\"agg\")\n\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport pandas as pd\n\nfrom msnoise.api import connect\nfrom msnoise.results import MSNoiseResult\nfrom msnoise.core.config import lineage_to_plot_tag, build_plot_outfile\n\n# \u2500\u2500 palette & style \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nplt.style.use(\"ggplot\")\n\nMETHOD_STYLE = {\n    \"mwcs_dtt_dvv\":    dict(color=\"#1f77b4\", label=\"MWCS\",       lw=1.8, zorder=3),\n    \"stretching_dvv\":  dict(color=\"#d62728\", label=\"Stretching\", lw=1.8, zorder=2),\n    \"wavelet_dtt_dvv\": dict(color=\"#2ca02c\", label=\"WCT\",        lw=1.8, zorder=1),\n}\n\n# \u2500\u2500 connect \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\ndb = connect()\n\n# \u2500\u2500 discover available DVV results \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nresults_by_method = {}\nfor category in (\"mwcs_dtt_dvv\", \"stretching_dvv\", \"wavelet_dtt_dvv\"):\n    found = MSNoiseResult.list(db, category)\n    if found:\n        results_by_method[category] = found[0]\n        tag = lineage_to_plot_tag(found[0].lineage_names)\n        print(f\"  {category}: {tag}\")\n    else:\n        print(f\"  {category}: not computed \u2014 skipping\")\n\nif not results_by_method:\n    print(\"No DVV results found \u2014 run the full pipeline first.\")\n    db.close()\n    raise SystemExit(0)\n\nPAIR_TYPE = \"CC\"\n\n# \u2500\u2500 collect all (components, mov_stack) combinations across methods \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nall_keys = set()\nraw = {}\nfor category, result in results_by_method.items():\n    data = result.get_dvv(pair_type=PAIR_TYPE)\n    filtered = {(comp, ms): ds\n                for (pt, comp, ms), ds in data.items()\n                if pt == PAIR_TYPE}\n    raw[category] = filtered\n    all_keys |= set(filtered.keys())\n\nif not all_keys:\n    print(\"No CC DVV data found.\")\n    db.close()\n    raise SystemExit(0)\n\n# \u2500\u2500 one figure per (components, mov_stack) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nfor comp, ms in sorted(all_keys):\n    ms_label = f\"{ms[0]}-{ms[1]}\"\n    title    = f\"dv/v comparison  |  comp={comp}  |  mov_stack={ms_label}\"\n\n    # Collect time series (converted to %)\n    series = {}\n    for category in (\"mwcs_dtt_dvv\", \"stretching_dvv\", \"wavelet_dtt_dvv\"):\n        if category not in raw or (comp, ms) not in raw[category]:\n            continue\n        ds = raw[category][(comp, ms)]\n        if \"mean\" not in ds:\n            continue\n        mean_da = ds[\"mean\"]\n        std_da  = ds.get(\"std\", ds.get(\"weighted_std\", None))\n        times   = pd.DatetimeIndex(mean_da.coords[\"times\"].values)\n        mean_s  = pd.Series(mean_da.values * 100.0, index=times, name=category)\n        std_s   = (pd.Series(std_da.values * 100.0, index=times)\n                   if std_da is not None else None)\n        series[category] = (mean_s, std_s)\n\n    if not series:\n        print(f\"  No data for comp={comp}, mov_stack={ms_label} \u2014 skipping\")\n        continue\n\n    # \u2500\u2500 figure: 2 rows \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    fig = plt.figure(figsize=(12, 7))\n    gs  = gridspec.GridSpec(2, 1, height_ratios=[3, 1], hspace=0.08)\n    ax_dvv = fig.add_subplot(gs[0])\n    ax_res = fig.add_subplot(gs[1], sharex=ax_dvv)\n\n    # Top panel\n    reference_s = None\n    for category, (mean_s, std_s) in series.items():\n        style = METHOD_STYLE[category]\n        ax_dvv.plot(mean_s.index, mean_s.values,\n                    color=style[\"color\"], lw=style[\"lw\"],\n                    label=style[\"label\"], zorder=style[\"zorder\"])\n        if std_s is not None:\n            ax_dvv.fill_between(mean_s.index,\n                                mean_s.values - std_s.values,\n                                mean_s.values + std_s.values,\n                                color=style[\"color\"], alpha=0.15,\n                                zorder=style[\"zorder\"] - 1)\n        if reference_s is None:\n            reference_s = mean_s\n\n    ax_dvv.axhline(0, color=\"k\", lw=0.8, ls=\"--\", alpha=0.5)\n    ax_dvv.set_ylabel(\"dv/v (%)\")\n    ax_dvv.set_title(title)\n    ax_dvv.legend(loc=\"upper right\", framealpha=0.9)\n    ax_dvv.tick_params(labelbottom=False)\n\n    # Bottom panel \u2014 residuals vs first method\n    ref_label = METHOD_STYLE[next(iter(series))][\"label\"]\n    has_residual = False\n    for category, (mean_s, _) in series.items():\n        if mean_s is reference_s:\n            continue\n        style = METHOD_STYLE[category]\n        common = reference_s.index.intersection(mean_s.index)\n        if len(common) == 0:\n            continue\n        resid = reference_s.loc[common] - mean_s.loc[common]\n        ax_res.plot(resid.index, resid.values,\n                    color=style[\"color\"], lw=1.4,\n                    label=f\"{ref_label} \\u2212 {style['label']}\",\n                    zorder=style[\"zorder\"])\n        has_residual = True\n\n    ax_res.axhline(0, color=\"k\", lw=0.8, ls=\"--\", alpha=0.5)\n    ax_res.set_ylabel(\"Residual (%)\")\n    ax_res.set_xlabel(\"Date\")\n    if has_residual:\n        ax_res.legend(loc=\"upper right\", framealpha=0.9, fontsize=8)\n\n    fig.autofmt_xdate()\n\n    # \u2500\u2500 save \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    first_result = next(iter(results_by_method.values()))\n    outfile = build_plot_outfile(\n        \"?.png\", \"dvv_comparison\",\n        first_result.lineage_names,\n        components=comp,\n        mov_stack=ms,\n    )\n    print(f\"  Saving: {outfile}\")\n    plt.savefig(outfile, dpi=150, bbox_inches=\"tight\")\n    plt.show()\n\ndb.close()\n\n# EOF"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.10.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}