opt-hist.py

from datetime import datetime, timezone, timedelta
import igrf12
import json
import numpy as np
import os
import re
from scipy.optimize import minimize_scalar
import subprocess
import urllib.request
import sys
import multiprocessing

now = datetime.utcnow()

data = None

def station_weight(station, tm):
    hour = float(tm.hour) + float(tm.minute) / 60. + float(tm.second) / 3600.
    local_time = (hour + float(station['longitude']) / 15.) % 24.
    abs_dip = abs(station['dip_angle'])

    time_weight = np.sqrt(2 + np.cos((local_time - 14.5) * np.pi / 12.)) - 0.73
    if abs_dip < 25.:
        return 0.75 * time_weight
    elif abs_dip < 50.:
        return 1. * time_weight
    elif abs_dip < 56.:
        return 0.75 * time_weight
    else:
        return 0.5 * time_weight

def recency_weight(tm):
    delta = now - tm
    return np.power(2.0, -(delta / timedelta(hours=6)))

def cs_to_stdev(cs):
    if cs == -1 or cs > 100:
        cs = 69
    if cs == 100:
        cs = 81
    return 0.191 - 0.00147 * cs

iri = None
pred_pool = None

def get_pred(station, ssn, tm):
    global iri
    if iri is None:
        iri = subprocess.Popen(
                ['/src/build/iri_opt'],
                stdin=subprocess.PIPE,
                stdout=subprocess.PIPE,
                universal_newlines=True,
                )

    iri.stdin.write("%f %f %f %s\n" % (station['latitude'], station['longitude'], ssn, tm.strftime("%Y %m %d %H %M %S")))
    iri.stdin.flush()
    data = [float(x) for x in iri.stdout.readline().split()]
    return (data[3], data[6], data[5]) # fof2, hmf2, mufd

def station_err(x):
    station, ssn = x

    station_total = 0.0
    station_tw = 0.0
    num_records = 0
    last_rw = 0.0

    if len(station['history']) < 3:
        return (0,0)

    last_tm = station['history'][-1][0]
    if (now - last_tm).total_seconds() > 7200:
        return (0,0)

    for record in station['history']:
        tm, cs, fof2, mufd, hmf2 = record

        if cs <= 25:
            continue

        fof2_pred, hmf2_pred, mufd_pred = get_pred(station, ssn, tm)
        stdev = cs_to_stdev(cs)

        err = 0.6 * np.square(np.log(fof2_pred) - np.log(fof2))
        err += 0.3 * np.square(np.log(mufd_pred) - np.log(mufd))
        err += 0.1 * np.square(np.log(hmf2_pred) - np.log(hmf2))
        err = np.sqrt(err) / stdev

        sw = station_weight(station, tm)
        rw = recency_weight(tm)
        weight = sw * rw
        last_rw = rw

        station_total += err * weight
        station_tw += weight
        num_records += 1

    if num_records >= 3:
        return (last_rw * station_total / float(num_records), last_rw * station_tw / float(num_records))
    else:
        return (0,0)


def err(ssn):
    total = 0.0
    total_weight = 0.0

#    print(ssn, file=sys.stderr)

    results = pred_pool.imap_unordered(station_err, [ (station, ssn) for station in data ])
    for result in results:
        total += result[0]
        total_weight += result[1]

    return total / total_weight

if __name__ == '__main__':
    with urllib.request.urlopen(os.getenv("HISTORY_URI")) as res:
        data = json.loads(res.read().decode())

    for station in data:
        mag = igrf12.igrf(
                '2019-12-31',
                glat=station['latitude'],
                glon=station['longitude'],
                alt_km=1.
                )

        station['dip_angle'] = mag['incl'].item()

        for record in station['history']:
            record[0] = re.sub(r"\.\d+$", "", record[0])
            record[0] = datetime.strptime(record[0], '%Y-%m-%d %H:%M:%S')


    pred_pool = multiprocessing.Pool(4)
    res = minimize_scalar(err, bounds=(-20.0, 200.0), method='Bounded', options={'xatol':0.1})
    ssn = res.x
    sfi = 63.75 + ssn * (0.725 + ssn*0.000089)
    print("%s\t%f\t%f\t%f\t%f" % (now.strftime("%s"), sfi, ssn, res.fun, err(-100.0)))