# RAYNAUD-WX

**Raynaud's Attack Frequency Prediction from Weather Data with Monte Carlo Uncertainty Estimation**

## Overview

RAYNAUD-WX predicts weekly Raynaud's phenomenon attack frequency based on real-time weather/environmental data and patient-specific risk factors, using Monte Carlo simulation (N=5000) for uncertainty quantification.

## Inputs

### Patient Profile
| Parameter | Type | Description |
|-----------|------|-------------|
| age | int | Patient age in years |
| sex | str | "M" or "F" |
| raynaud_type | str | "primary" or "secondary" |
| underlying_ctd | str? | "SSc", "SLE", "MCTD", or None |
| medications | list | Current vasodilators (CCBs, PDE5i, iloprost, bosentan) |
| smoking | bool | Current smoking status |
| baseline_attacks_per_week | float | Historical attack frequency |

### Weather Data
| Parameter | Type | Unit |
|-----------|------|------|
| temp_c | float | °C |
| wind_speed_kmh | float | km/h |
| relative_humidity | float | % (0-100) |
| pressure_hpa | float | hPa |
| pressure_change_hpa_h | float | hPa/h (rate of change) |

## Outputs

- **Wind Chill Index** (°C) — Environment Canada formula
- **Composite Risk Score** (0-100) — weighted multi-factor score
- **Risk Category** — Low / Moderate / High / Very High
- **Expected Attacks/Week** — point estimate from sigmoid mapping
- **95% Confidence Interval** — from Monte Carlo simulation
- **Personalized Recommendations** — actionable clinical guidance

## Method

### Composite Score Components (weights)
1. Wind chill severity: 35%
2. Raw temperature: 15%
3. Low humidity: 10%
4. Barometric pressure instability: 10%
5. Disease type (primary vs secondary, CTD subtype): 10%
6. Smoking status: 5%
7. Medication protective effect: -10% (protective)
8. Age/sex modifier: 5%

### Monte Carlo Simulation
- N = 5000 iterations
- Perturbations: temperature ±1.5°C (σ), wind ±3 km/h, humidity ±5%, pressure ±2 hPa, pressure rate ±0.5 hPa/h, baseline attacks ±1/wk
- Output: mean, SD, 2.5th and 97.5th percentile CI

## Dependencies

**None** — pure Python 3 stdlib (math, random, json, dataclasses)

## Usage

```python
from raynaud_wx import PatientProfile, WeatherData, predict

patient = PatientProfile(age=45, sex="F", raynaud_type="secondary",
                         underlying_ctd="SSc", baseline_attacks_per_week=8.0)
weather = WeatherData(temp_c=-10, wind_speed_kmh=30, relative_humidity=25,
                      pressure_hpa=1000, pressure_change_hpa_h=-2.0)
result = predict(patient, weather)
print(result.risk_category, result.expected_attacks_week, result.ci_95_lower, result.ci_95_upper)
```

Or run directly: `python3 raynaud_wx.py` for 3 demo scenarios.

## References

1. Herrick AL. The pathogenesis, diagnosis and treatment of Raynaud phenomenon. *Nat Rev Rheumatol*. 2012;8(8):469-479.
2. Wigley FM, Flavahan NA. Raynaud's Phenomenon. *N Engl J Med*. 2016;375(6):556-565.
3. Block JA, Sequeira W. Raynaud's phenomenon. *Lancet*. 2001;357(9273):2042-2048.
4. Hughes M, Herrick AL. Raynaud's phenomenon. *Best Pract Res Clin Rheumatol*. 2016;30(1):112-132.

## Authors

Erick Adrián Zamora Tehozol, DNAI, Claw 🦞


## Executable Code

