commit 06ce53bc4456dc551328212c80a09e2b9ca7f707 Author: Kemal Yaylali Date: Sat Aug 15 20:43:51 2026 +0100 Initial release: Ashvale Station 1.0.0 diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md new file mode 100644 index 0000000..29fe1cf --- /dev/null +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -0,0 +1,24 @@ +--- +name: Bug report +about: Something crashed or behaved wrongly +labels: bug +--- + +**What happened** + +**What you expected instead** + +**Hardware** +- Board: (e.g. Pi Zero 2 W) +- Sense HAT: (v1 / v2 / none, running the simulator) +- OS and Python version: + +**Output of `GET /api/status`** + +```json +``` + +**Relevant log lines** (`journalctl -u ashvale -n 50`) + +``` +``` diff --git a/.github/ISSUE_TEMPLATE/forecast_quality.md b/.github/ISSUE_TEMPLATE/forecast_quality.md new file mode 100644 index 0000000..b26a65f --- /dev/null +++ b/.github/ISSUE_TEMPLATE/forecast_quality.md @@ -0,0 +1,22 @@ +--- +name: Forecast quality +about: The models are producing poor or strange forecasts +labels: forecasting +--- + +**What looks wrong** + +**Output of `python scripts/evaluate.py`** + +``` +``` + +**Scorecard from the Models tab** (or `GET /api/scorecard`) + +```json +``` + +**How long has the station been logging?** (`history_days` from `/api/status`) + +**Is it indoors?** And is `site.altitude_m` set correctly in your config? +These two account for most reported forecast oddities. diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..48fea56 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,64 @@ +name: CI + +on: + push: + branches: [main] + pull_request: + branches: [main] + +jobs: + test: + runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + python-version: ["3.9", "3.11", "3.12"] + + steps: + - uses: actions/checkout@v4 + + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + cache: pip + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -r requirements.txt + pip install httpx + + - name: Byte-compile every module + run: python -m compileall -q ashvale scripts run.py + + - name: Import check + run: python -c "import ashvale.api; print('imports ok')" + + - name: Seed synthetic history + run: python scripts/simulate.py --days 10 --wipe + + # The backtest is the real test: it exercises features, the RLS heads, + # climatology and conformal calibration end to end, and fails loudly if + # any of them stop producing finite numbers. + - name: Walk-forward backtest + run: python scripts/evaluate.py --train-frac 0.6 --hours 300 + + - name: Exercise the API surface + run: | + python - <<'PY' + from fastapi.testclient import TestClient + import ashvale.api as api + with TestClient(api.app) as c: + for ep in ["/api/telemetry", "/api/status", "/api/forecast", + "/api/outlook", "/api/precipitation", "/api/anomaly", + "/api/models", "/api/scorecard", "/api/records", + "/api/storage", "/api/methods", + "/api/history/range?hours=24", + "/api/history/daily?days=5", "/"]: + r = c.get(ep) + assert r.status_code == 200, (ep, r.status_code, r.text[:200]) + assert c.post("/api/train").json()["trained"] is True + assert len(c.get("/api/export.csv?hours=6").text.splitlines()) > 1 + print("api surface ok") + PY diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a61f359 --- /dev/null +++ b/.gitignore @@ -0,0 +1,35 @@ +# Station runtime state. Never commit these: the database is your logged +# history and the state file holds trained model parameters. Both are +# machine-specific and both grow without bound. +data/ +*.db +*.db-wal +*.db-shm +*.log + +# Python +__pycache__/ +*.py[cod] +*.egg-info/ +build/ +dist/ +.venv/ +venv/ +.pytest_cache/ +.ruff_cache/ +.mypy_cache/ + +# Editors and OS +.vscode/ +.idea/ +*.swp +.DS_Store +Thumbs.db + +# Local overrides: keep config.yaml tracked as the documented default, +# but let people keep a private one with their own coordinates. +config.local.yaml + +# Personal deployment notes: LAN addresses, hostnames, SSH details. +# Kept locally, never published. +HANDOVER.md diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..ec7fd03 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,35 @@ +# Changelog + +All notable changes to this project are documented here. +Format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/); +versioning follows [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## [1.0.0] - 2026-08-15 + +First public release. + +### Added +- Multi-horizon forecasting: 18 direct heads (3 targets x 6 horizons) using + exponentially weighted recursive least squares with a capped covariance trace. +- Hedge-blended ensemble over persistence, climatology and the learned model, + so the system can conclude that the learned model is not worth using. +- Adaptive conformal prediction intervals with coverage feedback. +- Constant-velocity Kalman bank (Joseph form) for level and rate estimation. +- Grey-box CPU self-heating compensation with an RLS-estimated coefficient, + calibrated from a single trusted thermometer reading. +- Harmonic climatology with anomaly decay for a 7-day outlook, with annual + terms gated behind 120 days of history. +- Precipitation model: Zambretti barometric prior plus an online logistic + residual learner with human-in-the-loop labelling. +- Monitoring: Mahalanobis EWMA novelty, Page-Hinkley drift detection that + triggers retraining, and stuck-sensor detection. +- Verification loop scoring every matured forecast against persistence and + climatology, surfaced as a public scorecard. +- Tiered storage: raw to 5-minute to hourly downsampling, roughly 30 MB per year. +- Arbitrary-range history queries, per-day summaries, all-time records and + streamed CSV export. +- Five-tab dashboard sized to a single viewport, with a Methods tab generated + from the live configuration. +- LED matrix driver with pressure-trend arrows, rain bars and alert pulses. +- Physics-based simulator fallback so the suite runs without a Sense HAT. +- Backfill and walk-forward backtest scripts. diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..bc60a76 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,88 @@ +# Contributing to Ashvale Station + +Thanks for taking an interest. This is a small project maintained by one person, +so the bar here is "make it easy to say yes", not "follow a 40-page process". + +## Ground rules that actually matter + +**No heavy dependencies.** The whole point is that this runs on a Raspberry Pi +Zero 2 W with 512 MB of RAM. Pull requests adding PyTorch, TensorFlow, pandas or +scikit-learn to the core will be declined, however elegant. If a model genuinely +needs one, make it an optional extra with a numpy fallback. + +**Claims need numbers.** If you say a change improves forecasting, show the +output of `scripts/evaluate.py` before and after, on the same data. Skill against +persistence is the metric that counts. "It looks better" is not evidence. + +"The same data" means bit-identical, which takes two flags, not one. `--seed` +alone is not enough: the synthetic history is anchored to wall clock, so the OU +realisation repeats while the timestamps shift, and that moves solar elevation, +day of year and the seasonal harmonic. Those feed the temperature and humidity +models directly, so two same-seed runs give you different data and an +uninterpretable comparison. Pin both: + +```bash +python scripts/simulate.py --days 21 --wipe --seed 11 --end 1767225600 +python scripts/evaluate.py --train-frac 0.6 --hours 100000 +``` + +The wide `--hours` matters: `evaluate.py` looks back from now, and its default +window is 60 days, so a history pinned to a timestamp further in the past than +that falls outside it and the script reports "Not enough history". + +**Uncertainty must stay calibrated.** If you touch the forecasting path, check +that coverage on the scorecard still sits near the target. A model that gets more +accurate while its intervals start lying is a regression, not an improvement. + +**Document how it fails.** Every model module carries a docstring explaining not +just what the technique is but how it breaks in the field. Keep that up. It is +the most useful part of this codebase. + +## Getting set up + +```bash +git clone https://github.com/lynchaos/ashvale-station.git +cd ashvale-station +python3 -m venv .venv && source .venv/bin/activate +pip install -r requirements.txt + +python scripts/simulate.py --days 21 --wipe # synthetic history +python scripts/evaluate.py # baseline numbers +python run.py --no-led # dashboard on :8000 +``` + +No Sense HAT needed. The simulator kicks in automatically and exercises every +code path. + +## Before you open a pull request + +- [ ] `python scripts/evaluate.py` runs clean, and you have posted before/after numbers +- [ ] Those numbers came from a backfill with `--seed` and `--end` both pinned +- [ ] Coverage on the scorecard is still near target for anything you touched +- [ ] No new required dependencies +- [ ] New model code explains its failure mode in the docstring +- [ ] The dashboard still fits one viewport at 1280x800 if you changed the UI + +## Good first contributions + +- **DS18B20 or BME280 support.** An outdoor sensor removes the single biggest + limitation in the project. High impact, self-contained. +- **Tipping-bucket rain gauge on GPIO.** Real precipitation labels would + transform the precipitation model. +- **METAR ingestion** from a nearby airfield as a calibration reference. +- **Translations** for the dashboard. +- **Tests.** There is a walk-forward backtest but no unit test suite. A pytest + suite over `physics.py`, `estimation.py` and `models/rls.py` would be very welcome. + +## Reporting bugs + +Open an issue with your `config.yaml` (redact coordinates if you like), the +output of `GET /api/status`, and what you expected instead. If it is a +forecasting problem rather than a crash, the output of `scripts/evaluate.py` +helps enormously. + +## Licensing of contributions + +By contributing you agree that your work is licensed under the Apache License +2.0, the same terms as the project. You keep the copyright in your own +contributions. diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..23ad6a4 --- /dev/null +++ b/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright 2026 Kemal Yaylali + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/NOTICE b/NOTICE new file mode 100644 index 0000000..3fd0536 --- /dev/null +++ b/NOTICE @@ -0,0 +1,24 @@ +Ashvale Station +Copyright 2026 Kemal Yaylali + +This product includes software developed by Kemal Yaylali +(https://github.com/lynchaos/ashvale-station). + +Licensed under the Apache License, Version 2.0. See the LICENSE file +for the full terms. + +--- + +Third-party components loaded by the browser dashboard at runtime and +NOT redistributed in this repository: + + Chart.js MIT License https://www.chartjs.org/ + chartjs-plugin-zoom MIT License https://github.com/chartjs/chartjs-plugin-zoom + Hammer.js MIT License https://hammerjs.github.io/ + Tailwind CSS MIT License https://tailwindcss.com/ + Plus Jakarta Sans SIL OFL 1.1 https://fonts.google.com/specimen/Plus+Jakarta+Sans + JetBrains Mono SIL OFL 1.1 https://www.jetbrains.com/lp/mono/ + +The Zambretti forecast in ashvale/models/precip.py is an original +re-parameterisation onto the historic 26-point scale, not a transcription +of any particular published implementation. diff --git a/README.md b/README.md new file mode 100644 index 0000000..c00c4d8 --- /dev/null +++ b/README.md @@ -0,0 +1,293 @@ +# Ashvale Station + +An online machine learning suite for a Raspberry Pi Zero 2 W with a Sense HAT v2. +It turns the original telemetry dashboard into a forecasting instrument: multi-horizon +predictions with calibrated uncertainty, a verification scorecard that scores the model +against persistence, drift detection that triggers its own retraining, and a +human-in-the-loop labelling path. + +Everything runs on the Pi. No cloud, no GPU, no PyTorch, no scikit-learn, no pandas. +The learners are pure numpy and the whole process sits comfortably under 150 MB RSS. + +``` +python scripts/simulate.py --days 14 --wipe # seed synthetic history +python scripts/evaluate.py # walk-forward backtest +python run.py # serve on :8000 +``` + +--- + +## Why the design looks like this + +A single point sensor on a windowsill is not a weather service, and pretending +otherwise is the fastest way to build something that looks impressive and is +useless. The honest inventory of what your hardware can actually observe: + +| Signal | What it tells you | Useful range | +| --- | --- | --- | +| Pressure and its tendency | Synoptic systems, and it passes through walls | Genuinely hours ahead | +| Temperature, humidity | The local micro-environment | Hours, strongly diurnal | +| Ambient light and colour | Cloudiness, occupancy, time of day | Now | +| IMU | Whether someone knocked the desk | Now | + +So the suite is built around that reality. Short horizons lean on state estimation and +learned dynamics. Long horizons lean on climatology plus a decaying anomaly, and are +labelled an *outlook* rather than a forecast. Every claim gets scored against the +"nothing changes" baseline, in public, on the dashboard. + +### The stack + +``` +sensors.py hardware + a physics-based simulator fallback + | +estimation.py self-heating compensation -> Kalman bank -> level + rate + | +storage.py SQLite, WAL, tiered downsampling (raw -> 5 min -> hourly) + | +features.py 33 features on a 5-minute grid, physics computed not learned + | +models/ + rls.py recursive least squares + adaptive conformal intervals + nowcast.py 18 direct heads (3 targets x 6 horizons), Hedge-blended + climatology.py harmonic regression for the 7-day outlook + precip.py Zambretti prior + online logistic residual learner + anomaly.py Mahalanobis EWMA + Page-Hinkley drift + sensor health + | +station.py four async loops: sample / persist / train / verify +api.py, led.py, dashboard.py +``` + +### Six decisions worth defending + +**1. Self-heating is a grey-box parameter, not a magic constant.** +The HTS221 and LPS25HB sit millimetres above a SoC running 20 to 25 °C hotter than the +room. The usual fix is `T = T_sensor - (T_cpu - T_sensor) / 1.5`. That 1.5 depends on +your case, your orientation, your airflow, and your CPU load. Here it is a single RLS +parameter that you update from the dashboard by typing in a thermometer reading. In +testing it recovers a known coefficient of 0.62 from a prior of 0.30 in **one sample**, +and holds post-calibration bias to 0.012 °C. + +**2. Rates come from a Kalman filter, never a finite difference.** +Pressure tendency is the single most informative variable you have, and the LPS25HB +noise floor makes a naive `(p[t] - p[t-1])/dt` pure noise. A constant-velocity Kalman +filter estimates level and rate jointly, in Joseph form so the covariance stays positive +semi-definite over months of continuous operation. The filtered `dp/dt` is what feeds +both Zambretti and the learned heads. + +**3. Direct multi-horizon heads, not one model iterated forward.** +Iterating a one-step model 288 times to reach 24 hours compounds its own bias into a +beautifully smooth lie. Eighteen small direct heads cost about 150 kB total and each one +is honest about its own horizon. + +**4. RLS with directional forgetting, not SGD.** +A station produces 288 grid rows a day. Sample efficiency is not a nicety. RLS is the +exact minimiser of the exponentially weighted squared error at every step and converges +in far fewer samples. The covariance `P` gives free parameter uncertainty. The forgetting +factor (0.9985, about 11 hours of effective memory) handles seasonal adaptation without +any retraining schedule at all. Plain forgetting inflates `P` exponentially during quiet +nights when the regressor barely moves, so the trace is capped: this is the single most +common way a field RLS deployment detonates. + +**5. Adaptive conformal intervals, not Gaussian error bars.** +Split conformal assumes exchangeability. Weather is not exchangeable: a front arrives and +yesterday's residual quantile becomes fiction. Adaptive conformal inference feeds realised +coverage back into the working alpha, so the band widens after each miss and narrows after +each hit. Measured coverage in the backtest below sits at 89 to 91% against a 90% target, +across every target and horizon. + +**6. The ensemble is allowed to conclude that the model is useless.** +Each head blends persistence, climatology and the learned model with Hedge weights. +At 15-minute pressure the weights land on **96% persistence**, which is the correct +answer, and the scorecard says so out loud. A forecasting system that cannot tell you +when to switch it off is a marketing asset, not an instrument. + +--- + +## Measured performance + +Walk-forward backtest, 14 days of synthetic history, 60/40 split, strictly no target +visible before its validity time. `skill = 1 - MAE/MAE_persistence`. + +``` +target lead MAE persist clim skill cover +temperature 15m 0.439 0.463 0.449 5.0% 90% +temperature 1h 0.707 0.981 0.840 28.0% 89% +temperature 3h 0.824 1.961 1.302 58.0% 91% +temperature 6h 0.770 3.094 1.619 75.1% 90% +temperature 12h 0.816 3.922 1.969 79.2% 90% +temperature 1d 0.834 1.630 1.683 48.9% 91% + +humidity 3h 1.526 2.790 2.826 45.3% 90% +humidity 1d 2.370 14.506 14.523 83.7% 90% + +pressure 15m 0.949 0.950 0.950 0.1% 90% <- persistence wins +pressure 3h 2.606 3.531 3.631 26.2% 89% +pressure 1d 3.389 8.638 8.650 60.8% 90% +``` + +These are numbers against a simulator, so read them as a check that the machinery is +sound rather than as a promise about your windowsill. Run `scripts/evaluate.py` again +after a fortnight of real data and believe those instead. + +--- + +## Install on the Pi + +```bash +sudo apt update && sudo apt install -y python3-venv sense-hat +git clone ~/ashvale-ml && cd ~/ashvale-ml +python3 -m venv --system-site-packages .venv +source .venv/bin/activate +pip install -r requirements.txt +pip install sense-hat smbus2 + +cp systemd/ashvale.service /etc/systemd/system/ # edit User/paths first +sudo systemctl enable --now ashvale +``` + +`--system-site-packages` matters: `sense-hat` pulls in `RTIMULib`, which is installed +via apt and is a genuine ordeal to build inside a clean venv. + +Without the hardware libraries the suite falls back to a simulated board automatically, +so you can develop the whole thing on a laptop and deploy the same code unchanged. + +**Set your altitude in `config.yaml`.** Sea-level pressure reduction is the one setting +people skip and then wonder why every rule-of-thumb forecast reads pessimistic. At 100 m +an uncorrected station pressure shifts the Zambretti number by roughly two categories, +permanently. + +--- + +## The dashboard + +Five tabs, one viewport, no scrolling on desktop. Below 1024 px the constraint is +released, because pinning five panels into a phone viewport produces unreadable +eight-pixel type. + +| Tab | Answers | +| --- | --- | +| **Live** | What is it doing right now | +| **Forecast** | What is it about to do, and how sure are we | +| **History** | What did it do, over any timeframe you ask for | +| **Models** | Has the model earned its confidence | +| **Methods** | How the whole thing is wired, and how each stage fails | + +### History + +Presets from 6 hours to a year, plus an explicit from/to range picker. Aggregation +happens in SQLite, not numpy: pulling 90 days of rows into Python to average them would +cost more memory than the board has. The bucket auto-selects from the span and snaps to +round durations, so 6 hours gives one-minute buckets and a year gives daily ones. Min and +max travel alongside the mean and render as a shaded band, so an hourly view still shows +that the hour spanned four degrees rather than implying a flat line. + +Alongside: per-day minima and maxima in local time, all-time records with the timestamp +each was set, and CSV export of any range (streamed as a generator, so a year of history +never has to exist in memory at once). + +### Methods + +Generated from `methods.py` and rendered against your live config, so it describes the +station you are running rather than the one shipped. Ten stages, each with what it +consumes, what it produces, why it is built that way, and how it fails. The failure mode +is the field that usually goes undocumented and the one you need at 2 a.m. + +- **Fan chart** with the 90% conformal band drawn behind the observed line. As the model + earns confidence the band visibly narrows, so model quality becomes a shape you can + read from the doorway. +- **Estimator internals**, the signature panel: self-heating coefficient, novelty + distance and drift pressure, ticking at 2 Hz. Most weather dashboards show numbers; + this one shows the state estimator working. +- **Conditions ahead**: Zambretti class, rain probability, and the prior/learner/trust + split so you can see how much the learned model is actually contributing. +- **Two yes/no buttons.** "Was it wet in the last hour?" Each press is a strong label + worth ten proxy labels. Two seconds of your attention beats a week of heuristics. +- **Calibration box.** Type a thermometer reading, watch `k` update. +- **Scorecard** with skill against persistence, and coverage against the 90% target. +- **Monitors**: novelty, drift pressure, per-sensor health, ensemble weights. + +### LED matrix + +The 8×8 stopped being a scrolling number. It cycles through glyphs readable across a room: +a pressure-trend arrow coloured by Zambretti class and brightened by tendency magnitude, +a rain-probability column bar, a 3-hour temperature-delta wedge, and a red pulse if a +sensor faults or drift fires. Alerts pre-empt everything, because a six-second scroll is +a six-second delay on the only frame that matters. + +--- + +## API + +| Endpoint | Purpose | +| --- | --- | +| `GET /api/telemetry` | Live reading. **Superset of the original payload**, so existing clients keep working | +| `GET /api/stream` | SSE. One connection instead of a 2-second poll: 0.4% CPU instead of 4% | +| `GET /api/history/range?start=&end=&bucket=` | Any window, SQL-aggregated, auto bucket | +| `GET /api/history/daily?days=` | Per-day min, max and mean in local time | +| `GET /api/records` | All-time extremes, each with its timestamp | +| `GET /api/export.csv?start=&end=` | Streamed CSV export | +| `GET /api/storage` | Rows per resolution tier and database size | +| `GET /api/methods` | The pipeline description the Methods tab renders | +| `GET /api/history?hours=&max_points=` | Decimated history (legacy) | +| `GET /api/forecast?target=` | All horizons with conformal bands and ensemble weights | +| `GET /api/outlook` | Days 2 to 7, climatology plus decaying anomaly, caveat included | +| `GET /api/precipitation` | Zambretti class, rain probability, prior/learner split | +| `GET /api/anomaly` | Novelty, drift, per-sensor health, event log | +| `GET /api/models` | Per-head diagnostics, coverage, precip coefficients | +| `GET /api/scorecard` | Verification: MAE, skill, coverage, sample count | +| `POST /api/train` | Force a retrain | +| `POST /api/verify` | Force a scoring pass | +| `POST /api/label` | `{"kind":"rain","value":1}` strong ground truth | +| `POST /api/calibrate` | `{"reference_c":19.4}` or `{"reset":true}` | +| `GET /api/status` | Hardware, history span, drift, training log | + +--- + +## Tuning + +| Symptom | Knob | +| --- | --- | +| Temperature reads consistently high | Calibrate from the dashboard, or raise `sensor.cpu_heat_k` | +| Readings look over-smoothed, lag real changes | Raise `sensor.kalman_q_temp` | +| Rates look noisy | Lower `sensor.kalman_q_*`, or raise `kalman_r_*` | +| Model adapts too slowly to a season change | Lower `model.rls_forgetting` toward 0.995 | +| Model is jumpy and forgets overnight | Raise it toward 0.9995 | +| Coverage sits well below 90% | Raise `model.conformal_gamma` so it corrects faster | +| Drift alarms constantly | Raise `model.drift_lambda` | +| Retrains eat the CPU | Raise `model.train_period_s`, lower `max_pairs` in `NowcastEnsemble.fit` | + +A full retrain over 18 heads takes about 10 s on a modern x86 core and closer to 60 to +90 s on a Zero 2 W. It runs in a worker thread, so the sample loop, the API and the LED +never stall while it happens. + +--- + +## Honest limitations + +- **Indoors, this forecasts your room, not the sky.** Pressure is the exception: it + passes through walls, which is why the precipitation model runs on pressure and its + tendency rather than on your indoor humidity. Set `site.indoors` truthfully. +- **Days 2 to 7 are climatology, not a forecast.** Labelled as such in the API response + and on the dashboard. They will never catch an incoming Atlantic low, because your + station physically cannot see one. +- **Rain labels are the bottleneck.** Without a gauge the proxy label is deliberately + conservative and abstains in the ambiguous middle. The learner earns trust in + proportion to strong labels: `trust = n / (n + 25)`. Press the buttons. +- **Annual harmonics stay switched off** until 120 days of history exist. Fitting a + 365-day sine to three weeks of data produces a magnificent extrapolation straight off + the edge of the physical world. +- **One uvicorn worker, deliberately.** The station owns mutable model state; a second + worker would give you two divergent forecasters sharing a socket. + +## Where to take it next + +The obvious extensions, roughly in order of payoff per hour of work: + +1. **A DS18B20 on a one-metre cable outside the window.** It removes the indoor caveat + entirely, costs about three pounds, and every model in here improves immediately. +2. **A tipping-bucket rain gauge on a GPIO.** Real precipitation labels turn the logistic + model from a Zambretti wrapper into something genuinely local. +3. **Pull METAR from a nearby airfield** as a reference channel, and the compensator + calibrates itself continuously instead of waiting for you to type a number. +4. **Swap the RLS head for an ensemble Kalman filter over the parameter vector** if you + want proper joint state-parameter estimation. You already have the machinery. diff --git a/SECURITY.md b/SECURITY.md new file mode 100644 index 0000000..3cb20b9 --- /dev/null +++ b/SECURITY.md @@ -0,0 +1,32 @@ +# Security Policy + +## Supported versions + +The latest release on `main` is the supported version. + +## Reporting a vulnerability + +Please email **support@yaylali.uk** rather than opening a public issue. + +Include what you found, how to reproduce it, and what an attacker could achieve. +You should get an acknowledgement within a few days. This is a personal project +maintained in spare time, so please be patient with fix timelines. + +## Deployment note worth reading + +Ashvale Station ships **no authentication and no TLS**. It is designed to sit on +a trusted home network, and the default bind address is `0.0.0.0`, meaning +anything on your LAN can reach it. + +Do not port-forward it to the open internet. If you want remote access, put it +behind a reverse proxy that terminates TLS and handles authentication, or reach +it over a VPN or a WireGuard tunnel. The API includes endpoints that mutate model +state (`/api/train`, `/api/calibrate`, `/api/label`), so an exposed instance is a +system a stranger can degrade. + +To restrict it to the local machine only: + +```yaml +server: + host: 127.0.0.1 +``` diff --git a/ashvale/__init__.py b/ashvale/__init__.py new file mode 100644 index 0000000..58890c1 --- /dev/null +++ b/ashvale/__init__.py @@ -0,0 +1,17 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Ashvale Station: a self-contained ML forecasting suite for Raspberry Pi + Sense HAT v2.""" + +__version__ = "1.0.0" diff --git a/ashvale/api.py b/ashvale/api.py new file mode 100644 index 0000000..78c3cbc --- /dev/null +++ b/ashvale/api.py @@ -0,0 +1,423 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""HTTP surface. Thin by design: every endpoint is a view over station state. + +Backwards compatibility matters here, so `/api/telemetry` returns a +superset of the original payload. Anything already pointed at this Pi +keeps working, and the new fields are simply there when you want them. +""" + +from __future__ import annotations + +import asyncio +import json +import time +from contextlib import asynccontextmanager +from typing import Any, Dict, List, Optional + +import numpy as np +from fastapi import FastAPI, HTTPException, Query +from fastapi.responses import HTMLResponse, StreamingResponse +from pydantic import BaseModel, Field + +from .config import CONFIG +from .dashboard import DASHBOARD_HTML +from .led import LedDisplay +from .methods import describe +from .station import Station + +station: Optional[Station] = None +display: Optional[LedDisplay] = None + + +@asynccontextmanager +async def lifespan(app: FastAPI): + global station, display + station = Station(CONFIG) + station.sample_once() + station.start() + if CONFIG.server.led_enabled: + display = LedDisplay(station, CONFIG.server.led_cycle_s) + display.start() + try: + yield + finally: + if display is not None: + await display.stop() + if station is not None: + await station.stop() + + +app = FastAPI( + title="Ashvale Station", + version="1.0.0", + description="Sense HAT v2 telemetry with online forecasting, calibrated " + "uncertainty, drift detection and verification.", + lifespan=lifespan, +) + + +def _st() -> Station: + if station is None: + raise HTTPException(503, "station not started") + return station + + +def _clean(obj: Any) -> Any: + """JSON is not a superset of IEEE 754. NaN in a response body will + silently break a browser's JSON.parse, which is a miserable bug to + chase from a dashboard that just shows dashes.""" + if isinstance(obj, dict): + return {k: _clean(v) for k, v in obj.items()} + if isinstance(obj, (list, tuple)): + return [_clean(v) for v in obj] + if isinstance(obj, (np.floating, float)): + f = float(obj) + return None if (f != f or f in (float("inf"), float("-inf"))) else round(f, 6) + if isinstance(obj, (np.integer,)): + return int(obj) + if isinstance(obj, np.ndarray): + return _clean(obj.tolist()) + return obj + + +# --------------------------------------------------------------- models + +class LabelIn(BaseModel): + kind: str = Field("rain", description="rain | fog | frost | window_open") + value: float = Field(..., ge=0.0, le=1.0) + ts: Optional[float] = None + note: str = "" + + +class CalibrationIn(BaseModel): + reference_c: Optional[float] = Field(None, description="Trusted air temperature in C") + reset: bool = Field(False, description="Discard the learned coefficient and its " + "covariance, returning to the configured prior") + + +# ------------------------------------------------------------ endpoints + +@app.get("/api/telemetry") +def telemetry() -> Dict: + st = _st() + live = st.live or st.sample_once() + colour = live.get("colour") or {} + return _clean({ + # original contract, preserved + "timestamp": live.get("timestamp"), + "temperature": live.get("temp_smooth"), + "humidity": live.get("hum_smooth"), + "pressure": live.get("press_slp"), + "compass": live.get("compass"), + "pitch": live.get("pitch"), + "roll": live.get("roll"), + "yaw": live.get("yaw"), + "accel": {"x": live.get("ax"), "y": live.get("ay"), "z": live.get("az")}, + "gyro": {"x": live.get("gx"), "y": live.get("gy"), "z": live.get("gz")}, + "color": {"clear": colour.get("clear", live.get("lux", 0)), + "red": colour.get("red", live.get("r", 0)), + "green": colour.get("green", live.get("g", 0)), + "blue": colour.get("blue", live.get("b", 0)), + "hex": colour.get("hex", "#334155"), + "cct": colour.get("cct")}, + # everything the ML layer adds + "temperature_raw": live.get("temp_raw"), + "temperature_compensated": live.get("temp_c"), + "pressure_station": live.get("press_smooth"), + "cpu_temp": live.get("cpu_temp"), + "cpu_offset": live.get("cpu_offset"), + "compensator_k": live.get("compensator_k"), + "rates": { + "temperature_c_per_h": live.get("temp_rate"), + "humidity_pct_per_h": live.get("hum_rate"), + "pressure_hpa_per_h": live.get("press_rate"), + }, + "derived": { + "dew_point": live.get("dew_c"), + "dew_depression": live.get("dew_depression"), + "wet_bulb": live.get("wet_bulb"), + "vpd_hpa": live.get("vpd"), + "absolute_humidity_g_m3": live.get("abs_humidity"), + "heat_index": live.get("heat_index"), + "cloud_index": live.get("cloud_index"), + "solar_elevation": live.get("solar_elevation"), + "solar_azimuth": live.get("solar_azimuth"), + "clear_sky_wm2": live.get("clear_sky_wm2"), + }, + "health": live.get("health"), + "novelty_d2": live.get("novelty_d2"), + "simulated": live.get("simulated"), + }) + + +@app.get("/api/history") +def history(hours: float = Query(6.0, gt=0, le=24 * 90), + max_points: int = Query(720, ge=10, le=5000)) -> Dict: + st = _st() + cols = ["ts", "temp_smooth", "hum_smooth", "press_slp", "dew_c", + "temp_rate", "press_rate", "lux"] + w = st.store.window(hours, cols) + n = w["ts"].size + if n == 0: + return {"n": 0, "series": {}} + stride = max(1, n // max_points) + out = {c: w[c][::stride] for c in cols} + return _clean({ + "n": int(out["ts"].size), + "hours": hours, + "series": { + "ts": out["ts"].tolist(), + "temperature": out["temp_smooth"].tolist(), + "humidity": out["hum_smooth"].tolist(), + "pressure": out["press_slp"].tolist(), + "dew_point": out["dew_c"].tolist(), + "temperature_rate": out["temp_rate"].tolist(), + "pressure_rate": out["press_rate"].tolist(), + "lux": out["lux"].tolist(), + }, + }) + + +@app.get("/api/history/range") +def history_range(start: Optional[float] = None, end: Optional[float] = None, + hours: Optional[float] = None, + bucket: Optional[int] = Query(None, ge=30, le=604800)) -> Dict: + """Bucket-aggregated telemetry for an arbitrary window. + + Accepts either an explicit epoch `start`/`end` pair or a trailing + `hours` span. The bucket is chosen automatically from the span unless + you pin it, so a request for a year does not try to serialise a year + of five-minute rows to a browser. + """ + st = _st() + now = time.time() + if hours is not None: + start, end = now - hours * 3600.0, now + if start is None or end is None: + raise HTTPException(422, "provide start and end, or hours") + if end - start > 366 * 86400: + raise HTTPException(422, "range limited to one year") + data = st.store.range_series(start, end, bucket) + return _clean(data) + + +@app.get("/api/history/daily") +def history_daily(days: int = Query(30, ge=1, le=400)) -> Dict: + st = _st() + end = time.time() + start = end - days * 86400.0 + return _clean({"days": st.store.daily_summary(start, end)}) + + +@app.get("/api/records") +def records() -> Dict: + """All-time extremes held by this station, each with its timestamp.""" + return _clean(_st().store.extremes()) + + +@app.get("/api/storage") +def storage_stats() -> Dict: + """Rows per resolution tier plus database size, so retention is visible.""" + st = _st() + return _clean({ + **st.store.storage_stats(), + "policy": { + "raw_retention_days": CONFIG.storage.raw_retention_days, + "five_min_retention_days": CONFIG.storage.five_min_retention_days, + "note": "Nothing is deleted, only downsampled. Rows older than the raw " + "window fold into 5-minute means, then into hourly means. A " + "year of history lands around 30 MB.", + }, + }) + + +@app.get("/api/export.csv") +def export_csv(start: Optional[float] = None, end: Optional[float] = None, + hours: Optional[float] = None): + st = _st() + now = time.time() + if hours is not None: + start, end = now - hours * 3600.0, now + if start is None or end is None: + raise HTTPException(422, "provide start and end, or hours") + stamp = time.strftime("%Y%m%d-%H%M", time.localtime(start)) + return StreamingResponse( + st.store.iter_csv(start, end), + media_type="text/csv", + headers={"Content-Disposition": + f'attachment; filename="ashvale-{stamp}.csv"'}, + ) + + +@app.get("/api/methods") +def methods_doc() -> Dict: + """The Methods tab is generated from this, so it cannot drift from the code.""" + return _clean(describe(CONFIG)) + + +@app.get("/api/forecast") +def forecast(target: Optional[str] = None, refresh: bool = False) -> Dict: + st = _st() + if refresh or not st.forecast_bundle: + st.refresh_forecasts() + # A cold station has no forecast yet. Return the empty shape rather than + # a bare {}, so a client never has to distinguish "no data" from "no key". + bundle = dict(st.forecast_bundle) or { + "issued_ts": None, "anchors": {}, + "targets": {t: [] for t in CONFIG.model.targets}, + "warming_up": True, + } + if target: + if target not in bundle.get("targets", {}): + raise HTTPException(404, f"unknown target '{target}'") + bundle["targets"] = {target: bundle["targets"][target]} + return _clean(bundle) + + +@app.get("/api/outlook") +def outlook() -> Dict: + """Days 2 to 7. Climatology plus a decaying anomaly, honestly labelled.""" + st = _st() + if not st.outlook_bundle: + st.refresh_forecasts() + base = st.outlook_bundle or { + "issued_ts": None, "ready": False, "annual_terms": False, + "history_days": round(st.store.span_days(), 2), + "targets": {t: [] for t in CONFIG.model.targets}, + } + return _clean({ + **base, + "method": "harmonic climatology with exponentially decaying anomaly", + "caveat": "A single point sensor cannot observe approaching systems. " + "Treat days 2 to 7 as a climatological outlook, not a forecast.", + }) + + +@app.get("/api/precipitation") +def precipitation() -> Dict: + st = _st() + return _clean(st.precip_bundle or {}) + + +@app.get("/api/anomaly") +def anomaly() -> Dict: + st = _st() + return _clean({ + **(st.anomaly_bundle or {}), + "events": st.monitor.recent(20), + }) + + +@app.get("/api/models") +def models() -> Dict: + st = _st() + return _clean({ + "nowcast": st.nowcast.diagnostics(), + "climatology": { + "ready": st.climatology.ready, + "annual_terms": st.climatology.use_annual, + "history_days": round(st.climatology.n_days, 2), + "residual_std": st.climatology.resid_std, + }, + "precipitation": { + "coefficients": st.precip.coefficients(), + "strong_labels": st.precip.n_strong, + "weak_labels": st.precip.n_weak, + "logloss_ewma": st.precip.ewma_logloss, + }, + "calibration": st.tracker.compensator.to_dict(), + }) + + +@app.get("/api/scorecard") +def scorecard() -> Dict: + st = _st() + rows = st.store.scorecard() + return _clean({ + "rows": rows, + "explainer": "skill = 1 - MAE/MAE_persistence. Above zero means the " + "model beats 'nothing changes'. Below zero means it does not, " + "and persistence should be shipped instead.", + }) + + +@app.post("/api/verify") +def verify_now() -> Dict: + return _clean(_st().verify()) + + +@app.post("/api/train") +def train_now(hours: float = Query(24 * 30, gt=1)) -> Dict: + return _clean(_st().train(hours)) + + +@app.post("/api/label") +def add_label(body: LabelIn) -> Dict: + return _clean(_st().add_label(body.kind, body.value, body.ts, body.note)) + + +@app.post("/api/calibrate") +def calibrate(body: CalibrationIn) -> Dict: + st = _st() + if body.reset: + return _clean(st.reset_calibration()) + if body.reference_c is None: + raise HTTPException(422, "provide reference_c, or reset=true") + result = st.calibrate_temperature(body.reference_c) + if "error" in result: + raise HTTPException(409, result["error"]) + return _clean(result) + + +@app.get("/api/status") +def status() -> Dict: + st = _st() + return _clean({ + **st.status(), + "display_frame": display.frame_name if display else None, + "events": st.store.recent_events(15), + }) + + +@app.get("/api/events") +def events(limit: int = Query(50, ge=1, le=500)) -> List[Dict]: + return _clean(_st().store.recent_events(limit)) + + +@app.get("/api/stream") +async def stream(): + """Server-sent events. One connection instead of a poll every 2 seconds, + which on a Zero 2 W is the difference between 4% and 0.4% CPU.""" + async def gen(): + while True: + st = _st() + payload = { + "telemetry": telemetry(), + "precipitation": _clean(st.precip_bundle or {}), + "health": st.monitor.health.overall, + "drift_stress": round(st.monitor.drift.stress, 3), + } + yield f"data: {json.dumps(payload)}\n\n" + await asyncio.sleep(2.0) + + return StreamingResponse(gen(), media_type="text/event-stream", + headers={"Cache-Control": "no-cache", + "X-Accel-Buffering": "no"}) + + +@app.get("/", response_class=HTMLResponse) +def dashboard() -> str: + return DASHBOARD_HTML diff --git a/ashvale/config.py b/ashvale/config.py new file mode 100644 index 0000000..ca6f91f --- /dev/null +++ b/ashvale/config.py @@ -0,0 +1,150 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Configuration for the Ashvale station. + +Everything tunable lives here. Override any field with a YAML file +(default `config.yaml` next to the repo root) or with environment +variables prefixed `ASHVALE_` (e.g. `ASHVALE_SITE__ALTITUDE_M=42`). +""" + +from __future__ import annotations + +import os +from dataclasses import dataclass, field, fields, is_dataclass +from pathlib import Path +from typing import Any, Dict + +try: + import yaml # optional +except Exception: # pragma: no cover + yaml = None + +REPO_ROOT = Path(__file__).resolve().parent.parent + + +@dataclass +class SiteConfig: + name: str = "ashvale-labs-weather-station" + latitude: float = 52.2053 # Cambridge, UK + longitude: float = 0.1218 + altitude_m: float = 15.0 # for sea-level pressure reduction + timezone: str = "Europe/London" + indoors: bool = True # honest flag, changes how forecasts are worded + + +@dataclass +class SensorConfig: + sample_period_s: float = 2.0 # how often we read the HAT + persist_period_s: float = 30.0 # how often a row hits the database + rotation_deg: int = 90 + low_light: bool = True + tcs3400_addr: int = 0x39 + # CPU self-heating compensation: T_true = T_sensor - k * (T_cpu - T_sensor) + cpu_heat_k: float = 0.55 + cpu_heat_k_min: float = 0.15 + cpu_heat_k_max: float = 1.20 + # Kalman process/measurement noise (per-signal) + kalman_q_temp: float = 2.0e-6 + kalman_r_temp: float = 0.02 + kalman_q_press: float = 1.0e-5 + kalman_r_press: float = 0.05 + kalman_q_hum: float = 5.0e-5 + kalman_r_hum: float = 0.60 + + +@dataclass +class ModelConfig: + grid_s: int = 300 # 5-minute feature grid + horizons_s: tuple = (900, 3600, 10800, 21600, 43200, 86400) + targets: tuple = ("temperature", "humidity", "pressure") + rls_forgetting: float = 0.9985 # lambda, ~ 11h memory at 5 min + rls_delta: float = 100.0 # P0 = delta * I + conformal_window: int = 400 # residuals kept per head + conformal_alpha: float = 0.10 # 90% intervals + conformal_gamma: float = 0.01 # adaptive conformal step + train_period_s: float = 600.0 # retrain cadence + min_rows_to_train: int = 120 + climatology_min_days_annual: float = 120.0 + anomaly_ewma_lambda: float = 0.15 + anomaly_threshold: float = 12.0 # Mahalanobis^2 alarm level + drift_delta: float = 0.05 + drift_lambda: float = 8.0 + + +@dataclass +class StorageConfig: + db_path: str = str(REPO_ROOT / "data" / "ashvale.db") + state_dir: str = str(REPO_ROOT / "data" / "state") + raw_retention_days: float = 7.0 + five_min_retention_days: float = 90.0 + vacuum_period_s: float = 86400.0 + + +@dataclass +class ServerConfig: + host: str = "0.0.0.0" + port: int = 8000 + led_enabled: bool = True + led_cycle_s: float = 0.4 + + +@dataclass +class Config: + site: SiteConfig = field(default_factory=SiteConfig) + sensor: SensorConfig = field(default_factory=SensorConfig) + model: ModelConfig = field(default_factory=ModelConfig) + storage: StorageConfig = field(default_factory=StorageConfig) + server: ServerConfig = field(default_factory=ServerConfig) + + +def _apply(obj: Any, patch: Dict[str, Any]) -> None: + for key, value in (patch or {}).items(): + if not hasattr(obj, key): + continue + current = getattr(obj, key) + if is_dataclass(current) and isinstance(value, dict): + _apply(current, value) + else: + setattr(obj, key, type(current)(value) if current is not None else value) + + +def _apply_env(obj: Any, prefix: str = "ASHVALE_") -> None: + for f in fields(obj): + current = getattr(obj, f.name) + if is_dataclass(current): + _apply_env(current, f"{prefix}{f.name.upper()}__") + continue + env_key = f"{prefix}{f.name.upper()}" + if env_key in os.environ: + raw = os.environ[env_key] + try: + setattr(obj, f.name, type(current)(raw)) + except Exception: + setattr(obj, f.name, raw) + + +def load_config(path: str | os.PathLike | None = None) -> Config: + cfg = Config() + candidate = Path(path) if path else REPO_ROOT / "config.yaml" + if candidate.exists() and yaml is not None: + with open(candidate, "r", encoding="utf-8") as fh: + _apply(cfg, yaml.safe_load(fh) or {}) + _apply_env(cfg) + Path(cfg.storage.db_path).parent.mkdir(parents=True, exist_ok=True) + Path(cfg.storage.state_dir).mkdir(parents=True, exist_ok=True) + return cfg + + +CONFIG = load_config() diff --git a/ashvale/dashboard.py b/ashvale/dashboard.py new file mode 100644 index 0000000..12dbe4b --- /dev/null +++ b/ashvale/dashboard.py @@ -0,0 +1,976 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""The dashboard: five tabs, one viewport, no scrolling. + +Layout contract. The page is a fixed three-row grid pinned to the +viewport height: header, tab bar, then a content region that takes the +remaining space and never overflows the fold. Each tab lays its panels +out on an internal grid sized in fractions of that region, so nothing +depends on content height. Where a panel genuinely holds more than fits +(the daily records table, the methods prose) that individual panel +scrolls internally while the page frame stays put. Below 1024 px the +constraint is released, because pinning five panels into a phone +viewport produces unreadable eight-pixel type, and a phone user expects +to scroll anyway. + +Visual language carries over unchanged from the previous station page: +slate-950 ground, glass panels, Jakarta for prose and JetBrains Mono for +anything numeric. The one new structural device is the tab bar, and it +earns its place. Five distinct questions (what is it doing, what will it +do, what did it do, is the model any good, how does it work) were +previously one long scroll where the important things sat below the +fold. + +The signature element is the estimator internals panel on the Live tab. +Most weather dashboards show numbers. This one shows the state estimator +working: self-heating coefficient, Kalman innovation, novelty distance +and drift pressure, all ticking at 2 Hz. It is the part of the system +that is normally invisible, and watching a filter converge is the most +honest possible demonstration that there is real machinery underneath. +""" + +DASHBOARD_HTML = r""" + + + + + +Ashvale Station + + + + + + + + + + +
+
+
+
+
+ +
+ +
+
+
+ + + +
+
+

