"""Check the counting rule against USGS's own at-issue counts ('observations' in forecast.json) for 400 random versions:
count ComCat events from the mainshock time to the version's creation time inside its circle, mag >= M - eps."""
import sys; sys.path.insert(0, '/work/aftershock-forecasts')
from common import *
import random, collections
V = [v for v in csv.DictReader(open('versions.csv')) if v['rkm'] and v['obs3'] != '' and v['obs4'] != '']
random.seed(0); S = random.sample(V, 400)
print('event types in catalogs:', collections.Counter(t for v in S for t in catalog(v['eid'])['type']).most_common(6))
for eps in (0.05, 0.0):
    for eq_only in (False, True):
        ok3 = ok4 = 0; d = []
        for v in S:
            a = [float(v['rlat']), float(v['rlon']), float(v['rkm']), int(v['time_ms']), int(v['creation_ms'])]
            c3 = count(v['eid'], *a, 3, eps, eq_only); c4 = count(v['eid'], *a, 4, eps, eq_only)
            ok3 += c3 == int(v['obs3']); ok4 += c4 == int(v['obs4']); d.append(c3 - int(v['obs3']))
        d = np.array(d)
        print(f'eps={eps} earthquake-only={eq_only}: exact M3 {ok3}/400, M4 {ok4}/400; mine>USGS {np.mean(d > 0):.2f}, mine<USGS {np.mean(d < 0):.2f}')