#!/usr/bin/env python3
"""Independent numerical audit of ABLAQ-R05 (internal verification, not peer review).
Run: python3 -B verify-weather.v1.py --inputs inputs --results results --output audit.json
Uses the same pvlib astronomical library, separate polygon clipping and scalar
seasonal accumulation. Does not import the study implementation.
CC0-1.0; ABLAQ Studios.
"""
import argparse, csv, gzip, hashlib, json, math
from pathlib import Path
import numpy as np
import pandas as pd
from pvlib.solarposition import spa_python

def clipped_fraction(e,n,z,ratio,height,bearing):
    if z<=0:return 0.
    w,l=math.sqrt(100/ratio),math.sqrt(100*ratio)
    a=math.radians(bearing)
    shift=(height*(e*math.cos(a)-n*math.sin(a))/z,
           height*(e*math.sin(a)+n*math.cos(a))/z)
    poly=[(0.,0.),(w,0.),(w,l),(0.,l)]
    # Four half planes constrain the upward ray's intersection with the opening.
    for axis,edge,sign in ((0,0,1),(0,w,-1),(1,0,1),(1,l,-1)):
        if not poly:return 0.
        new=[]
        for a,b in zip(poly,poly[1:]+poly[:1]):
            sa=sign*(a[axis]+shift[axis]-edge)
            sb=sign*(b[axis]+shift[axis]-edge)
            if sa>=0:new.append(a)
            if (sa>=0)!=(sb>=0):
                t=sa/(sa-sb)
                new.append((a[0]+t*(b[0]-a[0]),a[1]+t*(b[1]-a[1])))
        poly=new
    return abs(math.fsum(a[0]*b[1]-a[1]*b[0] for a,b in zip(poly,poly[1:]+poly[:1])))/200

def audit(inputs, results):
    rows=list(csv.DictReader((results/"weather-seasonal.v1.csv").open()))
    aggregates=list(csv.DictReader((results/"weather-aggregate.v1.csv").open()))
    numeric_errors=[]
    for r in rows:
        b,d,t,f=map(float,(r["beam_horizontal_wh_m2"],r["diffuse_horizontal_wh_m2"],
                           r["transmitted_beam_wh_m2"],r["sky_isotropic"]))
        assert b>0 and d>=0 and 0<=t<=b
        numeric_errors.append(abs((t+d*f)/(b+d)-float(r["total_transmission"])))
        if r["period"]=="winter":
            assert int(r["records"]) in (2160,2184) and 2006<=int(r["year"])<=2023
    aggregate_errors=[]
    for a in aggregates:
        match=[r for r in rows if all(r[k]==a[k] for k in
               ("source","period","aspect_ratio","height_m","bearing_deg"))]
        assert len(match)==int(a["years"])
        b=math.fsum(float(r["beam_horizontal_wh_m2"]) for r in match)
        d=math.fsum(float(r["diffuse_horizontal_wh_m2"]) for r in match)
        t=math.fsum(float(r["transmitted_beam_wh_m2"]) for r in match)
        f=float(match[0]["sky_isotropic"])
        aggregate_errors.append(abs((t+d*f)/(b+d)-float(a["total_transmission"])))
    comparisons=[]
    for source in ("sarah3","era5"):
        p=inputs/f"pvgis-{source}-32.3-36.2-2005-2023.json.gz"
        receipt=json.loads((inputs/f"pvgis-{source}.receipt.json").read_text())
        assert hashlib.sha256(p.read_bytes()).hexdigest()==receipt["compressed_sha256"]
        raw=gzip.decompress(p.read_bytes())
        assert hashlib.sha256(raw).hexdigest()==receipt["raw_sha256"]
        doc=json.loads(raw)
        hourly=[r for r in doc["outputs"]["hourly"] if r["time"].startswith("2023")]
        t=pd.DatetimeIndex(pd.to_datetime([r["time"] for r in hourly],format="%Y%m%d:%H%M",utc=True))
        altitude=doc["inputs"]["location"]["elevation"]
        offsets=[0.] if source=="sarah3" else [(i+.5)*5-30 for i in range(12)]
        vectors=[]
        for shift in offsets:
            sun=spa_python(t+pd.to_timedelta(shift,unit="min"),32.3,36.2,
                           altitude=altitude,pressure=0,delta_t=None,how="numpy")
            el=np.radians(sun["elevation"].to_numpy())
            az=np.radians(sun["azimuth"].to_numpy())
            vectors.append(list(zip(np.cos(el)*np.sin(az),np.cos(el)*np.cos(az),np.sin(el))))
        b=math.fsum(r["Gb(i)"] for r in hourly);d=math.fsum(r["Gd(i)"] for r in hourly)
        for ratio,bearing in ((1,0),(4,0),(4,90)):
            source_row=next(r for r in rows if r["source"]==source and r["period"]=="annual"
                            and r["year"]=="2023" and r["aspect_ratio"]==str(ratio)
                            and r["height_m"]=="5" and r["bearing_deg"]==str(bearing))
            transmitted=[]
            for i,r in enumerate(hourly):
                if r["Gb(i)"]==0:continue
                if source=="sarah3":
                    f=clipped_fraction(*vectors[0][i],ratio,5,bearing)
                else:
                    numerator=[];denominator=[]
                    for v in vectors:
                        e,n,z=v[i]
                        if z>0:
                            numerator.append(z*clipped_fraction(e,n,z,ratio,5,bearing))
                            denominator.append(z)
                    f=math.fsum(numerator)/math.fsum(denominator) if denominator else 0.
                transmitted.append(r["Gb(i)"]*f)
            value=(math.fsum(transmitted)+d*float(source_row["sky_isotropic"]))/(b+d)
            reference=float(source_row["total_transmission"])
            comparisons.append(dict(source=source,year=2023,aspect_ratio=ratio,height_m=5,
                bearing_deg=bearing,clipping_accumulation=value,published=reference,
                absolute_difference=abs(value-reference)))
    assert max(numeric_errors)<2e-11
    assert max(aggregate_errors)<2e-11
    assert max(c["absolute_difference"] for c in comparisons)<2e-11
    return dict(status="passed",scope="Separately implemented internal audit, not independent researcher replication",
        shared_dependency="pvlib NREL SPA solar position; study code not imported",
        complete_table_checks=len(rows),aggregate_checks=len(aggregates),
        maximum_energy_identity_error=max(numeric_errors),
        maximum_aggregate_recalculation_error=max(aggregate_errors),
        end_to_end_cases=comparisons,tolerance=2e-11)

if __name__=="__main__":
    p=argparse.ArgumentParser(description=__doc__)
    p.add_argument("--inputs",type=Path,required=True);p.add_argument("--results",type=Path,required=True)
    p.add_argument("--output",type=Path,required=True);a=p.parse_args()
    record=audit(a.inputs,a.results)
    a.output.write_text(json.dumps(record,indent=2)+"\n")
    print(json.dumps({"status":"passed","end_to_end_cases":len(record["end_to_end_cases"])}))