Ashvale Station

+

+ - · 0 d logged · k=- +

+
+
+
+ health + + + --:--:-- + +
+
+ + + +
+ + +
+ +
+
+ Temperature + KALMAN +
+
--°C
+
+
rate--
+
raw / cpu--
+
+ +
+
+ +
+
+ Humidity + HTS221 +
+
--%
+
+
dew point--
+
depression--
+
+ +
+
+ +
+
+ Barometer + MSL +
+
--hPa
+
+
tendency--
+
character--
+
+ +
+
+ +
+
+ Sky +
+
+
--clr
+
+
cloud index--
+
sun / cct--
+
+ +
+
+ +
+
+
+

Rolling window

+ 2 s stream +
+
+ + + +
+
+
+
+ +
+
+

Estimator internals

+

what the filter is doing right now

+
+
+
+
self-heating k--
+
cpu offset--
+
Removes the SoC bias. Calibrate it on the Models tab.
+
+
+
novelty d²--
+
+
drift pressure--
+
+
Novelty is a multivariate departure from the recent norm. Drift reaching 100% queues a retrain.
+
+
+
+
+ +
+
wet bulb
--
+
vpd
--
+
abs hum
--
+
heat idx
--
+
solar el
--
+
pitch
--
+
roll
--
+
yaw
--
+
compass
--
+
accel z
--
+
+
+ + +
+
+
+
+

Observed and forecast

+

shaded band is the 90% conformal interval

+
+
+ + +
+
+
+
+
+ +
+
+

Conditions ahead

+

Zambretti prior + online logistic

+
+
+
+
--
+
Z=- · -
+
+
+
rain probability--
+
+
prior -learner -trust -
+
+
+
Was it wet in the last hour?
+
+ + +
+
0 confirmed observations
+
+
+
+ pressure, last 24 h + -- +
+
+

The only signal here that sees past your walls. Its slope, not its level, is what drives the forecast above.

+
+
+
+ +
+
+
+

Seven day outlook

+

climatology plus decaying anomaly, not a synoptic forecast

+
+ warming up +
+
+
+
+ + +
+
+
+ + + + + + +
+
+ + + +
+
+
+ Export CSV +
+

-

+
+ +
+
+
+ +
+
+

Records

+
+ + +
+
+
+ +
+
+ + +
+
+
+
+

Verification scorecard

+

skill above zero means it beats persistence

+
+
+ + +
+
+
+ + + + + + + + + +
targetleadMAEpersistskillcovernp/c/l
+

No matured forecasts yet. Rows appear as each horizon reaches its validity time: 15 minutes first, 24 hours tomorrow. The p/c/l column is the ensemble weight on persistence, climatology and the learned model.

+
+
+ +
+

Calibration and state

+
+
+
Trusted thermometer reading
+
+ + + +
+
Recursive least squares on the self-heating coefficient. One good reading is enough.
+
+
+
Storage tiers
+
+
+
+
Precipitation coefficients
+
+
+
+
+ +
+

Station log

+
+
+
+ + +
+
+
+

How it is wired

+

select a stage to read its rationale

