import numpy as np
rng=np.random.default_rng(42)
# Eén onafhankelijke modelwaarneming per plant; geen echte ziektegegevens.
y=np.r_[np.zeros(360,dtype=int),np.ones(90,dtype=int)]
deltaT=rng.normal(0.0+2.0*y,1.3,len(y))
vocht=np.clip(rng.normal(55-12*y,10,len(y)),5,95)
X=np.column_stack((deltaT,vocht))
delen=[[],[],[]]
for c in (0,1):
    ids=rng.permutation(np.flatnonzero(y==c))
    for j,a in enumerate(np.split(ids,[int(.6*len(ids)),int(.8*len(ids))])):
        delen[j].extend(a)
tr,va,te=[np.array(a) for a in delen]
gem=X[tr].mean(axis=0); sd=X[tr].std(axis=0)
Z=np.column_stack((np.ones(len(X)),(X-gem)/sd))
def sigmoid(a): return 1/(1+np.exp(-np.clip(a,-40,40)))
b=np.zeros(3)
for stap in range(3000):
    p=sigmoid(Z[tr]@b)
    grad=Z[tr].T@(p-y[tr])/len(tr)+.01*np.r_[0,b[1:]]
    b-=.1*grad
p=sigmoid(Z@b)
rooster=np.arange(.05,.96,.05)
def kosten(ids,q,drempel):
    pred=q>=drempel
    return int(5*np.sum((y[ids]==1)&~pred)+np.sum((y[ids]==0)&pred))
drempel=min(rooster,key=lambda t:kosten(va,p[va],t))
def maat(waar,pred):
    tn=np.sum((waar==0)&(pred==0)); fp=np.sum((waar==0)&(pred==1))
    fn=np.sum((waar==1)&(pred==0)); tp=np.sum((waar==1)&(pred==1))
    return [int(tn),int(fp),int(fn),int(tp)], tp/max(tp+fp,1),tp/max(tp+fn,1)
print("drempel validatie",round(drempel,2),"kosten",kosten(va,p[va],drempel))
for naam,q in (("logistiek0.5",p[te]>=.5),("logistiekgekozen",p[te]>=drempel)):
    cm,pr,re=maat(y[te],q)
    print(naam,"TN FP FN TP",cm,"precisie",round(pr,3),"recall",round(re,3))
# KNN: afstanden in de met training gestandaardiseerde ruimte.
afstand=((Z[te,None,1:]-Z[None,tr,1:])**2).sum(axis=2)
ids=np.argsort(afstand,axis=1)[:,:7]
knn=(y[tr][ids].mean(axis=1)>=.5)
print("KNN7",maat(y[te],knn))
# Kleine echte beslisboom: gewogen Gini, diepte3, minimaal12 per blad.
def gini(a):
    q=a.mean(); return 2*q*(1-q)
def boom(ids,diepte=0):
    blad=float(y[ids].mean())
    if diepte==3 or len(ids)<24 or blad in (0.,1.): return blad
    beste=(len(ids)*gini(y[ids]),None,None,None,None)
    for j in range(2):
        v=np.unique(X[ids,j]); grenzen=(v[:-1]+v[1:])/2
        for t in grenzen:
            links=ids[X[ids,j]<=t]; rechts=ids[X[ids,j]>t]
            if min(len(links),len(rechts))<12: continue
            verlies=len(links)*gini(y[links])+len(rechts)*gini(y[rechts])
            if verlies<beste[0]: beste=(verlies,j,t,links,rechts)
    if beste[1] is None: return blad
    _,j,t,l,r=beste
    return (j,float(t),boom(l,diepte+1),boom(r,diepte+1))
def boomkans(knoop,x):
    if isinstance(knoop,float): return knoop
    j,t,l,r=knoop; return boomkans(l if x[j]<=t else r,x)
tree=boom(tr)
pt=np.array([boomkans(tree,x) for x in X[te]])
print("boom3",maat(y[te],pt>=.5))
print("boomwortel feature/grens",tree[:2])
print("meerderheidsbasislijn accuracy",round(np.mean(y[te]==0),3))
assert len(set(tr)&set(te))==0