import numpy as np
rng=np.random.default_rng(43)
# Onafhankelijke synthetische omstandigheden; geen dag-vooruitprognoses.
G=rng.uniform(0,1000,500)
Tcel=rng.uniform(-5,55,500)+15*G/1000
X=np.column_stack((G,Tcel))
waar=1.40*(G/1000)*(1-.004*(Tcel-25))
y=waar+rng.normal(0,.015,len(G))
ids=rng.permutation(len(G)); tr,va,te=np.split(ids,[300,400])
gem=X[tr].mean(axis=0); sd=X[tr].std(axis=0)
Z=(X-gem)/sd
par=[rng.normal(0,.25,(2,12)),np.zeros(12),
     rng.normal(0,.25,(12,1)),np.zeros(1)]
lam=1e-4
def vooruit(z,p):
    W1,b1,W2,b2=p
    h=np.tanh(z@W1+b1)
    return (h@W2+b2).ravel(),h
def verlies(z,y,p):
    pred,h=vooruit(z,p)
    return np.mean((pred-y)**2)/2+lam*(np.sum(p[0]**2)+np.sum(p[2]**2))/2
def gradient(z,y,p):
    pred,h=vooruit(z,p)
    d2=(pred-y)[:,None]/len(y)
    d1=(d2@p[2].T)*(1-h*h)
    return [z.T@d1+lam*p[0],d1.sum(axis=0),
            h.T@d2+lam*p[2],d2.sum(axis=0)]
# Onafhankelijke centrale-differentiecontrole vóór training.
g=gradient(Z[tr],y[tr],par)[0][0,0]
eps=1e-6; oud=par[0][0,0]
par[0][0,0]=oud+eps; plus=verlies(Z[tr],y[tr],par)
par[0][0,0]=oud-eps; min_=verlies(Z[tr],y[tr],par)
par[0][0,0]=oud
fd=(plus-min_)/(2*eps)
print("gradient analytisch/numeriek",round(g,9),round(fd,9),"verschil",abs(g-fd))
assert abs(g-fd)<1e-7
# Adam: lopende eerste/tweede momenten, met biascorrectie.
m=[np.zeros_like(a) for a in par]; v=[np.zeros_like(a) for a in par]
best=np.inf; historie=[]; beste=None
for epoch in range(1,2001):
    gg=gradient(Z[tr],y[tr],par)
    for j in range(4):
        m[j]=.9*m[j]+.1*gg[j]; v[j]=.999*v[j]+.001*gg[j]**2
        par[j]-=.01*(m[j]/(1-.9**epoch))/(np.sqrt(v[j]/(1-.999**epoch))+1e-8)
    pt=vooruit(Z[tr],par)[0]; pv=vooruit(Z[va],par)[0]
    tv=np.mean((pv-y[va])**2)
    if tv<best: best=tv; beste=[a.copy() for a in par]; beste_epoch=epoch
    if epoch%50==0: historie.append((epoch,np.mean((pt-y[tr])**2),tv))
pred,h=vooruit(Z,beste)
lineair=np.column_stack((np.ones(len(Z)),Z))
coef=np.linalg.lstsq(lineair[tr],y[tr],rcond=None)[0]
rmse=lambda a,b: np.sqrt(np.mean((a-b)**2))
print("gekozen epoch",beste_epoch,"parameters",sum(a.size for a in par))
print("test RMSE netwerk",round(rmse(y[te],pred[te]),5),"kW")
print("test RMSE lineair",round(rmse(y[te],lineair[te]@coef),5),"kW")
print("model-noise RMSE",round(rmse(y[te],waar[te]),5),"kW")
assert np.all(np.isfinite(pred))