pray-calc-ml/.github/wiki/ML-Crunching.md

ML Crunching

This page explains how to run the machine learning analysis once you have a sufficient dataset.

---

Prerequisites

Software

python -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt

Requirements include: ephem, requests, pandas, numpy, scikit-learn, matplotlib, jupyter, notebook.

Data

You need the processed CSV files in data/processed/:

python -m src.pipeline

This produces:

• data/processed/fajr_angles.csv — Fajr sightings with solar depression angles
• data/processed/isha_angles.csv — Isha sightings with solar depression angles

Without these files, the notebook will fail immediately. See Data Collection for the full pipeline guide.

---

Step 1: Exploratory Analysis

Open the notebook:

jupyter notebook notebooks/01_exploratory_analysis.ipynb

Or run it headlessly and export:

jupyter nbconvert --to notebook --execute notebooks/01_exploratory_analysis.ipynb \
    --output notebooks/01_exploratory_analysis_executed.ipynb

The notebook covers nine analyses in sequence:

Cell	Analysis	What to look for
1	Load datasets	Record counts, column dtypes
2	Angle distributions	Histogram shape — should be roughly normal for Fajr
3	Latitude vs Fajr angle	The counter-intuitive equatorial-higher pattern
4	Birmingham seasonality	Sinusoidal pattern — confirms TOY effect
5	Latitude × Season interaction	Coloured scatter — should show lat × season interaction
6	Elevation vs Fajr angle	Weaker than lat/season but visible above 500m
7	Geographic coverage map	Reveals which regions are data-sparse
8	Linear regression baseline	R² and per-feature coefficients — sets the floor for ML
9	Isha analysis	Parallel analysis for Isha; currently sparse

A well-populated dataset produces:

• Fajr angle distribution: mean ~13.5°, std ~1.8°, range roughly 8°-20°
• Fajr linear regression R² ≥ 0.35 (lat + doy + elevation)
• Latitude coefficient: negative (higher lat = lower angle at mid-latitudes)

If you see a flat distribution or R² < 0.1, check the pipeline output for dropped records.

---

Step 2: Feature Engineering

The relevant features for predicting the solar depression angle at true dawn or dusk are:

Feature	Column	Notes
Latitude	lat	Decimal degrees
sin(day of year)	derived from day_of_year	Captures seasonality (365-day cycle)
cos(day of year)	derived from day_of_year	Paired with sin for full cycle encoding
Elevation	elevation_m	Metres above sea level
abs(lat)	derived	Symmetry across equator

Do not use longitude as a feature. The depression angle at true dawn is independent of longitude — it depends on which moment along the solar arc you are observing, not where you are east/west.

Do not use the observed time as a feature. The angle is the prediction target; the time is how you derived the angle. Using it as a feature would be data leakage.

Encode day of year as a unit circle pair:

import numpy as np
df["doy_sin"] = np.sin(2 * np.pi * df["day_of_year"] / 365.25)
df["doy_cos"] = np.cos(2 * np.pi * df["day_of_year"] / 365.25)

---

Step 3: Baseline Model

Before training any ML model, establish a linear baseline:

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
import numpy as np

features = ["lat", "doy_sin", "doy_cos", "elevation_m"]
X = df[features].values
y = df["fajr_angle"].values

lr = LinearRegression()
scores = cross_val_score(lr, X, y, cv=5, scoring="r2")
print(f"Linear baseline R²: {scores.mean():.3f} ± {scores.std():.3f}")

This gives the floor — any ML model should beat it. A linear model trained on the current data produces approximately R² = 0.38.

---

Step 4: Gradient Boosting (recommended)

Gradient boosting handles the non-linear lat × season interaction without explicit feature crosses. It is the recommended first ML model for this dataset.

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score, KFold
from sklearn.metrics import mean_absolute_error
import numpy as np

features = ["lat", "doy_sin", "doy_cos", "elevation_m"]
X = df[features].values
y = df["fajr_angle"].values

model = GradientBoostingRegressor(
    n_estimators=300,
    max_depth=4,
    learning_rate=0.05,
    subsample=0.8,
    random_state=42,
)

kf = KFold(n_splits=5, shuffle=True, random_state=42)
r2_scores = cross_val_score(model, X, y, cv=kf, scoring="r2")
mae_scores = -cross_val_score(model, X, y, cv=kf, scoring="neg_mean_absolute_error")

