HRV.py

#!/usr/bin/env python
# coding: utf-8
# SPDX-FileCopyrightText: 2023 Julian Andres Klode <jak@jak-linux.org>
#
# SPDX-License-Identifier: AGPL-3.0-or-later

# In[7]:


#!/usr/bin/python3
import json
import requests
import datetime
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sys
from matplotlib.ticker import FormatStrFormatter
from matplotlib.ticker import AutoLocator, ScalarFormatter


VIEW = int(sys.argv[1]) if len(sys.argv) > 1 else 60
SCORE = "noscore" not in sys.argv
WORKOUTS="noworkouts" not in sys.argv
today = datetime.date.today()
start = (today - datetime.timedelta(days=VIEW + 80))
first = start + datetime.timedelta(days=1)


url = "https://api.ouraring.com/v2/usercollection/sleep"
params = {
    "start_date": start.isoformat(),
    "end_date": today.isoformat(),
}
headers = {
    "Authorization": "Bearer " + open("/home/jak/Private/oura-token").read().strip()
}
try:
    with open("/home/jak/.cache/oura.json") as cached:
        data = json.load(cached)
    if datetime.date.fromisoformat(data["data"][-1]["day"]) < today:
        raise ValueError("Out of date")
    if datetime.date.fromisoformat(data["data"][0]["day"]) > first:
        raise ValueError("Out of date")
except (ValueError, FileNotFoundError) as e:
    print("Fetching due to error:", e)
    response = requests.request("GET", url, headers=headers, params=params)
    data = response.json()
    with open("/home/jak/.cache/oura.json", "w") as cached:
        json.dump(data, cached)


dates = []
hrs = []
hrvs = []
workouts = pd.read_json("/home/jak/.cache/strava.json", convert_dates=["start_date"])
for sleep in data["data"]:
    if sleep["type"] != "long_sleep":
        continue
    dates.append(pd.Timestamp(sleep["day"]))
    # hrs.append(sleep["average_heart_rate"])
    hrs.append(np.round(np.mean([x for x in sleep["heart_rate"]["items"] if x])))
    if SCORE:
        if sleep["average_hrv"]:
            #hrvs.append(np.round(np.log(sleep["average_hrv"]), 1))
            #hrvs.append(np.round(1.7 * np.log(sleep["average_hrv"]) + 1, 1))
            #hrvs.append(np.round(1.75202 * np.log(sleep["average_hrv"] * 1.57691), 1))
            hrvs.append(np.round(1.73677755 * np.log(sleep["average_hrv"] * 1.63747424), 1))
            #1.73677755, 1.63747424
        else:
            hrvs.append(None)
    else:
        hrvs.append(sleep["average_hrv"])


def normal_range(values):
    rolling = values.dropna().rolling(window="60d")
    mean = rolling.mean()
    std = rolling.std()
    return (mean - std), mean, (mean + std)


def plot_hrvs(hrvs, label):
    # some sample data
    ts = pd.Series(hrvs, index=dates)

    digits = 0
    if label == "HRV" and SCORE:
        digits = 1
    if "coeff" in label:
        digits = 2

    rolling = ts.dropna().rolling(window="7d")
    mean = rolling.mean()
    min_range, mean_range, max_range = (r.round(digits) for r in normal_range(ts))

    coeff = rolling.std() / rolling.mean() * 100
    min_coeff_range, mean_coeff_range, max_coeff_range = normal_range(coeff)

    table = (
        ts.to_frame(name="HRV")
        .join(min_range.round(digits).to_frame(name="MIN"))
        .join(max_range.round(digits).to_frame(name="MAX"))
        .join(mean.to_frame(name="MEAN"))
    )
    table['7day'] = mean
    table_coeff = (
        coeff.round(1).to_frame(name="Coeff")
        .join(min_coeff_range.round(1).to_frame(name="MIN"))
        .join(max_coeff_range.round(1).to_frame(name="MAX"))
    )

    if "coeff" in label:
        coeff.plot(label=label, style=".-")
        max_coeff_range.plot(label="60d mean + stddev", color="red")
        min_coeff_range.plot(label="60d mean - stddev", color="green")
    else:

        if VIEW > 160:
            ts.plot(style="o--", color="#" + ("d" * 6), zorder=1, label=label)
        else:
            ts.plot(style="o--", color="#" + ("a" * 6), zorder=1, label=label)
        mean.plot(style="-", label=f"{label} 7 day")
        #mean_range.plot(label=f"{label} 60 day", color="#333333", style="--", zorder=2)
        max_range.plot(label="60d mean + stddev", color="green" if "HRV" in label else "red")
        min_range.plot(label="60d mean - stddev", color="red" if "HRV" in label else "green")
    if "coeff" in label:
        return
    for index, entry in table.iterrows():
        if index.date() == today:
            print(
                    label + ":      ",
                index,
                entry["HRV"],
                "normal range",
                entry["MIN"],
                entry["MAX"],
                "" if entry["MIN"] < entry["HRV"] < entry["MAX"] else "(!)",
                sep="\t",
            )
            print(
                    "7 day " + label + ":",
                index,
                entry["7day"].round(digits),
                "normal range",
                entry["MIN"],
                entry["MAX"],
                "" if entry["MIN"] < entry["7day"].round(digits) < entry["MAX"] else "(!)",
                sep="\t",
            )
    for index, entry in table_coeff.iterrows():
        if index.date() == today:
            print(
                    label + " coeff:",
                index,
                entry["Coeff"],
                "normal range",
                entry["MIN"],
                entry["MAX"],
                "" if entry["MIN"] < entry["Coeff"] < entry["MAX"] else "(!)",
                sep="\t",
            )