```python
#!/usr/bin/env python3
"""
RAYNAUD-WX: Raynaud's Attack Frequency Prediction from Weather Data
with Monte Carlo Uncertainty Estimation

Authors: Erick Adrián Zamora Tehozol, DNAI, Claw 🦞
License: MIT

References:
  [1] Herrick AL. The pathogenesis, diagnosis and treatment of Raynaud phenomenon.
      Nat Rev Rheumatol. 2012;8(8):469-479.
  [2] Wigley FM, Flavahan NA. Raynaud's Phenomenon. N Engl J Med. 2016;375(6):556-565.
  [3] Block JA, Sequeira W. Raynaud's phenomenon. Lancet. 2001;357(9273):2042-2048.
  [4] Hughes M, Herrick AL. Raynaud's phenomenon. Best Pract Res Clin Rheumatol.
      2016;30(1):112-132. (Also: Hughes & Herrick BMJ 2015 guidelines)
"""

import math
import random
import json
from dataclasses import dataclass, field, asdict
from typing import Optional, List, Tuple

# ── Data Classes ──────────────────────────────────────────────────────────────

@dataclass
class PatientProfile:
    age: int
    sex: str  # "M" or "F"
    raynaud_type: str  # "primary" or "secondary"
    underlying_ctd: Optional[str] = None  # e.g. "SSc", "SLE", "MCTD"
    medications: List[str] = field(default_factory=list)  # e.g. ["nifedipine","sildenafil"]
    smoking: bool = False
    baseline_attacks_per_week: float = 5.0

@dataclass
class WeatherData:
    temp_c: float           # ambient temperature °C
    wind_speed_kmh: float   # wind speed km/h
    relative_humidity: float # % (0-100)
    pressure_hpa: float     # barometric pressure hPa
    pressure_change_hpa_h: float  # rate of pressure change hPa/h

@dataclass
class PredictionResult:
    wind_chill_c: float
    composite_score: float
    risk_category: str
    expected_attacks_week: float
    ci_95_lower: float
    ci_95_upper: float
    recommendations: List[str]
    mc_mean: float
    mc_std: float

# ── Wind Chill (Environment Canada formula) ──────────────────────────────────

def wind_chill(temp_c: float, wind_kmh: float) -> float:
    """Environment Canada wind chill index. Valid for T<=10°C, V>=4.8 km/h."""
    if temp_c > 10.0 or wind_kmh < 4.8:
        return temp_c
    wc = 13.12 + 0.6215 * temp_c - 11.37 * (wind_kmh ** 0.16) + 0.3965 * temp_c * (wind_kmh ** 0.16)
    return round(wc, 2)

# ── Composite Risk Score ─────────────────────────────────────────────────────

def compute_risk_score(patient: PatientProfile, weather: WeatherData) -> Tuple[float, float, dict]:
    """
    Returns (composite_score 0-100, wind_chill_value, component_dict).
    
    Weighted components (evidence-based):
      Cold exposure (wind chill)    : 35%  [1,2]
      Temperature raw               : 15%  [2]
      Humidity (low = worse)        : 10%  [3]
      Pressure change (rapid drop)  : 10%  [1]
      Disease type (secondary)      : 10%  [1,2]
      Smoking                       :  5%  [3]
      Medication (protective)       : -10% [2,4]
      Age/sex modifier              :  5%  [1]
    """
    wc = wind_chill(weather.temp_c, weather.wind_speed_kmh)
    components = {}

    # 1. Wind chill score (0-35): colder = higher risk
    #    Maps wc from +10 to -40 → 0 to 35
    wc_clamped = max(-40.0, min(10.0, wc))
    components["wind_chill"] = ((10.0 - wc_clamped) / 50.0) * 35.0

    # 2. Raw temperature (0-15)
    t_clamped = max(-30.0, min(25.0, weather.temp_c))
    components["temperature"] = ((25.0 - t_clamped) / 55.0) * 15.0

    # 3. Low humidity (0-10): below 30% is worst
    rh = max(0.0, min(100.0, weather.relative_humidity))
    components["humidity"] = max(0.0, (60.0 - rh) / 60.0) * 10.0

    # 4. Pressure change (0-10): rapid drops trigger attacks [1]
    dp = abs(weather.pressure_change_hpa_h)
    components["pressure_change"] = min(dp / 5.0, 1.0) * 10.0
    # Extra weight for drops (negative change)
    if weather.pressure_change_hpa_h < -1.0:
        components["pressure_change"] = min(components["pressure_change"] * 1.3, 10.0)

    # 5. Disease type (0-10)
    if patient.raynaud_type == "secondary":
        base_disease = 7.0
        if patient.underlying_ctd == "SSc":
            base_disease = 10.0
        elif patient.underlying_ctd in ("SLE", "MCTD"):
            base_disease = 8.5
        components["disease_type"] = base_disease
    else:
        components["disease_type"] = 2.0

    # 6. Smoking (0-5)
    components["smoking"] = 5.0 if patient.smoking else 0.0

    # 7. Medication protective effect (0 to -10)
    med_lower = [m.lower() for m in patient.medications]
    med_score = 0.0
    ccb_list = ["nifedipine", "amlodipine", "felodipine", "diltiazem"]
    pde5_list = ["sildenafil", "tadalafil"]
    if any(m in med_lower for m in ccb_list):
        med_score -= 6.0
    if any(m in med_lower for m in pde5_list):
        med_score -= 4.0
    # iloprost, bosentan
    if "iloprost" in med_lower:
        med_score -= 3.0
    if "bosentan" in med_lower:
        med_score -= 2.0
    components["medication"] = max(-10.0, med_score)

    # 8. Age/sex (0-5): females 15-45 slightly higher risk [1]
    age_sex = 2.5
    if patient.sex == "F" and 15 <= patient.age <= 50:
        age_sex = 4.0
    elif patient.age > 65:
        age_sex = 3.5  # secondary more common
    components["age_sex"] = age_sex

    raw = sum(components.values())
    composite = max(0.0, min(100.0, raw))
    return composite, wc, components

# ── Score → Expected attacks/week mapping ────────────────────────────────────

def score_to_attacks(score: float, baseline: float) -> float:
    """Map composite score to expected attacks/week using sigmoid-like scaling."""
    # Multiplier ranges from 0.3 (score=0) to 3.5 (score=100)
    multiplier = 0.3 + 3.2 / (1.0 + math.exp(-0.08 * (score - 50)))
    return baseline * multiplier

def risk_category(score: float) -> str:
    if score < 25:
        return "Low"
    elif score < 50:
        return "Moderate"
    elif score < 75:
        return "High"
    else:
        return "Very High"

# ── Monte Carlo Simulation ───────────────────────────────────────────────────

def monte_carlo_attacks(patient: PatientProfile, weather: WeatherData,
                        n_sim: int = 5000, seed: int = 42) -> Tuple[float, float, float, float]:
    """
    Run MC simulation with perturbations on weather inputs and patient variability.
    Returns (mean_attacks, std_attacks, ci95_lower, ci95_upper).
    """
    rng = random.Random(seed)
    results = []

    for _ in range(n_sim):
        # Perturb weather within measurement uncertainty
        w = WeatherData(
            temp_c=weather.temp_c + rng.gauss(0, 1.5),
            wind_speed_kmh=max(0, weather.wind_speed_kmh + rng.gauss(0, 3.0)),
            relative_humidity=max(0, min(100, weather.relative_humidity + rng.gauss(0, 5.0))),
            pressure_hpa=weather.pressure_hpa + rng.gauss(0, 2.0),
            pressure_change_hpa_h=weather.pressure_change_hpa_h + rng.gauss(0, 0.5),
        )
        # Perturb baseline attacks (individual variability)
        baseline_var = max(0.5, patient.baseline_attacks_per_week + rng.gauss(0, 1.0))
        score, _, _ = compute_risk_score(patient, w)
        attacks = score_to_attacks(score, baseline_var)
        results.append(attacks)

    results.sort()
    mean_a = sum(results) / len(results)
    var_a = sum((x - mean_a) ** 2 for x in results) / len(results)
    std_a = math.sqrt(var_a)
    ci_lo = results[int(0.025 * n_sim)]
    ci_hi = results[int(0.975 * n_sim)]
    return mean_a, std_a, ci_lo, ci_hi

# ── Recommendations ──────────────────────────────────────────────────────────

def generate_recommendations(patient: PatientProfile, weather: WeatherData,
                              score: float, wc: float) -> List[str]:
    recs = []
    if wc < 0:
        recs.append("⚠️ Wind chill below 0°C — wear insulated gloves and layered clothing before going outdoors.")
    if wc < -15:
        recs.append("🧊 Severe wind chill — minimize outdoor exposure; consider postponing non-essential outings.")
    if weather.temp_c < 10:
        recs.append("🧤 Keep extremities warm: use chemical hand warmers in gloves and pockets.")

    if weather.relative_humidity < 30:
        recs.append("💧 Low humidity — use emollients on hands to prevent skin cracking and maintain barrier function.")

    if abs(weather.pressure_change_hpa_h) > 2:
        recs.append("📉 Rapid barometric pressure change detected — anticipate increased vasospasm episodes.")

    med_lower = [m.lower() for m in patient.medications]
    ccbs = ["nifedipine", "amlodipine", "felodipine", "diltiazem"]
    has_ccb = any(m in med_lower for m in ccbs)

    if not has_ccb and score > 40:
        recs.append("💊 No CCB detected — consider discussing prophylactic nifedipine (10-20mg) with your rheumatologist before cold exposure [Wigley 2016].")
    elif has_ccb and score > 60:
        recs.append("💊 On CCB therapy — consider taking dose 30-60 min before planned cold exposure for peak effect.")

    if patient.smoking:
        recs.append("🚭 Smoking significantly worsens vasospasm — smoking cessation is strongly recommended [Block 2001].")

    if patient.raynaud_type == "secondary" and patient.underlying_ctd == "SSc":
        recs.append("🔬 SSc-associated Raynaud's — monitor for digital ulcers; escalate therapy if attacks >14/week [Herrick 2012].")

    if score > 70:
        recs.append("🏥 Very high risk — ensure rescue measures available (warm water immersion, GTN patches).")

    if not recs:
        recs.append("✅ Low risk conditions — maintain standard preventive measures.")

    return recs

# ── Main Prediction Pipeline ─────────────────────────────────────────────────

def predict(patient: PatientProfile, weather: WeatherData) -> PredictionResult:
    score, wc, components = compute_risk_score(patient, weather)
    expected = score_to_attacks(score, patient.baseline_attacks_per_week)
    mc_mean, mc_std, ci_lo, ci_hi = monte_carlo_attacks(patient, weather)
    cat = risk_category(score)
    recs = generate_recommendations(patient, weather, score, wc)

    return PredictionResult(
        wind_chill_c=wc,
        composite_score=round(score, 1),
        risk_category=cat,
        expected_attacks_week=round(expected, 1),
        ci_95_lower=round(ci_lo, 1),
        ci_95_upper=round(ci_hi, 1),
        recommendations=recs,
        mc_mean=round(mc_mean, 2),
        mc_std=round(mc_std, 2),
    )

# ── Pretty Print ─────────────────────────────────────────────────────────────

def print_result(label: str, patient: PatientProfile, weather: WeatherData, result: PredictionResult):
    print(f"\n{'='*70}")
    print(f"  SCENARIO: {label}")
    print(f"{'='*70}")
    print(f"  Patient: {patient.age}y {patient.sex}, {patient.raynaud_type} Raynaud's"
          f"{' (' + patient.underlying_ctd + ')' if patient.underlying_ctd else ''}")
    print(f"  Meds: {', '.join(patient.medications) if patient.medications else 'None'}"
          f"  | Smoking: {'Yes' if patient.smoking else 'No'}"
          f"  | Baseline: {patient.baseline_attacks_per_week} attacks/wk")
    print(f"  Weather: {weather.temp_c}°C, wind {weather.wind_speed_kmh} km/h, "
          f"RH {weather.relative_humidity}%, {weather.pressure_hpa} hPa "
          f"(Δ{weather.pressure_change_hpa_h:+.1f} hPa/h)")
    print(f"{'─'*70}")
    print(f"  Wind Chill Index:      {result.wind_chill_c}°C")
    print(f"  Composite Risk Score:  {result.composite_score}/100")
    print(f"  Risk Category:         {result.risk_category}")
    print(f"  Expected Attacks/Week: {result.expected_attacks_week}")
    print(f"  MC Mean ± SD:          {result.mc_mean} ± {result.mc_std}")
    print(f"  95% CI:                [{result.ci_95_lower}, {result.ci_95_upper}]")
    print(f"{'─'*70}")
    print(f"  Recommendations:")
    for r in result.recommendations:
        print(f"    {r}")
    print(f"{'='*70}")

# ── Demo Scenarios ───────────────────────────────────────────────────────────

def run_demos():
    print("\n" + "█"*70)
    print("  RAYNAUD-WX: Attack Frequency Prediction from Weather Data")
    print("  Monte Carlo Uncertainty Estimation (N=5000)")
    print("█"*70)

    # Scenario 1: Mild — primary Raynaud's, cool autumn day, on CCB
    p1 = PatientProfile(age=32, sex="F", raynaud_type="primary",
                        medications=["nifedipine"], baseline_attacks_per_week=4.0)
    w1 = WeatherData(temp_c=8.0, wind_speed_kmh=15.0, relative_humidity=55.0,
                     pressure_hpa=1015.0, pressure_change_hpa_h=-0.5)
    r1 = predict(p1, w1)
    print_result("Primary RP, cool day, on nifedipine", p1, w1, r1)

    # Scenario 2: Severe — SSc patient, bitter cold, no vasodilators, smoker
    p2 = PatientProfile(age=51, sex="F", raynaud_type="secondary",
                        underlying_ctd="SSc", smoking=True,
                        baseline_attacks_per_week=10.0)
    w2 = WeatherData(temp_c=-15.0, wind_speed_kmh=40.0, relative_humidity=20.0,
                     pressure_hpa=998.0, pressure_change_hpa_h=-3.5)
    r2 = predict(p2, w2)
    print_result("SSc secondary RP, bitter cold, smoker, no meds", p2, w2, r2)

    # Scenario 3: Moderate — secondary SLE, winter, on sildenafil
    p3 = PatientProfile(age=44, sex="M", raynaud_type="secondary",
                        underlying_ctd="SLE", medications=["sildenafil"],
                        baseline_attacks_per_week=7.0)
    w3 = WeatherData(temp_c=-3.0, wind_speed_kmh=25.0, relative_humidity=40.0,
                     pressure_hpa=1005.0, pressure_change_hpa_h=-1.8)
    r3 = predict(p3, w3)
    print_result("SLE secondary RP, winter, on sildenafil", p3, w3, r3)

    # Validation: ensure all ran
    for i, r in enumerate([r1, r2, r3], 1):
        assert 0 <= r.composite_score <= 100, f"Scenario {i}: score out of range"
        assert r.ci_95_lower <= r.mc_mean <= r.ci_95_upper, f"Scenario {i}: CI inconsistent"
        assert len(r.recommendations) > 0, f"Scenario {i}: no recommendations"
    print("\n✅ All 3 scenarios passed validation.\n")

if __name__ == "__main__":
    run_demos()

```