print(f"R²:  {r2_scores.mean():.3f} ± {r2_scores.std():.3f}")
print(f"MAE: {mae_scores.mean():.3f}° ± {mae_scores.std():.3f}°")

Target performance with a well-populated dataset (10k+ records):

• R² ≥ 0.55
• MAE ≤ 0.9°

---

Step 5: Evaluating the Model

Residual analysis

from sklearn.model_selection import cross_val_predict
import matplotlib.pyplot as plt

model.fit(X, y)
y_pred = cross_val_predict(model, X, y, cv=5)
residuals = y - y_pred

plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.scatter(y_pred, residuals, alpha=0.3, s=10)
plt.axhline(0, color="red")
plt.xlabel("Predicted angle (°)")
plt.ylabel("Residual (°)")
plt.title("Residuals vs Predicted")

plt.subplot(1, 2, 2)
plt.scatter(df["lat"], residuals, alpha=0.3, s=10)
plt.axhline(0, color="red")
plt.xlabel("Latitude")
plt.ylabel("Residual (°)")
plt.title("Residuals vs Latitude")
plt.tight_layout()
plt.show()

Watch for:

• Systematic residuals at high latitudes (55°N+) — the model may underfit
• Residuals correlated with season at a single location — the model may underfit seasonality
• Outliers > 3° from the line — these may be data entry errors or unusual atmospheric events

Leave-location-out cross-validation

Standard k-fold mixes records from the same location across train/test splits, making the model look better than it generalises to new locations. For this dataset, location-aware CV is more informative:

from sklearn.model_selection import LeaveOneGroupOut
import numpy as np

# Group by location (round lat/lng to 1 decimal for grouping)
groups = (df["lat"].round(1).astype(str) + "," + df["lng"].round(1).astype(str))

logo = LeaveOneGroupOut()
scores = cross_val_score(model, X, y, cv=logo, groups=groups, scoring="r2")
print(f"Leave-location-out R²: {scores.mean():.3f} ± {scores.std():.3f}")

This tests whether the model generalises to locations it has never seen.

---

Step 6: Feature Importance

model.fit(X, y)
importances = model.feature_importances_

for name, imp in zip(features, importances):
    print(f"  {name}: {imp:.3f}")

Expected order: doy_sin or doy_cos highest, then lat, then elevation_m lowest. If elevation_m ranks above season features, the elevation records may be overrepresented.

---

Step 7: Exporting the Model

Once satisfied with validation performance:

import joblib
import json
import numpy as np

model.fit(X, y)

joblib.dump(model, "models/fajr_gbm.pkl")

# Export feature ranges for the pray-calc DPC algorithm
meta = {
    "features": features,
    "lat_range": [float(df["lat"].min()), float(df["lat"].max())],
    "elevation_range": [float(df["elevation_m"].min()), float(df["elevation_m"].max())],
    "angle_mean": float(y.mean()),
    "angle_std": float(y.std()),
    "n_records": int(len(df)),
    "r2_cv": float(r2_scores.mean()),
    "mae_cv": float(mae_scores.mean()),
}
with open("models/fajr_gbm_meta.json", "w") as f:
    json.dump(meta, f, indent=2)

print(f"Saved fajr_gbm.pkl ({len(df)} training records, R²={r2_scores.mean():.3f})")

---

Current Model Status

The current dataset has:

• Fajr: ~4,100 records, but 98% are from Birmingham, UK. The model heavily reflects one location.
• Isha: ~43 records. Not enough to train a reliable ML model.

The priority is data collection before further ML work. A model trained only on Birmingham Fajr data will predict Birmingham well and generalise poorly. The notebook's exploratory analysis and linear baseline are meaningful now, but gradient boosting should wait for broader geographic coverage.

Target before training a production model:

• Fajr: 10,000+ records from 100+ locations across all latitude bands
• Isha: 500+ records from 30+ locations

See Data Collection for how to contribute new sighting records.

---

Connecting to pray-calc

The output of the ML model feeds the DPC (Dynamic Prayer Calc) algorithm in pray-calc. The DPC algorithm takes:

• Latitude
• Day of year
• Elevation

And returns a recommended depression angle for that location and date.

The current DPC implementation uses a simplified physics model. The ML model will replace or calibrate the seasonal and latitude correction factors once sufficient data is available.

---

← Data Collection · Architecture →