# In[9]:


#fig=plt.figure(figsize=(40, 3*25))
if VIEW > 30:
    fig=plt.figure(figsize=(8.27*2,11.69*2)) # for landscape
else:
    fig=plt.figure(figsize=(8.27,11.69)) # for landscape


plot_id = 410 if WORKOUTS else 310

# plt.show()
if WORKOUTS:
    plot_id += 1
    fig.add_subplot(plot_id, title="Workouts")
    ts = pd.Series(hrvs, index=dates)
    efforts = ts.to_frame("HRV").join(workouts.groupby([workouts["start_date"].dt.date])["suffer_score"].sum()).fillna(0)["suffer_score"]
    #plt.gca().set_yscale('log')
    #plt.gca().yaxis.set_major_locator(AutoLocator())
    #plt.gca().yaxis.set_major_formatter(ScalarFormatter())
    if "notrend" in sys.argv or (VIEW < 120 and "trend" not in sys.argv):
        efforts.plot(style=".--", label="Relative Effort")
        plt.ylim([0, efforts.rolling(VIEW).max()[-1]*1.1])
    else:
        efforts.ewm(span=60).mean().plot(style="-", label="Relative Effort 60 day EWM (fitness)")
        efforts.ewm(span=7).mean().plot(style="-", label="Relative Effort 7 day EWM (fatigue)", zorder=1, color="#dddddd")
        plt.ylim([0, efforts.ewm(span=7).mean().rolling(VIEW).max()[-1]])
        #plt.ylim([0, efforts.ewm(span=60).mean().rolling(VIEW).max()[-1]])
    #(efforts.rolling(7).mean().max() / efforts.max() * efforts).plot(style="-", label="Relative Effort 7 day mean", zorder=1, color="#000000")
    plt.xlim([today - datetime.timedelta(days=VIEW), today + datetime.timedelta(days=0)])
    plt.gca().minorticks_off()
    plt.gca().legend()
    plt.grid()

plot_id += 1
fig.add_subplot(plot_id, title="HRV")
plot_hrvs(hrvs, "HRV")
plt.xlim([today - datetime.timedelta(days=VIEW), today + datetime.timedelta(days=0)])
plt.gca().legend()

if not SCORE:
    plt.gca().set_yscale('log')
    plt.gca().yaxis.set_major_locator(AutoLocator())
    plt.gca().yaxis.set_major_formatter(ScalarFormatter())
plt.gca().minorticks_off()

plt.grid()

# In[4]:
if 1:
    plot_id += 1
    fig.add_subplot(plot_id, title="HRV coeff")
    plot_hrvs(hrvs, "HRV coeff")
    plt.xlim([today - datetime.timedelta(days=VIEW), today + datetime.timedelta(days=0)])
    plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%d%%'))
    plt.gca().legend()
    plt.grid(axis="x")


#plt.figure(figsize=(38, 24))
plot_id += 1
fig.add_subplot(plot_id, title="HR")
plot_hrvs(hrs, "HR")
plt.gca().legend()
plt.xlim([today - datetime.timedelta(days=VIEW), today + datetime.timedelta(days=0)])
plt.gca().yaxis.set_major_locator(plt.MultipleLocator(1))
plt.grid()

if VIEW >= 210:
    from matplotlib.dates import MO, MonthLocator
    for ax in plt.gcf().get_axes():
        ax.xaxis.set_major_locator(MonthLocator())
else:
    from matplotlib.dates import MO, WeekdayLocator
    for ax in plt.gcf().get_axes():
        ax.xaxis.set_major_locator(WeekdayLocator(byweekday=MO, interval=max(VIEW//70, 1)) )
if VIEW <= 30:
    fig.autofmt_xdate()
plt.text(0.05,0.95, f"{VIEW} day HRV analysis", transform=fig.transFigure, size=24)
plt.savefig("Figure.png",  orientation = 'portrait', format = 'png', dpi=300)
plt.savefig("Figure.pdf",  orientation = 'portrait', format = 'pdf', dpi=300)
plt.savefig("Figure.svg",  orientation = 'portrait', format = 'svg', dpi=300)
#plt.show()

This used to be a Jupyter notebook so it is a bit whacky.