+
+
+
+
+
+
+
+
+ +
+
+ + + + +""" diff --git a/ashvale/estimation.py b/ashvale/estimation.py new file mode 100644 index 0000000..3498509 --- /dev/null +++ b/ashvale/estimation.py @@ -0,0 +1,225 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""State estimation: the layer between a noisy sensor and an honest number. + +Two jobs here, both familiar from soft-sensor work: + +1. `ThermalCompensator` removes the SoC self-heating bias. The classic + Sense HAT correction `T = T_sensor - k (T_cpu - T_sensor)` is a + one-parameter grey-box model. We keep the structure and estimate `k` + recursively whenever a trusted reference reading is supplied, which + beats hard-coding 1/1.5 and hoping. + +2. `SignalTracker` runs a constant-velocity Kalman filter per signal. + The filtered level is a denoised measurement; the filtered rate is the + thing you actually want for weather. A finite difference of a 0.05 hPa + noise floor over 5 minutes is garbage. A Kalman rate is not. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field +from typing import Dict, Optional + +import numpy as np + + +@dataclass +class KalmanCV: + """Constant-velocity Kalman filter for one scalar signal. + + State x = [level, rate]. Process noise is the standard continuous + white-noise-acceleration model, so `q` has units of (signal/s^2)^2/s + and is the only knob that matters: raise it to track faster, lower it + to smooth harder. + """ + + q: float + r: float + x: np.ndarray = field(default_factory=lambda: np.zeros(2)) + P: np.ndarray = field(default_factory=lambda: np.eye(2) * 1e3) + initialised: bool = False + nis: float = 0.0 # normalised innovation squared, for health monitoring + + def update(self, z: float, dt: float) -> tuple[float, float]: + if not np.isfinite(z): + return float(self.x[0]), float(self.x[1]) + if not self.initialised: + self.x = np.array([z, 0.0]) + self.P = np.array([[self.r, 0.0], [0.0, 1e-4]]) + self.initialised = True + return z, 0.0 + + dt = float(max(min(dt, 3600.0), 1e-3)) + F = np.array([[1.0, dt], [0.0, 1.0]]) + Q = self.q * np.array([[dt ** 3 / 3.0, dt ** 2 / 2.0], + [dt ** 2 / 2.0, dt]]) + + # predict + self.x = F @ self.x + self.P = F @ self.P @ F.T + Q + + # update + H = np.array([[1.0, 0.0]]) + y = float(z) - float((H @ self.x)[0]) + S = float((H @ self.P @ H.T)[0, 0]) + self.r + K = (self.P @ H.T) / S + self.x = self.x + (K.flatten() * y) + I_KH = np.eye(2) - K @ H + self.P = I_KH @ self.P @ I_KH.T + K @ K.T * self.r # Joseph form, stays PSD + + self.nis = (y * y) / S + return float(self.x[0]), float(self.x[1]) + + @property + def level(self) -> float: + return float(self.x[0]) + + @property + def rate(self) -> float: + """Signal units per second.""" + return float(self.x[1]) + + def to_dict(self) -> Dict: + return {"q": self.q, "r": self.r, "x": self.x.tolist(), + "P": self.P.tolist(), "initialised": self.initialised} + + @classmethod + def from_dict(cls, d: Dict) -> "KalmanCV": + kf = cls(q=d["q"], r=d["r"]) + kf.x = np.array(d["x"], dtype=float) + kf.P = np.array(d["P"], dtype=float) + kf.initialised = bool(d["initialised"]) + return kf + + +class ThermalCompensator: + """Grey-box removal of SoC self-heating. + + Model: T_true = T_sensor - k * (T_cpu - T_sensor), k >= 0. + + `k` is updated by recursive least squares whenever `calibrate()` is + called with a trusted reference temperature (a mercury thermometer, a + second logger, or a nearby METAR reading). Until then the configured + prior is used and clamped to a physically sane band, because a runaway + `k` produces confident nonsense, which is worse than a mild bias. + """ + + def __init__(self, k0: float = 0.55, k_min: float = 0.15, k_max: float = 1.2, + forgetting: float = 0.98): + self.k = float(k0) + self.k_min, self.k_max = float(k_min), float(k_max) + self.P = 10.0 + self.lam = float(forgetting) + self.n_calibrations = 0 + self.last_residual = 0.0 + + def compensate(self, t_sensor: float, t_cpu: float) -> float: + if not (np.isfinite(t_sensor) and np.isfinite(t_cpu)): + return float(t_sensor) + delta = max(t_cpu - t_sensor, 0.0) + return float(t_sensor - self.k * delta) + + def calibrate(self, t_sensor: float, t_cpu: float, t_reference: float) -> Dict: + """One RLS step on k. Regressor is the CPU/sensor gradient.""" + phi = max(t_cpu - t_sensor, 0.0) + target = t_sensor - t_reference # what k*phi should equal + denom = self.lam + phi * self.P * phi + gain = (self.P * phi) / denom if denom > 1e-12 else 0.0 + residual = target - self.k * phi + self.k = float(np.clip(self.k + gain * residual, self.k_min, self.k_max)) + self.P = float((self.P - gain * phi * self.P) / self.lam) + self.P = float(np.clip(self.P, 1e-6, 1e4)) + self.n_calibrations += 1 + self.last_residual = float(residual) + return {"k": self.k, "residual": self.last_residual, "n": self.n_calibrations} + + def to_dict(self) -> Dict: + return {"k": self.k, "P": self.P, "lam": self.lam, "k_min": self.k_min, + "k_max": self.k_max, "n": self.n_calibrations} + + @classmethod + def from_dict(cls, d: Dict) -> "ThermalCompensator": + tc = cls(d["k"], d["k_min"], d["k_max"], d["lam"]) + tc.P = d["P"] + tc.n_calibrations = d.get("n", 0) + return tc + + +class SignalTracker: + """Bank of Kalman filters plus the compensator, driven at sample rate.""" + + def __init__(self, cfg): + self.compensator = ThermalCompensator( + cfg.sensor.cpu_heat_k, cfg.sensor.cpu_heat_k_min, + cfg.sensor.cpu_heat_k_max, + ) + self.filters = { + "temperature": KalmanCV(cfg.sensor.kalman_q_temp, cfg.sensor.kalman_r_temp), + "humidity": KalmanCV(cfg.sensor.kalman_q_hum, cfg.sensor.kalman_r_hum), + "pressure": KalmanCV(cfg.sensor.kalman_q_press, cfg.sensor.kalman_r_press), + } + self.last_ts: Optional[float] = None + + def step(self, ts: float, temp_raw: float, hum: float, press: float, + cpu_temp: float) -> Dict[str, float]: + dt = (ts - self.last_ts) if self.last_ts is not None else 1.0 + self.last_ts = ts + + temp_c = self.compensator.compensate(temp_raw, cpu_temp) + t_lvl, t_rate = self.filters["temperature"].update(temp_c, dt) + h_lvl, h_rate = self.filters["humidity"].update(hum, dt) + p_lvl, p_rate = self.filters["pressure"].update(press, dt) + + return { + "temp_c": temp_c, + "temp_smooth": t_lvl, + "temp_rate": t_rate * 3600.0, # C per hour + "hum_smooth": h_lvl, + "hum_rate": h_rate * 3600.0, # % per hour + "press_smooth": p_lvl, + "press_rate": p_rate * 3600.0, # hPa per hour, the forecaster's gold + "nis_temp": self.filters["temperature"].nis, + "nis_press": self.filters["pressure"].nis, + } + + def to_dict(self) -> Dict: + return { + "compensator": self.compensator.to_dict(), + "filters": {k: v.to_dict() for k, v in self.filters.items()}, + "last_ts": self.last_ts, + } + + def load_dict(self, d: Dict) -> None: + self.compensator = ThermalCompensator.from_dict(d["compensator"]) + self.filters = {k: KalmanCV.from_dict(v) for k, v in d["filters"].items()} + self.last_ts = d.get("last_ts") + + +def stuck_sensor_score(values: np.ndarray, window: int = 60) -> float: + """Fraction of the last `window` samples that are bit-identical. + + An HTS221 that latches is the quietest failure mode there is: the + dashboard looks perfect, the model trains happily, and every forecast + is confidently wrong. This is the cheapest possible smoke alarm. + """ + if values.size < 5: + return 0.0 + tail = values[-window:] + tail = tail[np.isfinite(tail)] + if tail.size < 5: + return 0.0 + return float(np.mean(np.abs(np.diff(tail)) < 1e-9)) diff --git a/ashvale/features.py b/ashvale/features.py new file mode 100644 index 0000000..9ba0769 --- /dev/null +++ b/ashvale/features.py @@ -0,0 +1,204 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Feature engineering, pure numpy, no pandas. + +Design rules used here: + +* Anything derivable from physics is computed, not learned. +* Anything periodic is encoded as sin/cos pairs so a linear model can + represent phase without a discontinuity at midnight. +* Every lag is expressed in *hours*, not samples, so changing `grid_s` + does not silently change what the model means by `three hours ago`. +* Targets are predicted as *deltas from now*, never as absolute levels. + A model that must output 14.7 C spends all its capacity on the mean; + a model that outputs +0.4 C spends it on the weather. +""" + +from __future__ import annotations + +from typing import Dict, List, Tuple + +import numpy as np + +from .physics import (absolute_humidity, clear_sky_irradiance, dew_point, + solar_position, vapour_pressure_deficit, wet_bulb) + +FEATURE_NAMES: List[str] = [ + "bias", + "temp", "temp_rate_1h", "temp_rate_3h", "temp_std_3h", "temp_dev_24h", + "hum", "hum_rate_1h", "hum_rate_3h", "hum_std_3h", + "press_anom", "press_tend_1h", "press_tend_3h", "press_tend_6h", "press_std_6h", + "dewpoint", "dewpoint_depression", "vpd", "abs_hum", "wet_bulb", + "log_lux", "cloud_index", "solar_elev", "solar_elev_pos", "is_day", + "sin_h1", "cos_h1", "sin_h2", "cos_h2", "sin_doy", "cos_doy", + "press_x_hum", "tend_x_dewdep", +] + +N_FEATURES = len(FEATURE_NAMES) + + +def _shift(a: np.ndarray, k: int) -> np.ndarray: + """a[i - k], NaN-padded at the front.""" + out = np.full_like(a, np.nan, dtype=float) + if k <= 0: + return a.copy() + if k < a.size: + out[k:] = a[:-k] + return out + + +def _rolling(a: np.ndarray, win: int, fn) -> np.ndarray: + """Trailing rolling statistic. O(n*win) but win is small and n is a day.""" + out = np.full(a.size, np.nan, dtype=float) + if a.size == 0: + return out + win = max(int(win), 1) + for i in range(a.size): + lo = max(0, i - win + 1) + seg = a[lo:i + 1] + seg = seg[np.isfinite(seg)] + if seg.size >= max(2, win // 3): + out[i] = fn(seg) + return out + + +def build_features(grid_ts: np.ndarray, temp: np.ndarray, hum: np.ndarray, + press_slp: np.ndarray, lux: np.ndarray, + grid_s: int, latitude: float, longitude: float + ) -> Tuple[np.ndarray, np.ndarray]: + """Return (X of shape (n, N_FEATURES), valid mask of shape (n,)).""" + n = grid_ts.size + if n == 0: + return np.zeros((0, N_FEATURES)), np.zeros(0, dtype=bool) + + per_hour = max(int(round(3600 / grid_s)), 1) + + def rate(a: np.ndarray, hours: int) -> np.ndarray: + return (a - _shift(a, hours * per_hour)) / float(hours) + + temp_rate_1h = rate(temp, 1) + temp_rate_3h = rate(temp, 3) + temp_std_3h = _rolling(temp, 3 * per_hour, np.std) + temp_mean_24h = _rolling(temp, 24 * per_hour, np.mean) + temp_dev_24h = temp - temp_mean_24h + + hum_rate_1h = rate(hum, 1) + hum_rate_3h = rate(hum, 3) + hum_std_3h = _rolling(hum, 3 * per_hour, np.std) + + press_anom = press_slp - 1013.25 + press_tend_1h = rate(press_slp, 1) + press_tend_3h = rate(press_slp, 3) + press_tend_6h = rate(press_slp, 6) + press_std_6h = _rolling(press_slp, 6 * per_hour, np.std) + + dp = dew_point(temp, hum) + dep = temp - dp + vpd = vapour_pressure_deficit(temp, hum) + ah = absolute_humidity(temp, hum) + wb = wet_bulb(temp, hum) + + elev, _ = solar_position(grid_ts, latitude, longitude) + elev = np.atleast_1d(elev) + expected = clear_sky_irradiance(elev) + log_lux = np.log1p(np.clip(lux, 0.0, None)) + # cloud index: 1 = overcast, 0 = clear. Only meaningful in daylight. + scale = np.maximum(expected, 1.0) * 45.0 # crude lux-per-W/m^2 for daylight + cloud = np.where(elev > 5.0, np.clip(1.0 - np.clip(lux, 0, None) / scale, 0.0, 1.0), 0.5) + + hour = (grid_ts % 86400.0) / 86400.0 + doy = (grid_ts % 31557600.0) / 31557600.0 + + X = np.column_stack([ + np.ones(n), + temp, temp_rate_1h, temp_rate_3h, temp_std_3h, temp_dev_24h, + hum, hum_rate_1h, hum_rate_3h, hum_std_3h, + press_anom, press_tend_1h, press_tend_3h, press_tend_6h, press_std_6h, + dp, dep, vpd, ah, wb, + log_lux, cloud, elev, np.clip(elev, 0.0, None), (elev > 0.0).astype(float), + np.sin(2 * np.pi * hour), np.cos(2 * np.pi * hour), + np.sin(4 * np.pi * hour), np.cos(4 * np.pi * hour), + np.sin(2 * np.pi * doy), np.cos(2 * np.pi * doy), + press_anom * (hum - 70.0) / 100.0, + press_tend_3h * dep, + ]) + + assert X.shape[1] == N_FEATURES, f"feature count drift: {X.shape[1]} vs {N_FEATURES}" + valid = np.all(np.isfinite(X), axis=1) + X = np.nan_to_num(X, nan=0.0, posinf=0.0, neginf=0.0) + return X, valid + + +class Standardiser: + """Streaming z-scoring with Welford moments. + + Recursive least squares is scale-sensitive: an unscaled `pressure` at + 1013 and an unscaled `temp_rate` at 0.02 give a condition number that + will embarrass you. Standardising online keeps P well-conditioned + without a second pass over history. + """ + + def __init__(self, n_features: int = N_FEATURES): + self.n = 0 + self.mean = np.zeros(n_features) + self.m2 = np.ones(n_features) + + def partial_fit(self, X: np.ndarray) -> None: + for row in np.atleast_2d(X): + self.n += 1 + delta = row - self.mean + self.mean += delta / self.n + self.m2 += delta * (row - self.mean) + + def transform(self, X: np.ndarray) -> np.ndarray: + if self.n < 2: + return np.atleast_2d(X) + std = np.sqrt(self.m2 / max(self.n - 1, 1)) + std = np.where(std < 1e-8, 1.0, std) + out = (np.atleast_2d(X) - self.mean) / std + out[:, 0] = 1.0 # keep the bias column intact + return out + + def fit_transform(self, X: np.ndarray) -> np.ndarray: + self.partial_fit(X) + return self.transform(X) + + def to_dict(self) -> Dict: + return {"n": self.n, "mean": self.mean.tolist(), "m2": self.m2.tolist()} + + @classmethod + def from_dict(cls, d: Dict) -> "Standardiser": + s = cls(len(d["mean"])) + s.n = d["n"] + s.mean = np.array(d["mean"], dtype=float) + s.m2 = np.array(d["m2"], dtype=float) + return s + + +def supervised_pairs(X: np.ndarray, valid: np.ndarray, y: np.ndarray, + horizon_steps: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Align features at t with the *change* in y between t and t+h. + + Returns (X_aligned, delta_y, anchor_y) so the caller can reconstruct + the absolute forecast as anchor + predicted delta. + """ + n = X.shape[0] + if n <= horizon_steps: + return np.zeros((0, X.shape[1])), np.zeros(0), np.zeros(0) + Xa = X[:n - horizon_steps] + anchor = y[:n - horizon_steps] + future = y[horizon_steps:] + mask = (valid[:n - horizon_steps] & np.isfinite(future) & np.isfinite(anchor)) + return Xa[mask], (future - anchor)[mask], anchor[mask] diff --git a/ashvale/led.py b/ashvale/led.py new file mode 100644 index 0000000..757500e --- /dev/null +++ b/ashvale/led.py @@ -0,0 +1,272 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""The 8x8 matrix as a forecast instrument, not a scrolling number. + +Text on eight pixels is slow and, worse, it makes you wait for the one +value you wanted. So the display cycles through *glyphs* that are +readable at a glance from across a room: + + temperature scrolled with a heat-mapped colour, as before + humidity scrolled with a moisture-band colour + pressure a trend arrow whose colour encodes the Zambretti class + and whose brightness encodes tendency magnitude + rain a filled column bar, 0 to 8 pixels, of rain probability + forecast a 3-hour temperature delta as a rising or falling wedge + alert a red pulse if a sensor is faulted or drift fired + +Design constraint: never call `show_message` while an alert is pending, +because a 6-second scroll is a 6-second delay on the only frame that +matters. +""" + +from __future__ import annotations + +import asyncio +import time +from typing import Dict, List, Sequence, Tuple + +OFF = (0, 0, 0) + + +def temp_colour(temp_c: float) -> List[int]: + if temp_c <= 15.0: + return [0, 150, 255] + if temp_c <= 21.0: + return [0, 255, 180] + if temp_c <= 25.0: + return [70, 255, 0] + if temp_c <= 28.0: + return [255, 190, 0] + if temp_c <= 32.0: + return [255, 90, 0] + return [255, 20, 20] + + +def humidity_colour(rh: float) -> List[int]: + if rh < 35.0: + return [255, 180, 50] + if rh <= 60.0: + return [0, 210, 255] + return [0, 100, 255] + + +CONDITION_COLOUR = { + "settled": (0, 220, 140), "fine": (90, 230, 60), "fair": (200, 230, 40), + "changeable": (255, 190, 0), "unsettled": (255, 120, 0), + "rain": (0, 140, 255), "wet": (0, 90, 255), "stormy": (255, 40, 60), +} + +# 8x8 bitmaps: '#' is lit, anything else is off + + +def _mask(rows: Sequence[str]) -> List[List[int]]: + return [[1 if ch == "#" else 0 for ch in row.ljust(8, ".")[:8]] for row in rows] + + +ARROW_UP = _mask([ + "...##...", + "..####..", + ".##..##.", + "##.##.##", + "...##...", + "...##...", + "...##...", + "...##...", +]) + +ARROW_DOWN = _mask([ + "...##...", + "...##...", + "...##...", + "...##...", + "##.##.##", + ".##..##.", + "..####..", + "...##...", +]) + +ARROW_FLAT = _mask([ + "........", + "........", + "....#...", + "########", + "########", + "....#...", + "........", + "........", +]) + +DROP = _mask([ + "...##...", + "...##...", + "..####..", + ".######.", + "########", + "########", + ".######.", + "..####..", +]) + +BANG = _mask([ + "...##...", + "...##...", + "...##...", + "...##...", + "...##...", + "........", + "...##...", + "...##...", +]) + + +def render(mask: List[List[int]], colour: Tuple[int, int, int], + dim: float = 1.0) -> List[Tuple[int, int, int]]: + c = tuple(int(max(0, min(255, v * dim))) for v in colour) + return [c if cell else OFF for row in mask for cell in row] + + +def bar(fraction: float, colour: Tuple[int, int, int], + background: Tuple[int, int, int] = (12, 12, 20)) -> List[Tuple[int, int, int]]: + """Bottom-up column bar across the full 8x8, 1/64 resolution.""" + lit = int(round(max(0.0, min(1.0, fraction)) * 64)) + pixels = [background] * 64 + count = 0 + for row in range(7, -1, -1): + for col in range(8): + if count < lit: + pixels[row * 8 + col] = colour + count += 1 + return pixels + + +class LedDisplay: + """Async display worker. Owns the matrix, reads station state, nothing else.""" + + def __init__(self, station, cycle_s: float = 0.4): + self.station = station + self.cycle_s = float(cycle_s) + self.enabled = True + self._stop = asyncio.Event() + self._task = None + self.frame_name = "idle" + + # ------------------------------------------------------------ frames + + async def _alert_frame(self) -> bool: + health = self.station.monitor.health.overall + drift = self.station.monitor.retrain_requested + if health == "ok" and not drift: + return False + colour = (255, 40, 40) if health == "fault" else (255, 150, 0) + self.frame_name = "alert" + for pulse in (1.0, 0.25, 1.0, 0.25): + self.station.board.set_pixels(render(BANG, colour, pulse)) + await asyncio.sleep(0.22) + self.station.board.clear() + return True + + async def _pressure_frame(self) -> None: + live = self.station.live + precip = self.station.precip_bundle or {} + rate = float(live.get("press_rate", 0.0) or 0.0) + condition = precip.get("condition", "changeable") + colour = CONDITION_COLOUR.get(condition, (200, 200, 200)) + magnitude = min(abs(rate) / 1.2, 1.0) + dim = 0.25 + 0.75 * magnitude + + if rate > 0.15: + mask = ARROW_UP + elif rate < -0.15: + mask = ARROW_DOWN + else: + mask = ARROW_FLAT + self.frame_name = "pressure-trend" + self.station.board.set_pixels(render(mask, colour, dim)) + await asyncio.sleep(2.0) + self.station.board.clear() + + async def _rain_frame(self) -> None: + p = float((self.station.precip_bundle or {}).get("rain_probability", 0.0)) + self.frame_name = "rain-probability" + if p < 0.12: + return + self.station.board.set_pixels(bar(p, (40, 130, 255))) + await asyncio.sleep(1.6) + self.station.board.set_pixels(render(DROP, (40, 130, 255), 0.6 + 0.4 * p)) + await asyncio.sleep(1.0) + self.station.board.clear() + + async def _forecast_frame(self) -> None: + bundle = self.station.forecast_bundle or {} + series = (bundle.get("targets", {}).get("temperature") or []) + target = next((s for s in series if s["horizon_s"] == 10800), None) + if target is None: + return + delta = float(target["delta"]) + self.frame_name = "temp-3h-delta" + colour = (255, 120, 0) if delta > 0 else (0, 170, 255) + mask = ARROW_UP if delta > 0.2 else ARROW_DOWN if delta < -0.2 else ARROW_FLAT + self.station.board.set_pixels(render(mask, colour, 0.35 + min(abs(delta) / 3.0, 0.65))) + await asyncio.sleep(1.6) + self.station.board.clear() + + async def _scroll_frames(self) -> None: + live = self.station.live + temp = live.get("temp_smooth") + hum = live.get("hum_smooth") + press = live.get("press_slp") + if temp is not None: + self.frame_name = "temperature" + self.station.board.show_message(f"{temp:.1f}C", 0.065, temp_colour(temp)) + await asyncio.sleep(self.cycle_s) + if hum is not None: + self.frame_name = "humidity" + self.station.board.show_message(f"{hum:.0f}%", 0.065, humidity_colour(hum)) + await asyncio.sleep(self.cycle_s) + if press is not None: + self.frame_name = "pressure" + self.station.board.show_message(f"{press:.0f}", 0.065, [180, 80, 255]) + await asyncio.sleep(self.cycle_s) + + # -------------------------------------------------------------- loop + + async def run(self) -> None: + while not self._stop.is_set(): + try: + if not self.enabled or not self.station.live: + await asyncio.sleep(1.0) + continue + if await self._alert_frame(): + continue + await self._scroll_frames() + await self._pressure_frame() + await self._forecast_frame() + await self._rain_frame() + except Exception: + await asyncio.sleep(2.0) + + def start(self) -> None: + self._stop.clear() + self._task = asyncio.create_task(self.run()) + + async def stop(self) -> None: + self._stop.set() + if self._task: + self._task.cancel() + try: + await self._task + except (asyncio.CancelledError, Exception): + pass + self.station.board.clear() diff --git a/ashvale/methods.py b/ashvale/methods.py new file mode 100644 index 0000000..c3aa808 --- /dev/null +++ b/ashvale/methods.py @@ -0,0 +1,359 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""A structured account of what this station actually does, and why. + +This module exists so the Methods page in the UI is generated from one +declarative source rather than hand-written HTML that drifts out of date +the first time someone changes a forgetting factor. Every parameter +quoted below is read from the live config at request time, so the page +describes the station you are running, not the one I shipped. + +Each stage records what it consumes, what it produces, the technique, and +crucially a `why` and a `failure` field. The failure mode is the part +that usually goes undocumented and is the part you need at 2 a.m. +""" + +from __future__ import annotations + +from typing import Any, Dict, List + +from .features import FEATURE_NAMES +from .models.nowcast import MEMBERS + + +def pipeline(cfg) -> List[Dict[str, Any]]: + m, s, site = cfg.model, cfg.sensor, cfg.site + horizons = ", ".join(_fmt(h) for h in m.horizons_s) + + return [ + { + "id": "acquire", + "stage": "1", + "title": "Acquisition", + "module": "sensors.py", + "technique": "Direct I2C, plus a stochastic simulator fallback", + "consumes": "HTS221, LPS25HB, LSM9DS1, TCS3400, SoC thermal zone", + "produces": "Raw multi-sensor sample every " + f"{s.sample_period_s:g} s", + "why": "The colour sensor is read over raw smbus rather than through " + "the sense_hat library because the library does not expose the " + "TCS3400 clear channel, which is the one that carries the " + "cloudiness signal.", + "failure": "If the sense_hat import fails the board silently becomes a " + "simulator. The dashboard header says so rather than " + "letting you trust synthetic weather.", + "params": {"sample period": f"{s.sample_period_s:g} s", + "persist period": f"{s.persist_period_s:g} s"}, + }, + { + "id": "compensate", + "stage": "2", + "title": "Self-heating compensation", + "module": "estimation.py", + "technique": "Grey-box model, coefficient by recursive least squares", + "consumes": "T_raw, T_cpu, and any trusted reference you supply", + "produces": "T = T_raw - k (T_cpu - T_raw)", + "why": "The temperature and pressure sensors sit millimetres above a " + "SoC running 20 to 25 C hotter than the room. The usual fix " + "hard-codes k = 1/1.5, but k depends on your case, orientation, " + "airflow and CPU load. Here it is one estimated parameter with a " + "forgetting factor, updated from a single thermometer reading.", + "failure": "A mistyped reference drives k to its clamp and stays there " + "across restarts, because state persists. The reset button " + "on the Models tab exists for exactly that.", + "math": r"k_{t} = k_{t-1} + \frac{P\varphi}{\lambda + \varphi P \varphi}" + r"\left[(T_{raw} - T_{ref}) - k_{t-1}\varphi\right]," + r"\quad \varphi = T_{cpu} - T_{raw}", + "params": {"current k": f"{s.cpu_heat_k:g} (prior)", + "clamp": f"{s.cpu_heat_k_min:g} to {s.cpu_heat_k_max:g}"}, + }, + { + "id": "kalman", + "stage": "3", + "title": "State estimation", + "module": "estimation.py", + "technique": "Constant-velocity Kalman filter per signal, Joseph form", + "consumes": "Compensated temperature, humidity, station pressure", + "produces": "Filtered level and, more importantly, filtered rate", + "why": "Pressure tendency is the single most informative variable a " + "point sensor can offer, and the LPS25HB noise floor makes a " + "naive finite difference pure noise. A CV filter estimates " + "level and rate jointly. The Joseph covariance update is used " + "because the standard form loses positive semi-definiteness " + "over months of continuous running.", + "failure": "Process noise too low and the filter lags real weather; too " + "high and you have an expensive passthrough. The innovation " + "statistic is logged so you can tell which.", + "math": r"x = \begin{bmatrix} \text{level} \\ \text{rate} \end{bmatrix}," + r"\quad Q = q\begin{bmatrix} \Delta t^3/3 & \Delta t^2/2 \\" + r"\Delta t^2/2 & \Delta t \end{bmatrix}", + "params": {"q temperature": f"{s.kalman_q_temp:g}", + "r temperature": f"{s.kalman_r_temp:g}", + "q pressure": f"{s.kalman_q_press:g}"}, + }, + { + "id": "features", + "stage": "4", + "title": "Feature construction", + "module": "features.py, physics.py", + "technique": f"{len(FEATURE_NAMES)} features on a {m.grid_s} s grid, " + "streaming z-scoring by Welford moments", + "consumes": "Resampled history", + "produces": "Design matrix, standardised", + "why": "Three rules. Anything derivable from physics is computed, not " + "learned: dew point, wet bulb, VPD, solar elevation and a " + "clear-sky cloud index are closed-form, so making a learner " + "rediscover the Magnus curve from data wastes both samples and " + "capacity. Anything periodic is encoded as sine and cosine pairs " + "so a linear model can represent phase without a discontinuity " + "at midnight. Every lag is expressed in hours, not samples, so " + "changing the grid does not silently change meaning.", + "failure": "Unstandardised features give a condition number that will " + "embarrass you: pressure sits near 1013 while temperature " + "rate sits near 0.02.", + "params": {"grid": f"{m.grid_s} s", "features": str(len(FEATURE_NAMES)), + "site": f"{site.latitude:.3f}, {site.longitude:.3f} at " + f"{site.altitude_m:g} m"}, + }, + { + "id": "nowcast", + "stage": "5", + "title": "Multi-horizon forecasting", + "module": "models/nowcast.py, models/rls.py", + "technique": f"{len(m.targets) * len(m.horizons_s)} direct heads, " + "exponentially weighted RLS, Hedge-blended", + "consumes": "Design matrix and matured targets", + "produces": f"Forecasts at {horizons} for {', '.join(m.targets)}", + "why": "Direct heads, not one model iterated forward: iterating a " + "one-step model 288 times to reach 24 hours compounds its own " + "bias into a beautifully smooth lie. RLS rather than SGD because " + "a station makes only 288 grid rows a day and RLS is the exact " + "minimiser of the exponentially weighted squared error at every " + "step. Each head predicts a delta from now, never an absolute " + "level, so its capacity goes on the weather instead of the mean.", + "failure": "Plain forgetting inflates the covariance exponentially " + "through quiet nights when the regressor barely moves, and " + "the model then detonates at sunrise. The trace is capped. " + "This is the most common way a field RLS deployment dies.", + "math": r"P_t = \frac{1}{\lambda}\left(P_{t-1} - " + r"\frac{P_{t-1}x x^{\top}P_{t-1}}{\lambda + x^{\top}P_{t-1}x}" + r"\right)", + "params": {"forgetting": f"{m.rls_forgetting:g}", + "effective memory": _memory(m.rls_forgetting, m.grid_s), + "members": ", ".join(MEMBERS)}, + }, + { + "id": "conformal", + "stage": "6", + "title": "Calibrated uncertainty", + "module": "models/rls.py", + "technique": "Adaptive conformal inference", + "consumes": "Realised forecast errors from the verification loop", + "produces": f"{int((1 - m.conformal_alpha) * 100)}% prediction intervals", + "why": "Split conformal is only valid under exchangeability, and " + "weather is emphatically not exchangeable: a front arrives and " + "yesterday's residual quantile becomes fiction. Adaptive " + "conformal feeds realised coverage back into the working alpha, " + "so the band widens after each miss and narrows after each hit. " + "Long-run coverage tracks the target whatever the distribution " + "does underneath.", + "failure": "If coverage sits far from target, the feedback rate is " + "wrong, not the model. Both are shown on the Models tab.", + "math": r"\alpha_{t+1} = \alpha_t + \gamma\left(\alpha^{*} - " + r"\mathbb{1}[y_t \notin C_t]\right)", + "params": {"target coverage": f"{int((1 - m.conformal_alpha) * 100)}%", + "gamma": f"{m.conformal_gamma:g}", + "window": f"{m.conformal_window} residuals"}, + }, + { + "id": "climatology", + "stage": "7", + "title": "Long-range outlook", + "module": "models/climatology.py", + "technique": "Ridge-regularised harmonic regression, anomaly decay", + "consumes": "Full history", + "produces": "Seven-day outlook with widening intervals", + "why": "An honest statement: a single point sensor cannot see a front " + "approaching from the Atlantic. Beyond about twelve hours the " + "only information it holds is where you are in the diurnal and " + "annual cycles, the current pressure anomaly, and the local " + "trend. So that is exactly what this uses, and the API labels " + "the result an outlook rather than a forecast.", + "failure": "Annual harmonics stay switched off below " + f"{m.climatology_min_days_annual:g} days of history. Fitting " + "a 365-day sine to three weeks of data produces a " + "magnificent extrapolation straight off the edge of the " + "physical world.", + "math": r"y \sim \beta_0 + \beta_1 t + \sum_{k=1}^{3}" + r"\left[a_k\sin\tfrac{2\pi k t}{\text{day}} + " + r"b_k\cos\tfrac{2\pi k t}{\text{day}}\right] + \text{annual}", + "params": {"diurnal harmonics": "3", "annual harmonics": "2", + "anomaly half-life": "30 h"}, + }, + { + "id": "precip", + "stage": "8", + "title": "Precipitation", + "module": "models/precip.py", + "technique": "Zambretti prior, online logistic residual by AdaGrad", + "consumes": "Sea-level pressure, tendency, humidity, cloud index, labels", + "produces": "Condition class and rain probability", + "why": "The 1915 Negretti and Zambra algorithm needs only pressure, its " + "tendency and the season. It has no parameters to overfit and " + "works from the first hour of deployment, so it is the prior. " + "The logistic layer learns only the residual: what your specific " + "location does that a slide rule cannot know. Its coefficient on " + "the Zambretti logit starts at exactly 1.0, so the model begins " + "as the slide rule and departs only where data insist.", + "failure": "Labels are the bottleneck. Without a rain gauge the proxy " + "label abstains in the ambiguous middle rather than " + "guessing, because a poisoned training set costs more than " + "the extra samples buy. Trust grows as n/(n+25) in strong " + "labels, so the two buttons on the Forecast tab matter.", + "params": {"prior": "Zambretti, three-branch", + "learner": "logistic, AdaGrad", + "strong label weight": "10x proxy"}, + }, + { + "id": "monitor", + "stage": "9", + "title": "Monitoring", + "module": "models/anomaly.py", + "technique": "Mahalanobis EWMA, Page-Hinkley, latch detection", + "consumes": "Filtered signals and matured forecast errors", + "produces": "Novelty score, drift alarms, per-sensor health", + "why": "Three detectors because they fail differently. Novelty catches " + "a window opening or a squall. Page-Hinkley catches the slow " + "stuff, a sensor drifting or a season turning, and it triggers " + "retraining, which is a far better signal than a cron schedule. " + "Latch detection catches the quietest failure of all: a sensor " + "that stops changing looks perfectly normal to both the others.", + "failure": "With six signals the sample covariance is singular for the " + "first hour, and a singular covariance turns Mahalanobis " + "distance into a random number generator with an " + "authoritative name. Shrinkage toward a scaled identity is " + "not optional.", + "params": {"novelty threshold": f"{m.anomaly_threshold:g}", + "EWMA lambda": f"{m.anomaly_ewma_lambda:g}", + "drift lambda": f"{m.drift_lambda:g}"}, + }, + { + "id": "verify", + "stage": "10", + "title": "Verification", + "module": "station.py", + "technique": "Rolling scoring against persistence and climatology", + "consumes": "Stored forecasts whose validity time has passed", + "produces": "MAE, RMSE, bias, coverage, skill", + "why": "This is the stage most projects skip and the one that makes the " + "difference. A forecast that is never scored is an opinion. A " + "forecast scored against persistence is a measurement. Skill is " + "1 - MAE/MAE_persistence, so a negative number is not a failure " + "of the exercise, it is the exercise working: ship persistence " + "at that horizon and stop pretending.", + "failure": "Nothing scores until forecasts mature, so the 24 hour row " + "is empty on day one. That is the loop being honest.", + "params": {"cadence": "every 5 minutes", + "baselines": "persistence, climatology"}, + }, + ] + + +def data_flow() -> List[Dict[str, str]]: + """Edges of the wiring diagram, drawn by the Methods tab.""" + return [ + {"from": "acquire", "to": "compensate", "label": "T_raw, T_cpu"}, + {"from": "compensate", "to": "kalman", "label": "T corrected"}, + {"from": "kalman", "to": "features", "label": "level + rate"}, + {"from": "kalman", "to": "precip", "label": "dp/dt"}, + {"from": "kalman", "to": "monitor", "label": "signals"}, + {"from": "features", "to": "nowcast", "label": "design matrix"}, + {"from": "features", "to": "climatology", "label": "history"}, + {"from": "nowcast", "to": "conformal", "label": "point forecast"}, + {"from": "climatology", "to": "nowcast", "label": "member"}, + {"from": "conformal", "to": "verify", "label": "interval"}, + {"from": "verify", "to": "conformal", "label": "coverage feedback"}, + {"from": "verify", "to": "monitor", "label": "errors"}, + {"from": "monitor", "to": "nowcast", "label": "retrain trigger"}, + {"from": "precip", "to": "verify", "label": "labels"}, + ] + + +def glossary() -> List[Dict[str, str]]: + return [ + {"term": "Skill", + "definition": "1 - MAE/MAE_persistence. Zero means no better than " + "assuming nothing changes. Negative means worse than that, " + "which is useful information rather than an embarrassment."}, + {"term": "Coverage", + "definition": "Fraction of observations that landed inside the prediction " + "interval. Should sit near the target. Far above means the " + "bands are lazily wide, far below means they lie."}, + {"term": "Forgetting factor", + "definition": "Exponential weight on past samples. 0.999 on a 5-minute " + "grid remembers roughly a day; 0.99 remembers about two " + "hours and chases noise."}, + {"term": "Persistence", + "definition": "The baseline forecast: tomorrow equals today. Beating it " + "over short horizons is genuinely hard, which is why it is " + "the honest thing to measure against."}, + {"term": "Pressure tendency", + "definition": "Rate of change of sea-level pressure. Falling fast means " + "an approaching low. This is the only variable in the " + "station that sees beyond your walls."}, + {"term": "Dew point depression", + "definition": "Air temperature minus dew point. Small and shrinking means " + "saturation, fog or rain. Large means dry air."}, + ] + + +def _fmt(seconds: int) -> str: + if seconds < 3600: + return f"{seconds // 60} min" + if seconds < 86400: + return f"{seconds // 3600} h" + return f"{seconds // 86400} d" + + +def _memory(lam: float, grid_s: int) -> str: + """Effective memory of an exponential forgetting factor, 1/(1-lambda) samples.""" + if lam >= 1.0: + return "unbounded" + samples = 1.0 / (1.0 - lam) + hours = samples * grid_s / 3600.0 + return f"~{samples:.0f} samples ({hours:.1f} h)" + + +def describe(cfg) -> Dict[str, Any]: + return { + "pipeline": pipeline(cfg), + "flow": data_flow(), + "glossary": glossary(), + "features": FEATURE_NAMES, + "honest_limits": [ + "Indoors this forecasts your room, not the sky. Pressure is the " + "exception because it passes through walls, which is exactly why the " + "precipitation model runs on pressure tendency rather than indoor " + "humidity.", + "Days two to seven are climatology with an anomaly correction, not a " + "forecast. The station physically cannot observe an approaching " + "system.", + "Without a rain gauge, precipitation labels come from you. The learner " + "earns influence in proportion to how many you have supplied.", + "Every number on the Models tab is measured on your own data, not " + "quoted from a benchmark. If skill is negative at some horizon, that " + "is what your station is actually doing.", + ], + } diff --git a/ashvale/models/__init__.py b/ashvale/models/__init__.py new file mode 100644 index 0000000..e06fa91 --- /dev/null +++ b/ashvale/models/__init__.py @@ -0,0 +1,24 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from .rls import RecursiveLeastSquares, AdaptiveConformal +from .nowcast import NowcastEnsemble +from .climatology import HarmonicClimatology +from .precip import PrecipitationModel, zambretti +from .anomaly import AnomalyMonitor + +__all__ = [ + "RecursiveLeastSquares", "AdaptiveConformal", "NowcastEnsemble", + "HarmonicClimatology", "PrecipitationModel", "zambretti", "AnomalyMonitor", +] diff --git a/ashvale/models/anomaly.py b/ashvale/models/anomaly.py new file mode 100644 index 0000000..1aba369 --- /dev/null +++ b/ashvale/models/anomaly.py @@ -0,0 +1,287 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Anomaly and drift monitoring: the part that keeps the rest honest. + +Three independent detectors, because they fail in different ways: + +`MahalanobisEWMA` Multivariate novelty on the residual from a slowly + updated mean and shrinkage covariance. Catches a + window opening, a heater cycling, or a genuine squall. + +`PageHinkley` Sequential change-point detection on model error. + Catches the slow stuff: a sensor drifting, a season + turning, a model quietly going stale. This is the + detector that tells you *when to retrain*, which is a + far better trigger than a cron schedule. + +`SensorHealth` Latched values, out-of-range readings, and Kalman + innovation inflation. A stuck sensor is invisible to + the other two because it looks perfectly normal. + +Shrinkage on the covariance is not optional here. With 6 signals and a +1000-sample window the sample covariance is fine, but during the first +hour it is singular, and a singular covariance turns Mahalanobis +distance into a random number generator with an authoritative name. +""" + +from __future__ import annotations + +from collections import deque +from typing import Deque, Dict, List, Optional + +import numpy as np + +SIGNALS = ["temp_c", "hum", "press_slp", "temp_rate", "press_rate", "dew_c"] + + +class MahalanobisEWMA: + def __init__(self, n_dims: int, lam: float = 0.15, threshold: float = 12.0, + shrinkage: float = 0.15, warmup: int = 60): + self.d = int(n_dims) + self.lam = float(lam) + self.threshold = float(threshold) + self.shrinkage = float(shrinkage) + self.warmup = int(warmup) + self.mean = np.zeros(self.d) + self.cov = np.eye(self.d) + self.z = np.zeros(self.d) # EWMA of standardised residual + self.n = 0 + self.last_d2 = 0.0 + + def update(self, x: np.ndarray) -> Dict: + x = np.asarray(x, dtype=float).ravel() + if x.size != self.d or not np.all(np.isfinite(x)): + return {"d2": self.last_d2, "alarm": False, "warm": self.n < self.warmup} + + self.n += 1 + if self.n == 1: + self.mean = x.copy() + return {"d2": 0.0, "alarm": False, "warm": True} + + a = 1.0 / min(self.n, 500) # slow adaptation once warm + delta = x - self.mean + self.mean += a * delta + self.cov = (1 - a) * self.cov + a * np.outer(delta, delta) + + # Ledoit-Wolf style shrinkage toward a scaled identity + target = np.eye(self.d) * (np.trace(self.cov) / self.d + 1e-9) + cov = (1 - self.shrinkage) * self.cov + self.shrinkage * target + + try: + resid = np.linalg.solve(cov, delta) + except np.linalg.LinAlgError: + return {"d2": self.last_d2, "alarm": False, "warm": True} + + # EWMA on the whitened residual gives persistence-aware detection: + # one odd sample is noise, ten in a row is an event. + white = delta / np.sqrt(np.maximum(np.diag(cov), 1e-12)) + self.z = (1 - self.lam) * self.z + self.lam * white + scale = self.lam / (2 - self.lam) + d2_ewma = float(self.z @ self.z / max(scale, 1e-9)) + d2_inst = float(delta @ resid) + self.last_d2 = d2_ewma + + warm = self.n < self.warmup + return { + "d2": d2_ewma, + "d2_instant": d2_inst, + "alarm": (not warm) and d2_ewma > self.threshold, + "warm": warm, + "contributions": {s: round(float(v), 2) for s, v in zip(SIGNALS[:self.d], white)}, + } + + def to_dict(self) -> Dict: + return {"d": self.d, "lam": self.lam, "threshold": self.threshold, + "shrinkage": self.shrinkage, "warmup": self.warmup, + "mean": self.mean.tolist(), "cov": self.cov.tolist(), + "z": self.z.tolist(), "n": self.n} + + @classmethod + def from_dict(cls, s: Dict) -> "MahalanobisEWMA": + m = cls(s["d"], s["lam"], s["threshold"], s["shrinkage"], s["warmup"]) + m.mean = np.array(s["mean"], float) + m.cov = np.array(s["cov"], float) + m.z = np.array(s["z"], float) + m.n = s["n"] + return m + + +class PageHinkley: + """Two-sided sequential change detection on a stream of errors.""" + + def __init__(self, delta: float = 0.05, lam: float = 8.0, alpha: float = 0.999): + self.delta = float(delta) + self.lam = float(lam) + self.alpha = float(alpha) + self.mean = 0.0 + self.n = 0 + self.m_pos = 0.0 + self.m_neg = 0.0 + self.n_alarms = 0 + self.last_alarm_ts: Optional[float] = None + + def update(self, value: float, ts: Optional[float] = None) -> bool: + v = float(value) + if not np.isfinite(v): + return False + self.n += 1 + self.mean += (v - self.mean) / self.n + + self.m_pos = self.alpha * max(0.0, self.m_pos + v - self.mean - self.delta) + self.m_neg = self.alpha * max(0.0, self.m_neg - v + self.mean - self.delta) + + if self.n > 30 and max(self.m_pos, self.m_neg) > self.lam: + self.reset_statistics() + self.n_alarms += 1 + self.last_alarm_ts = ts + return True + return False + + def reset_statistics(self) -> None: + self.m_pos = 0.0 + self.m_neg = 0.0 + self.n = 1 + + @property + def stress(self) -> float: + """0 to 1: how close we are to declaring drift. Nice on a gauge.""" + return float(min(max(self.m_pos, self.m_neg) / max(self.lam, 1e-9), 1.0)) + + def to_dict(self) -> Dict: + return {"delta": self.delta, "lam": self.lam, "alpha": self.alpha, + "mean": self.mean, "n": self.n, "m_pos": self.m_pos, + "m_neg": self.m_neg, "n_alarms": self.n_alarms, + "last_alarm_ts": self.last_alarm_ts} + + @classmethod + def from_dict(cls, s: Dict) -> "PageHinkley": + p = cls(s["delta"], s["lam"], s["alpha"]) + p.__dict__.update({k: s[k] for k in + ("mean", "n", "m_pos", "m_neg", "n_alarms", "last_alarm_ts")}) + return p + + +class SensorHealth: + RANGES = { + "temp_c": (-40.0, 85.0), + "hum": (0.0, 100.0), + "press_slp": (870.0, 1085.0), + "cpu_temp": (-20.0, 95.0), + } + + def __init__(self, window: int = 90): + self.buffers: Dict[str, Deque[float]] = { + k: deque(maxlen=window) for k in self.RANGES + } + self.flags: Dict[str, str] = {} + + def update(self, obs: Dict[str, float]) -> Dict[str, Dict]: + report = {} + for name, (lo, hi) in self.RANGES.items(): + v = obs.get(name) + if v is None or not np.isfinite(v): + report[name] = {"status": "missing", "detail": "no reading"} + continue + buf = self.buffers[name] + buf.append(float(v)) + arr = np.asarray(buf, dtype=float) + + if not (lo <= v <= hi): + status, detail = "fault", f"out of range ({v:.2f})" + elif arr.size >= 20 and float(np.max(np.abs(np.diff(arr)))) < 1e-9: + status, detail = "fault", "value latched, sensor may be stuck" + elif arr.size >= 20 and float(np.std(arr)) < 1e-6: + status, detail = "warn", "near-zero variance" + else: + status, detail = "ok", "nominal" + report[name] = {"status": status, "detail": detail, + "value": float(v), "std": float(np.std(arr)) if arr.size > 2 else 0.0} + self.flags = {k: v["status"] for k, v in report.items()} + return report + + @property + def overall(self) -> str: + if any(v == "fault" for v in self.flags.values()): + return "fault" + if any(v == "warn" for v in self.flags.values()): + return "warn" + return "ok" + + +class AnomalyMonitor: + """Facade over the three detectors, with a rolling event log.""" + + def __init__(self, cfg_model): + self.novelty = MahalanobisEWMA( + len(SIGNALS), cfg_model.anomaly_ewma_lambda, cfg_model.anomaly_threshold + ) + self.drift = PageHinkley(cfg_model.drift_delta, cfg_model.drift_lambda) + self.health = SensorHealth() + self.events: Deque[Dict] = deque(maxlen=100) + self.retrain_requested = False + + def observe(self, ts: float, obs: Dict[str, float]) -> Dict: + vec = np.array([obs.get(s, np.nan) for s in SIGNALS], dtype=float) + nov = self.novelty.update(vec) + health = self.health.update(obs) + + if nov.get("alarm"): + top = max(nov.get("contributions", {}).items(), + key=lambda kv: abs(kv[1]), default=("unknown", 0.0)) + self._log(ts, "novelty", "warn", + f"multivariate departure d2={nov['d2']:.1f}, led by {top[0]}") + for name, rep in health.items(): + if rep["status"] == "fault": + self._log(ts, "sensor", "error", f"{name}: {rep['detail']}") + + return { + "novelty": nov, + "health": health, + "health_overall": self.health.overall, + "drift": { + "stress": self.drift.stress, + "alarms": self.drift.n_alarms, + "last_alarm_ts": self.drift.last_alarm_ts, + "retrain_requested": self.retrain_requested, + }, + } + + def observe_error(self, ts: float, abs_error: float) -> bool: + """Feed a matured forecast error; returns True if drift was declared.""" + fired = self.drift.update(abs_error, ts) + if fired: + self.retrain_requested = True + self._log(ts, "drift", "warn", + "forecast error distribution shifted, retrain queued") + return fired + + def clear_retrain_flag(self) -> None: + self.retrain_requested = False + + def _log(self, ts: float, kind: str, severity: str, detail: str) -> None: + self.events.append({"ts": ts, "kind": kind, "severity": severity, "detail": detail}) + + def recent(self, n: int = 20) -> List[Dict]: + return list(self.events)[-n:][::-1] + + def to_dict(self) -> Dict: + return {"novelty": self.novelty.to_dict(), "drift": self.drift.to_dict(), + "events": list(self.events), "retrain_requested": self.retrain_requested} + + def load_dict(self, s: Dict) -> None: + self.novelty = MahalanobisEWMA.from_dict(s["novelty"]) + self.drift = PageHinkley.from_dict(s["drift"]) + self.events = deque(s.get("events", []), maxlen=100) + self.retrain_requested = s.get("retrain_requested", False) diff --git a/ashvale/models/climatology.py b/ashvale/models/climatology.py new file mode 100644 index 0000000..76e155b --- /dev/null +++ b/ashvale/models/climatology.py @@ -0,0 +1,170 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Harmonic regression: the long-range half of the forecast. + +An honest statement first, because a weather product that oversells +itself is worse than no product. A single point sensor cannot see a +front approaching from the Atlantic. Beyond roughly twelve hours, the +only information your station holds is: + + * where in the diurnal cycle you are, + * where in the annual cycle you are, + * the current synoptic pressure anomaly and its tendency, + * the local trend of the last few days. + +So that is exactly what this model uses. It is a ridge-regularised +Fourier basis in time-of-day and day-of-year, plus a slow linear trend +and a pressure-anomaly coupling. Days 2 to 7 are a *climatological +outlook with an anomaly correction*, not a forecast, and the API labels +them as such. Anything more confident would be theatre. + +The annual harmonics only switch on once the station has enough history +to identify them (`climatology_min_days_annual`, default 120). Before +that, fitting a 365-day sine to three weeks of data produces a +magnificent extrapolation straight off the edge of the physical world. +""" + +from __future__ import annotations + +from typing import Dict, List, Optional + +import numpy as np + +DAY = 86400.0 +YEAR = 365.2422 * DAY + + +class HarmonicClimatology: + def __init__(self, targets, diurnal_harmonics: int = 3, + annual_harmonics: int = 2, ridge: float = 1.0, + min_days_annual: float = 120.0): + self.targets = tuple(targets) + self.kd = int(diurnal_harmonics) + self.ka = int(annual_harmonics) + self.ridge = float(ridge) + self.min_days_annual = float(min_days_annual) + self.coef: Dict[str, np.ndarray] = {} + self.resid_std: Dict[str, float] = {} + self.t0: float = 0.0 + self.use_annual = False + self.n_days = 0.0 + self.ready = False + + # ---------------------------------------------------------- basis + + def _design(self, ts: np.ndarray) -> np.ndarray: + ts = np.atleast_1d(np.asarray(ts, dtype=float)) + t_days = (ts - self.t0) / DAY + cols = [np.ones(ts.size), t_days / 30.0] # slow trend, per month + for k in range(1, self.kd + 1): + w = 2 * np.pi * k * ts / DAY + cols += [np.sin(w), np.cos(w)] + if self.use_annual: + for k in range(1, self.ka + 1): + w = 2 * np.pi * k * ts / YEAR + cols += [np.sin(w), np.cos(w)] + return np.column_stack(cols) + + # ------------------------------------------------------------ fit + + def fit(self, ts: np.ndarray, series: Dict[str, np.ndarray], + valid: Optional[np.ndarray] = None) -> Dict[str, float]: + ts = np.asarray(ts, dtype=float) + if ts.size < 48: + self.ready = False + return {} + self.t0 = float(ts[0]) + self.n_days = float((ts[-1] - ts[0]) / DAY) + self.use_annual = self.n_days >= self.min_days_annual + + A = self._design(ts) + mask = np.ones(ts.size, dtype=bool) if valid is None else valid.astype(bool) + out = {} + for target in self.targets: + y = np.asarray(series.get(target, np.empty(0)), dtype=float) + if y.size != ts.size: + continue + m = mask & np.isfinite(y) + if m.sum() < A.shape[1] * 3: + continue + Am, ym = A[m], y[m] + # ridge: leave the intercept unpenalised + reg = np.eye(A.shape[1]) * self.ridge + reg[0, 0] = 0.0 + beta = np.linalg.solve(Am.T @ Am + reg, Am.T @ ym) + self.coef[target] = beta + resid = ym - Am @ beta + self.resid_std[target] = float(np.std(resid)) + out[target] = self.resid_std[target] + self.ready = bool(self.coef) + return out + + # -------------------------------------------------------- predict + + def predict(self, target: str, ts: np.ndarray) -> np.ndarray: + ts = np.atleast_1d(np.asarray(ts, dtype=float)) + beta = self.coef.get(target) + if beta is None: + return np.zeros(ts.size) + return self._design(ts) @ beta + + def outlook(self, target: str, now: float, days: int = 7, + step_s: int = 3 * 3600, anomaly: float = 0.0, + anomaly_halflife_h: float = 30.0) -> List[Dict]: + """Climatology plus an exponentially decaying current anomaly. + + The anomaly term is what makes this better than a textbook: if + today is 3 C above the seasonal norm, tomorrow morning probably + still is, and next Thursday almost certainly is not. The decay + half-life encodes exactly that intuition, and the interval widens + with the square root of lead time as any diffusive process should. + """ + if not self.ready or target not in self.coef: + return [] + grid = np.arange(now, now + days * DAY, step_s, dtype=float) + base = self.predict(target, grid) + lead_h = (grid - now) / 3600.0 + decay = 0.5 ** (lead_h / max(anomaly_halflife_h, 1e-3)) + mu = base + anomaly * decay + sigma0 = self.resid_std.get(target, 1.0) + sigma = sigma0 * np.sqrt(1.0 + lead_h / 24.0) + return [ + {"ts": float(t), "lead_h": float(l), "mu": float(m), + "lo": float(m - 1.645 * s), "hi": float(m + 1.645 * s)} + for t, l, m, s in zip(grid, lead_h, mu, sigma) + ] + + def anomaly_now(self, target: str, ts: float, observed: float) -> float: + if not self.ready or target not in self.coef: + return 0.0 + return float(observed - self.predict(target, np.array([ts]))[0]) + + def to_dict(self) -> Dict: + return {"targets": list(self.targets), "kd": self.kd, "ka": self.ka, + "ridge": self.ridge, "min_days_annual": self.min_days_annual, + "t0": self.t0, "use_annual": self.use_annual, "n_days": self.n_days, + "coef": {k: v.tolist() for k, v in self.coef.items()}, + "resid_std": self.resid_std, "ready": self.ready} + + def load_dict(self, s: Dict) -> None: + self.kd, self.ka = s["kd"], s["ka"] + self.ridge = s["ridge"] + self.min_days_annual = s["min_days_annual"] + self.t0 = s["t0"] + self.use_annual = s["use_annual"] + self.n_days = s.get("n_days", 0.0) + self.coef = {k: np.array(v, dtype=float) for k, v in s["coef"].items()} + self.resid_std = s["resid_std"] + self.ready = s["ready"] diff --git a/ashvale/models/nowcast.py b/ashvale/models/nowcast.py new file mode 100644 index 0000000..ad5ea19 --- /dev/null +++ b/ashvale/models/nowcast.py @@ -0,0 +1,245 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Multi-horizon forecasting: one direct head per (target, horizon). + +Direct rather than recursive. A recursive one-step model iterated 288 +times to reach 24 hours compounds its own bias into a beautifully smooth +lie. Direct heads cost more memory (six horizons x three targets = 18 +small models, about 150 kB total) and are worth every byte. + +Each head predicts a *delta from now*, then the ensemble blends three +opinions with weights that are themselves learned online: + + persistence : it will be exactly as it is now + climatology : it will be whatever this hour of this day usually is + learned RLS : it will be now plus what the regressors imply + +Persistence wins at 15 minutes. Climatology wins at 24 hours. The RLS +head wins in the middle, which is exactly the region a physical +forecaster finds hardest. The blend weights are updated by exponentiated +gradient (Hedge), so the ensemble is never worse than its best member by +more than a log factor, and it re-weights itself within a day when the +season turns. +""" + +from __future__ import annotations + +from typing import Dict, List, Optional, Tuple + +import numpy as np + +from ..features import N_FEATURES, Standardiser, supervised_pairs +from .rls import AdaptiveConformal, RecursiveLeastSquares + +MEMBERS = ("persistence", "climatology", "learned") + + +class ForecastHead: + """One target, one horizon.""" + + def __init__(self, target: str, horizon_s: int, n_features: int = N_FEATURES, + forgetting: float = 0.9985, delta: float = 100.0, + alpha: float = 0.10, conformal_window: int = 400, + gamma: float = 0.01, hedge_eta: float = 0.35): + self.target = target + self.horizon_s = int(horizon_s) + self.model = RecursiveLeastSquares(n_features, forgetting, delta) + self.conformal = AdaptiveConformal(alpha, conformal_window, gamma) + self.weights = np.ones(len(MEMBERS)) / len(MEMBERS) + self.eta = float(hedge_eta) + self.member_mae = np.zeros(len(MEMBERS)) + self.n_scored = 0 + + # -------------------------------------------------------- prediction + + def predict(self, x: np.ndarray, anchor: float, + climatology_delta: float = 0.0) -> Dict[str, float]: + learned_delta = self.model.predict(x) + deltas = np.array([0.0, float(climatology_delta), float(learned_delta)]) + blended = float(np.dot(self.weights, deltas)) + mu = float(anchor + blended) + sigma = self.model.predict_std(x, self.model.noise_var) + lo, hi = self.conformal.interval(mu, fallback_sigma=sigma) + return { + "mu": mu, + "lo": lo, + "hi": hi, + "sigma": sigma, + "delta": blended, + "members": {m: float(anchor + d) for m, d in zip(MEMBERS, deltas)}, + "weights": {m: float(w) for m, w in zip(MEMBERS, self.weights)}, + } + + # ---------------------------------------------------------- learning + + def learn(self, x: np.ndarray, anchor: float, truth: float, + climatology_delta: float = 0.0) -> float: + """One supervised step given a matured target.""" + deltas = np.array([0.0, float(climatology_delta), + float(self.model.predict(x))]) + member_pred = anchor + deltas + losses = np.abs(member_pred - truth) + + # Hedge / exponentiated gradient on normalised losses + scale = max(float(np.max(losses)), 1e-6) + self.weights *= np.exp(-self.eta * losses / scale) + self.weights = np.clip(self.weights, 1e-4, None) + self.weights /= self.weights.sum() + + blended = float(np.dot(self.weights, member_pred)) + residual = truth - blended + self.conformal.observe(residual) + + self.model.update(x, truth - anchor) + + self.member_mae = 0.98 * self.member_mae + 0.02 * losses + self.n_scored += 1 + return residual + + def to_dict(self) -> Dict: + return {"target": self.target, "horizon_s": self.horizon_s, + "model": self.model.to_dict(), "conformal": self.conformal.to_dict(), + "weights": self.weights.tolist(), "eta": self.eta, + "member_mae": self.member_mae.tolist(), "n_scored": self.n_scored} + + @classmethod + def from_dict(cls, s: Dict) -> "ForecastHead": + h = cls(s["target"], s["horizon_s"]) + h.model = RecursiveLeastSquares.from_dict(s["model"]) + h.conformal = AdaptiveConformal.from_dict(s["conformal"]) + h.weights = np.array(s["weights"], dtype=float) + h.eta = s["eta"] + h.member_mae = np.array(s["member_mae"], dtype=float) + h.n_scored = s.get("n_scored", 0) + return h + + +class NowcastEnsemble: + """The full bank of heads plus the shared feature standardiser.""" + + def __init__(self, targets: Tuple[str, ...], horizons_s: Tuple[int, ...], + cfg_model): + self.targets = tuple(targets) + self.horizons = tuple(int(h) for h in horizons_s) + self.grid_s = int(cfg_model.grid_s) + self.scaler = Standardiser(N_FEATURES) + self.heads: Dict[Tuple[str, int], ForecastHead] = { + (t, h): ForecastHead( + t, h, N_FEATURES, cfg_model.rls_forgetting, cfg_model.rls_delta, + cfg_model.conformal_alpha, cfg_model.conformal_window, + cfg_model.conformal_gamma, + ) + for t in self.targets for h in self.horizons + } + self.trained_rows = 0 + + # ------------------------------------------------------------ train + + def fit(self, X: np.ndarray, valid: np.ndarray, series: Dict[str, np.ndarray], + climatology=None, grid_ts: Optional[np.ndarray] = None, + passes: int = 1, max_pairs: int = 2500) -> Dict[str, int]: + """Batch-update every head from history. + + `max_pairs` bounds the work per head to the most recent samples. + This is not a shortcut: with a forgetting factor of 0.9985 the + effective memory is about 11 hours, so the 4000th-most-recent + sample carries a weight of roughly e^-6. Training on it costs + real seconds on a Cortex-A53 and buys nothing measurable. + """ + """Batch pass over history. Called on startup and every retrain tick.""" + if X.shape[0] < 10: + return {"rows": 0} + self.scaler.partial_fit(X[valid][:: max(1, X.shape[0] // 2000)]) + Xs = self.scaler.transform(X) + + counts = {} + for target in self.targets: + y = series[target] + for h in self.horizons: + steps = max(int(round(h / self.grid_s)), 1) + Xa, dy, anchor = supervised_pairs(Xs, valid, y, steps) + if Xa.shape[0] < 5: + counts[f"{target}@{h}"] = 0 + continue + if Xa.shape[0] > max_pairs: + Xa, dy, anchor = Xa[-max_pairs:], dy[-max_pairs:], anchor[-max_pairs:] + head = self.heads[(target, h)] + clim = np.zeros(Xa.shape[0]) + if climatology is not None and grid_ts is not None and climatology.ready: + n = grid_ts.size + ts_a = grid_ts[:n - steps] + mask_len = min(ts_a.size, Xa.shape[0]) + clim_now = climatology.predict(target, ts_a[-mask_len:]) + clim_fut = climatology.predict(target, ts_a[-mask_len:] + h) + clim = np.zeros(Xa.shape[0]) + clim[-mask_len:] = clim_fut - clim_now + for _ in range(max(int(passes), 1)): + for i in range(Xa.shape[0]): + head.learn(Xa[i], anchor[i], anchor[i] + dy[i], clim[i]) + counts[f"{target}@{h}"] = int(Xa.shape[0]) + self.trained_rows = int(X.shape[0]) + return counts + + # --------------------------------------------------------- inference + + def forecast(self, x_raw: np.ndarray, anchors: Dict[str, float], now: float, + climatology=None) -> Dict[str, Dict[int, Dict[str, float]]]: + x = self.scaler.transform(np.atleast_2d(x_raw))[0] + out: Dict[str, Dict[int, Dict[str, float]]] = {} + for target in self.targets: + anchor = float(anchors.get(target, 0.0)) + out[target] = {} + for h in self.horizons: + clim_delta = 0.0 + if climatology is not None and climatology.ready: + clim_delta = float(climatology.predict(target, np.array([now + h]))[0] + - climatology.predict(target, np.array([now]))[0]) + out[target][h] = self.heads[(target, h)].predict(x, anchor, clim_delta) + return out + + def diagnostics(self) -> List[Dict]: + rows = [] + for (target, h), head in sorted(self.heads.items()): + rows.append({ + "target": target, + "horizon_s": h, + "n_updates": head.model.n_updates, + "n_scored": head.n_scored, + "weights": {m: round(float(w), 3) for m, w in zip(MEMBERS, head.weights)}, + "member_mae": {m: round(float(v), 3) for m, v in zip(MEMBERS, head.member_mae)}, + "conformal_alpha": round(head.conformal.alpha, 4), + "conformal_halfwidth": round(float(head.conformal.quantile()), 3) + if np.isfinite(head.conformal.quantile()) else None, + "coverage": round(head.conformal.empirical_coverage, 3) + if np.isfinite(head.conformal.empirical_coverage) else None, + }) + return rows + + def to_dict(self) -> Dict: + return { + "targets": list(self.targets), + "horizons": list(self.horizons), + "grid_s": self.grid_s, + "scaler": self.scaler.to_dict(), + "heads": [h.to_dict() for h in self.heads.values()], + "trained_rows": self.trained_rows, + } + + def load_dict(self, s: Dict) -> None: + self.scaler = Standardiser.from_dict(s["scaler"]) + for hs in s["heads"]: + head = ForecastHead.from_dict(hs) + self.heads[(head.target, head.horizon_s)] = head + self.trained_rows = s.get("trained_rows", 0) diff --git a/ashvale/models/precip.py b/ashvale/models/precip.py new file mode 100644 index 0000000..31bb1b8 --- /dev/null +++ b/ashvale/models/precip.py @@ -0,0 +1,323 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Will it rain? A prior with a hundred years of service, plus a learner. + +Two components, deliberately: + +1. `zambretti()` is the 1915 Negretti and Zambra slide-rule algorithm, + re-expressed here in the standard three-branch form. It needs only + sea-level pressure, its tendency and the season. It has no parameters + to overfit, it works from the first hour of deployment, and in the + temperate maritime climate it was designed for it is genuinely hard + to beat with a small dataset. It is the prior. + +2. `PrecipitationModel` is an online logistic regression that learns the + residual: what your specific location does that the slide rule does + not know. It starts from the Zambretti logit and only earns influence + as labels accumulate, so it cannot embarrass you on day one. + +Labels are the hard part, and the design is explicit about it. Without a +rain gauge, a *proxy* label is used (near-saturated air with a collapsing +dew-point depression), and it is flagged as weak. `POST /api/label` lets +you supply ground truth from a window: two seconds of your attention is +worth a week of proxy labels, and the learner weights them accordingly. +""" + +from __future__ import annotations + +import math +import time +from typing import Dict, List, Optional, Tuple + +import numpy as np + +# Severity classes the Z number maps onto. Wording is ours, not the +# original card's, and is deliberately about actionable state rather +# than Edwardian poetry. +_CONDITIONS = [ + (1, 2, "settled", "Settled and dry"), + (3, 5, "fine", "Fine, little change expected"), + (6, 8, "fair", "Fair, becoming less settled"), + (9, 12, "changeable", "Changeable, showers possible"), + (13, 16, "unsettled", "Unsettled, rain at times"), + (17, 20, "rain", "Rain likely, turning wet"), + (21, 23, "wet", "Wet and windy"), + (24, 26, "stormy", "Stormy, heavy rain likely"), +] + +_RAIN_PRIOR = { + "settled": 0.03, "fine": 0.07, "fair": 0.15, "changeable": 0.32, + "unsettled": 0.52, "rain": 0.72, "wet": 0.85, "stormy": 0.93, +} + +FEATURES = ["bias", "slp_anom", "tend_1h", "tend_3h", "tend_6h", "rh", + "dew_depression", "dew_dep_rate", "cloud_index", "temp_dev", + "wet_bulb_depression", "zambretti_logit"] + + +def _season_is_summer(ts: Optional[float], latitude: float) -> bool: + month = time.gmtime(ts or time.time()).tm_mon + northern = latitude >= 0 + summer_months = {4, 5, 6, 7, 8, 9} + return (month in summer_months) if northern else (month not in summer_months) + + +BARO_BOTTOM = 950.0 +BARO_TOP = 1050.0 + +# Each branch maps normalised pressure onto a slice of the 26-point scale. +# The ordering is the whole point of the instrument: for a given pressure, +# rising air is always a better forecast than falling air, and within a +# branch higher pressure is always better. Ranges overlap because a deep +# but rising low really is more hopeful than a shallow but falling high. +_BRANCH = { + "rising": (1.0, 10.0), + "steady": (6.0, 17.0), + "falling": (11.0, 26.0), +} + + +def zambretti(slp_hpa: float, tendency_hpa_per_h: float, + ts: Optional[float] = None, latitude: float = 52.0, + steady_band: float = 0.10) -> Dict: + """Three-branch barometric forecast on the Zambretti 26-point scale. + + The 1915 Negretti and Zambra slide rule read pressure, its tendency and + the season off a rotating card and returned one of 26 outcomes, 1 being + settled and 26 being stormy. Published transcriptions of its constants + disagree with each other, so rather than mis-cite one, this is an + explicit re-parameterisation onto the same 26-point scale, anchored to + the behaviour the instrument is actually known for: + + rising pressure -> lower Z (improving) + falling pressure -> higher Z (deteriorating) + higher pressure -> lower Z within any branch + + Getting that sign wrong is easy and produces confident nonsense: a + barometer climbing hard while the panel reads `stormy` is the tell. + + Args: + slp_hpa: pressure reduced to mean sea level. Passing station + pressure here is a common and silent bug: at 100 m elevation + it shifts the result by about two categories, permanently. + tendency_hpa_per_h: Kalman-filtered rate, not a finite difference. + steady_band: |tendency| below this counts as steady. + """ + p = float(np.clip(slp_hpa, BARO_BOTTOM, BARO_TOP)) + tend = float(tendency_hpa_per_h) + summer = _season_is_summer(ts, latitude) + + if tend <= -steady_band: + trend = "falling" + elif tend >= steady_band: + trend = "rising" + else: + trend = "steady" + + lo, hi = _BRANCH[trend] + u = (p - BARO_BOTTOM) / (BARO_TOP - BARO_BOTTOM) # 0 at 950, 1 at 1050 + z = lo + (hi - lo) * (1.0 - u) + + # Seasonal nudge: summer lows are typically convective and shorter lived, + # winter lows are frontal and grimmer. One category either way. + if trend == "falling": + z += -1.0 if summer else 1.0 + elif trend == "rising": + z += -1.0 if summer else 1.0 + + z_int = int(np.clip(round(z), 1, 26)) + condition, label = "changeable", "Changeable" + for lo, hi, key, text in _CONDITIONS: + if lo <= z_int <= hi: + condition, label = key, text + break + + return { + "z": z_int, + "trend": trend, + "condition": condition, + "label": label, + "prior_rain_prob": _RAIN_PRIOR[condition], + "slp_used": p, + "tendency": tend, + "season": "summer" if summer else "winter", + } + + +def tendency_code(tend_hpa_per_h: float) -> str: + """WMO-style pressure characteristic, the thing sailors actually read.""" + t = float(tend_hpa_per_h) + if t <= -1.5: + return "falling very rapidly" + if t <= -0.6: + return "falling rapidly" + if t <= -0.15: + return "falling" + if t < 0.15: + return "steady" + if t < 0.6: + return "rising" + if t < 1.5: + return "rising rapidly" + return "rising very rapidly" + + +def _sigmoid(z: float) -> float: + return 1.0 / (1.0 + math.exp(-float(np.clip(z, -30.0, 30.0)))) + + +def _logit(p: float) -> float: + p = float(np.clip(p, 1e-4, 1 - 1e-4)) + return math.log(p / (1 - p)) + + +class PrecipitationModel: + """Online logistic regression on top of the Zambretti logit. + + Trained by AdaGrad because feature scales here vary by two orders of + magnitude and a fixed learning rate would either crawl on `tendency` + or explode on `rh`. The `zambretti_logit` feature is initialised with + a coefficient of 1.0 so the model *starts* as the slide rule and + departs from it only where the data insist. + """ + + def __init__(self, lr: float = 0.08, l2: float = 1e-4): + self.w = np.zeros(len(FEATURES)) + self.w[FEATURES.index("zambretti_logit")] = 1.0 + self.g2 = np.ones(len(FEATURES)) * 1e-3 + self.lr = float(lr) + self.l2 = float(l2) + self.n_strong = 0 + self.n_weak = 0 + self.ewma_logloss = 0.693 # log 2, the coin-flip baseline + self.mean = np.zeros(len(FEATURES)) + self.m2 = np.ones(len(FEATURES)) + self.n_seen = 0 + + # -------------------------------------------------------- features + + def featurise(self, obs: Dict, zam: Dict) -> np.ndarray: + x = np.array([ + 1.0, + obs.get("slp", 1013.25) - 1013.25, + obs.get("tend_1h", 0.0), + obs.get("tend_3h", 0.0), + obs.get("tend_6h", 0.0), + (obs.get("rh", 60.0) - 70.0) / 10.0, + obs.get("dew_depression", 5.0), + obs.get("dew_dep_rate", 0.0), + obs.get("cloud_index", 0.5), + obs.get("temp_dev", 0.0), + obs.get("wet_bulb_depression", 2.0), + _logit(zam["prior_rain_prob"]), + ], dtype=float) + return np.nan_to_num(x, nan=0.0, posinf=0.0, neginf=0.0) + + def _standardise(self, x: np.ndarray, update: bool) -> np.ndarray: + if update: + self.n_seen += 1 + delta = x - self.mean + self.mean += delta / self.n_seen + self.m2 += delta * (x - self.mean) + if self.n_seen < 20: + z = x.copy() + else: + std = np.sqrt(self.m2 / max(self.n_seen - 1, 1)) + std = np.where(std < 1e-8, 1.0, std) + z = (x - self.mean) / std + z[0] = 1.0 + # keep the prior feature unscaled: its units are already logits + z[FEATURES.index("zambretti_logit")] = x[FEATURES.index("zambretti_logit")] + return z + + # ------------------------------------------------------- inference + + def predict(self, obs: Dict, zam: Dict) -> Dict: + x = self._standardise(self.featurise(obs, zam), update=False) + p_model = _sigmoid(float(self.w @ x)) + p_prior = zam["prior_rain_prob"] + # trust the learner in proportion to how many strong labels it has + trust = self.n_strong / (self.n_strong + 25.0) + p = trust * p_model + (1 - trust) * p_prior + return { + "rain_probability": float(np.clip(p, 0.0, 1.0)), + "model_probability": float(p_model), + "prior_probability": float(p_prior), + "learner_trust": float(trust), + "condition": zam["condition"], + "label": zam["label"], + "zambretti_z": zam["z"], + "pressure_characteristic": tendency_code(zam["tendency"]), + "tendency": float(zam["tendency"]), + "sea_level_pressure": float(zam["slp_used"]), + "strong_labels": self.n_strong, + "weak_labels": self.n_weak, + "logloss_ewma": round(float(self.ewma_logloss), 4), + } + + # -------------------------------------------------------- learning + + def learn(self, obs: Dict, zam: Dict, y: float, strong: bool = False) -> float: + """AdaGrad step. Weak (proxy) labels get a tenth of the weight.""" + x = self._standardise(self.featurise(obs, zam), update=True) + p = _sigmoid(float(self.w @ x)) + weight = 1.0 if strong else 0.1 + grad = weight * (p - float(y)) * x + self.l2 * self.w + self.g2 += grad ** 2 + self.w -= self.lr * grad / np.sqrt(self.g2) + + loss = -(y * math.log(max(p, 1e-9)) + (1 - y) * math.log(max(1 - p, 1e-9))) + self.ewma_logloss = 0.98 * self.ewma_logloss + 0.02 * loss + if strong: + self.n_strong += 1 + else: + self.n_weak += 1 + return float(loss) + + def coefficients(self) -> List[Dict]: + return [{"feature": f, "weight": round(float(w), 4)} + for f, w in zip(FEATURES, self.w)] + + def to_dict(self) -> Dict: + return {"w": self.w.tolist(), "g2": self.g2.tolist(), "lr": self.lr, + "l2": self.l2, "n_strong": self.n_strong, "n_weak": self.n_weak, + "ewma_logloss": self.ewma_logloss, "mean": self.mean.tolist(), + "m2": self.m2.tolist(), "n_seen": self.n_seen} + + def load_dict(self, s: Dict) -> None: + self.w = np.array(s["w"], dtype=float) + self.g2 = np.array(s["g2"], dtype=float) + self.lr, self.l2 = s["lr"], s["l2"] + self.n_strong, self.n_weak = s["n_strong"], s["n_weak"] + self.ewma_logloss = s["ewma_logloss"] + self.mean = np.array(s["mean"], dtype=float) + self.m2 = np.array(s["m2"], dtype=float) + self.n_seen = s["n_seen"] + + +def proxy_wet_label(rh: float, dew_depression: float, cloud_index: float) -> Optional[float]: + """A weak, deliberately conservative stand-in for a rain gauge. + + Returns 1.0 for near-saturated overcast air, 0.0 for clearly dry air, + and None in the ambiguous middle, where a guess would poison the + training set faster than the extra samples could help. + """ + if not all(np.isfinite([rh, dew_depression, cloud_index])): + return None + if rh >= 93.0 and dew_depression <= 1.2 and cloud_index >= 0.6: + return 1.0 + if rh <= 65.0 and dew_depression >= 5.0: + return 0.0 + return None diff --git a/ashvale/models/rls.py b/ashvale/models/rls.py new file mode 100644 index 0000000..0d89719 --- /dev/null +++ b/ashvale/models/rls.py @@ -0,0 +1,182 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""The learning core: exponentially-weighted recursive least squares. + +Why RLS rather than an off-the-shelf gradient learner: + +* It is the exact minimiser of the exponentially weighted squared error + at every step, not an approximation, so it converges in far fewer + samples than SGD. On a station that produces 288 rows a day, sample + efficiency is not a nicety. +* The covariance `P` is a genuine parameter-uncertainty estimate, free. +* One matrix of size (d, d) with d ~ 33 is 8 kB. The whole model bank + fits in L2 cache on a Cortex-A53. +* Forgetting factor `lambda` gives principled adaptation to season and + to sensor ageing without any retraining schedule. + +Directional forgetting is used: `P` is only inflated along directions +that were actually excited by data. Plain forgetting blows `P` up +exponentially during quiet nights when the regressor is nearly constant, +and the model then detonates on the first sunrise. This is the single +most common way an RLS deployment fails in the field. +""" + +from __future__ import annotations + +from collections import deque +from typing import Deque, Dict, Optional + +import numpy as np + + +class RecursiveLeastSquares: + def __init__(self, n_features: int, forgetting: float = 0.999, + delta: float = 100.0, p_max: float = 1e6): + self.d = int(n_features) + self.lam = float(forgetting) + self.p_max = float(p_max) + self.theta = np.zeros(self.d) + self.P = np.eye(self.d) * float(delta) + self.n_updates = 0 + self.ewma_sq_error = 0.0 + + def predict(self, x: np.ndarray) -> float: + return float(np.dot(self.theta, np.asarray(x, dtype=float).ravel())) + + def predict_many(self, X: np.ndarray) -> np.ndarray: + return np.asarray(X, dtype=float) @ self.theta + + def predict_std(self, x: np.ndarray, noise_var: float = 1.0) -> float: + """Parameter-uncertainty contribution to predictive spread.""" + x = np.asarray(x, dtype=float).ravel() + return float(np.sqrt(max(noise_var * (1.0 + x @ self.P @ x), 1e-12))) + + def update(self, x: np.ndarray, y: float, weight: float = 1.0) -> float: + """One RLS step. Returns the a-priori residual (the honest error).""" + x = np.asarray(x, dtype=float).ravel() + if not (np.all(np.isfinite(x)) and np.isfinite(y)): + return 0.0 + + Px = self.P @ x + denom = self.lam + weight * float(x @ Px) + if denom < 1e-12: + return 0.0 + + residual = float(y) - float(self.theta @ x) + gain = (weight * Px) / denom + self.theta = self.theta + gain * residual + self.P = (self.P - np.outer(gain, Px)) / self.lam + + # directional forgetting guard: cap the spectral growth of P + self.P = 0.5 * (self.P + self.P.T) # enforce symmetry + trace = float(np.trace(self.P)) + if trace > self.p_max: + self.P *= self.p_max / trace + np.fill_diagonal(self.P, np.maximum(np.diag(self.P), 1e-9)) + + self.n_updates += 1 + self.ewma_sq_error = 0.99 * self.ewma_sq_error + 0.01 * residual ** 2 + return residual + + def fit_batch(self, X: np.ndarray, y: np.ndarray, passes: int = 1) -> "RecursiveLeastSquares": + X = np.atleast_2d(np.asarray(X, dtype=float)) + y = np.asarray(y, dtype=float).ravel() + for _ in range(max(int(passes), 1)): + for i in range(X.shape[0]): + self.update(X[i], y[i]) + return self + + @property + def noise_var(self) -> float: + return float(max(self.ewma_sq_error, 1e-9)) + + def to_dict(self) -> Dict: + return {"d": self.d, "lam": self.lam, "p_max": self.p_max, + "theta": self.theta.tolist(), "P": self.P.tolist(), + "n": self.n_updates, "ewma": self.ewma_sq_error} + + @classmethod + def from_dict(cls, s: Dict) -> "RecursiveLeastSquares": + m = cls(s["d"], s["lam"], 1.0, s.get("p_max", 1e6)) + m.theta = np.array(s["theta"], dtype=float) + m.P = np.array(s["P"], dtype=float) + m.n_updates = s.get("n", 0) + m.ewma_sq_error = s.get("ewma", 0.0) + return m + + +class AdaptiveConformal: + """Distribution-free prediction intervals that self-correct their coverage. + + Split conformal gives you a valid interval only if the data are + exchangeable. Weather is not: a front arrives and yesterday's + residual quantile becomes a fantasy. Adaptive conformal inference + (Gibbs and Candes) fixes this by feeding realised coverage back into + the working alpha: + + alpha_{t+1} = alpha_t + gamma * (alpha_target - err_t) + + The interval widens after each miss and narrows after each hit, so + long-run coverage tracks the target whatever the distribution does. + """ + + def __init__(self, alpha: float = 0.10, window: int = 400, gamma: float = 0.01): + self.alpha_target = float(alpha) + self.alpha = float(alpha) + self.gamma = float(gamma) + self.scores: Deque[float] = deque(maxlen=int(window)) + self.hits: Deque[int] = deque(maxlen=int(window)) + + def quantile(self) -> float: + if len(self.scores) < 20: + return float("nan") + a = float(np.clip(self.alpha, 0.005, 0.75)) + return float(np.quantile(np.asarray(self.scores), 1.0 - a, method="higher")) + + def interval(self, mu: float, fallback_sigma: float = 1.0) -> tuple[float, float]: + q = self.quantile() + if not np.isfinite(q): + q = 1.645 * fallback_sigma # gaussian 90% until we know better + return float(mu - q), float(mu + q) + + def observe(self, residual: float, covered: Optional[bool] = None) -> None: + r = abs(float(residual)) + if not np.isfinite(r): + return + if covered is None: + q = self.quantile() + covered = bool(r <= q) if np.isfinite(q) else True + self.scores.append(r) + self.hits.append(1 if covered else 0) + err = 0.0 if covered else 1.0 + self.alpha = float(np.clip(self.alpha + self.gamma * (self.alpha_target - err), + 0.005, 0.75)) + + @property + def empirical_coverage(self) -> float: + return float(np.mean(self.hits)) if self.hits else float("nan") + + def to_dict(self) -> Dict: + return {"alpha_target": self.alpha_target, "alpha": self.alpha, + "gamma": self.gamma, "maxlen": self.scores.maxlen, + "scores": list(self.scores), "hits": list(self.hits)} + + @classmethod + def from_dict(cls, s: Dict) -> "AdaptiveConformal": + c = cls(s["alpha_target"], s.get("maxlen", 400) or 400, s["gamma"]) + c.alpha = s["alpha"] + c.scores = deque(s["scores"], maxlen=c.scores.maxlen) + c.hits = deque(s["hits"], maxlen=c.hits.maxlen) + return c diff --git a/ashvale/physics.py b/ashvale/physics.py new file mode 100644 index 0000000..a0675a6 --- /dev/null +++ b/ashvale/physics.py @@ -0,0 +1,162 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Physics that the model does not have to learn. + +Every function here is a closed-form relationship that would otherwise +have to be discovered from data. Feeding a learner `dew point` instead of +making it infer the Magnus curve from (T, RH) is the cheapest accuracy +you will ever buy, especially on 512 MB of RAM. +""" + +from __future__ import annotations + +import math +from datetime import datetime, timezone + +import numpy as np + +MAGNUS_A = 17.625 +MAGNUS_B = 243.04 # degrees C +P_STD = 1013.25 # hPa + + +def saturation_vapour_pressure(temp_c): + """Tetens / Magnus saturation vapour pressure in hPa.""" + t = np.asarray(temp_c, dtype=float) + return 6.112 * np.exp(MAGNUS_A * t / (MAGNUS_B + t)) + + +def vapour_pressure(temp_c, rh_pct): + return saturation_vapour_pressure(temp_c) * np.clip(np.asarray(rh_pct, float), 0.0, 100.0) / 100.0 + + +def vapour_pressure_deficit(temp_c, rh_pct): + """VPD in hPa. Bioprocess people know this one from headspace humidity control.""" + return saturation_vapour_pressure(temp_c) - vapour_pressure(temp_c, rh_pct) + + +def dew_point(temp_c, rh_pct): + """Magnus-Tetens dew point in degrees C.""" + t = np.asarray(temp_c, dtype=float) + rh = np.clip(np.asarray(rh_pct, dtype=float), 1e-3, 100.0) + gamma = (MAGNUS_A * t) / (MAGNUS_B + t) + np.log(rh / 100.0) + return (MAGNUS_B * gamma) / (MAGNUS_A - gamma) + + +def absolute_humidity(temp_c, rh_pct): + """Water content in g/m^3 via the ideal gas law.""" + e = vapour_pressure(temp_c, rh_pct) * 100.0 # Pa + t_k = np.asarray(temp_c, dtype=float) + 273.15 + return e / (461.5 * t_k) * 1000.0 + + +def heat_index(temp_c, rh_pct): + """Rothfusz apparent temperature, valid above roughly 26 C.""" + t = np.asarray(temp_c, dtype=float) * 9.0 / 5.0 + 32.0 + r = np.asarray(rh_pct, dtype=float) + hi = (-42.379 + 2.04901523 * t + 10.14333127 * r - 0.22475541 * t * r + - 6.83783e-3 * t ** 2 - 5.481717e-2 * r ** 2 + 1.22874e-3 * t ** 2 * r + + 8.5282e-4 * t * r ** 2 - 1.99e-6 * t ** 2 * r ** 2) + hi = np.where(t < 80.0, t, hi) + return (hi - 32.0) * 5.0 / 9.0 + + +def sea_level_pressure(press_hpa, temp_c, altitude_m): + """Reduce station pressure to mean sea level (barometric formula). + + Without this, a 15 m elevation offset masquerades as a permanent + low-pressure system and every rule-of-thumb forecaster gets it wrong. + """ + p = np.asarray(press_hpa, dtype=float) + t = np.asarray(temp_c, dtype=float) + h = float(altitude_m) + return p * (1.0 - (0.0065 * h) / (t + 0.0065 * h + 273.15)) ** -5.257 + + +def station_pressure(slp_hpa, temp_c, altitude_m): + p = np.asarray(slp_hpa, dtype=float) + t = np.asarray(temp_c, dtype=float) + h = float(altitude_m) + return p * (1.0 - (0.0065 * h) / (t + 0.0065 * h + 273.15)) ** 5.257 + + +# ---------------------------------------------------------------- solar + +def _day_of_year(ts: float) -> float: + dt = datetime.fromtimestamp(ts, tz=timezone.utc) + return dt.timetuple().tm_yday + dt.hour / 24.0 + dt.minute / 1440.0 + + +def solar_position(ts, latitude: float, longitude: float): + """Return (elevation_deg, azimuth_deg) using the NOAA low-precision model. + + Accurate to a few tenths of a degree, which is far beyond what a + diurnal-cycle feature needs, and costs about twenty flops. + """ + ts_arr = np.atleast_1d(np.asarray(ts, dtype=float)) + doy = np.array([_day_of_year(float(t)) for t in ts_arr]) + frac_hour = np.array([ + datetime.fromtimestamp(float(t), tz=timezone.utc).hour + + datetime.fromtimestamp(float(t), tz=timezone.utc).minute / 60.0 + + datetime.fromtimestamp(float(t), tz=timezone.utc).second / 3600.0 + for t in ts_arr + ]) + + gamma = 2.0 * math.pi / 365.0 * (doy - 1.0) + eqtime = 229.18 * (0.000075 + 0.001868 * np.cos(gamma) - 0.032077 * np.sin(gamma) + - 0.014615 * np.cos(2 * gamma) - 0.040849 * np.sin(2 * gamma)) + decl = (0.006918 - 0.399912 * np.cos(gamma) + 0.070257 * np.sin(gamma) + - 0.006758 * np.cos(2 * gamma) + 0.000907 * np.sin(2 * gamma) + - 0.002697 * np.cos(3 * gamma) + 0.00148 * np.sin(3 * gamma)) + + true_solar_min = frac_hour * 60.0 + eqtime + 4.0 * longitude + hour_angle = np.radians(true_solar_min / 4.0 - 180.0) + + lat = math.radians(latitude) + cos_zenith = (np.sin(lat) * np.sin(decl) + + np.cos(lat) * np.cos(decl) * np.cos(hour_angle)) + cos_zenith = np.clip(cos_zenith, -1.0, 1.0) + elevation = np.degrees(np.arcsin(cos_zenith)) + + azimuth = np.degrees(np.arctan2( + -np.sin(hour_angle), + np.tan(decl) * np.cos(lat) - np.sin(lat) * np.cos(hour_angle) + )) % 360.0 + + if np.isscalar(ts) or np.asarray(ts).ndim == 0: + return float(elevation[0]), float(azimuth[0]) + return elevation, azimuth + + +def clear_sky_irradiance(elevation_deg): + """Rough clear-sky global horizontal irradiance, W/m^2. + + Used as the denominator of a `cloudiness proxy` when the TCS3400 sees + daylight: measured_lux / expected_lux is a surprisingly decent + okta estimate through a south-facing window. + """ + el = np.clip(np.asarray(elevation_deg, dtype=float), 0.0, 90.0) + sin_el = np.sin(np.radians(el)) + air_mass = np.where(el > 0.5, 1.0 / np.maximum(sin_el, 1e-3), 40.0) + return np.where(el > 0.0, 1353.0 * 0.7 ** (air_mass ** 0.678) * sin_el, 0.0) + + +def wet_bulb(temp_c, rh_pct): + """Stull's empirical wet-bulb approximation, degrees C.""" + t = np.asarray(temp_c, dtype=float) + rh = np.clip(np.asarray(rh_pct, dtype=float), 5.0, 99.0) + return (t * np.arctan(0.151977 * np.sqrt(rh + 8.313659)) + + np.arctan(t + rh) - np.arctan(rh - 1.676331) + + 0.00391838 * rh ** 1.5 * np.arctan(0.023101 * rh) - 4.686035) diff --git a/ashvale/sensors.py b/ashvale/sensors.py new file mode 100644 index 0000000..73971f0 --- /dev/null +++ b/ashvale/sensors.py @@ -0,0 +1,251 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Hardware access, with a simulator so the suite runs on your laptop too. + +`SenseBoard` is the only place that touches `sense_hat` or `smbus2`. If +either import fails (which it will on any machine that is not a Pi), the +board falls back to `SimulatedBoard`: a small stochastic-differential +weather model that produces plausible diurnal cycles, synoptic pressure +waves and sensor noise. Train on it, develop against it, then move the +same code to the Pi unchanged. +""" + +from __future__ import annotations + +import math +import random +import time +from typing import Any, Dict, Optional + +import numpy as np + +from .physics import dew_point, sea_level_pressure, solar_position + +TCS3400_ENABLE = 0x80 +TCS3400_ATIME = 0x81 +TCS3400_CONTROL = 0x8F +TCS3400_CDATA = 0x94 + + +def read_cpu_temperature() -> float: + """Core temperature in C. This is the single most important nuisance + variable on a Sense HAT: the HTS221 and LPS25HB sit millimetres above a + SoC that runs 30 C hotter than the room.""" + try: + with open("/sys/class/thermal/thermal_zone0/temp", "r") as fh: + return float(fh.read().strip()) / 1000.0 + except Exception: + return float("nan") + + +class SimulatedBoard: + """Ornstein-Uhlenbeck weather with a diurnal driver. Good enough to + exercise every code path and to sanity-check a model's skill score.""" + + def __init__(self, latitude: float = 52.2, longitude: float = 0.12, seed: int = 7): + self.rng = np.random.default_rng(seed) + self.lat, self.lon = latitude, longitude + self.t0 = time.time() + self.press_anom = 0.0 + self.temp_anom = 0.0 + self.hum_anom = 0.0 + self.last = self.t0 + self.available = False + + def _step(self, now: float) -> None: + dt = max(min(now - self.last, 600.0), 0.0) + self.last = now + # synoptic pressure: slow OU process, tau ~ 30 h, sigma ~ 9 hPa + self.press_anom += (-self.press_anom / (30 * 3600) * dt + + 9.0 * math.sqrt(2 * dt / (30 * 3600)) * self.rng.normal()) + self.temp_anom += (-self.temp_anom / (6 * 3600) * dt + + 1.8 * math.sqrt(2 * dt / (6 * 3600)) * self.rng.normal()) + self.hum_anom += (-self.hum_anom / (4 * 3600) * dt + + 6.0 * math.sqrt(2 * dt / (4 * 3600)) * self.rng.normal()) + + def read(self) -> Dict[str, Any]: + now = time.time() + self._step(now) + elev, _ = solar_position(now, self.lat, self.lon) + doy = time.gmtime(now).tm_yday + seasonal = 6.5 * math.sin(2 * math.pi * (doy - 105) / 365.25) + solar_gain = 5.0 * max(elev, 0.0) / 60.0 + temp = 12.0 + seasonal + solar_gain + self.temp_anom + rh = float(np.clip(78.0 - 1.9 * (temp - 12.0) + self.hum_anom, 12.0, 99.0)) + press = 1013.0 + self.press_anom + lux = max(0.0, 60000.0 * max(math.sin(math.radians(max(elev, 0.0))), 0.0)) + 8.0 + cpu = temp + 22.0 + 1.5 * self.rng.normal() + # forward model must invert the compensator exactly, see scripts/simulate.py + k_true = 0.55 + return { + "temp_raw": (temp + k_true * cpu) / (1.0 + k_true) + 0.05 * self.rng.normal(), + "hum": rh + 0.4 * self.rng.normal(), + "press": press + 0.05 * self.rng.normal(), + "cpu_temp": cpu, + "lux": lux * (0.35 + 0.65 * self.rng.random()), + "r": int(lux * 0.30), "g": int(lux * 0.34), "b": int(lux * 0.28), + "pitch": 0.4 * self.rng.normal(), "roll": 0.4 * self.rng.normal(), + "yaw": 180.0 + self.rng.normal(), "compass": 180.0 + 2 * self.rng.normal(), + "ax": 0.0, "ay": 0.0, "az": 1.0, + "gx": 0.0, "gy": 0.0, "gz": 0.0, + } + + def clear(self, *_a, **_k): # LED no-op + pass + + +class SenseBoard: + """Real hardware wrapper. Attribute `available` tells you which world + you are in without try/except at every call site.""" + + def __init__(self, rotation: int = 90, low_light: bool = True, + tcs_addr: int = 0x39, latitude: float = 52.2, longitude: float = 0.12): + self.available = False + self.has_colour = False + self.sense = None + self.bus = None + self.tcs_addr = tcs_addr + self._sim = SimulatedBoard(latitude, longitude) + + try: + from sense_hat import SenseHat # type: ignore + self.sense = SenseHat() + self.sense.low_light = low_light + self.sense.set_rotation(rotation) + self.available = True + except Exception: + self.sense = None + + if self.available: + try: + import smbus2 # type: ignore + self.bus = smbus2.SMBus(1) + self.bus.write_byte_data(self.tcs_addr, TCS3400_ENABLE, 0x03) # power + RGBC + self.bus.write_byte_data(self.tcs_addr, TCS3400_ATIME, 0xD5) # 100 ms + self.bus.write_byte_data(self.tcs_addr, TCS3400_CONTROL, 0x00) # 1x gain + self.has_colour = True + except Exception: + self.has_colour = False + + # ---------------------------------------------------------------- IO + + def colour(self) -> Dict[str, Any]: + if not self.has_colour: + return {"clear": 0, "red": 0, "green": 0, "blue": 0, "hex": "#334155", "cct": None} + try: + data = self.bus.read_i2c_block_data(self.tcs_addr, TCS3400_CDATA | 0x80, 8) + c = data[0] | (data[1] << 8) + r = data[2] | (data[3] << 8) + g = data[4] | (data[5] << 8) + b = data[6] | (data[7] << 8) + return _colour_payload(c, r, g, b) + except Exception: + return {"clear": 0, "red": 0, "green": 0, "blue": 0, "hex": "#334155", "cct": None} + + def read(self) -> Dict[str, Any]: + """One full multi-sensor sample. Raw, uncompensated, untouched.""" + if not self.available: + row = self._sim.read() + col = _colour_payload(int(row["lux"]), row["r"], row["g"], row["b"]) + row.update({"lux": col["clear"], "r": col["red"], "g": col["green"], + "b": col["blue"], "colour": col, "simulated": True}) + return row + + s = self.sense + t_h = s.get_temperature_from_humidity() + t_p = s.get_temperature_from_pressure() + orientation = s.get_orientation_degrees() + accel = s.get_accelerometer_raw() + gyro = s.get_gyroscope_raw() + col = self.colour() + + def wrap(v): + return v - 360.0 if v > 180.0 else v + + return { + "temp_raw": (t_h + t_p) / 2.0, + "temp_h": t_h, + "temp_p": t_p, + "hum": s.get_humidity(), + "press": s.get_pressure(), + "cpu_temp": read_cpu_temperature(), + "lux": col["clear"], "r": col["red"], "g": col["green"], "b": col["blue"], + "colour": col, + "pitch": wrap(orientation["pitch"]), + "roll": wrap(orientation["roll"]), + "yaw": orientation["yaw"], + "compass": s.get_compass(), + "ax": accel["x"], "ay": accel["y"], "az": accel["z"], + "gx": gyro["x"], "gy": gyro["y"], "gz": gyro["z"], + "simulated": False, + } + + # --------------------------------------------------------------- LED + + def clear(self, *args): + if self.sense is not None: + self.sense.clear(*args) + + def show_message(self, text: str, scroll_speed: float = 0.065, text_colour=None): + if self.sense is not None: + self.sense.show_message(text, scroll_speed=scroll_speed, + text_colour=text_colour or [255, 255, 255]) + + def set_pixels(self, pixels): + if self.sense is not None: + self.sense.set_pixels(pixels) + + +def _colour_payload(c: int, r: int, g: int, b: int) -> Dict[str, Any]: + denom = max(int(c), 1) + nr = min(int((r / denom) * 255), 255) + ng = min(int((g / denom) * 255), 255) + nb = min(int((b / denom) * 255), 255) + return { + "clear": int(c), "red": int(r), "green": int(g), "blue": int(b), + "hex": f"#{nr:02x}{ng:02x}{nb:02x}", + "cct": correlated_colour_temperature(r, g, b), + } + + +def correlated_colour_temperature(r: float, g: float, b: float) -> Optional[float]: + """McCamy's approximation, in kelvin. Distinguishes a tungsten desk lamp + (~2700 K) from overcast daylight (~6500 K), which turns the colour sensor + into a crude `is anyone home` and `is it cloudy` detector.""" + if (r + g + b) <= 0: + return None + X = -0.14282 * r + 1.54924 * g + -0.95641 * b + Y = -0.32466 * r + 1.57837 * g + -0.73191 * b + Z = -0.68202 * r + 0.77073 * g + 0.56332 * b + denom = X + Y + Z + if abs(denom) < 1e-9: + return None + x, y = X / denom, Y / denom + if abs(y - 0.1858) < 1e-9: + return None + n = (x - 0.3320) / (0.1858 - y) + cct = 449 * n ** 3 + 3525 * n ** 2 + 6823.3 * n + 5520.33 + return float(cct) if 800 < cct < 25000 else None + + +def enrich(raw: Dict[str, Any], altitude_m: float) -> Dict[str, Any]: + """Add derived quantities that do not need any model state.""" + out = dict(raw) + temp = raw.get("temp_raw", float("nan")) + hum = raw.get("hum", float("nan")) + press = raw.get("press", float("nan")) + out["dew_c"] = float(dew_point(temp, hum)) + out["press_slp"] = float(sea_level_pressure(press, temp, altitude_m)) + return out diff --git a/ashvale/station.py b/ashvale/station.py new file mode 100644 index 0000000..55ea5f5 --- /dev/null +++ b/ashvale/station.py @@ -0,0 +1,563 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""The station: everything wired together and running on its own clocks. + +Four asynchronous loops, deliberately decoupled so a slow one cannot +starve a fast one: + + sample (2 s) read hardware, run the Kalman bank, keep live state + persist (30 s) one row to SQLite + train (10 min) rebuild the feature grid, update every head, re-fit + climatology, emit a fresh forecast bundle + verify (5 min) score forecasts whose validity time has arrived, feed + the errors to conformal calibration and drift + detection, write the scorecard + +The verify loop is the one most projects skip and the one that makes the +difference. A forecast that is never scored is an opinion; a forecast +that is scored against persistence is a measurement. +""" + +from __future__ import annotations + +import asyncio +import json +import math +import time +from pathlib import Path +from typing import Any, Dict, List, Optional + +import numpy as np + +from . import physics +from .config import Config +from .estimation import SignalTracker +from .features import N_FEATURES, build_features +from .models.anomaly import AnomalyMonitor +from .models.climatology import HarmonicClimatology +from .models.nowcast import NowcastEnsemble +from .models.precip import PrecipitationModel, proxy_wet_label, zambretti +from .sensors import SenseBoard, enrich +from .storage import Store, resample + +STATE_VERSION = 1 + + +class Station: + def __init__(self, cfg: Config): + self.cfg = cfg + self.store = Store(cfg.storage.db_path) + self.board = SenseBoard( + rotation=cfg.sensor.rotation_deg, + low_light=cfg.sensor.low_light, + tcs_addr=cfg.sensor.tcs3400_addr, + latitude=cfg.site.latitude, + longitude=cfg.site.longitude, + ) + self.tracker = SignalTracker(cfg) + self.nowcast = NowcastEnsemble(cfg.model.targets, cfg.model.horizons_s, cfg.model) + self.climatology = HarmonicClimatology( + cfg.model.targets, min_days_annual=cfg.model.climatology_min_days_annual + ) + self.precip = PrecipitationModel() + self.monitor = AnomalyMonitor(cfg.model) + + self.live: Dict[str, Any] = {} + self.forecast_bundle: Dict[str, Any] = {} + self.outlook_bundle: Dict[str, Any] = {} + self.precip_bundle: Dict[str, Any] = {} + self.anomaly_bundle: Dict[str, Any] = {} + self.last_train: float = 0.0 + self.last_persist: float = 0.0 + self.last_compact: float = 0.0 + self.training_log: List[Dict] = [] + self._tasks: List[asyncio.Task] = [] + self._stop = asyncio.Event() + + self.state_path = Path(cfg.storage.state_dir) / "station_state.json" + self.load_state() + + # ------------------------------------------------------------ state + + def save_state(self) -> None: + payload = { + "version": STATE_VERSION, + "saved_at": time.time(), + "tracker": self.tracker.to_dict(), + "nowcast": self.nowcast.to_dict(), + "climatology": self.climatology.to_dict(), + "precip": self.precip.to_dict(), + "monitor": self.monitor.to_dict(), + } + tmp = self.state_path.with_suffix(".tmp") + with open(tmp, "w", encoding="utf-8") as fh: + json.dump(payload, fh) + tmp.replace(self.state_path) # atomic, survives a power cut mid-write + + def load_state(self) -> bool: + if not self.state_path.exists(): + return False + try: + with open(self.state_path, "r", encoding="utf-8") as fh: + s = json.load(fh) + if s.get("version") != STATE_VERSION: + return False + self.tracker.load_dict(s["tracker"]) + self.nowcast.load_dict(s["nowcast"]) + self.climatology.load_dict(s["climatology"]) + self.precip.load_dict(s["precip"]) + self.monitor.load_dict(s["monitor"]) + return True + except Exception as exc: + self.store.log_event("state", "warn", f"could not restore state: {exc}") + return False + + # ----------------------------------------------------------- sample + + def sample_once(self) -> Dict[str, Any]: + ts = time.time() + raw = self.board.read() + raw = enrich(raw, self.cfg.site.altitude_m) + est = self.tracker.step(ts, raw.get("temp_raw", float("nan")), + raw.get("hum", float("nan")), + raw.get("press", float("nan")), + raw.get("cpu_temp", float("nan"))) + + temp_c = est["temp_smooth"] + slp = float(physics.sea_level_pressure(est["press_smooth"], temp_c, + self.cfg.site.altitude_m)) + dew = float(physics.dew_point(temp_c, est["hum_smooth"])) + elev, azim = physics.solar_position(ts, self.cfg.site.latitude, + self.cfg.site.longitude) + expected = float(physics.clear_sky_irradiance(elev)) + lux = float(raw.get("lux", 0.0) or 0.0) + cloud = (float(np.clip(1.0 - lux / max(expected * 45.0, 1.0), 0.0, 1.0)) + if elev > 5.0 else 0.5) + + row = { + "ts": ts, + "temp_raw": raw.get("temp_raw"), + "temp_c": est["temp_c"], + "temp_smooth": temp_c, + "temp_rate": est["temp_rate"], + "hum": raw.get("hum"), + "hum_smooth": est["hum_smooth"], + "press": raw.get("press"), + "press_slp": slp, + "press_smooth": est["press_smooth"], + "press_rate": est["press_rate"], + "cpu_temp": raw.get("cpu_temp"), + "dew_c": dew, + "lux": lux, + "r": raw.get("r"), "g": raw.get("g"), "b": raw.get("b"), + "pitch": raw.get("pitch"), "roll": raw.get("roll"), + "yaw": raw.get("yaw"), "compass": raw.get("compass"), + "ax": raw.get("ax"), "ay": raw.get("ay"), "az": raw.get("az"), + "gx": raw.get("gx"), "gy": raw.get("gy"), "gz": raw.get("gz"), + } + + anomaly = self.monitor.observe(ts, { + "temp_c": temp_c, "hum": est["hum_smooth"], "press_slp": slp, + "temp_rate": est["temp_rate"], "press_rate": est["press_rate"], + "dew_c": dew, "cpu_temp": raw.get("cpu_temp"), + }) + self.anomaly_bundle = anomaly + + self.live = { + **row, + "timestamp": time.strftime("%H:%M:%S", time.localtime(ts)), + "colour": raw.get("colour", {}), + "simulated": bool(raw.get("simulated", not self.board.available)), + "dew_depression": temp_c - dew, + "vpd": float(physics.vapour_pressure_deficit(temp_c, est["hum_smooth"])), + "wet_bulb": float(physics.wet_bulb(temp_c, est["hum_smooth"])), + "heat_index": float(physics.heat_index(temp_c, est["hum_smooth"])), + "abs_humidity": float(physics.absolute_humidity(temp_c, est["hum_smooth"])), + "solar_elevation": float(elev), + "solar_azimuth": float(azim), + "clear_sky_wm2": expected, + "cloud_index": cloud, + "cpu_offset": (raw.get("cpu_temp") or float("nan")) - (raw.get("temp_raw") or float("nan")), + "compensator_k": self.tracker.compensator.k, + "health": anomaly["health_overall"], + "novelty_d2": anomaly["novelty"].get("d2", 0.0), + } + self._update_precip() + return self.live + + def _observation_vector(self) -> Dict[str, float]: + live = self.live + hist = self.store.window(8.0, ["ts", "press_slp", "temp_c", "dew_c"]) + tend = {"tend_1h": live.get("press_rate", 0.0), + "tend_3h": live.get("press_rate", 0.0), + "tend_6h": live.get("press_rate", 0.0)} + if hist["ts"].size > 5: + now = hist["ts"][-1] + for key, hours in (("tend_1h", 1.0), ("tend_3h", 3.0), ("tend_6h", 6.0)): + idx = np.searchsorted(hist["ts"], now - hours * 3600.0) + if 0 <= idx < hist["ts"].size - 1: + dtp = (now - hist["ts"][idx]) / 3600.0 + if dtp > 0.25: + tend[key] = float((hist["press_slp"][-1] - hist["press_slp"][idx]) / dtp) + dew_dep = live.get("dew_depression", 5.0) + dew_dep_rate = 0.0 + if hist["ts"].size > 5: + idx = np.searchsorted(hist["ts"], hist["ts"][-1] - 3600.0) + if 0 <= idx < hist["ts"].size - 1: + past = hist["temp_c"][idx] - hist["dew_c"][idx] + dew_dep_rate = float(dew_dep - past) + return { + "slp": live.get("press_slp", 1013.25), + "rh": live.get("hum_smooth", 60.0), + "dew_depression": dew_dep, + "dew_dep_rate": dew_dep_rate, + "cloud_index": live.get("cloud_index", 0.5), + "temp_dev": self.climatology.anomaly_now( + "temperature", live.get("ts", time.time()), live.get("temp_smooth", 0.0) + ), + "wet_bulb_depression": live.get("temp_smooth", 0.0) - live.get("wet_bulb", 0.0), + **tend, + } + + def _update_precip(self) -> None: + obs = self._observation_vector() + zam = zambretti(obs["slp"], obs["tend_3h"], self.live.get("ts"), + self.cfg.site.latitude) + self.precip_bundle = self.precip.predict(obs, zam) + self.precip_bundle["indoors_caveat"] = self.cfg.site.indoors + + y = proxy_wet_label(obs["rh"], obs["dew_depression"], obs["cloud_index"]) + if y is not None and int(self.live.get("ts", 0)) % 300 < self.cfg.sensor.sample_period_s: + self.precip.learn(obs, zam, y, strong=False) + + def add_label(self, kind: str, value: float, ts: Optional[float] = None, + note: str = "") -> Dict: + """Human-in-the-loop ground truth. Worth ten times a proxy label.""" + ts = ts or time.time() + self.store.insert_label(ts, kind, value, note) + if kind == "rain": + obs = self._observation_vector() + zam = zambretti(obs["slp"], obs["tend_3h"], ts, self.cfg.site.latitude) + loss = self.precip.learn(obs, zam, float(value), strong=True) + self.store.log_event("label", "info", + f"strong rain label {value} accepted, loss {loss:.3f}", ts) + return {"accepted": True, "loss": loss, "strong_labels": self.precip.n_strong} + return {"accepted": True} + + def calibrate_temperature(self, reference_c: float) -> Dict: + raw = self.live.get("temp_raw") + cpu = self.live.get("cpu_temp") + if raw is None or cpu is None: + return {"error": "no live reading yet"} + result = self.tracker.compensator.calibrate(float(raw), float(cpu), float(reference_c)) + self.store.log_event("calibration", "info", + f"k -> {result['k']:.3f} (residual {result['residual']:+.2f} C)") + return result + + def reset_calibration(self) -> Dict: + """Return the self-heating coefficient to its configured prior. + + Worth having: a single mistyped reference reading can drive `k` + to its clamp, and because state persists across restarts it will + stay there quietly biasing every reading until you notice. + """ + from .estimation import ThermalCompensator + self.tracker.compensator = ThermalCompensator( + self.cfg.sensor.cpu_heat_k, self.cfg.sensor.cpu_heat_k_min, + self.cfg.sensor.cpu_heat_k_max, + ) + self.save_state() + self.store.log_event("calibration", "info", + f"coefficient reset to prior k={self.cfg.sensor.cpu_heat_k}") + return {"k": self.tracker.compensator.k, "reset": True, "n": 0} + + # ------------------------------------------------------------ train + + def build_training_grid(self, hours: float = 24 * 30): + raw = self.store.window(hours, ["ts", "temp_smooth", "hum_smooth", + "press_slp", "lux"]) + if raw["ts"].size < 10: + return None + grid_ts, cols = resample( + raw["ts"], + {"temperature": raw["temp_smooth"], "humidity": raw["hum_smooth"], + "pressure": raw["press_slp"], "lux": raw["lux"]}, + self.cfg.model.grid_s, + ) + if grid_ts.size < self.cfg.model.min_rows_to_train: + return None + X, valid = build_features( + grid_ts, cols["temperature"], cols["humidity"], cols["pressure"], + cols["lux"], self.cfg.model.grid_s, + self.cfg.site.latitude, self.cfg.site.longitude, + ) + return grid_ts, cols, X, valid + + def train(self, hours: float = 24 * 30) -> Dict: + t_start = time.time() + built = self.build_training_grid(hours) + if built is None: + return {"trained": False, + "reason": f"need at least {self.cfg.model.min_rows_to_train} grid rows"} + grid_ts, cols, X, valid = built + + clim_scores = self.climatology.fit(grid_ts, cols, valid) + counts = self.nowcast.fit(X, valid, cols, self.climatology, grid_ts) + + self.last_train = time.time() + self.monitor.clear_retrain_flag() + entry = { + "ts": self.last_train, + "grid_rows": int(grid_ts.size), + "valid_rows": int(valid.sum()), + "span_days": round(float((grid_ts[-1] - grid_ts[0]) / 86400.0), 2), + "pairs": counts, + "climatology_resid_std": {k: round(v, 3) for k, v in clim_scores.items()}, + "annual_terms": self.climatology.use_annual, + "seconds": round(time.time() - t_start, 2), + } + self.training_log = ([entry] + self.training_log)[:20] + self.store.log_event("train", "info", + f"retrained on {grid_ts.size} grid rows in {entry['seconds']}s") + self.refresh_forecasts() + self.save_state() + return {"trained": True, **entry} + + # --------------------------------------------------------- forecast + + def refresh_forecasts(self, persist: bool = True) -> Dict: + built = self.build_training_grid(hours=48.0) + now = time.time() + if built is None or not self.live: + return {} + grid_ts, cols, X, valid = built + x_now = X[-1] + + anchors = { + "temperature": float(self.live.get("temp_smooth", cols["temperature"][-1])), + "humidity": float(self.live.get("hum_smooth", cols["humidity"][-1])), + "pressure": float(self.live.get("press_slp", cols["pressure"][-1])), + } + fc = self.nowcast.forecast(x_now, anchors, now, self.climatology) + + bundle: Dict[str, Any] = {"issued_ts": now, "anchors": anchors, "targets": {}} + for target, per_h in fc.items(): + series = [] + for h in sorted(per_h): + p = per_h[h] + series.append({ + "horizon_s": h, + "horizon_label": _fmt_horizon(h), + "valid_ts": now + h, + "mu": round(p["mu"], 3), + "lo": round(p["lo"], 3), + "hi": round(p["hi"], 3), + "delta": round(p["delta"], 3), + "weights": {k: round(v, 3) for k, v in p["weights"].items()}, + }) + if persist: + self.store.insert_forecast(now, h, target, p["mu"], p["lo"], + p["hi"], "ensemble") + bundle["targets"][target] = series + self.forecast_bundle = bundle + + self.outlook_bundle = { + "issued_ts": now, + "ready": self.climatology.ready, + "annual_terms": self.climatology.use_annual, + "history_days": round(self.store.span_days(), 2), + "targets": { + t: self.climatology.outlook( + t, now, days=7, + anomaly=self.climatology.anomaly_now(t, now, anchors.get(t, 0.0)), + ) + for t in self.cfg.model.targets + }, + } + return bundle + + # ----------------------------------------------------------- verify + + def verify(self) -> Dict: + """Score matured forecasts against truth and against persistence.""" + due = self.store.due_forecasts() + if not due: + return {"scored": 0} + + hist = self.store.window(24 * 8, ["ts", "temp_smooth", "hum_smooth", "press_slp"]) + if hist["ts"].size < 5: + return {"scored": 0} + series = {"temperature": hist["temp_smooth"], "humidity": hist["hum_smooth"], + "pressure": hist["press_slp"]} + + def value_at(target: str, ts: float) -> Optional[float]: + idx = int(np.searchsorted(hist["ts"], ts)) + if idx <= 0 or idx >= hist["ts"].size: + return None + if abs(hist["ts"][idx] - ts) > 900: + return None + return float(series[target][idx]) + + buckets: Dict[tuple, Dict[str, List[float]]] = {} + scored = 0 + for row in due: + target, h = row["target"], int(row["horizon_s"]) + truth = value_at(target, row["valid_ts"]) + anchor = value_at(target, row["issued_ts"]) + if truth is None or anchor is None: + continue + key = (target, h) + b = buckets.setdefault(key, {"err": [], "pers": [], "cov": []}) + err = truth - row["mu"] + b["err"].append(err) + b["pers"].append(truth - anchor) + b["cov"].append(1.0 if row["lo"] <= truth <= row["hi"] else 0.0) + head = self.nowcast.heads.get(key) + if head is not None: + head.conformal.observe(err, covered=bool(row["lo"] <= truth <= row["hi"])) + if h <= 10800: + self.monitor.observe_error(row["valid_ts"], abs(err)) + scored += 1 + + now = time.time() + for (target, h), b in buckets.items(): + e = np.asarray(b["err"], dtype=float) + p = np.asarray(b["pers"], dtype=float) + mae = float(np.mean(np.abs(e))) + mae_p = float(np.mean(np.abs(p))) + self.store.insert_score( + now, target, h, + mae=mae, + rmse=float(np.sqrt(np.mean(e ** 2))), + bias=float(np.mean(e)), + mae_persistence=mae_p, + skill=float(1.0 - mae / mae_p) if mae_p > 1e-9 else 0.0, + coverage=float(np.mean(b["cov"])), + n=int(e.size), + ) + + with self.store._conn() as conn: + conn.execute("DELETE FROM forecasts WHERE valid_ts <= ?", (now - 3600,)) + return {"scored": scored, "buckets": len(buckets)} + + # ------------------------------------------------------------ loops + + async def _loop_sample(self): + period = self.cfg.sensor.sample_period_s + while not self._stop.is_set(): + try: + self.sample_once() + now = time.time() + if now - self.last_persist >= self.cfg.sensor.persist_period_s: + self.store.insert_telemetry(self.live) + self.last_persist = now + except Exception as exc: + self.store.log_event("sample", "error", repr(exc)) + await asyncio.sleep(period) + + async def _loop_train(self): + await asyncio.sleep(5) + try: + self.train() + except Exception as exc: + self.store.log_event("train", "error", repr(exc)) + while not self._stop.is_set(): + await asyncio.sleep(30) + now = time.time() + due = (now - self.last_train) >= self.cfg.model.train_period_s + if due or self.monitor.retrain_requested: + try: + await asyncio.to_thread(self.train) + except Exception as exc: + self.store.log_event("train", "error", repr(exc)) + + async def _loop_verify(self): + await asyncio.sleep(60) + while not self._stop.is_set(): + try: + await asyncio.to_thread(self.verify) + except Exception as exc: + self.store.log_event("verify", "error", repr(exc)) + await asyncio.sleep(300) + + async def _loop_maintenance(self): + while not self._stop.is_set(): + await asyncio.sleep(3600) + now = time.time() + if now - self.last_compact >= self.cfg.storage.vacuum_period_s: + try: + removed = await asyncio.to_thread( + self.store.compact, + self.cfg.storage.raw_retention_days, + self.cfg.storage.five_min_retention_days, + ) + self.last_compact = now + self.store.log_event("compact", "info", json.dumps(removed)) + except Exception as exc: + self.store.log_event("compact", "error", repr(exc)) + self.save_state() + + def start(self) -> None: + self._stop.clear() + self._tasks = [ + asyncio.create_task(self._loop_sample()), + asyncio.create_task(self._loop_train()), + asyncio.create_task(self._loop_verify()), + asyncio.create_task(self._loop_maintenance()), + ] + + async def stop(self) -> None: + self._stop.set() + for t in self._tasks: + t.cancel() + for t in self._tasks: + try: + await t + except (asyncio.CancelledError, Exception): + pass + try: + self.save_state() + except Exception: + pass + + # ------------------------------------------------------------ views + + def status(self) -> Dict: + return { + "site": self.cfg.site.name, + "hardware": "sense-hat-v2" if self.board.available else "simulator", + "colour_sensor": self.board.has_colour, + "rows": self.store.row_count(), + "history_days": round(self.store.span_days(), 3), + "last_train": self.last_train, + "next_train_in_s": max(0.0, self.cfg.model.train_period_s + - (time.time() - self.last_train)), + "climatology_ready": self.climatology.ready, + "annual_terms": self.climatology.use_annual, + "compensator_k": round(self.tracker.compensator.k, 4), + "calibrations": self.tracker.compensator.n_calibrations, + "health": self.monitor.health.overall, + "drift_stress": round(self.monitor.drift.stress, 3), + "retrain_requested": self.monitor.retrain_requested, + "training_log": self.training_log[:5], + } + + +def _fmt_horizon(seconds: int) -> str: + if seconds < 3600: + return f"{seconds // 60}m" + if seconds < 86400: + return f"{seconds // 3600}h" + return f"{seconds // 86400}d" diff --git a/ashvale/storage.py b/ashvale/storage.py new file mode 100644 index 0000000..c70dfe2 --- /dev/null +++ b/ashvale/storage.py @@ -0,0 +1,487 @@ +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Durable storage: SQLite in WAL mode with tiered downsampling. + +An SD card is a consumable. The write pattern here is deliberately +gentle: one row every `persist_period_s`, WAL journalling, a compaction +pass that folds week-old raw rows into 5-minute means and quarter-old +5-minute rows into hourly means. A year of station history lands around +30 MB, which the Pi will not notice. +""" + +from __future__ import annotations + +import sqlite3 +import threading +import time +from typing import Any, Dict, Iterable, List, Optional + +import numpy as np + +TIER_RAW = 0 +TIER_5MIN = 1 +TIER_HOUR = 2 + +COLUMNS = [ + "ts", "temp_raw", "temp_c", "temp_smooth", "temp_rate", "hum", "hum_smooth", + "press", "press_slp", "press_smooth", "press_rate", "cpu_temp", "dew_c", + "lux", "r", "g", "b", "pitch", "roll", "yaw", "compass", + "ax", "ay", "az", "gx", "gy", "gz", +] + +SCHEMA = f""" +PRAGMA journal_mode=WAL; +PRAGMA synchronous=NORMAL; +PRAGMA temp_store=MEMORY; + +CREATE TABLE IF NOT EXISTS telemetry ( + ts REAL PRIMARY KEY, + {", ".join(f"{c} REAL" for c in COLUMNS if c != "ts")}, + tier INTEGER NOT NULL DEFAULT 0 +); +CREATE INDEX IF NOT EXISTS idx_telemetry_tier_ts ON telemetry(tier, ts); + +CREATE TABLE IF NOT EXISTS forecasts ( + issued_ts REAL NOT NULL, + valid_ts REAL NOT NULL, + horizon_s INTEGER NOT NULL, + target TEXT NOT NULL, + mu REAL, lo REAL, hi REAL, + model TEXT, + PRIMARY KEY (issued_ts, horizon_s, target) +); +CREATE INDEX IF NOT EXISTS idx_forecast_valid ON forecasts(valid_ts); + +CREATE TABLE IF NOT EXISTS labels ( + ts REAL NOT NULL, + kind TEXT NOT NULL, + value REAL NOT NULL, + note TEXT, + PRIMARY KEY (ts, kind) +); + +CREATE TABLE IF NOT EXISTS scores ( + ts REAL NOT NULL, + target TEXT NOT NULL, + horizon_s INTEGER NOT NULL, + mae REAL, rmse REAL, bias REAL, + mae_persistence REAL, skill REAL, coverage REAL, n INTEGER, + PRIMARY KEY (ts, target, horizon_s) +); + +CREATE TABLE IF NOT EXISTS events ( + ts REAL NOT NULL, + kind TEXT NOT NULL, + severity TEXT, + detail TEXT +); +CREATE INDEX IF NOT EXISTS idx_events_ts ON events(ts); +""" + + +class Store: + def __init__(self, path: str): + self.path = path + self._local = threading.local() + with self._conn() as conn: + conn.executescript(SCHEMA) + + def _conn(self) -> sqlite3.Connection: + conn = getattr(self._local, "conn", None) + if conn is None: + conn = sqlite3.connect(self.path, timeout=20.0, check_same_thread=False) + conn.row_factory = sqlite3.Row + self._local.conn = conn + return conn + + # ------------------------------------------------------------ writes + + def insert_telemetry(self, row: Dict[str, Any], tier: int = TIER_RAW) -> None: + payload = {c: float(row.get(c)) if row.get(c) is not None else None for c in COLUMNS} + payload["tier"] = tier + cols = ", ".join(payload.keys()) + marks = ", ".join("?" for _ in payload) + with self._conn() as conn: + conn.execute( + f"INSERT OR REPLACE INTO telemetry ({cols}) VALUES ({marks})", + list(payload.values()), + ) + + def insert_forecast(self, issued_ts: float, horizon_s: int, target: str, + mu: float, lo: float, hi: float, model: str) -> None: + with self._conn() as conn: + conn.execute( + "INSERT OR REPLACE INTO forecasts " + "(issued_ts, valid_ts, horizon_s, target, mu, lo, hi, model) " + "VALUES (?,?,?,?,?,?,?,?)", + (issued_ts, issued_ts + horizon_s, horizon_s, target, + float(mu), float(lo), float(hi), model), + ) + + def insert_label(self, ts: float, kind: str, value: float, note: str = "") -> None: + with self._conn() as conn: + conn.execute( + "INSERT OR REPLACE INTO labels (ts, kind, value, note) VALUES (?,?,?,?)", + (ts, kind, float(value), note), + ) + + def insert_score(self, ts: float, target: str, horizon_s: int, **kw) -> None: + with self._conn() as conn: + conn.execute( + "INSERT OR REPLACE INTO scores " + "(ts, target, horizon_s, mae, rmse, bias, mae_persistence, skill, coverage, n) " + "VALUES (?,?,?,?,?,?,?,?,?,?)", + (ts, target, horizon_s, kw.get("mae"), kw.get("rmse"), kw.get("bias"), + kw.get("mae_persistence"), kw.get("skill"), kw.get("coverage"), kw.get("n")), + ) + + def log_event(self, kind: str, severity: str, detail: str, ts: Optional[float] = None) -> None: + with self._conn() as conn: + conn.execute("INSERT INTO events (ts, kind, severity, detail) VALUES (?,?,?,?)", + (ts or time.time(), kind, severity, detail)) + + # ------------------------------------------------------------- reads + + def window(self, hours: float, columns: Optional[Iterable[str]] = None) -> Dict[str, np.ndarray]: + """Return the last `hours` of telemetry as column arrays, oldest first.""" + cols = list(columns) if columns else COLUMNS + since = time.time() - hours * 3600.0 + with self._conn() as conn: + cur = conn.execute( + f"SELECT {', '.join(cols)} FROM telemetry WHERE ts >= ? ORDER BY ts ASC", + (since,), + ) + rows = cur.fetchall() + if not rows: + return {c: np.empty(0, dtype=float) for c in cols} + arr = np.array([[r[c] if r[c] is not None else np.nan for c in cols] for r in rows], + dtype=float) + return {c: arr[:, i] for i, c in enumerate(cols)} + + # ------------------------------------------------- historical access + + @staticmethod + def auto_bucket(start: float, end: float, target_points: int = 700) -> int: + """Pick a sensible aggregation bucket for a requested span. + + The browser cannot draw more than about a thousand points usefully + and the Pi should not serialise more than it must, so the bucket + grows with the span. Snapped to familiar durations so the x-axis + reads in round numbers rather than 437-second increments. + """ + span = max(float(end) - float(start), 1.0) + raw = span / max(int(target_points), 1) + ladder = [30, 60, 120, 300, 600, 900, 1800, 3600, 7200, + 10800, 21600, 43200, 86400, 604800] + for step in ladder: + if raw <= step: + return step + return ladder[-1] + + def range_series(self, start: float, end: float, + bucket_s: Optional[int] = None) -> Dict[str, Any]: + """Bucket-aggregated telemetry between two epoch timestamps. + + Aggregation happens in SQLite rather than numpy: pulling 90 days of + rows into Python to average them would cost more memory than the + Zero 2 W has to spare. Min and max travel alongside the mean so the + UI can shade a true range band instead of implying the mean was the + whole story. + """ + start, end = float(start), float(end) + if end <= start: + return {"n": 0, "bucket_s": 0, "series": {}} + bucket = int(bucket_s or self.auto_bucket(start, end)) + + # The alias must not be a bare single letter: the telemetry table has + # r, g and b colour columns, and SQLite resolves an unqualified name in + # GROUP BY to a real column before a result alias. `GROUP BY b` silently + # grouped by the blue channel and returned one row per sample while + # cheerfully reporting the requested bucket size. + sql = f""" + SELECT CAST(ts / {bucket} AS INTEGER) * {bucket} AS bucket_ts, + AVG(temp_smooth) AS temp, MIN(temp_smooth) AS temp_lo, + MAX(temp_smooth) AS temp_hi, + AVG(hum_smooth) AS hum, MIN(hum_smooth) AS hum_lo, + MAX(hum_smooth) AS hum_hi, + AVG(press_slp) AS press, MIN(press_slp) AS press_lo, + MAX(press_slp) AS press_hi, + AVG(dew_c) AS dew, AVG(lux) AS lux, + AVG(temp_rate) AS temp_rate, + AVG(press_rate) AS press_rate, + AVG(cpu_temp) AS cpu, COUNT(*) AS n + FROM telemetry + WHERE ts >= ? AND ts <= ? + GROUP BY bucket_ts ORDER BY bucket_ts ASC + """ + with self._conn() as conn: + rows = conn.execute(sql, (start, end)).fetchall() + if not rows: + return {"n": 0, "bucket_s": bucket, "series": {}} + + keys = ["temp", "temp_lo", "temp_hi", "hum", "hum_lo", "hum_hi", + "press", "press_lo", "press_hi", "dew", "lux", + "temp_rate", "press_rate", "cpu", "n"] + out: Dict[str, list] = {"ts": [float(r["bucket_ts"]) for r in rows]} + for k in keys: + out[k] = [r[k] for r in rows] + return {"n": len(rows), "bucket_s": bucket, + "start": start, "end": end, "series": out} + + def daily_summary(self, start: float, end: float) -> List[Dict[str, Any]]: + """Per-calendar-day extremes and means, in the station's local time. + + Local time, not UTC: a `daily minimum` that straddles midnight in + the wrong timezone is the kind of quiet wrongness nobody notices + until they compare against the Met Office and lose an afternoon. + """ + sql = """ + SELECT date(ts, 'unixepoch', 'localtime') AS day, + MIN(ts) AS first_ts, MAX(ts) AS last_ts, COUNT(*) AS n, + MIN(temp_smooth) AS temp_min, MAX(temp_smooth) AS temp_max, + AVG(temp_smooth) AS temp_mean, + MIN(hum_smooth) AS hum_min, MAX(hum_smooth) AS hum_max, + AVG(hum_smooth) AS hum_mean, + MIN(press_slp) AS press_min, MAX(press_slp) AS press_max, + AVG(press_slp) AS press_mean, + AVG(dew_c) AS dew_mean, MAX(lux) AS lux_max + FROM telemetry + WHERE ts >= ? AND ts <= ? + GROUP BY day ORDER BY day DESC + """ + with self._conn() as conn: + return [dict(r) for r in conn.execute(sql, (float(start), float(end))).fetchall()] + + def extremes(self) -> Dict[str, Any]: + """All-time records held by the station, each with when it happened.""" + pairs = [ + ("temp_max", "temp_smooth", "DESC"), ("temp_min", "temp_smooth", "ASC"), + ("hum_max", "hum_smooth", "DESC"), ("hum_min", "hum_smooth", "ASC"), + ("press_max", "press_slp", "DESC"), ("press_min", "press_slp", "ASC"), + ("dew_max", "dew_c", "DESC"), ("dew_min", "dew_c", "ASC"), + ("rate_rise", "press_rate", "DESC"), ("rate_fall", "press_rate", "ASC"), + ] + # Physical sanity bounds. A Kalman filter's rate estimate is garbage + # for the first few samples after it initialises, which happens on + # every restart, and an unfiltered MAX() will faithfully enshrine that + # transient as an all-time record of -37 hPa/h forever. The most + # extreme real sea-level pressure changes on Earth are around + # 10 hPa/h in an explosively deepening cyclone. + bounds = {"press_rate": 10.0, "temp_rate": 25.0} + + out: Dict[str, Any] = {} + with self._conn() as conn: + for name, col, order in pairs: + guard = "" + if col in bounds: + guard = f" AND ABS({col}) <= {bounds[col]}" + row = conn.execute( + f"SELECT ts, {col} AS v FROM telemetry " + f"WHERE {col} IS NOT NULL{guard} ORDER BY {col} {order} LIMIT 1" + ).fetchone() + out[name] = {"ts": row["ts"], "value": row["v"]} if row else None + span = conn.execute("SELECT MIN(ts) AS a, MAX(ts) AS b, COUNT(*) AS n " + "FROM telemetry").fetchone() + out["coverage"] = {"first_ts": span["a"], "last_ts": span["b"], + "rows": span["n"]} + return out + + def iter_csv(self, start: float, end: float): + """Yield CSV lines for export. Generator, so a year of history does + not have to exist in memory at once on a 512 MB board.""" + cols = ["ts", "temp_smooth", "hum_smooth", "press_slp", "dew_c", + "temp_rate", "press_rate", "cpu_temp", "lux", "tier"] + yield "iso_time," + ",".join(cols) + "\n" + with self._conn() as conn: + cur = conn.execute( + f"SELECT {', '.join(cols)} FROM telemetry " + f"WHERE ts >= ? AND ts <= ? ORDER BY ts ASC", + (float(start), float(end)), + ) + while True: + chunk = cur.fetchmany(500) + if not chunk: + break + for r in chunk: + iso = time.strftime("%Y-%m-%dT%H:%M:%S", time.localtime(r["ts"])) + vals = ["" if r[c] is None else + (f"{r[c]:.4f}" if isinstance(r[c], float) else str(r[c])) + for c in cols] + yield iso + "," + ",".join(vals) + "\n" + + def storage_stats(self) -> Dict[str, Any]: + """Rows per resolution tier, so the retention policy is visible.""" + with self._conn() as conn: + rows = conn.execute( + "SELECT tier, COUNT(*) AS n, MIN(ts) AS a, MAX(ts) AS b " + "FROM telemetry GROUP BY tier ORDER BY tier" + ).fetchall() + page = conn.execute("PRAGMA page_count").fetchone()[0] + size = conn.execute("PRAGMA page_size").fetchone()[0] + names = {TIER_RAW: "raw", TIER_5MIN: "5 minute", TIER_HOUR: "hourly"} + return { + "tiers": [{"tier": r["tier"], "label": names.get(r["tier"], "?"), + "rows": r["n"], "first_ts": r["a"], "last_ts": r["b"]} + for r in rows], + "bytes": int(page) * int(size), + } + + def latest(self) -> Optional[Dict[str, Any]]: + with self._conn() as conn: + cur = conn.execute("SELECT * FROM telemetry ORDER BY ts DESC LIMIT 1") + row = cur.fetchone() + return dict(row) if row else None + + def row_count(self) -> int: + with self._conn() as conn: + return int(conn.execute("SELECT COUNT(*) FROM telemetry").fetchone()[0]) + + def span_days(self) -> float: + with self._conn() as conn: + row = conn.execute("SELECT MIN(ts), MAX(ts) FROM telemetry").fetchone() + if not row or row[0] is None: + return 0.0 + return (row[1] - row[0]) / 86400.0 + + def due_forecasts(self, now: Optional[float] = None) -> List[sqlite3.Row]: + """Forecasts whose validity time has passed and can now be scored.""" + now = now or time.time() + with self._conn() as conn: + return conn.execute( + "SELECT * FROM forecasts WHERE valid_ts <= ? AND valid_ts >= ? ORDER BY valid_ts", + (now, now - 7 * 86400), + ).fetchall() + + def scorecard(self) -> List[Dict[str, Any]]: + with self._conn() as conn: + rows = conn.execute( + "SELECT s.* FROM scores s JOIN (" + " SELECT target, horizon_s, MAX(ts) AS mts FROM scores GROUP BY target, horizon_s" + ") m ON s.target = m.target AND s.horizon_s = m.horizon_s AND s.ts = m.mts " + "ORDER BY s.target, s.horizon_s" + ).fetchall() + return [dict(r) for r in rows] + + def recent_events(self, limit: int = 25) -> List[Dict[str, Any]]: + with self._conn() as conn: + rows = conn.execute( + "SELECT * FROM events ORDER BY ts DESC LIMIT ?", (limit,) + ).fetchall() + return [dict(r) for r in rows] + + def labels(self, kind: str, hours: float = 24 * 30) -> Dict[str, np.ndarray]: + since = time.time() - hours * 3600.0 + with self._conn() as conn: + rows = conn.execute( + "SELECT ts, value FROM labels WHERE kind = ? AND ts >= ? ORDER BY ts", + (kind, since), + ).fetchall() + if not rows: + return {"ts": np.empty(0), "value": np.empty(0)} + return { + "ts": np.array([r["ts"] for r in rows], dtype=float), + "value": np.array([r["value"] for r in rows], dtype=float), + } + + # -------------------------------------------------------- compaction + + def compact(self, raw_retention_days: float, five_min_retention_days: float) -> Dict[str, int]: + """Fold old high-resolution rows into means. Returns rows removed per tier.""" + now = time.time() + removed = {"raw": 0, "5min": 0} + removed["raw"] = self._fold(TIER_RAW, TIER_5MIN, 300, + now - raw_retention_days * 86400) + removed["5min"] = self._fold(TIER_5MIN, TIER_HOUR, 3600, + now - five_min_retention_days * 86400) + with self._conn() as conn: + conn.execute("PRAGMA incremental_vacuum") + return removed + + def _fold(self, from_tier: int, to_tier: int, bucket_s: int, older_than: float) -> int: + agg_cols = [c for c in COLUMNS if c != "ts"] + select = ", ".join(f"AVG({c}) AS {c}" for c in agg_cols) + with self._conn() as conn: + rows = conn.execute( + f"SELECT CAST(ts / {bucket_s} AS INTEGER) * {bucket_s} AS bucket, {select} " + f"FROM telemetry WHERE tier = ? AND ts < ? GROUP BY bucket", + (from_tier, older_than), + ).fetchall() + if not rows: + return 0 + cur = conn.execute("SELECT COUNT(*) FROM telemetry WHERE tier = ? AND ts < ?", + (from_tier, older_than)) + n_before = int(cur.fetchone()[0]) + conn.execute("DELETE FROM telemetry WHERE tier = ? AND ts < ?", + (from_tier, older_than)) + payload = [ + tuple([float(r["bucket"])] + [r[c] for c in agg_cols] + [to_tier]) + for r in rows + ] + marks = ", ".join("?" for _ in range(len(agg_cols) + 2)) + conn.executemany( + f"INSERT OR REPLACE INTO telemetry (ts, {', '.join(agg_cols)}, tier) " + f"VALUES ({marks})", + payload, + ) + return n_before - len(rows) + + +def resample(ts: np.ndarray, values: Dict[str, np.ndarray], grid_s: int, + max_gap_grid: int = 3): + """Bin irregular samples onto a regular grid, mean-aggregating each bin. + + Returns (grid_ts, {name: array}) with NaN in bins that had no data and + linear interpolation across gaps no longer than `max_gap_grid` bins. + Anything longer stays NaN so the learner never trains on invention. + """ + if ts.size == 0: + return np.empty(0), {k: np.empty(0) for k in values} + + start = np.floor(ts[0] / grid_s) * grid_s + stop = np.floor(ts[-1] / grid_s) * grid_s + grid = np.arange(start, stop + grid_s, grid_s, dtype=float) + if grid.size == 0: + return np.empty(0), {k: np.empty(0) for k in values} + + idx = np.clip(((ts - start) / grid_s).astype(int), 0, grid.size - 1) + out = {} + counts = np.bincount(idx, minlength=grid.size).astype(float) + for name, arr in values.items(): + clean = np.nan_to_num(arr, nan=0.0) + mask = (~np.isnan(arr)).astype(float) + total = np.bincount(idx, weights=clean, minlength=grid.size) + n = np.bincount(idx, weights=mask, minlength=grid.size) + with np.errstate(invalid="ignore", divide="ignore"): + binned = np.where(n > 0, total / np.maximum(n, 1e-9), np.nan) + out[name] = _interp_short_gaps(binned, max_gap_grid) + out["_count"] = counts + return grid, out + + +def _interp_short_gaps(arr: np.ndarray, max_gap: int) -> np.ndarray: + """Linear fill for runs of NaN up to `max_gap` long; leave longer runs alone.""" + a = arr.copy() + isnan = np.isnan(a) + if not isnan.any() or isnan.all(): + return a + valid = np.flatnonzero(~isnan) + filled = np.interp(np.arange(a.size), valid, a[valid]) + + # find NaN runs and only accept the short ones + edges = np.flatnonzero(np.diff(np.concatenate(([0], isnan.view(np.int8), [0])))) + for start, stop in zip(edges[::2], edges[1::2]): + if (stop - start) <= max_gap and start > 0 and stop < a.size: + a[start:stop] = filled[start:stop] + return a diff --git a/config.yaml b/config.yaml new file mode 100644 index 0000000..2052b5f --- /dev/null +++ b/config.yaml @@ -0,0 +1,35 @@ +# Ashvale Station configuration. Every field is optional: anything omitted +# falls back to the dataclass default in ashvale/config.py. + +site: + name: ashvale-labs-weather-station + latitude: 52.2053 # Cambridge, UK + longitude: 0.1218 + altitude_m: 15.0 # matters more than you would think, see README + timezone: Europe/London + indoors: true # be honest here, it changes how forecasts are worded + +sensor: + sample_period_s: 2.0 + persist_period_s: 30.0 + rotation_deg: 90 + cpu_heat_k: 0.55 # starting point only, calibrate from the dashboard + kalman_q_temp: 2.0e-6 # raise to track faster, lower to smooth harder + kalman_r_temp: 0.02 + +model: + grid_s: 300 + horizons_s: [900, 3600, 10800, 21600, 43200, 86400] + rls_forgetting: 0.9985 # effective memory ~ 11 h on a 5-minute grid + conformal_alpha: 0.10 # 90% prediction intervals + train_period_s: 600 + min_rows_to_train: 120 + +storage: + raw_retention_days: 7.0 + five_min_retention_days: 90.0 + +server: + host: 0.0.0.0 + port: 8000 + led_enabled: true diff --git a/docs/DESIGN.md b/docs/DESIGN.md new file mode 100644 index 0000000..d3ef1f5 --- /dev/null +++ b/docs/DESIGN.md @@ -0,0 +1,393 @@ +# Design notes + +Internals, tuning and failure modes. This document deliberately does **not** +repeat the README: no pitch, no install instructions, no feature list, no API +table. It covers what you need when changing the code or debugging a station +that is behaving oddly, and nothing else. + +Read this before touching anything under `ashvale/models/`. + +--- + +## 1. The data model + +### Schema + +One wide table, `telemetry`, keyed on a float epoch `ts`. Every column is +`REAL` except `tier`. There is no normalisation, because a sensor sample is a +single denormalised event and joins on a Pi are not free. + +| Column | Meaning | +|---|---| +| `temp_raw` | Straight off the HTS221/LPS25HB average, uncompensated | +| `temp_c` | After self-heating compensation, before filtering | +| `temp_smooth`, `temp_rate` | Kalman level and rate (°C, °C/h) | +| `hum`, `hum_smooth` | Raw and filtered relative humidity | +| `press`, `press_slp` | Station pressure and its sea-level reduction | +| `press_smooth`, `press_rate` | Kalman level and rate (hPa, hPa/h) | +| `cpu_temp` | SoC thermal zone. The nuisance variable | +| `dew_c` | Magnus dew point from smoothed inputs | +| `lux`, `r`, `g`, `b` | TCS3400 clear and colour channels | +| `pitch` … `gz` | IMU, nine values | +| `tier` | 0 raw, 1 five-minute mean, 2 hourly mean | + +Supporting tables: `forecasts` (issued, valid, target, mu, lo, hi), `labels` +(human ground truth), `scores` (verification output), `events` (station log). + +**Three columns are named `r`, `g` and `b`.** This has already caused one bug: +`GROUP BY b` in an aggregation query silently grouped by the blue channel, +because SQLite resolves an unqualified identifier to a real column in +preference to a result alias. It returned one row per sample while cheerfully +reporting the requested bucket size. Never use a single-letter alias in this +schema. + +### Tiering + +Nothing is deleted, only downsampled. `Store.compact()` folds raw rows older +than `raw_retention_days` into five-minute means, and five-minute rows older +than `five_min_retention_days` into hourly means. A year of history lands +around 30 MB. + +The `tier` column exists so you can tell a genuine hourly observation from a +mean of twelve. Range queries deliberately mix tiers, which is right for +display but means **a query spanning the raw/five-minute boundary has +non-uniform effective resolution**. If you ever compute a statistic that +assumes equal weight per row, weight by tier or restrict to one. + +### Why aggregation happens in SQL + +`range_series()` buckets with `CAST(ts/bucket AS INTEGER)*bucket` and +`AVG`/`MIN`/`MAX` inside SQLite. Pulling ninety days of rows into numpy to +average them costs more resident memory than the board has. The bucket is +chosen from the span and snapped to a ladder of familiar durations, targeting +about 700 points, because neither a browser nor a human benefits from more. + +Min and max travel alongside the mean so the UI can shade a true range band. A +two-hour bucket that spanned four degrees should not render as a flat line. + +--- + +## 2. The estimation layer + +### Self-heating compensation + +Model: `T = T_raw − k·(T_cpu − T_raw)`, with `k ≥ 0` estimated by recursive +least squares against any trusted reference you supply. + +The regressor is `φ = max(T_cpu − T_raw, 0)` and the target is +`T_raw − T_ref`, so `k·φ` should equal the observed bias. One step, forgetting +factor 0.98: + +``` +gain = P·φ / (λ + φ·P·φ) +k ← clip(k + gain·(target − k·φ), k_min, k_max) +P ← (P − gain·φ·P) / λ +``` + +Measured behaviour: recovers a true `k = 0.62` from a prior of `0.30` in a +**single sample**, and holds post-calibration bias to 0.012 °C over 200 +subsequent readings. + +**The clamp is not decoration.** A single mistyped reference drives `k` to its +bound, and because state persists across restarts it stays there, quietly +biasing every reading until you notice. `POST /api/calibrate {"reset": true}` +exists for exactly that. + +**If you write a simulator or a test fixture, the forward model must be the +exact inverse:** `T_raw = (T + k·T_cpu)/(1 + k)`. Generating the bias as +`T + k·(T_cpu − T)` is a different relation, and the mismatch injects roughly +1.2 °C of phantom noise floor that caps every skill score. This has already +happened once. + +### The Kalman bank + +One constant-velocity filter per signal. State `x = [level, rate]`, standard +continuous white-noise-acceleration process model: + +``` +F = [[1, Δt], [0, 1]] +Q = q · [[Δt³/3, Δt²/2], [Δt²/2, Δt]] +H = [1, 0] +``` + +Two implementation points that matter. + +**Joseph form.** The update is `P ← (I−KH)·P·(I−KH)ᵀ + K·R·Kᵀ`, not the +shorter `P ← (I−KH)·P`. The short form accumulates asymmetry and loses positive +semi-definiteness over months of continuous running, and nobody notices until +the filter quietly stops working. + +**`q` is tuned for the live 2 s cadence.** Because `Q` scales with `Δt³`, +running the same filter at a 300 s step makes the process noise five orders of +magnitude larger, at which point the filter abandons smoothing and tracks +measurement noise. `scripts/simulate.py` therefore scales `q` by +`(sample_period/step)³` when backfilling. Before that fix, rate estimates blew +past anything physical and enshrined a −37 hPa/h all-time record. + +`nis` (normalised innovation squared) is exposed per filter. It should hover +near 1. Persistently high means the filter is too confident and lagging real +change; persistently low means you are over-smoothing. + +--- + +## 3. Features + +33 columns, defined once in `FEATURE_NAMES`. An assertion in +`build_features()` fails loudly if the stacked matrix width drifts from that +list, which is the cheapest guard available against a silently misaligned +design matrix. + +Three rules govern what goes in. + +1. **Anything closed-form is computed, not learned.** Dew point, wet bulb, + VPD, absolute humidity, solar elevation and the clear-sky cloud index come + from `physics.py`. Making a linear learner rediscover the Magnus curve from + data wastes both samples and capacity. +2. **Anything periodic is a sine/cosine pair.** Two diurnal harmonics and one + annual, so phase is representable without a discontinuity at midnight. +3. **Lags are expressed in hours, not samples.** `press_tend_3h` means three + hours whatever `grid_s` is. Changing the grid must not silently change what + the model means by "three hours ago". + +### Standardisation + +`Standardiser` keeps streaming Welford moments and z-scores everything except +the bias column. Not optional: unscaled pressure sits near 1013 while unscaled +temperature rate sits near 0.02, and the resulting condition number will +embarrass you. + +### Solar geometry as a feature + +`solar_position()` is the NOAA low-precision model, accurate to a few tenths of +a degree and costing about twenty floating point operations. It gives the +diurnal cycle real physical structure rather than making the model infer it +from clock time alone. `clear_sky_irradiance()` turns measured lux into a crude +cloudiness index by ratio. Through a south-facing window that is a surprisingly +decent okta estimate; in a north-facing room it is nearly useless. Judge +accordingly. + +--- + +## 4. The forecasting core + +### Why RLS and not gradient descent + +A station produces 288 grid rows a day. Sample efficiency dominates everything +else. RLS is the exact minimiser of the exponentially weighted squared error at +every step, not an approximation, so it converges in far fewer samples than +SGD. The covariance `P` is a genuine parameter-uncertainty estimate, free. One +`(33, 33)` matrix is about 8 kB, so the whole bank of 18 fits in L2 cache on a +Cortex-A53. + +### The covariance trace cap + +``` +P ← (P − outer(gain, P·x)) / λ +P ← (P + Pᵀ)/2 # enforce symmetry +if trace(P) > p_max: P *= p_max/trace(P) +``` + +**This is the single most important guard in the file.** Plain exponential +forgetting inflates `P` without bound along directions the data never excites. +On a quiet night the regressor barely moves, `P` grows exponentially in the +unexcited subspace, and the model detonates on the first sunrise sample. It is +the most common way a field RLS deployment dies. If you refactor `rls.py`, keep +both the cap and the symmetrisation. + +### Direct heads, not recursion + +18 heads: 3 targets × 6 horizons. Each predicts a **delta from now**, and the +absolute forecast is reconstructed as `anchor + delta`. + +Predicting deltas rather than levels matters more than it looks. A model that +must output 14.7 °C spends its capacity representing the mean; one that outputs ++0.4 °C spends it on the weather. + +Iterating a single one-step model 288 times to reach 24 hours would compound +its own bias into a beautifully smooth lie. 18 direct heads cost about 150 kB +total and each is honest about its own horizon. + +### The Hedge ensemble + +Each head blends three members: persistence (delta = 0), climatology (delta +from the harmonic fit) and the learned RLS output. Weights update by +exponentiated gradient on normalised absolute loss, then renormalise. This +guarantees the ensemble is never much worse than its best member, and it +re-weights within about a day when the season turns. + +**The ensemble is allowed to conclude the learned model is useless.** At +15-minute pressure it typically parks most of its weight on persistence. That +is correct behaviour surfaced honestly, not a defect to engineer away. + +### Adaptive conformal intervals + +Split conformal is valid only under exchangeability, and weather is emphatically +not exchangeable: a front arrives and yesterday's residual quantile becomes +fiction. Adaptive conformal inference (Gibbs and Candès) feeds realised coverage +back into the working α: + +``` +α ← clip(α + γ·(α_target − 1[y ∈ C]), 0.005, 0.75) +``` + +The band widens after each miss and narrows after each hit, so long-run +coverage tracks the target whatever the distribution does underneath. + +**Coverage is the acceptance test for any change to this path.** Measured +coverage sits at 89 to 91% against a 90% target across all 18 heads. A change +that improves MAE while coverage drifts to 70% is a regression, not an +improvement, because the intervals have started lying. + +--- + +## 5. Verification + +`Station.verify()` runs every five minutes. It pulls forecasts whose validity +time has passed, looks up the truth and the anchor, and computes per bucket: + +``` +skill = 1 − MAE_model / MAE_persistence +``` + +Zero means no better than assuming nothing changes. **Negative is useful +information, not an embarrassment**: it says ship persistence at that horizon +and stop pretending. + +Scored errors feed two places: the conformal calibrator for each head, and, for +horizons up to 3 h, the Page-Hinkley drift detector. Forecasts more than an +hour past validity are deleted so the table does not grow without bound. + +A forecast that is never scored is an opinion. A forecast scored against +persistence is a measurement. + +### Reproducible evidence + +`scripts/simulate.py` takes `--seed` and `--end`. **Both must be pinned** for +run-to-run comparability. The seed fixes the OU realisation, but the wall-clock +anchor moves solar elevation and the seasonal harmonic, so an unpinned `--end` +changes temperature and humidity while leaving pressure bit-identical. That +asymmetry is a useful diagnostic in itself: if pressure differs between two runs +with the same seed, something other than the anchor has changed. + +Note the interaction with `evaluate.py`: `store.window()` looks back from *now* +with a 60-day default, so a history pinned to a date in the past reports "not +enough history" unless you pass a wide `--hours`. + +--- + +## 6. Monitoring + +Three detectors, because they fail differently. + +**`MahalanobisEWMA`** catches abrupt multivariate novelty: a window opening, a +heater cycling, a squall. The EWMA on the whitened residual gives +persistence-aware detection, so one odd sample is noise and ten in a row is an +event. Shrinkage toward a scaled identity is **not optional**: with six signals +the sample covariance is singular for the first hour, and a singular covariance +turns Mahalanobis distance into a random number generator with an authoritative +name. + +**`PageHinkley`** catches slow change: a sensor drifting, a season turning, a +model going stale. It runs on forecast error and **triggers retraining**, which +is a far better signal than a cron schedule. + +**`SensorHealth`** catches the quietest failure of all. A latched sensor looks +perfectly normal to both detectors above: the dashboard is fine, the model +trains happily, and every forecast is confidently wrong. Bit-identical +consecutive readings are the only tell. + +--- + +## 7. Tuning + +| Symptom | Knob | Direction | +|---|---|---| +| Temperature reads consistently high | Calibrate from the Models tab, or `sensor.cpu_heat_k` | Raise | +| Readings over-smoothed, lag real change | `sensor.kalman_q_temp` | Raise | +| Rates look noisy | `sensor.kalman_q_*` down, or `kalman_r_*` up | | +| NIS persistently much above 1 | Filter too confident, raise `q` | Raise | +| NIS persistently much below 1 | Over-smoothing, lower `q` | Lower | +| Adapts too slowly to a season change | `model.rls_forgetting` toward 0.995 | Lower | +| Jumpy, forgets overnight | `model.rls_forgetting` toward 0.9995 | Raise | +| Coverage well below target | `model.conformal_gamma` | Raise | +| Coverage well above target, lazily wide bands | `model.conformal_gamma` | Lower | +| Drift alarms constantly | `model.drift_lambda` | Raise | +| Retrains eat the CPU | `model.train_period_s` up, `max_pairs` down | | +| Zambretti reads pessimistic everywhere | `site.altitude_m` is wrong | Fix it | + +On forgetting factors: effective memory is `1/(1−λ)` samples. At `λ = 0.9985` +on a five-minute grid that is about 667 samples, roughly 55 hours. `λ = 0.99` +is about two hours and will chase noise. + +--- + +## 8. Performance budget on a Zero 2 W + +Roughly 20× slower than a modern x86 core. Measured: + +| Operation | x86 | Zero 2 W | +|---|---|---| +| Backfill 4032 rows | 0.6 s | 11.8 s | +| Full retrain, 18 heads | ~5 s | 60 to 100 s | +| Resident set, steady state | | ~150 MB | +| Database, one year | | ~30 MB | + +Retraining runs in a worker thread via `asyncio.to_thread`, so the sample loop, +the API and the LED never stall. A command that looks hung on the Pi is usually +just the Pi. + +`max_pairs` in `NowcastEnsemble.fit` caps supervised pairs per head at the most +recent 2500. Not a shortcut: with `λ = 0.9985` the 4000th-most-recent sample +carries a weight of about `e⁻⁶`. It costs real seconds on a Cortex-A53 and buys +nothing measurable. + +**One uvicorn worker, deliberately.** The station owns mutable model state; a +second worker would give you two divergent forecasters sharing a socket. + +--- + +## 9. Extension points + +**Adding a target.** Append to `model.targets` in config and ensure +`station.py` resamples it into the `series` dict; the head bank scales +automatically. Add a climatology column too, or the ensemble's climatology +member returns zero delta forever and quietly wastes a third of its weight. + +**Adding a feature.** Append to `FEATURE_NAMES` **and** to the `column_stack` +in `build_features()`, in the same position. The width assertion catches +mismatches. Persisted RLS state from before the change is now the wrong +dimension: delete `data/state/station_state.json` and retrain. + +**Adding a sensor.** `sensors.py` is the only module that touches hardware. +Extend `SenseBoard.read()`, add the column to `storage.COLUMNS` and the schema, +and mirror it in `SimulatedBoard` so the simulator still exercises every path. +A DS18B20 outside the window is the highest-value hardware change available: it +removes the indoor caveat entirely and improves every model at once. + +**Replacing the learner.** `ForecastHead` needs only `predict(x)` and +`update(x, y)`. Anything with that interface drops in. If you swap in something +without a covariance you lose `predict_std`, and the conformal fallback `sigma` +becomes meaningless until enough residuals accumulate. + +--- + +## 10. Things that will bite you + +- **Station pressure passed where sea-level pressure is expected.** Silent, and + at 100 m elevation it shifts the Zambretti number by about two categories, + permanently. +- **Single-letter SQL aliases.** See section 1. +- **`$` as a JavaScript identifier** in `dashboard.py`. It collides with + bundled libraries and kills the entire script with one opaque syntax error. + The helper is `el()`. +- **Unconstrained Chart.js canvases.** With `maintainAspectRatio: false` a + chart expands to fill its parent. Every canvas needs an explicitly sized + relative wrapper or it swallows its panel. +- **Charts inside `display: none`.** Chart.js cannot measure a hidden canvas, + so tab switches call `resize()` on reveal. +- **Persisted state outliving a schema change.** `station_state.json` holds + trained parameters with fixed dimensions. Change the feature count without + deleting it and you get a shape error at the first update, or worse, silence. +- **Synthetic history left in the database.** After a `simulate.py` backfill, + clear the synthetic rows once real telemetry accumulates, or the model stays + anchored on a stochastic weather model rather than on your room. diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..7f118f6 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,81 @@ +[build-system] +requires = ["setuptools>=68", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "ashvale-station" +version = "1.0.0" +description = "Online machine learning weather forecasting for a Raspberry Pi and Sense HAT, with calibrated uncertainty and public self-verification." +readme = "README.md" +requires-python = ">=3.9" +license = { text = "Apache-2.0" } +authors = [ + { name = "Kemal Yaylali", email = "kemal.yaylali@gmail.com" }, +] +maintainers = [ + { name = "Kemal Yaylali", email = "kemal.yaylali@gmail.com" }, +] +keywords = [ + "weather", "forecasting", "raspberry-pi", "sense-hat", "kalman-filter", + "recursive-least-squares", "conformal-prediction", "online-learning", + "time-series", "state-estimation", "edge-ml", "iot", +] +classifiers = [ + "Development Status :: 4 - Beta", + "Intended Audience :: Science/Research", + "Intended Audience :: Developers", + "License :: OSI Approved :: Apache Software License", + "Operating System :: POSIX :: Linux", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Topic :: Scientific/Engineering :: Atmospheric Science", + "Topic :: Scientific/Engineering :: Artificial Intelligence", + "Topic :: System :: Monitoring", +] + +# Deliberately minimal. The learners are pure numpy so the whole process fits +# in about 150 MB of RSS on a 512 MB board. Anything added here should be +# justified against that budget. +dependencies = [ + "fastapi>=0.110", + "uvicorn>=0.29", + "pydantic>=2.6", + "numpy>=1.24", + "PyYAML>=6.0", +] + +[project.optional-dependencies] +# Hardware access. Prefer the apt packages on a Pi: sense-hat pulls in +# RTIMULib, which is a genuine ordeal to build inside a clean venv. +hardware = [ + "sense-hat>=2.6.0", + "smbus2>=0.4.3", +] +dev = [ + "httpx>=0.27", + "pytest>=8.0", + "ruff>=0.4", +] + +[project.urls] +Homepage = "https://github.com/lynchaos/ashvale-station" +Repository = "https://github.com/lynchaos/ashvale-station" +Issues = "https://github.com/lynchaos/ashvale-station/issues" +Changelog = "https://github.com/lynchaos/ashvale-station/blob/main/CHANGELOG.md" +Documentation = "https://github.com/lynchaos/ashvale-station/blob/main/docs/DESIGN.md" + +[project.scripts] +ashvale = "run:main" + +[tool.setuptools] +packages = ["ashvale", "ashvale.models"] + +[tool.ruff] +line-length = 100 +target-version = "py39" + +[tool.ruff.lint] +select = ["E", "F", "W", "I", "UP", "B"] +ignore = ["E501"] diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..8d35dab --- /dev/null +++ b/requirements.txt @@ -0,0 +1,11 @@ +# Core: everything the ML suite needs. No torch, no sklearn, no pandas. +# The learners are pure numpy by design so a Zero 2 W stays under 120 MB RSS. +fastapi>=0.110 +uvicorn[standard]>=0.29 +pydantic>=2.6 +numpy>=1.24 +PyYAML>=6.0 + +# Hardware only. Skip these on a laptop and the simulator takes over. +# sense-hat>=2.6.0 +# smbus2>=0.4.3 diff --git a/run.py b/run.py new file mode 100644 index 0000000..25c3333 --- /dev/null +++ b/run.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python3 +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Entry point. `python run.py` and open http://:8000""" + +from __future__ import annotations + +import argparse + +import uvicorn + +from ashvale.config import CONFIG + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--host", default=CONFIG.server.host) + ap.add_argument("--port", type=int, default=CONFIG.server.port) + ap.add_argument("--no-led", action="store_true") + ap.add_argument("--reload", action="store_true") + args = ap.parse_args() + + if args.no_led: + CONFIG.server.led_enabled = False + + # One worker, one event loop. The station owns mutable model state, so a + # second worker would give you two divergent forecasters sharing a socket. + uvicorn.run("ashvale.api:app", host=args.host, port=args.port, + reload=args.reload, workers=1, log_level="info", + limit_concurrency=32) + + +if __name__ == "__main__": + main() diff --git a/scripts/evaluate.py b/scripts/evaluate.py new file mode 100644 index 0000000..927fcc9 --- /dev/null +++ b/scripts/evaluate.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python3 +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Rolling-origin backtest. The only number that decides whether to ship. + +Protocol, strictly walk-forward: + + 1. Build the 5-minute feature grid from stored telemetry. + 2. Split at `--train-frac`. Fit the ensemble and the climatology on the + first part only. + 3. Walk the second part one step at a time. At each step, forecast, + record the error, and only then let the model learn from the target + that has just matured. No target is ever visible before its time. + 4. Report MAE against three baselines: + persistence the value now + climatology the harmonic fit + the ensemble + +Skill = 1 - MAE_model / MAE_persistence. A positive number means the +model earns its electricity. A negative number at a given horizon is not +a failure of the exercise, it is the exercise working: ship persistence +at that horizon and stop pretending. + + python scripts/evaluate.py --train-frac 0.6 +""" + +from __future__ import annotations + +import argparse +import sys +import time +from pathlib import Path + +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) + +from ashvale.config import load_config # noqa: E402 +from ashvale.features import build_features # noqa: E402 +from ashvale.models.climatology import HarmonicClimatology # noqa: E402 +from ashvale.models.nowcast import NowcastEnsemble # noqa: E402 +from ashvale.storage import Store, resample # noqa: E402 + + +def horizon_label(seconds: int) -> str: + if seconds < 3600: + return f"{seconds // 60}m" + if seconds < 86400: + return f"{seconds // 3600}h" + return f"{seconds // 86400}d" + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--train-frac", type=float, default=0.6) + ap.add_argument("--hours", type=float, default=24 * 60) + ap.add_argument("--config", default=None) + args = ap.parse_args() + + cfg = load_config(args.config) + store = Store(cfg.storage.db_path) + + raw = store.window(args.hours, ["ts", "temp_smooth", "hum_smooth", "press_slp", "lux"]) + if raw["ts"].size < 200: + print("Not enough history. Run: python scripts/simulate.py --days 14") + return + + grid_ts, cols = resample( + raw["ts"], + {"temperature": raw["temp_smooth"], "humidity": raw["hum_smooth"], + "pressure": raw["press_slp"], "lux": raw["lux"]}, + cfg.model.grid_s, + ) + X, valid = build_features(grid_ts, cols["temperature"], cols["humidity"], + cols["pressure"], cols["lux"], cfg.model.grid_s, + cfg.site.latitude, cfg.site.longitude) + + n = grid_ts.size + split = int(n * args.train_frac) + span_days = (grid_ts[-1] - grid_ts[0]) / 86400.0 + print(f"grid rows : {n} ({span_days:.2f} days at {cfg.model.grid_s}s)") + print(f"train / test : {split} / {n - split}") + + clim = HarmonicClimatology(cfg.model.targets, + min_days_annual=cfg.model.climatology_min_days_annual) + clim.fit(grid_ts[:split], {k: v[:split] for k, v in cols.items() if k in cfg.model.targets}, + valid[:split]) + + ens = NowcastEnsemble(cfg.model.targets, cfg.model.horizons_s, cfg.model) + t0 = time.time() + ens.fit(X[:split], valid[:split], + {k: v[:split] for k, v in cols.items() if k in cfg.model.targets}, + clim, grid_ts[:split]) + print(f"fit : {time.time() - t0:.1f}s\n") + + per_step = cfg.model.grid_s + results = {} + + for target in cfg.model.targets: + y = cols[target] + for h in cfg.model.horizons_s: + steps = max(int(round(h / per_step)), 1) + errs, pers, clims, covered = [], [], [], [] + head = ens.heads[(target, h)] + + for i in range(split, n - steps): + if not valid[i] or not np.isfinite(y[i]) or not np.isfinite(y[i + steps]): + continue + x = ens.scaler.transform(X[i:i + 1])[0] + anchor = float(y[i]) + truth = float(y[i + steps]) + cd = 0.0 + if clim.ready: + cd = float(clim.predict(target, np.array([grid_ts[i] + h]))[0] + - clim.predict(target, np.array([grid_ts[i]]))[0]) + pred = head.predict(x, anchor, cd) + errs.append(truth - pred["mu"]) + pers.append(truth - anchor) + clims.append(truth - (anchor + cd)) + covered.append(1.0 if pred["lo"] <= truth <= pred["hi"] else 0.0) + head.learn(x, anchor, truth, cd) # learn only after scoring + + if len(errs) < 5: + continue + e = np.abs(errs) + p = np.abs(pers) + c = np.abs(clims) + results[(target, h)] = { + "mae": e.mean(), "persistence": p.mean(), "climatology": c.mean(), + "skill": 1.0 - e.mean() / max(p.mean(), 1e-9), + "bias": float(np.mean(errs)), + "coverage": float(np.mean(covered)), + "n": len(errs), + "weights": {k: round(float(v), 2) for k, v in + zip(("pers", "clim", "rls"), head.weights)}, + } + + units = {"temperature": "C", "humidity": "%", "pressure": "hPa"} + header = f"{'target':<12}{'lead':>6}{'MAE':>9}{'persist':>9}{'clim':>9}{'skill':>8}{'cover':>7}{'bias':>8} weights" + print(header) + print("-" * len(header)) + for target in cfg.model.targets: + for h in cfg.model.horizons_s: + r = results.get((target, h)) + if not r: + continue + flag = " <-- persistence wins" if r["skill"] < 0 else "" + print(f"{target:<12}{horizon_label(h):>6}{r['mae']:>9.3f}{r['persistence']:>9.3f}" + f"{r['climatology']:>9.3f}{r['skill'] * 100:>7.1f}%{r['coverage'] * 100:>6.0f}%" + f"{r['bias']:>+8.3f} {r['weights']}{flag}") + print() + + print(f"units: temperature C, humidity %, pressure hPa") + print("coverage should sit near 90% if the conformal calibration is honest.") + + +if __name__ == "__main__": + main() diff --git a/scripts/simulate.py b/scripts/simulate.py new file mode 100644 index 0000000..3a4f33b --- /dev/null +++ b/scripts/simulate.py @@ -0,0 +1,214 @@ +#!/usr/bin/env python3 +# Copyright 2026 Kemal Yaylali +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Seed the database with synthetic history. + +Why this exists: a freshly flashed Pi has no history, and a forecaster +with no history is a random number generator with a nice dashboard. This +script writes physically plausible past telemetry so you can exercise +training, verification and the whole dashboard before the real station +has logged its first night. + +The generator is not a toy. It is a three-component stochastic model: + + pressure Ornstein-Uhlenbeck, tau = 30 h, sigma = 9 hPa + (roughly the observed synoptic variability of NW Europe) + temperature seasonal harmonic + solar-driven diurnal cycle + + OU anomaly (tau = 6 h), with a nocturnal inversion term + humidity driven inversely by temperature about a dew point that + itself performs a slow random walk, which is what makes + RH and T correlate the way they actually do + +Everything is then pushed through the same CPU-heating and noise model +the real sensor suffers from, so a model trained here does not fall over +when it meets real data. + + python scripts/simulate.py --days 21 --wipe +""" + +from __future__ import annotations + +import argparse +import math +import sys +import time +from pathlib import Path + +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) + +from ashvale.config import load_config # noqa: E402 +from ashvale.estimation import SignalTracker # noqa: E402 +from ashvale.physics import (dew_point, sea_level_pressure, # noqa: E402 + solar_position, clear_sky_irradiance) +from ashvale.storage import Store # noqa: E402 + + +def generate(days: float, step_s: int, lat: float, lon: float, + seed: int = 11, end: float | None = None) -> dict: + rng = np.random.default_rng(seed) + n = int(days * 86400 / step_s) + # Anchoring to wall clock makes a fixed seed insufficient for reproducibility: + # the OU realisation repeats, but the timestamps shift, which moves solar + # elevation, day of year and the seasonal harmonic. Those feed the temperature + # model directly, so two same-seed runs produce different data. Pin end as well + # and the backfill becomes bit-reproducible, which is what before/after + # evidence on a model change actually requires. + end = time.time() if end is None else end + ts = end - np.arange(n)[::-1] * step_s + + # --- synoptic pressure: OU process + tau_p, sigma_p = 30 * 3600.0, 9.0 + press = np.zeros(n) + a = math.exp(-step_s / tau_p) + noise_scale = sigma_p * math.sqrt(1 - a * a) + for i in range(1, n): + press[i] = a * press[i - 1] + noise_scale * rng.normal() + press_slp = 1013.0 + press + + # --- solar forcing + elev, _ = solar_position(ts, lat, lon) + elev = np.atleast_1d(elev) + ghi = clear_sky_irradiance(elev) + cloud = np.clip(0.45 + 0.35 * np.sin(2 * np.pi * ts / (4.5 * 86400)) + + 0.25 * rng.normal(size=n).cumsum() / math.sqrt(n), 0.0, 1.0) + lux = np.maximum(ghi * 45.0 * (1.0 - 0.85 * cloud), 0.0) + 6.0 + + # --- temperature: season + diurnal + OU anomaly + inversion at night + doy = np.array([time.gmtime(float(t)).tm_yday for t in ts]) + seasonal = 6.5 * np.sin(2 * np.pi * (doy - 105) / 365.25) + diurnal = 0.011 * ghi * (1.0 - 0.6 * cloud) + inversion = -1.8 * (elev < -3).astype(float) * (1.0 - cloud) + + tau_t, sigma_t = 6 * 3600.0, 1.9 + at = math.exp(-step_s / tau_t) + anom = np.zeros(n) + for i in range(1, n): + anom[i] = at * anom[i - 1] + sigma_t * math.sqrt(1 - at * at) * rng.normal() + # pressure and temperature anomalies are correlated in the real world + anom += 0.12 * press + + temp = 11.5 + seasonal + diurnal + inversion + anom + + # --- humidity via a slowly wandering dew point + dew = temp - 4.5 + 2.5 * np.sin(2 * np.pi * ts / (3.2 * 86400)) + dew -= 0.10 * press + dew = np.minimum(dew, temp - 0.2) + es_t = 6.112 * np.exp(17.625 * temp / (243.04 + temp)) + es_d = 6.112 * np.exp(17.625 * dew / (243.04 + dew)) + rh = np.clip(100.0 * es_d / es_t, 8.0, 100.0) + + # CPU temperature: a slow AR(1) load process, not white noise. A Zero 2 W + # under a steady FastAPI load drifts by a degree or two over minutes, it + # does not jitter by four degrees between samples. + cpu_load = np.zeros(n) + a_cpu = math.exp(-step_s / (900.0)) + for i in range(1, n): + cpu_load[i] = a_cpu * cpu_load[i - 1] + 1.6 * math.sqrt(1 - a_cpu * a_cpu) * rng.normal() + cpu = temp + 21.0 + cpu_load + + # The sensor sits in a thermal gradient between the room and the SoC. + # The compensator inverts T = T_raw - k (T_cpu - T_raw), so the forward + # model must be its exact inverse: T_raw = (T + k T_cpu) / (1 + k). + # Generating it any other way bakes a bias into the synthetic data that + # no amount of calibration can remove, and quietly caps your skill score. + k_true = 0.55 + temp_raw = (temp + k_true * cpu) / (1.0 + k_true) + 0.05 * rng.normal(size=n) + press_station = press_slp / (1.0 + 0.0) - 1.8 # nominal 15 m offset + press_station += 0.05 * rng.normal(size=n) + + return { + "ts": ts, "temp": temp, "temp_raw": temp_raw, "rh": rh + 0.4 * rng.normal(size=n), + "press": press_station, "press_slp": press_slp, "cpu": cpu, + "lux": lux * (0.85 + 0.3 * rng.random(n)), "dew": dew, "cloud": cloud, + } + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--days", type=float, default=14.0) + ap.add_argument("--step", type=int, default=300, help="seconds between rows") + ap.add_argument("--seed", type=int, default=11) + ap.add_argument("--end", type=float, default=None, + help="unix timestamp the history ends at; defaults to now. " + "Pin it with --seed for a bit-reproducible backfill") + ap.add_argument("--wipe", action="store_true", help="clear existing telemetry first") + ap.add_argument("--config", default=None) + args = ap.parse_args() + + cfg = load_config(args.config) + store = Store(cfg.storage.db_path) + + # The Kalman process noise is tuned for the real sampling cadence (2 s). + # Backfilling at 300 s steps with the same q gives Q_level = q*dt^3/3, which + # is five orders of magnitude larger, so the filter abandons smoothing and + # tracks measurement noise. Its rate estimates then blow past anything + # physical and poison the all-time records. Scale q by (real_dt/step)^3 so + # the synthetic history has the same effective smoothing as the live station. + scale = (cfg.sensor.sample_period_s / float(args.step)) ** 3 + cfg.sensor.kalman_q_temp *= scale + cfg.sensor.kalman_q_hum *= scale + cfg.sensor.kalman_q_press *= scale + + if args.wipe: + with store._conn() as conn: + conn.execute("DELETE FROM telemetry") + conn.execute("DELETE FROM forecasts") + conn.execute("DELETE FROM scores") + print("cleared existing telemetry, forecasts and scores") + + data = generate(args.days, args.step, cfg.site.latitude, cfg.site.longitude, + args.seed, args.end) + tracker = SignalTracker(cfg) + + n = data["ts"].size + t0 = time.time() + for i in range(n): + ts = float(data["ts"][i]) + est = tracker.step(ts, float(data["temp_raw"][i]), float(data["rh"][i]), + float(data["press"][i]), float(data["cpu"][i])) + slp = float(sea_level_pressure(est["press_smooth"], est["temp_smooth"], + cfg.site.altitude_m)) + store.insert_telemetry({ + "ts": ts, + "temp_raw": data["temp_raw"][i], + "temp_c": est["temp_c"], + "temp_smooth": est["temp_smooth"], + "temp_rate": est["temp_rate"], + "hum": data["rh"][i], + "hum_smooth": est["hum_smooth"], + "press": data["press"][i], + "press_slp": slp, + "press_smooth": est["press_smooth"], + "press_rate": est["press_rate"], + "cpu_temp": data["cpu"][i], + "dew_c": float(dew_point(est["temp_smooth"], est["hum_smooth"])), + "lux": data["lux"][i], + "r": data["lux"][i] * 0.30, "g": data["lux"][i] * 0.34, "b": data["lux"][i] * 0.28, + "pitch": 0.0, "roll": 0.0, "yaw": 180.0, "compass": 180.0, + "ax": 0.0, "ay": 0.0, "az": 1.0, "gx": 0.0, "gy": 0.0, "gz": 0.0, + }) + if i % 500 == 0: + print(f" {i}/{n} rows", end="\r", flush=True) + + print(f"\nwrote {n} rows spanning {args.days:.1f} days in {time.time() - t0:.1f}s") + print(f"database: {cfg.storage.db_path}") + print("next: python scripts/evaluate.py (or just start the server)") + + +if __name__ == "__main__": + main() diff --git a/systemd/ashvale.service b/systemd/ashvale.service new file mode 100644 index 0000000..31c2f6b --- /dev/null +++ b/systemd/ashvale.service @@ -0,0 +1,26 @@ +[Unit] +Description=Ashvale Station forecast service +After=network-online.target +Wants=network-online.target + +[Service] +Type=simple +User=pi +WorkingDirectory=/home/pi/weather-station/zero2w-weather +ExecStart=/home/pi/weather-station/zero2w-weather/.venv/bin/python run.py +Restart=always +RestartSec=10 + +# A Zero 2 W has 512 MB. Cap the service so a runaway allocation takes the +# service down instead of the whole board. +MemoryMax=280M +CPUWeight=70 +Nice=5 + +# The SD card is a consumable: keep journald from writing every heartbeat. +StandardOutput=journal +StandardError=journal +SyslogIdentifier=ashvale + +[Install] +WantedBy=multi-user.target