kemal/ashvale-station
1
0
Fork 0
mirror of https://github.com/lynchaos/ashvale-station.git synced 2026-09-12 20:52:23 +00:00
Code Issues Packages Projects Releases Wiki Activity

Initial release: Ashvale Station 1.0.0

Browse Source
This commit is contained in:
kemal
2026-08-15 20:43:51 +01:00
commit cbaf3f6c48
36 changed files with 7116 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.github/ISSUE_TEMPLATE/bug_report.md
+24
View File
@@ -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`)
```
```
.github/ISSUE_TEMPLATE/forecast_quality.md
+22
View File
@@ -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.
.github/workflows/ci.yml
+64
View File
@@ -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
.gitignore
+35
View File
@@ -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
CHANGELOG.md
+35
View File
@@ -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.
CONTRIBUTING.md
+88
View File
@@ -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.
LICENSE
+202
View File
@@ -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.
NOTICE
+24
View File
@@ -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.
README.md
+293
View File
@@ -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 <your-repo> ~/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.
SECURITY.md
+32
View File
@@ -0,0 +1,32 @@
# Security Policy
## Supported versions
The latest release on `main` is the supported version.
## Reporting a vulnerability
Please email **[email protected]** 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
```
ashvale/__init__.py
+17
View File
@@ -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"
ashvale/api.py
+423
View File
@@ -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
ashvale/config.py
+150
View File
@@ -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()
ashvale/dashboard.py
+976
View File
@@ -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"""
<!DOCTYPE html>
<html lang="en" class="dark">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Ashvale Station</title>
<script src="https://cdn.tailwindcss.com"></script>
<script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/chart.umd.min.js"></script>
<script src="https://cdn.jsdelivr.net/npm/[email protected]/hammer.min.js"></script>
<script src="https://cdn.jsdelivr.net/npm/[email protected]/dist/chartjs-plugin-zoom.min.js"></script>
<link rel="preconnect" href="https://fonts.googleapis.com">
<link rel="stylesheet" href="https://fonts.googleapis.com/css2?family=Plus+Jakarta+Sans:wght@400;500;600;700;800&family=JetBrains+Mono:wght@400;500;600;700&display=swap">
<style>
body { font-family:'Plus Jakarta Sans',sans-serif; }
.font-mono { font-family:'JetBrains Mono',monospace; }
.glass {
background: radial-gradient(130% 130% at 50% 0%, rgba(30,41,59,.5) 0%, rgba(15,23,42,.75) 100%);
backdrop-filter: blur(20px); -webkit-backdrop-filter: blur(20px);
border: 1px solid rgba(255,255,255,.08);
box-shadow: 0 10px 30px -10px rgba(0,0,0,.5);
}
.tick { transition: width .5s cubic-bezier(.4,0,.2,1); }
.tabbtn { transition: all .18s ease; }
.tabbtn[aria-selected="true"] {
background: rgba(99,102,241,.16); color:#c7d2fe; border-color: rgba(99,102,241,.4);
}
.pane { display:none; }
.pane.active { display:grid; }
.scroller { overflow-y:auto; scrollbar-width:thin; }
.scroller::-webkit-scrollbar { width:7px; }
.scroller::-webkit-scrollbar-thumb { background:rgba(148,163,184,.28); border-radius:8px; }
.flash { animation: flash .5s ease-out; }
@keyframes flash { from { color:#a5b4fc; } to { color:inherit; } }
:focus-visible { outline:2px solid #818cf8; outline-offset:2px; border-radius:6px; }
@media (prefers-reduced-motion: reduce) { *{animation:none!important;transition:none!important} }
@media (min-width:1024px) {
html,body { height:100%; overflow:hidden; }
#shell { height:100dvh; }
}
</style>
</head>
<body class="bg-slate-950 text-slate-100 antialiased selection:bg-indigo-500 selection:text-white">
<div class="fixed inset-0 pointer-events-none overflow-hidden -z-10">
<div class="absolute -top-32 left-1/4 w-[500px] h-[500px] bg-indigo-600/15 rounded-full blur-[120px]"></div>
<div class="absolute top-1/3 -right-32 w-[500px] h-[500px] bg-emerald-600/10 rounded-full blur-[120px]"></div>
<div class="absolute bottom-10 left-10 w-[400px] h-[400px] bg-amber-600/10 rounded-full blur-[100px]"></div>
</div>
<div id="shell" class="max-w-[1600px] mx-auto px-3 sm:px-5 py-3 grid grid-rows-[auto_auto_1fr] gap-3 min-h-0">
<header class="glass rounded-2xl px-4 py-2.5 flex flex-wrap items-center justify-between gap-3">
<div class="flex items-center gap-3">
<div class="p-2 bg-gradient-to-tr from-indigo-500/20 to-emerald-500/20 border border-white/10 rounded-xl">
<svg class="w-5 h-5 text-indigo-400" fill="none" stroke="currentColor" viewBox="0 0 24 24">
<path stroke-linecap="round" stroke-linejoin="round" stroke-width="1.8"
d="M3 15a4 4 0 004 4h9a5 5 0 10-.1-9.999 5.002 5.002 0 00-9.78 2.096A4.001 4.001 0 003 15z"/>
</svg>
</div>
<div>
<h1 class="text-lg font-extrabold tracking-tight leading-none bg-gradient-to-r from-white to-slate-400 bg-clip-text text-transparent">Ashvale Station</h1>
<p class="text-[10px] text-slate-500 font-mono mt-0.5">
<span id="hd-hw">-</span> &middot; <span id="hd-days">0</span> d logged &middot; k=<span id="hd-k">-</span>
</p>
</div>
</div>
<div class="flex items-center gap-2 font-mono text-[11px]">
<span id="hd-health" class="px-2 py-1 rounded-lg border bg-slate-500/15 text-slate-300 border-slate-500/20 uppercase font-semibold">health</span>
<span class="flex items-center gap-1.5 bg-slate-900/80 px-2.5 py-1 rounded-lg border border-slate-800">
<span id="hd-pulse" class="w-1.5 h-1.5 rounded-full bg-emerald-400 animate-pulse"></span>
<span id="hd-time" class="text-white font-semibold">--:--:--</span>
</span>
</div>
</header>
<nav role="tablist" class="glass rounded-2xl p-1.5 flex gap-1.5 overflow-x-auto">
<button role="tab" data-tab="live" aria-selected="true" class="tabbtn shrink-0 px-4 py-2 rounded-xl text-xs font-semibold text-slate-400 border border-transparent hover:text-slate-200">Live</button>
<button role="tab" data-tab="forecast" aria-selected="false" class="tabbtn shrink-0 px-4 py-2 rounded-xl text-xs font-semibold text-slate-400 border border-transparent hover:text-slate-200">Forecast</button>
<button role="tab" data-tab="history" aria-selected="false" class="tabbtn shrink-0 px-4 py-2 rounded-xl text-xs font-semibold text-slate-400 border border-transparent hover:text-slate-200">History</button>
<button role="tab" data-tab="models" aria-selected="false" class="tabbtn shrink-0 px-4 py-2 rounded-xl text-xs font-semibold text-slate-400 border border-transparent hover:text-slate-200">Models</button>
<button role="tab" data-tab="methods" aria-selected="false" class="tabbtn shrink-0 px-4 py-2 rounded-xl text-xs font-semibold text-slate-400 border border-transparent hover:text-slate-200">Methods</button>
</nav>
<main class="min-h-0">
<!-- ---------------- LIVE ---------------- -->
<section id="pane-live" class="pane active h-full min-h-0 gap-3 grid-cols-1 lg:grid-cols-4 lg:grid-rows-[auto_1fr_auto]">
<div class="glass rounded-2xl p-4 flex flex-col justify-between">
<div class="flex items-center justify-between text-[10px] font-semibold uppercase tracking-wider text-amber-400">
<span class="flex items-center gap-1.5"><span class="w-1.5 h-1.5 rounded-full bg-amber-400"></span>Temperature</span>
<span class="text-slate-600 font-mono">KALMAN</span>
</div>
<div class="flex items-baseline gap-1 my-1"><span id="l-temp" class="text-5xl font-extrabold tracking-tight">--</span><span class="text-lg text-amber-400/70 font-semibold">&deg;C</span></div>
<div class="font-mono text-[10px] text-slate-500 space-y-0.5">
<div class="flex justify-between"><span>rate</span><span id="l-temp-rate" class="text-amber-300">--</span></div>
<div class="flex justify-between"><span>raw / cpu</span><span id="l-temp-raw" class="text-slate-400">--</span></div>
</div>
<!-- Chart.js with maintainAspectRatio:false fills its parent, so a
sparkline needs an explicitly sized relative wrapper or it eats the card. -->
<div class="h-9 mt-1.5 shrink-0 relative"><canvas id="spark-temp"></canvas></div>
</div>
<div class="glass rounded-2xl p-4 flex flex-col justify-between">
<div class="flex items-center justify-between text-[10px] font-semibold uppercase tracking-wider text-cyan-400">
<span class="flex items-center gap-1.5"><span class="w-1.5 h-1.5 rounded-full bg-cyan-400"></span>Humidity</span>
<span class="text-slate-600 font-mono">HTS221</span>
</div>
<div class="flex items-baseline gap-1 my-1"><span id="l-hum" class="text-5xl font-extrabold tracking-tight">--</span><span class="text-lg text-cyan-400/70 font-semibold">%</span></div>
<div class="font-mono text-[10px] text-slate-500 space-y-0.5">
<div class="flex justify-between"><span>dew point</span><span id="l-dew" class="text-cyan-300">--</span></div>
<div class="flex justify-between"><span>depression</span><span id="l-dep" class="text-slate-400">--</span></div>
</div>
<!-- Chart.js with maintainAspectRatio:false fills its parent, so a
sparkline needs an explicitly sized relative wrapper or it eats the card. -->
<div class="h-9 mt-1.5 shrink-0 relative"><canvas id="spark-hum"></canvas></div>
</div>
<div class="glass rounded-2xl p-4 flex flex-col justify-between">
<div class="flex items-center justify-between text-[10px] font-semibold uppercase tracking-wider text-violet-400">
<span class="flex items-center gap-1.5"><span class="w-1.5 h-1.5 rounded-full bg-violet-400"></span>Barometer</span>
<span class="text-slate-600 font-mono">MSL</span>
</div>
<div class="flex items-baseline gap-1 my-1"><span id="l-press" class="text-5xl font-extrabold tracking-tight">--</span><span class="text-lg text-violet-400/70 font-semibold">hPa</span></div>
<div class="font-mono text-[10px] text-slate-500 space-y-0.5">
<div class="flex justify-between"><span>tendency</span><span id="l-press-rate" class="text-violet-300">--</span></div>
<div class="flex justify-between"><span>character</span><span id="l-press-char" class="text-slate-400 truncate ml-2">--</span></div>
</div>
<!-- Chart.js with maintainAspectRatio:false fills its parent, so a
sparkline needs an explicitly sized relative wrapper or it eats the card. -->
<div class="h-9 mt-1.5 shrink-0 relative"><canvas id="spark-press"></canvas></div>
</div>
<div class="glass rounded-2xl p-4 flex flex-col justify-between">
<div class="flex items-center justify-between text-[10px] font-semibold uppercase tracking-wider text-emerald-400">
<span class="flex items-center gap-1.5"><span class="w-1.5 h-1.5 rounded-full bg-emerald-400"></span>Sky</span>
<div id="l-swatch" class="w-3 h-3 rounded-full border border-white/30"></div>
</div>
<div class="flex items-baseline gap-1 my-1"><span id="l-lux" class="text-5xl font-extrabold tracking-tight">--</span><span class="text-lg text-emerald-400/70 font-semibold">clr</span></div>
<div class="font-mono text-[10px] text-slate-500 space-y-0.5">
<div class="flex justify-between"><span>cloud index</span><span id="l-cloud" class="text-emerald-300">--</span></div>
<div class="flex justify-between"><span>sun / cct</span><span id="l-sun" class="text-slate-400">--</span></div>
</div>
<!-- Chart.js with maintainAspectRatio:false fills its parent, so a
sparkline needs an explicitly sized relative wrapper or it eats the card. -->
<div class="h-9 mt-1.5 shrink-0 relative"><canvas id="spark-lux"></canvas></div>
</div>
<div class="glass rounded-2xl p-4 lg:col-span-3 flex flex-col min-h-0">
<div class="flex items-center justify-between pb-2 mb-2 border-b border-slate-800 shrink-0">
<div class="flex items-center gap-2">
<h2 class="text-sm font-bold">Rolling window</h2>
<span class="px-1.5 py-0.5 text-[9px] font-mono rounded bg-emerald-500/10 text-emerald-300 border border-emerald-500/20 uppercase">2 s stream</span>
</div>
<div class="flex gap-1" id="live-span">
<button data-min="60" class="px-2 py-1 rounded-lg text-[10px] font-mono border border-slate-800 bg-indigo-600/20 text-indigo-300">1 h</button>
<button data-min="360" class="px-2 py-1 rounded-lg text-[10px] font-mono border border-slate-800 text-slate-400 hover:text-slate-200">6 h</button>
<button data-min="1440" class="px-2 py-1 rounded-lg text-[10px] font-mono border border-slate-800 text-slate-400 hover:text-slate-200">24 h</button>
</div>
</div>
<div class="flex-1 min-h-0 relative"><canvas id="liveChart"></canvas></div>
</div>
<div class="glass rounded-2xl p-4 flex flex-col min-h-0">
<div class="pb-2 mb-2 border-b border-slate-800 shrink-0">
<h2 class="text-sm font-bold">Estimator internals</h2>
<p class="text-[10px] text-indigo-400 font-mono">what the filter is doing right now</p>
</div>
<div class="flex-1 min-h-0 scroller space-y-2 font-mono text-[10px] pr-1">
<div class="bg-slate-900/70 rounded-xl border border-slate-800/80 p-2.5">
<div class="flex justify-between text-slate-400"><span>self-heating k</span><span id="e-k" class="text-emerald-300 font-bold">--</span></div>
<div class="flex justify-between text-slate-500 mt-1"><span>cpu offset</span><span id="e-off">--</span></div>
<div class="text-[9px] text-slate-600 mt-1.5 leading-snug">Removes the SoC bias. Calibrate it on the Models tab.</div>
</div>
<div class="bg-slate-900/70 rounded-xl border border-slate-800/80 p-2.5">
<div class="flex justify-between text-slate-400 mb-1"><span>novelty d&sup2;</span><span id="e-nov" class="text-slate-200 font-bold">--</span></div>
<div class="h-1.5 bg-slate-950 rounded-full overflow-hidden border border-slate-800"><div id="e-nov-bar" class="tick h-full bg-gradient-to-r from-emerald-500 to-amber-500" style="width:0%"></div></div>
<div class="flex justify-between text-slate-400 mt-2 mb-1"><span>drift pressure</span><span id="e-drift" class="text-slate-200 font-bold">--</span></div>
<div class="h-1.5 bg-slate-950 rounded-full overflow-hidden border border-slate-800"><div id="e-drift-bar" class="tick h-full bg-gradient-to-r from-indigo-500 to-rose-500" style="width:0%"></div></div>
<div class="text-[9px] text-slate-600 mt-1.5 leading-snug">Novelty is a multivariate departure from the recent norm. Drift reaching 100% queues a retrain.</div>
</div>
<div id="e-health" class="bg-slate-900/70 rounded-xl border border-slate-800/80 p-2.5 space-y-1"></div>
</div>
</div>
<div class="glass rounded-2xl px-4 py-2.5 lg:col-span-4 grid grid-cols-3 sm:grid-cols-5 lg:grid-cols-10 gap-x-4 gap-y-1.5 font-mono text-[10px]">
<div><div class="text-slate-600 uppercase">wet bulb</div><div id="d-wb" class="text-slate-200 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">vpd</div><div id="d-vpd" class="text-slate-200 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">abs hum</div><div id="d-ah" class="text-slate-200 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">heat idx</div><div id="d-hi" class="text-slate-200 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">solar el</div><div id="d-el" class="text-slate-200 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">pitch</div><div id="d-pitch" class="text-amber-300 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">roll</div><div id="d-roll" class="text-cyan-300 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">yaw</div><div id="d-yaw" class="text-indigo-300 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">compass</div><div id="d-comp" class="text-emerald-300 font-semibold">--</div></div>
<div><div class="text-slate-600 uppercase">accel z</div><div id="d-az" class="text-slate-200 font-semibold">--</div></div>
</div>
</section>
<!-- ---------------- FORECAST ---------------- -->
<section id="pane-forecast" class="pane h-full min-h-0 gap-3 grid-cols-1 lg:grid-cols-3 lg:grid-rows-[1fr_auto]">
<div class="glass rounded-2xl p-4 lg:col-span-2 flex flex-col min-h-0">
<div class="flex flex-wrap items-center justify-between gap-2 pb-2 mb-2 border-b border-slate-800 shrink-0">
<div>
<h2 class="text-sm font-bold">Observed and forecast</h2>
<p class="text-[10px] text-slate-500 font-mono">shaded band is the 90% conformal interval</p>
</div>
<div class="flex items-center gap-1.5">
<select id="fc-target" class="bg-slate-900/90 border border-slate-800 text-slate-300 text-[11px] font-mono rounded-lg px-2 py-1">
<option value="temperature">temperature</option><option value="humidity">humidity</option><option value="pressure">pressure</option>
</select>
<button id="fc-reset" class="px-2.5 py-1 rounded-lg text-[11px] font-mono bg-slate-800/60 text-slate-300 border border-slate-700 hover:bg-slate-800">Reset</button>
</div>
</div>
<div class="flex-1 min-h-0 relative"><canvas id="fanChart"></canvas></div>
<div id="fc-table" class="shrink-0 mt-2 pt-2 border-t border-slate-800 grid grid-cols-3 sm:grid-cols-6 gap-2 font-mono text-[10px]"></div>
</div>
<div class="glass rounded-2xl p-4 flex flex-col min-h-0">
<div class="pb-2 mb-2 border-b border-slate-800 shrink-0">
<h2 class="text-sm font-bold">Conditions ahead</h2>
<p class="text-[10px] text-indigo-400 font-mono">Zambretti prior + online logistic</p>
</div>
<div class="flex-1 min-h-0 scroller space-y-3 pr-1">
<div class="text-center py-1">
<div id="c-label" class="text-lg font-extrabold leading-tight">--</div>
<div class="text-[10px] text-slate-500 font-mono mt-0.5">Z=<span id="c-z">-</span> &middot; <span id="c-trend">-</span></div>
</div>
<div>
<div class="flex justify-between items-baseline mb-1"><span class="text-[10px] text-slate-500 font-mono">rain probability</span><span id="c-pct" class="text-base font-bold text-cyan-300 font-mono">--</span></div>
<div class="h-2 bg-slate-900 rounded-full overflow-hidden border border-slate-800"><div id="c-bar" class="tick h-full bg-gradient-to-r from-cyan-500 to-blue-600" style="width:0%"></div></div>
<div class="flex justify-between text-[9px] font-mono text-slate-600 mt-1"><span>prior <span id="c-prior">-</span></span><span>learner <span id="c-model">-</span></span><span>trust <span id="c-trust">-</span></span></div>
</div>
<div class="bg-slate-900/70 rounded-xl border border-slate-800/80 p-2.5 space-y-2">
<div class="text-[10px] text-slate-400 font-mono">Was it wet in the last hour?</div>
<div class="flex gap-2">
<button data-label="1" class="rain-label flex-1 px-2 py-1.5 rounded-lg bg-blue-600/20 hover:bg-blue-600/30 text-blue-300 border border-blue-500/30 text-[11px] font-medium">Yes, rain</button>
<button data-label="0" class="rain-label flex-1 px-2 py-1.5 rounded-lg bg-slate-800/60 hover:bg-slate-800 text-slate-300 border border-slate-700 text-[11px] font-medium">No, dry</button>
</div>
<div id="c-labstat" class="text-[9px] text-slate-600 font-mono"><span id="c-labn">0</span> confirmed observations</div>
</div>
<div>
<div class="flex justify-between items-baseline mb-1">
<span class="text-[10px] text-slate-500 font-mono">pressure, last 24 h</span>
<span id="c-tend" class="text-[10px] text-violet-300 font-mono">--</span>
</div>
<div class="h-20 relative"><canvas id="tendChart"></canvas></div>
<p class="text-[9px] text-slate-600 font-mono mt-1 leading-snug">The only signal here that sees past your walls. Its slope, not its level, is what drives the forecast above.</p>
</div>
</div>
</div>
<div class="glass rounded-2xl p-4 lg:col-span-3 shrink-0">
<div class="flex items-center justify-between mb-2">
<div>
<h2 class="text-sm font-bold">Seven day outlook</h2>
<p class="text-[10px] text-slate-500 font-mono">climatology plus decaying anomaly, not a synoptic forecast</p>
</div>
<span id="o-badge" class="px-2 py-0.5 rounded-lg bg-amber-500/10 text-amber-300 border border-amber-500/20 text-[9px] font-mono uppercase font-semibold">warming up</span>
</div>
<div id="o-strip" class="grid grid-cols-4 sm:grid-cols-7 gap-2"></div>
</div>
</section>
<!-- ---------------- HISTORY ---------------- -->
<section id="pane-history" class="pane h-full min-h-0 gap-3 grid-cols-1 lg:grid-cols-4 lg:grid-rows-[auto_1fr]">
<div class="glass rounded-2xl px-4 py-3 lg:col-span-4 flex flex-wrap items-end gap-3">
<div class="flex gap-1 flex-wrap" id="h-presets">
<button data-h="6" class="px-2.5 py-1.5 rounded-lg text-[11px] font-mono border border-slate-800 text-slate-400 hover:text-slate-200">6 h</button>
<button data-h="24" class="px-2.5 py-1.5 rounded-lg text-[11px] font-mono border border-slate-800 bg-indigo-600/20 text-indigo-300">24 h</button>
<button data-h="168" class="px-2.5 py-1.5 rounded-lg text-[11px] font-mono border border-slate-800 text-slate-400 hover:text-slate-200">7 d</button>
<button data-h="720" class="px-2.5 py-1.5 rounded-lg text-[11px] font-mono border border-slate-800 text-slate-400 hover:text-slate-200">30 d</button>
<button data-h="2160" class="px-2.5 py-1.5 rounded-lg text-[11px] font-mono border border-slate-800 text-slate-400 hover:text-slate-200">90 d</button>
<button data-h="8760" class="px-2.5 py-1.5 rounded-lg text-[11px] font-mono border border-slate-800 text-slate-400 hover:text-slate-200">1 y</button>
</div>
<div class="flex items-end gap-2">
<label class="block"><span class="block text-[9px] text-slate-600 font-mono uppercase mb-0.5">from</span>
<input id="h-from" type="datetime-local" class="bg-slate-950/70 border border-slate-800 rounded-lg px-2 py-1.5 text-[11px] font-mono text-white"></label>
<label class="block"><span class="block text-[9px] text-slate-600 font-mono uppercase mb-0.5">to</span>
<input id="h-to" type="datetime-local" class="bg-slate-950/70 border border-slate-800 rounded-lg px-2 py-1.5 text-[11px] font-mono text-white"></label>
<button id="h-apply" class="px-3 py-1.5 rounded-lg text-[11px] font-medium bg-emerald-600/20 hover:bg-emerald-600/30 text-emerald-300 border border-emerald-500/30">Apply range</button>
</div>
<div class="flex items-center gap-2 ml-auto">
<div id="h-series" class="flex gap-1"></div>
<a id="h-csv" href="#" class="px-3 py-1.5 rounded-lg text-[11px] font-medium bg-slate-800/60 hover:bg-slate-800 text-slate-300 border border-slate-700">Export CSV</a>
</div>
<p id="h-meta" class="w-full text-[10px] text-slate-600 font-mono">-</p>
</div>
<div class="glass rounded-2xl p-4 lg:col-span-3 flex flex-col min-h-0">
<div class="flex-1 min-h-0 relative"><canvas id="histChart"></canvas></div>
</div>
<div class="glass rounded-2xl p-4 flex flex-col min-h-0">
<div class="pb-2 mb-2 border-b border-slate-800 shrink-0 flex items-center justify-between">
<h2 class="text-sm font-bold">Records</h2>
<div class="flex gap-1">
<button id="rec-tab-daily" class="px-2 py-0.5 rounded text-[10px] font-mono border bg-indigo-600/20 text-indigo-300 border-indigo-500/30">daily</button>
<button id="rec-tab-all" class="px-2 py-0.5 rounded text-[10px] font-mono border border-slate-800 text-slate-500">all time</button>
</div>
</div>
<div id="rec-daily" class="flex-1 min-h-0 scroller pr-1"></div>
<div id="rec-all" class="flex-1 min-h-0 scroller pr-1 hidden space-y-1.5"></div>
</div>
</section>
<!-- ---------------- MODELS ---------------- -->
<section id="pane-models" class="pane h-full min-h-0 gap-3 grid-cols-1 lg:grid-cols-3 lg:grid-rows-[1fr_auto]">
<div class="glass rounded-2xl p-4 lg:col-span-2 flex flex-col min-h-0">
<div class="flex items-center justify-between pb-2 mb-2 border-b border-slate-800 shrink-0">
<div>
<h2 class="text-sm font-bold">Verification scorecard</h2>
<p class="text-[10px] text-slate-500 font-mono">skill above zero means it beats persistence</p>
</div>
<div class="flex gap-1.5">
<button id="m-verify" class="px-2.5 py-1 rounded-lg text-[11px] font-mono bg-slate-800/60 hover:bg-slate-800 text-slate-300 border border-slate-700">Score now</button>
<button id="m-train" class="px-2.5 py-1 rounded-lg text-[11px] font-mono bg-indigo-600/20 hover:bg-indigo-600/30 text-indigo-300 border border-indigo-500/30">Retrain</button>
</div>
</div>
<div class="flex-1 min-h-0 scroller pr-1">
<table class="w-full text-[11px] font-mono">
<thead class="text-slate-600 uppercase text-[9px] sticky top-0 bg-slate-950/90 backdrop-blur">
<tr class="border-b border-slate-800">
<th class="text-left py-1.5">target</th><th class="text-right">lead</th><th class="text-right">MAE</th>
<th class="text-right">persist</th><th class="text-right">skill</th><th class="text-right">cover</th>
<th class="text-right">n</th><th class="text-right pl-3">p/c/l</th>
</tr>
</thead>
<tbody id="m-score" class="text-slate-300"></tbody>
</table>
<p id="m-empty" class="text-[10px] text-slate-600 font-mono mt-3 leading-relaxed">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.</p>
</div>
</div>
<div class="glass rounded-2xl p-4 flex flex-col min-h-0">
<div class="pb-2 mb-2 border-b border-slate-800 shrink-0"><h2 class="text-sm font-bold">Calibration and state</h2></div>
<div class="flex-1 min-h-0 scroller space-y-2.5 pr-1 font-mono text-[10px]">
<div class="bg-slate-900/70 rounded-xl border border-slate-800/80 p-2.5 space-y-2">
<div class="text-slate-400">Trusted thermometer reading</div>
<div class="flex gap-1.5">
<input id="m-calin" type="number" step="0.1" placeholder="20.5" class="flex-1 min-w-0 bg-slate-950/70 border border-slate-800 rounded-lg px-2 py-1.5 text-white">
<button id="m-calgo" class="px-2.5 py-1.5 rounded-lg bg-emerald-600/20 hover:bg-emerald-600/30 text-emerald-300 border border-emerald-500/30">Set</button>
<button id="m-calrst" class="px-2.5 py-1.5 rounded-lg bg-slate-800/60 hover:bg-slate-800 text-slate-400 border border-slate-700">Reset</button>
</div>
<div id="m-calstat" class="text-[9px] text-slate-600 leading-snug">Recursive least squares on the self-heating coefficient. One good reading is enough.</div>
</div>
<div class="bg-slate-900/70 rounded-xl border border-slate-800/80 p-2.5">
<div class="text-slate-400 mb-1.5">Storage tiers</div>
<div id="m-storage" class="space-y-1"></div>
</div>
<div class="bg-slate-900/70 rounded-xl border border-slate-800/80 p-2.5">
<div class="text-slate-400 mb-1.5">Precipitation coefficients</div>
<div id="m-coef" class="space-y-0.5"></div>
</div>
</div>
</div>
<div class="glass rounded-2xl p-4 lg:col-span-3 shrink-0">
<h2 class="text-sm font-bold mb-2">Station log</h2>
<div id="m-log" class="max-h-28 scroller space-y-1 font-mono text-[10px] pr-1"></div>
</div>
</section>
<!-- ---------------- METHODS ---------------- -->
<section id="pane-methods" class="pane h-full min-h-0 gap-3 grid-cols-1 lg:grid-cols-5">
<div class="glass rounded-2xl p-4 lg:col-span-2 flex flex-col min-h-0">
<div class="pb-2 mb-2 border-b border-slate-800 shrink-0">
<h2 class="text-sm font-bold">How it is wired</h2>
<p class="text-[10px] text-slate-500 font-mono">select a stage to read its rationale</p>
</div>
<div class="flex-1 min-h-0 scroller pr-1"><div id="me-diagram"></div></div>
</div>
<div class="glass rounded-2xl p-4 lg:col-span-3 flex flex-col min-h-0">
<div id="me-head" class="pb-2 mb-2 border-b border-slate-800 shrink-0"></div>
<div id="me-body" class="flex-1 min-h-0 scroller pr-1 space-y-3"></div>
</div>
</section>
</main>
</div>
<script>
const el = (id) => document.getElementById(id);
const fmt = (v,d=1) => (v===null||v===undefined||Number.isNaN(v)) ? '--' : Number(v).toFixed(d);
const SEV = { info:'text-slate-400', warn:'text-amber-300', error:'text-rose-300' };
const tsFmt = (t) => new Date(t*1000).toLocaleString([], {month:'short',day:'numeric',hour:'2-digit',minute:'2-digit'});
const GRID = 'rgba(255,255,255,.04)';
const MONO = { family:'JetBrains Mono', size:10 };
const charts = {}, loaders = {};
let activeTab = 'live';
document.querySelectorAll('[role=tab]').forEach(b => b.addEventListener('click', () => selectTab(b.dataset.tab)));
function selectTab(name) {
activeTab = name;
document.querySelectorAll('[role=tab]').forEach(b => b.setAttribute('aria-selected', String(b.dataset.tab===name)));
document.querySelectorAll('.pane').forEach(p => p.classList.toggle('active', p.id==='pane-'+name));
if (loaders[name]) loaders[name]();
// Chart.js cannot measure a canvas inside display:none, so resize on reveal.
setTimeout(() => Object.values(charts).forEach(c => c && c.resize()), 40);
}
/* ---------------- LIVE ---------------- */
const sparks = {};
function makeSpark(id, colour) {
return new Chart(el(id).getContext('2d'), {
type:'line',
data:{ labels:[], datasets:[{ data:[], borderColor:colour, borderWidth:1.5, pointRadius:0, tension:.35, fill:false }] },
options:{ responsive:true, maintainAspectRatio:false, animation:false,
scales:{ x:{display:false}, y:{display:false} },
plugins:{ legend:{display:false}, tooltip:{enabled:false} } }
});
}
['temp','hum','press','lux'].forEach((k,i) =>
sparks[k] = makeSpark('spark-'+k, ['#f59e0b','#06b6d4','#a78bfa','#34d399'][i]));
function pushSpark(k,v) {
if (v===null || v===undefined) return;
const d = sparks[k].data;
d.labels.push(''); d.datasets[0].data.push(v);
if (d.labels.length > 90) { d.labels.shift(); d.datasets[0].data.shift(); }
sparks[k].update('none');
}
charts.live = new Chart(el('liveChart').getContext('2d'), {
type:'line',
data:{ datasets:[
{ label:'temperature', data:[], borderColor:'#f59e0b', backgroundColor:'rgba(245,158,11,.10)',
fill:true, borderWidth:2, pointRadius:0, tension:.3, parsing:false, yAxisID:'y' },
{ label:'humidity', data:[], borderColor:'#06b6d4', borderWidth:2, pointRadius:0,
tension:.3, parsing:false, yAxisID:'y1' } ]},
options:{ responsive:true, maintainAspectRatio:false, animation:false,
interaction:{ mode:'index', intersect:false },
scales:{
x:{ type:'linear', grid:{color:GRID}, ticks:{ color:'#64748b', font:MONO, maxTicksLimit:7,
callback:v=>new Date(v*1000).toLocaleTimeString([], {hour:'2-digit',minute:'2-digit'}) } },
y:{ position:'left', grid:{color:GRID}, ticks:{ color:'#f59e0b', font:MONO, callback:v=>v.toFixed(1)+'\u00b0' } },
y1:{ position:'right', grid:{drawOnChartArea:false}, ticks:{ color:'#06b6d4', font:MONO, callback:v=>v.toFixed(0)+'%' } } },
plugins:{ legend:{display:false},
tooltip:{ backgroundColor:'rgba(15,23,42,.95)', borderColor:'rgba(255,255,255,.1)', borderWidth:1,
titleFont:MONO, bodyFont:{family:'Plus Jakarta Sans',size:11},
callbacks:{ title:i=>new Date(i[0].parsed.x*1000).toLocaleTimeString() } } } }
});
let liveMin = 60;
el('live-span').addEventListener('click', e => {
const b = e.target.closest('button'); if (!b) return;
liveMin = Number(b.dataset.min);
[...el('live-span').children].forEach(x => x.className =
'px-2 py-1 rounded-lg text-[10px] font-mono border border-slate-800 ' +
(x===b ? 'bg-indigo-600/20 text-indigo-300' : 'text-slate-400 hover:text-slate-200'));
loadLiveChart();
});
async function loadLiveChart() {
const r = await fetch('/api/history/range?hours='+(liveMin/60)).then(r=>r.json());
const s = r.series||{};
charts.live.data.datasets[0].data = (s.ts||[]).map((t,i)=>({x:t,y:s.temp[i]})).filter(p=>p.y!=null);
charts.live.data.datasets[1].data = (s.ts||[]).map((t,i)=>({x:t,y:s.hum[i]})).filter(p=>p.y!=null);
charts.live.update('none');
}
loaders.live = loadLiveChart;
function setFlash(id,val) {
const node = el(id); if (!node || node.innerText===val) return;
node.innerText = val; node.classList.remove('flash'); void node.offsetWidth; node.classList.add('flash');
}
function applyTelemetry(d) {
if (!d) return;
el('hd-time').innerText = d.timestamp || '--:--:--';
el('hd-hw').innerText = d.simulated ? 'simulator' : 'sense hat v2';
el('hd-k').innerText = fmt(d.compensator_k,3);
setFlash('l-temp', fmt(d.temperature,2));
setFlash('l-hum', fmt(d.humidity,1));
setFlash('l-press', fmt(d.pressure,1));
setFlash('l-lux', d.color ? String(d.color.clear) : '--');
const rt = d.rates||{}, dv = d.derived||{};
el('l-temp-rate').innerText = (rt.temperature_c_per_h>=0?'+':'')+fmt(rt.temperature_c_per_h,2)+' \u00b0C/h';
el('l-press-rate').innerText = (rt.pressure_hpa_per_h>=0?'+':'')+fmt(rt.pressure_hpa_per_h,2)+' hPa/h';
el('l-temp-raw').innerText = fmt(d.temperature_raw,1)+' / '+fmt(d.cpu_temp,0)+'\u00b0';
el('l-dew').innerText = fmt(dv.dew_point,1)+' \u00b0C';
el('l-dep').innerText = fmt(dv.dew_depression,1)+' K';
el('l-cloud').innerText = fmt(dv.cloud_index,2);
el('l-sun').innerText = fmt(dv.solar_elevation,0)+'\u00b0 / '+(d.color&&d.color.cct?Math.round(d.color.cct)+'K':'n/a');
if (d.color && d.color.hex) el('l-swatch').style.backgroundColor = d.color.hex;
el('d-wb').innerText = fmt(dv.wet_bulb,1)+'\u00b0';
el('d-vpd').innerText = fmt(dv.vpd_hpa,2);
el('d-ah').innerText = fmt(dv.absolute_humidity_g_m3,1);
el('d-hi').innerText = fmt(dv.heat_index,1)+'\u00b0';
el('d-el').innerText = fmt(dv.solar_elevation,0)+'\u00b0';
el('d-pitch').innerText = fmt(d.pitch,1)+'\u00b0';
el('d-roll').innerText = fmt(d.roll,1)+'\u00b0';
el('d-yaw').innerText = fmt(d.yaw,1)+'\u00b0';
el('d-comp').innerText = fmt(d.compass,1)+'\u00b0';
el('d-az').innerText = fmt(d.accel && d.accel.z,2);
el('e-k').innerText = fmt(d.compensator_k,4);
el('e-off').innerText = fmt(d.cpu_offset,1)+' K';
el('e-nov').innerText = fmt(d.novelty_d2,1);
el('e-nov-bar').style.width = Math.min((d.novelty_d2||0)/24,1)*100+'%';
pushSpark('temp', d.temperature); pushSpark('hum', d.humidity);
pushSpark('press', d.pressure); pushSpark('lux', d.color && d.color.clear);
}
function applyPrecip(p) {
if (!p || p.rain_probability===undefined) return;
el('c-label').innerText = p.label||'--';
el('c-z').innerText = p.zambretti_z!==undefined ? p.zambretti_z : '-';
el('c-trend').innerText = p.pressure_characteristic||'-';
if (p.tendency!==undefined) el('c-tend').innerText = (p.tendency>=0?'+':'')+fmt(p.tendency,2)+' hPa/h';
el('l-press-char').innerText = p.pressure_characteristic||'--';
const pct = Math.round(p.rain_probability*100);
el('c-pct').innerText = pct+'%'; el('c-bar').style.width = pct+'%';
el('c-prior').innerText = Math.round((p.prior_probability||0)*100)+'%';
el('c-model').innerText = Math.round((p.model_probability||0)*100)+'%';
el('c-trust').innerText = Math.round((p.learner_trust||0)*100)+'%';
el('c-labn').innerText = p.strong_labels||0;
}
function connectStream() {
const es = new EventSource('/api/stream');
es.onmessage = ev => {
const d = JSON.parse(ev.data);
applyTelemetry(d.telemetry); applyPrecip(d.precipitation);
const dr = Math.round((d.drift_stress||0)*100);
el('e-drift').innerText = dr+'%'; el('e-drift-bar').style.width = dr+'%';
const b = el('hd-health'); b.innerText = d.health||'unknown';
b.className = 'px-2 py-1 rounded-lg border uppercase font-semibold ' + (
d.health==='ok' ? 'bg-emerald-500/15 text-emerald-300 border-emerald-500/20'
: d.health==='warn' ? 'bg-amber-500/15 text-amber-300 border-amber-500/20'
: 'bg-rose-500/15 text-rose-300 border-rose-500/20');
el('hd-pulse').className = 'w-1.5 h-1.5 rounded-full bg-emerald-400 animate-pulse';
};
es.onerror = () => { el('hd-pulse').className='w-1.5 h-1.5 rounded-full bg-rose-500'; es.close(); setTimeout(connectStream,5000); };
}
/* ---------------- FORECAST ---------------- */
charts.fan = new Chart(el('fanChart').getContext('2d'), {
type:'line',
data:{ datasets:[
{ label:'observed', data:[], borderColor:'#f59e0b', backgroundColor:'rgba(245,158,11,.12)',
fill:true, borderWidth:2, pointRadius:0, tension:.3, parsing:false },
{ label:'upper', data:[], borderColor:'rgba(99,102,241,.25)', backgroundColor:'rgba(99,102,241,.14)',
borderWidth:1, pointRadius:0, tension:.3, fill:'+1', parsing:false },
{ label:'lower', data:[], borderColor:'rgba(99,102,241,.25)', borderWidth:1,
pointRadius:0, tension:.3, fill:false, parsing:false },
{ label:'forecast', data:[], borderColor:'#818cf8', borderDash:[6,4], borderWidth:2,
pointRadius:3, pointBackgroundColor:'#818cf8', tension:.3, parsing:false } ]},
options:{ responsive:true, maintainAspectRatio:false, animation:false,
interaction:{ mode:'nearest', axis:'x', intersect:false },
scales:{ x:{ type:'linear', grid:{color:GRID}, ticks:{ color:'#64748b', font:MONO, maxTicksLimit:7,
callback:v=>new Date(v*1000).toLocaleTimeString([], {hour:'2-digit',minute:'2-digit'}) } },
y:{ grid:{color:GRID}, ticks:{ color:'#cbd5e1', font:MONO } } },
plugins:{ legend:{display:false},
tooltip:{ backgroundColor:'rgba(15,23,42,.95)', borderColor:'rgba(255,255,255,.1)', borderWidth:1,
filter:i=>i.dataset.label!=='lower', titleFont:MONO, callbacks:{ title:i=>tsFmt(i[0].parsed.x) } },
zoom:{ pan:{enabled:true,mode:'xy'}, zoom:{ wheel:{enabled:true,speed:.08}, pinch:{enabled:true}, mode:'xy' } } } }
});
charts.tend = new Chart(el('tendChart').getContext('2d'), {
type:'line',
data:{ datasets:[{ data:[], borderColor:'#a78bfa', backgroundColor:'rgba(167,139,250,.12)',
fill:true, borderWidth:1.6, pointRadius:0, tension:.3, parsing:false }] },
options:{ responsive:true, maintainAspectRatio:false, animation:false,
scales:{ x:{ type:'linear', grid:{display:false}, ticks:{ color:'#475569', font:{family:'JetBrains Mono',size:8}, maxTicksLimit:4,
callback:v=>new Date(v*1000).toLocaleTimeString([], {hour:'2-digit'}) } },
y:{ grid:{color:GRID}, ticks:{ color:'#a78bfa', font:{family:'JetBrains Mono',size:8}, maxTicksLimit:4, callback:v=>v.toFixed(0) } } },
plugins:{ legend:{display:false}, tooltip:{ backgroundColor:'rgba(15,23,42,.95)', titleFont:MONO,
callbacks:{ title:i=>tsFmt(i[0].parsed.x) } } } }
});
el('fc-reset').addEventListener('click', ()=>charts.fan.resetZoom());
el('fc-target').addEventListener('change', loadForecast);
async function loadForecast() {
const target = el('fc-target').value;
const key = {temperature:'temp', humidity:'hum', pressure:'press'}[target];
const res = await Promise.all([
fetch('/api/history/range?hours=24').then(r=>r.json()),
fetch('/api/forecast').then(r=>r.json()) ]);
const hist = res[0], fc = res[1], s = hist.series||{};
charts.fan.data.datasets[0].data = (s.ts||[]).map((t,i)=>({x:t,y:s[key][i]})).filter(p=>p.y!=null);
// the pressure trace is the same fetch, reused: one request, two panels
charts.tend.data.datasets[0].data = (s.ts||[]).map((t,i)=>({x:t,y:s.press[i]})).filter(p=>p.y!=null);
charts.tend.update('none');
const series = (fc.targets&&fc.targets[target])||[];
const anchor = fc.anchors ? fc.anchors[target] : undefined;
const head = (anchor!==undefined && fc.issued_ts) ? [{x:fc.issued_ts,y:anchor}] : [];
charts.fan.data.datasets[3].data = head.concat(series.map(p=>({x:p.valid_ts,y:p.mu})));
charts.fan.data.datasets[1].data = head.concat(series.map(p=>({x:p.valid_ts,y:p.hi})));
charts.fan.data.datasets[2].data = head.concat(series.map(p=>({x:p.valid_ts,y:p.lo})));
charts.fan.update('none');
el('fc-table').innerHTML = series.length ? series.map(p=>
'<div class="bg-slate-900/60 rounded-lg border border-slate-800/70 px-2 py-1.5">'+
'<div class="text-slate-600 uppercase text-[9px]">'+p.horizon_label+'</div>'+
'<div class="text-slate-100 font-bold text-xs">'+p.mu.toFixed(1)+'</div>'+
'<div class="text-indigo-400/80 text-[9px]">&plusmn;'+((p.hi-p.lo)/2).toFixed(2)+'</div></div>').join('')
: '<p class="col-span-full text-[10px] text-slate-600">Awaiting the first training pass.</p>';
}
async function loadOutlook() {
const d = await fetch('/api/outlook').then(r=>r.json());
el('o-badge').innerText = d.ready ? (d.annual_terms?'seasonal terms on':'diurnal only') : 'warming up';
const rows = (d.targets&&d.targets.temperature)||[];
if (!rows.length) { el('o-strip').innerHTML = '<p class="col-span-full text-[10px] text-slate-600 font-mono">Needs about two days of history before the harmonic fit means anything.</p>'; return; }
const byDay = {};
rows.forEach(r => { const k = new Date(r.ts*1000).toLocaleDateString([], {weekday:'short'}); (byDay[k]=byDay[k]||[]).push(r); });
el('o-strip').innerHTML = Object.keys(byDay).slice(0,7).map(day=>{
const v = byDay[day];
const hi = Math.max.apply(null, v.map(x=>x.mu)), lo = Math.min.apply(null, v.map(x=>x.mu));
const sp = Math.max.apply(null, v.map(x=>x.hi-x.lo))/2;
return '<div class="bg-slate-900/70 rounded-xl border border-slate-800/80 p-2 text-center">'+
'<div class="text-[9px] uppercase text-slate-600 font-mono">'+day+'</div>'+
'<div class="text-base font-bold font-mono mt-0.5">'+hi.toFixed(1)+'\u00b0</div>'+
'<div class="text-[10px] text-slate-500 font-mono">'+lo.toFixed(1)+'\u00b0</div>'+
'<div class="text-[9px] text-indigo-400/80 font-mono">&plusmn;'+sp.toFixed(1)+'</div></div>';
}).join('');
}
loaders.forecast = () => { loadForecast(); loadOutlook(); };
document.querySelectorAll('.rain-label').forEach(b => b.addEventListener('click', async () => {
const res = await fetch('/api/label', { method:'POST', headers:{'Content-Type':'application/json'},
body: JSON.stringify({kind:'rain', value:Number(b.dataset.label)}) }).then(r=>r.json());
el('c-labstat').innerHTML = '<span class="text-emerald-300">Recorded.</span> '+(res.strong_labels||0)+' confirmed, loss '+fmt(res.loss,3);
}));
/* ---------------- HISTORY ---------------- */
let hRange = { hours:24, start:null, end:null }, histData = null;
const SERIES_META = {
temp:{ label:'temperature', colour:'#f59e0b', axis:'y', on:true },
hum:{ label:'humidity', colour:'#06b6d4', axis:'y1', on:false },
press:{ label:'pressure', colour:'#a78bfa', axis:'y2', on:true },
dew:{ label:'dew point', colour:'#34d399', axis:'y', on:false }
};
el('h-series').innerHTML = Object.keys(SERIES_META).map(k=>
'<button data-s="'+k+'" class="px-2 py-1 rounded-lg text-[10px] font-mono border"></button>').join('');
function paintSeriesButtons() {
document.querySelectorAll('#h-series button').forEach(b => {
const m = SERIES_META[b.dataset.s];
b.innerText = b.dataset.s;
b.style.borderColor = m.on ? m.colour+'66' : 'rgb(30,41,59)';
b.style.backgroundColor = m.on ? m.colour+'22' : 'transparent';
b.style.color = m.on ? m.colour : 'rgb(100,116,139)';
});
}
paintSeriesButtons();
charts.hist = new Chart(el('histChart').getContext('2d'), {
type:'line', data:{ datasets:[] },
options:{ responsive:true, maintainAspectRatio:false, animation:false,
interaction:{ mode:'index', intersect:false },
scales:{
x:{ type:'linear', grid:{color:GRID}, ticks:{ color:'#64748b', font:MONO, maxTicksLimit:9,
callback:v=>{ const d=new Date(v*1000); const span=(hRange.end||0)-(hRange.start||0);
return span > 3*86400 ? d.toLocaleDateString([], {month:'short',day:'numeric'})
: d.toLocaleTimeString([], {hour:'2-digit',minute:'2-digit'}); } } },
y:{ position:'left', grid:{color:GRID}, ticks:{ color:'#f59e0b', font:MONO, callback:v=>v.toFixed(1)+'\u00b0' } },
y1:{ position:'right', display:false, grid:{drawOnChartArea:false}, ticks:{ color:'#06b6d4', font:MONO } },
y2:{ position:'right', grid:{drawOnChartArea:false}, ticks:{ color:'#a78bfa', font:MONO, callback:v=>v.toFixed(0) } } },
plugins:{ legend:{display:false},
tooltip:{ backgroundColor:'rgba(15,23,42,.95)', borderColor:'rgba(255,255,255,.1)', borderWidth:1,
titleFont:MONO, bodyFont:{family:'Plus Jakarta Sans',size:11},
filter:i=>i.dataset.label.charAt(0)!=='_',
callbacks:{ title:i=>tsFmt(i[0].parsed.x) } },
zoom:{ pan:{enabled:true,mode:'x'}, zoom:{ wheel:{enabled:true,speed:.08}, pinch:{enabled:true}, mode:'x' } } } }
});
el('h-presets').addEventListener('click', e => {
const b = e.target.closest('button'); if (!b) return;
[...el('h-presets').children].forEach(x => x.className =
'px-2.5 py-1.5 rounded-lg text-[11px] font-mono border border-slate-800 ' +
(x===b ? 'bg-indigo-600/20 text-indigo-300' : 'text-slate-400 hover:text-slate-200'));
hRange = { hours:Number(b.dataset.h), start:null, end:null };
loadHistory();
});
el('h-series').addEventListener('click', e => {
const b = e.target.closest('button'); if (!b) return;
SERIES_META[b.dataset.s].on = !SERIES_META[b.dataset.s].on;
paintSeriesButtons(); drawHistory();
});
el('h-apply').addEventListener('click', () => {
const a = el('h-from').value, z = el('h-to').value;
if (!a || !z) { el('h-meta').innerText = 'Pick both a from and a to date.'; return; }
const s = new Date(a).getTime()/1000, e = new Date(z).getTime()/1000;
if (e <= s) { el('h-meta').innerText = 'The to date must be after the from date.'; return; }
hRange = { hours:null, start:s, end:e };
[...el('h-presets').children].forEach(x => x.className =
'px-2.5 py-1.5 rounded-lg text-[11px] font-mono border border-slate-800 text-slate-400 hover:text-slate-200');
loadHistory();
});
function rangeQuery() {
return hRange.hours!=null ? 'hours='+hRange.hours : 'start='+hRange.start+'&end='+hRange.end;
}
async function loadHistory() {
el('h-meta').innerText = 'Loading...';
const r = await fetch('/api/history/range?'+rangeQuery()).then(r=>r.json());
histData = r;
if (r.start) hRange.start = r.start;
if (r.end) hRange.end = r.end;
el('h-csv').href = '/api/export.csv?'+rangeQuery();
const bs = r.bucket_s||0;
const bl = bs>=86400 ? (bs/86400)+' d' : bs>=3600 ? (bs/3600)+' h' : (bs/60)+' min';
el('h-meta').innerText = r.n
? r.n+' points at '+bl+' resolution, '+tsFmt(r.start)+' to '+tsFmt(r.end)
: 'No data in that range. The station may not have been running then.';
drawHistory(); loadRecords();
}
function drawHistory() {
if (!histData || !histData.n) { charts.hist.data.datasets = []; charts.hist.update('none'); return; }
const s = histData.series, ds = [];
// Range band: min and max within each bucket, so an hourly view still shows
// that the hour spanned four degrees rather than implying a flat mean.
if (SERIES_META.temp.on && s.temp_hi) {
ds.push({ label:'_hi', data:s.ts.map((t,i)=>({x:t,y:s.temp_hi[i]})), borderColor:'transparent',
backgroundColor:'rgba(245,158,11,.10)', fill:'+1', pointRadius:0, parsing:false, yAxisID:'y', order:9 });
ds.push({ label:'_lo', data:s.ts.map((t,i)=>({x:t,y:s.temp_lo[i]})), borderColor:'transparent',
fill:false, pointRadius:0, parsing:false, yAxisID:'y', order:9 });
}
Object.keys(SERIES_META).forEach(k => {
const m = SERIES_META[k];
if (!m.on || !s[k]) return;
ds.push({ label:m.label, data:s.ts.map((t,i)=>({x:t,y:s[k][i]})).filter(p=>p.y!=null),
borderColor:m.colour, borderWidth:1.8, pointRadius:0, tension:.25, parsing:false, yAxisID:m.axis, order:1 });
});
charts.hist.options.scales.y1.display = SERIES_META.hum.on;
charts.hist.options.scales.y2.display = SERIES_META.press.on;
charts.hist.data.datasets = ds;
charts.hist.update('none');
}
async function loadRecords() {
const span = (hRange.end-hRange.start)||86400;
const days = Math.max(2, Math.min(400, Math.ceil(span/86400)));
const res = await Promise.all([
fetch('/api/history/daily?days='+days).then(r=>r.json()),
fetch('/api/records').then(r=>r.json()) ]);
const daily = res[0].days||[], all = res[1];
el('rec-daily').innerHTML = daily.length ?
'<table class="w-full text-[10px] font-mono"><thead class="text-slate-600 uppercase text-[9px] sticky top-0 bg-slate-950/90 backdrop-blur">'+
'<tr class="border-b border-slate-800"><th class="text-left py-1">day</th><th class="text-right">min</th>'+
'<th class="text-right">max</th><th class="text-right">hPa</th></tr></thead><tbody>'+
daily.map(d=>'<tr class="border-b border-slate-800/40"><td class="py-1 text-slate-400">'+d.day.slice(5)+'</td>'+
'<td class="text-right text-cyan-300">'+fmt(d.temp_min,1)+'</td>'+
'<td class="text-right text-amber-300">'+fmt(d.temp_max,1)+'</td>'+
'<td class="text-right text-slate-500">'+fmt(d.press_mean,0)+'</td></tr>').join('')+
'</tbody></table>'
: '<p class="text-[10px] text-slate-600 font-mono">No completed days yet.</p>';
const R = (k,label,unit,dec) => all[k] ?
'<div class="flex justify-between items-baseline bg-slate-900/60 rounded-lg border border-slate-800/70 px-2 py-1.5">'+
'<span class="text-slate-500 text-[10px]">'+label+'</span><span class="text-right">'+
'<span class="text-slate-100 font-bold text-[11px]">'+fmt(all[k].value,dec===undefined?1:dec)+unit+'</span>'+
'<span class="block text-slate-600 text-[9px]">'+tsFmt(all[k].ts)+'</span></span></div>' : '';
el('rec-all').innerHTML = [
R('temp_max','warmest','\u00b0'), R('temp_min','coldest','\u00b0'),
R('press_max','highest pressure',''), R('press_min','lowest pressure',''),
R('hum_max','most humid','%'), R('hum_min','driest','%'),
R('rate_fall','fastest fall','/h',2), R('rate_rise','fastest rise','/h',2)
].join('') || '<p class="text-[10px] text-slate-600 font-mono">No records yet.</p>';
}
function toggleRec(daily) {
el('rec-daily').classList.toggle('hidden', !daily);
el('rec-all').classList.toggle('hidden', daily);
el('rec-tab-daily').className = 'px-2 py-0.5 rounded text-[10px] font-mono border '+(daily?'bg-indigo-600/20 text-indigo-300 border-indigo-500/30':'border-slate-800 text-slate-500');
el('rec-tab-all').className = 'px-2 py-0.5 rounded text-[10px] font-mono border '+(!daily?'bg-indigo-600/20 text-indigo-300 border-indigo-500/30':'border-slate-800 text-slate-500');
}
el('rec-tab-daily').addEventListener('click', ()=>toggleRec(true));
el('rec-tab-all').addEventListener('click', ()=>toggleRec(false));
loaders.history = () => { if (!histData) loadHistory(); };
/* ---------------- MODELS ---------------- */
async function loadModels() {
const res = await Promise.all([
fetch('/api/scorecard').then(r=>r.json()), fetch('/api/models').then(r=>r.json()),
fetch('/api/status').then(r=>r.json()), fetch('/api/storage').then(r=>r.json()) ]);
const sc = res[0], md = res[1], st = res[2], sg = res[3];
el('hd-days').innerText = fmt(st.history_days,2);
const wmap = {};
(md.nowcast||[]).forEach(h => { wmap[h.target+'@'+h.horizon_s] = h.weights; });
const rows = sc.rows||[];
el('m-empty').style.display = rows.length ? 'none' : 'block';
el('m-score').innerHTML = rows.map(r=>{
const sk = r.skill||0;
const cls = sk>0.05?'text-emerald-300':sk<-0.05?'text-rose-300':'text-slate-400';
const lead = r.horizon_s<3600 ? (r.horizon_s/60)+'m' : r.horizon_s<86400 ? (r.horizon_s/3600)+'h' : (r.horizon_s/86400)+'d';
const w = wmap[r.target+'@'+r.horizon_s];
const ws = w ? Math.round(w.persistence*100)+'/'+Math.round(w.climatology*100)+'/'+Math.round(w.learned*100) : '-';
return '<tr class="border-b border-slate-800/40"><td class="py-1.5 text-slate-300">'+r.target+'</td>'+
'<td class="text-right text-slate-500">'+lead+'</td><td class="text-right">'+fmt(r.mae,3)+'</td>'+
'<td class="text-right text-slate-600">'+fmt(r.mae_persistence,3)+'</td>'+
'<td class="text-right font-bold '+cls+'">'+Math.round(sk*100)+'%</td>'+
'<td class="text-right text-slate-400">'+Math.round((r.coverage||0)*100)+'%</td>'+
'<td class="text-right text-slate-700">'+r.n+'</td>'+
'<td class="text-right text-slate-500 pl-3">'+ws+'</td></tr>';
}).join('');
el('m-storage').innerHTML = (sg.tiers||[]).map(t=>
'<div class="flex justify-between text-slate-500"><span>'+t.label+'</span>'+
'<span class="text-slate-300">'+t.rows.toLocaleString()+' rows</span></div>').join('')+
'<div class="flex justify-between text-slate-600 pt-1 mt-1 border-t border-slate-800"><span>database</span>'+
'<span>'+(sg.bytes/1e6).toFixed(2)+' MB</span></div>';
const coef = ((md.precipitation||{}).coefficients||[]).slice()
.sort((a,b)=>Math.abs(b.weight)-Math.abs(a.weight)).slice(0,6);
const mx = Math.max.apply(null, coef.map(c=>Math.abs(c.weight)).concat([1e-6]));
el('m-coef').innerHTML = coef.map(c=>
'<div class="flex items-center gap-1.5"><span class="w-24 truncate text-slate-600 text-[9px]">'+c.feature+'</span>'+
'<div class="flex-1 h-1 bg-slate-950 rounded-full overflow-hidden"><div class="h-full '+
(c.weight>=0?'bg-emerald-500':'bg-rose-500')+'" style="width:'+(Math.abs(c.weight)/mx*100)+'%"></div></div>'+
'<span class="w-10 text-right text-slate-500 text-[9px]">'+c.weight.toFixed(2)+'</span></div>').join('');
el('m-log').innerHTML = (st.events||[]).map(e=>
'<div class="flex gap-2 border-b border-slate-800/40 pb-1">'+
'<span class="text-slate-700 shrink-0">'+new Date(e.ts*1000).toLocaleTimeString()+'</span>'+
'<span class="text-slate-600 shrink-0 w-16">'+e.kind+'</span>'+
'<span class="'+(SEV[e.severity]||'text-slate-400')+'">'+e.detail+'</span></div>').join('')
|| '<span class="text-slate-700">Nothing logged yet.</span>';
const an = await fetch('/api/anomaly').then(r=>r.json());
el('e-health').innerHTML = Object.keys(an.health||{}).map(k=>{
const v = an.health[k];
const dot = v.status==='ok'?'bg-emerald-400':v.status==='warn'?'bg-amber-400':'bg-rose-400';
return '<div class="flex justify-between items-center text-slate-500">'+
'<span class="flex items-center gap-1.5"><span class="w-1.5 h-1.5 rounded-full '+dot+'"></span>'+k+'</span>'+
'<span class="text-slate-600 text-[9px] truncate ml-2">'+v.detail+'</span></div>';
}).join('');
}
loaders.models = loadModels;
el('m-train').addEventListener('click', async () => {
el('m-train').innerText = 'Training...';
const r = await fetch('/api/train', {method:'POST'}).then(r=>r.json());
el('m-train').innerText = r.trained ? 'Trained '+r.grid_rows : 'Not enough data';
setTimeout(()=>{ el('m-train').innerText='Retrain'; }, 4000);
loadModels();
});
el('m-verify').addEventListener('click', async () => { await fetch('/api/verify',{method:'POST'}); loadModels(); });
el('m-calgo').addEventListener('click', async () => {
const v = parseFloat(el('m-calin').value);
if (Number.isNaN(v)) { el('m-calstat').innerText = 'Enter a temperature in degrees Celsius.'; return; }
const r = await fetch('/api/calibrate', {method:'POST', headers:{'Content-Type':'application/json'},
body: JSON.stringify({reference_c:v})}).then(r=>r.json());
el('m-calstat').innerHTML = r.k!==undefined
? 'k is now <span class="text-emerald-300">'+r.k.toFixed(4)+'</span>, residual '+r.residual.toFixed(2)+' \u00b0C'
: 'Rejected: no live reading yet.';
});
el('m-calrst').addEventListener('click', async () => {
const r = await fetch('/api/calibrate', {method:'POST', headers:{'Content-Type':'application/json'},
body: JSON.stringify({reset:true})}).then(r=>r.json());
el('m-calstat').innerHTML = 'Reset to prior k = <span class="text-emerald-300">'+r.k+'</span>';
});
/* ---------------- METHODS ---------------- */
let methodsDoc = null, methodSel = 'acquire';
const STAGE_COLOUR = { acquire:'#94a3b8', compensate:'#34d399', kalman:'#f59e0b', features:'#06b6d4',
nowcast:'#818cf8', conformal:'#a78bfa', climatology:'#38bdf8', precip:'#60a5fa',
monitor:'#fb7185', verify:'#4ade80' };
async function loadMethods() {
if (!methodsDoc) methodsDoc = await fetch('/api/methods').then(r=>r.json());
drawDiagram(); drawStage();
}
function drawDiagram() {
const p = methodsDoc.pipeline;
el('me-diagram').innerHTML = p.map((s,i)=>{
const c = STAGE_COLOUR[s.id]||'#94a3b8', on = s.id===methodSel;
const edge = methodsDoc.flow.filter(f=>f.from===s.id)[0];
return '<button data-stage="'+s.id+'" class="stagebtn w-full text-left rounded-xl border px-3 py-2 transition" '+
'style="border-color:'+(on?c+'88':'rgba(255,255,255,.07)')+';background:'+(on?c+'1a':'rgba(15,23,42,.5)')+'">'+
'<div class="flex items-center gap-2.5"><span class="font-mono text-[10px] w-5 shrink-0" style="color:'+c+'">'+s.stage+'</span>'+
'<div class="min-w-0 flex-1"><div class="text-xs font-semibold '+(on?'text-white':'text-slate-300')+'">'+s.title+'</div>'+
'<div class="text-[9px] font-mono text-slate-600 truncate">'+s.technique+'</div></div></div></button>'+
(i<p.length-1 ? '<div class="flex items-center gap-1.5 pl-[26px] h-4"><div class="w-px h-full" style="background:'+c+'55"></div>'+
'<span class="text-[8px] font-mono text-slate-700">'+(edge?edge.label:'')+'</span></div>' : '');
}).join('')+
'<div class="mt-3 pt-3 border-t border-slate-800"><div class="text-[9px] font-mono text-slate-600 uppercase mb-1.5">feedback edges</div>'+
methodsDoc.flow.filter(f=>['verify','monitor','climatology'].indexOf(f.from)>=0).map(f=>
'<div class="text-[9px] font-mono text-slate-600 flex items-center gap-1.5 mb-0.5">'+
'<span class="text-slate-500">'+f.from+'</span><span class="text-indigo-500">&rarr;</span>'+
'<span class="text-slate-500">'+f.to+'</span><span class="text-slate-700">'+f.label+'</span></div>').join('')+'</div>';
document.querySelectorAll('.stagebtn').forEach(b => b.addEventListener('click', () => {
methodSel = b.dataset.stage; drawDiagram(); drawStage();
}));
}
function drawStage() {
const s = methodsDoc.pipeline.filter(x=>x.id===methodSel)[0];
if (!s) return;
const c = STAGE_COLOUR[s.id]||'#94a3b8';
el('me-head').innerHTML =
'<div class="flex flex-wrap items-baseline gap-2.5"><span class="font-mono text-xs" style="color:'+c+'">stage '+s.stage+'</span>'+
'<h2 class="text-base font-bold">'+s.title+'</h2>'+
'<span class="ml-auto font-mono text-[10px] text-slate-600">'+s.module+'</span></div>'+
'<p class="text-[11px] text-slate-500 font-mono mt-0.5">'+s.technique+'</p>';
const params = Object.keys(s.params||{}).map(k=>
'<div class="bg-slate-900/60 rounded-lg border border-slate-800/70 px-2.5 py-1.5">'+
'<div class="text-[9px] text-slate-600 font-mono uppercase">'+k+'</div>'+
'<div class="text-[11px] text-slate-200 font-mono font-semibold">'+s.params[k]+'</div></div>').join('');
el('me-body').innerHTML =
'<div class="grid grid-cols-2 gap-2 font-mono text-[10px]">'+
'<div class="bg-slate-900/40 rounded-lg border border-slate-800/60 px-2.5 py-1.5">'+
'<div class="text-slate-600 uppercase text-[9px]">consumes</div><div class="text-slate-300 mt-0.5">'+s.consumes+'</div></div>'+
'<div class="bg-slate-900/40 rounded-lg border border-slate-800/60 px-2.5 py-1.5">'+
'<div class="text-slate-600 uppercase text-[9px]">produces</div><div class="text-slate-300 mt-0.5">'+s.produces+'</div></div></div>'+
(s.math ? '<div class="bg-slate-950/60 rounded-lg border border-slate-800/70 px-3 py-2.5 overflow-x-auto">'+
'<div class="text-[9px] text-slate-600 font-mono uppercase mb-1">core relation</div>'+
'<div class="text-[11px] font-mono text-indigo-300">'+s.math.replace(/[{}\\]/g,' ').replace(/\s+/g,' ')+'</div></div>' : '')+
'<div><div class="text-[9px] text-slate-600 font-mono uppercase mb-1">why it is done this way</div>'+
'<p class="text-[12px] text-slate-300 leading-relaxed">'+s.why+'</p></div>'+
'<div class="border-l-2 pl-3" style="border-color:'+c+'66">'+
'<div class="text-[9px] font-mono uppercase mb-1" style="color:'+c+'">how it fails</div>'+
'<p class="text-[12px] text-slate-400 leading-relaxed">'+s.failure+'</p></div>'+
(params ? '<div><div class="text-[9px] text-slate-600 font-mono uppercase mb-1">live parameters</div>'+
'<div class="grid grid-cols-2 sm:grid-cols-3 gap-2">'+params+'</div></div>' : '')+
(methodSel==='verify' ? '<div class="pt-2 border-t border-slate-800">'+
'<div class="text-[9px] text-slate-600 font-mono uppercase mb-1.5">honest limits of this station</div><ul class="space-y-1.5">'+
methodsDoc.honest_limits.map(l=>'<li class="text-[11px] text-slate-400 leading-relaxed flex gap-2">'+
'<span class="text-amber-500 shrink-0">&middot;</span><span>'+l+'</span></li>').join('')+'</ul></div>' : '')+
(methodSel==='features' ? '<div class="pt-2 border-t border-slate-800">'+
'<div class="text-[9px] text-slate-600 font-mono uppercase mb-1.5">the '+methodsDoc.features.length+' features</div>'+
'<div class="flex flex-wrap gap-1">'+methodsDoc.features.map(f=>
'<span class="px-1.5 py-0.5 rounded bg-slate-900/70 border border-slate-800 text-[9px] font-mono text-slate-500">'+f+'</span>').join('')+
'</div></div>' : '')+
(methodSel==='acquire' ? '<div class="pt-2 border-t border-slate-800">'+
'<div class="text-[9px] text-slate-600 font-mono uppercase mb-1.5">glossary</div>'+
methodsDoc.glossary.map(g=>'<div class="mb-2"><span class="text-[11px] font-semibold text-slate-300">'+g.term+'</span>'+
'<p class="text-[11px] text-slate-500 leading-relaxed">'+g.definition+'</p></div>').join('')+'</div>' : '');
}
loaders.methods = loadMethods;
/* ---------------- BOOT ---------------- */
connectStream();
loadLiveChart();
fetch('/api/status').then(r=>r.json()).then(s => { el('hd-days').innerText = fmt(s.history_days,2); });
setInterval(() => { if (activeTab==='live') loadLiveChart(); }, 60000);
setInterval(() => { if (activeTab==='forecast') { loadForecast(); loadOutlook(); } }, 120000);
setInterval(() => { if (activeTab==='models') loadModels(); }, 60000);
setInterval(() => { if (activeTab==='history') loadHistory(); }, 300000);
</script>
</body>
</html>
"""
ashvale/estimation.py
+225
View File
@@ -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))
ashvale/features.py
+204
View File
@@ -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]
ashvale/led.py
+272
View File
@@ -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()
ashvale/methods.py
+359
View File
@@ -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.",
],
}
ashvale/models/__init__.py
+24
View File
@@ -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",
]
ashvale/models/anomaly.py
+287
View File
@@ -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)
ashvale/models/climatology.py
+170
View File
@@ -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"]
ashvale/models/nowcast.py
+245
View File
@@ -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)
ashvale/models/precip.py
+323
View File
@@ -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
ashvale/models/rls.py
+182
View File
@@ -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
ashvale/physics.py
+162
View File
@@ -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)
ashvale/sensors.py
+251
View File
@@ -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
ashvale/station.py
+563
View File
@@ -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"
ashvale/storage.py
+487
View File
@@ -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
config.yaml
+35
View File
@@ -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
docs/DESIGN.md
+393
View File
@@ -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.
pyproject.toml
+81
View File
@@ -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 = "[email protected]" },
]
maintainers = [
{ name = "Kemal Yaylali", email = "[email protected]" },
]
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"]
requirements.txt
+11
View File
@@ -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
run.py
+46
View File
@@ -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://<pi>: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()
scripts/evaluate.py
+171
View File
@@ -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()
scripts/simulate.py
+214
View File
@@ -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()
systemd/ashvale.service
+26
View File
@@ -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