airfoil-self-noise

cat /app/train_and_predict.py

cd /app && ls -la data/ && echo "---TRAIN---" && head -5 data/train.csv && echo "..." && wc -l data/train.csv && echo "---VALIDATION---" && head -5 data/validation.csv && wc -l data/validation.csv && echo "---TEST---" && head -5 data/test_features.csv && wc -l data/test_features.csv

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv')
va = pd.read_csv('data/validation.csv')
print('train shape', tr.shape, 'val shape', va.shape)
print()
print('TRAIN describe:')
print(tr.describe().T)
print()
for c in ['frequency','attack_angle','chord_length','free_stream_velocity','suction_side_displacement_thickness','scaled_sound_pressure_level']:
    print(c, 'n_unique train', tr[c].nunique(), 'val', va[c].nunique())
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv')
va = pd.read_csv('data/validation.csv')
alld = pd.concat([tr,va])
for c in ['chord_length','free_stream_velocity','attack_angle']:
    print('===',c,'===')
    v = np.sort(alld[c].unique())
    print('min',v.min(),'max',v.max())
    # histogram of sorted values to detect clusters
    print(np.round(np.sort(alld[c].values),4)[:20])
"

cd /app && python3 -c "
import sklearn; print('sklearn', sklearn.__version__)
import numpy, pandas, scipy; print('numpy', numpy.__version__, 'pandas', pandas.__version__, 'scipy', scipy.__version__)
for m in ['xgboost','lightgbm','catboost','torch','tensorflow']:
    try:
        mod=__import__(m); print(m, getattr(mod,'__version__','?'))
    except Exception as e:
        print(m, 'NOT available')
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv'); va = pd.read_csv('data/validation.csv')
alld = pd.concat([tr,va], ignore_index=True)
# chord canonical
print('chord unique rounded 3:', sorted(alld['chord_length'].round(3).unique()))
print('velocity rounded 0:', sorted(alld['free_stream_velocity'].round(0).unique()))
# attack angle: cluster
aa = np.sort(alld['attack_angle'].unique())
# find gaps
gaps = np.where(np.diff(aa) > 0.05)[0]
clusters=[]
start=0
for g in gaps:
    clusters.append((aa[start], aa[g]))
    start=g+1
clusters.append((aa[start], aa[-1]))
print('n attack clusters', len(clusters))
for c in clusters: print(round(c[0],3), round(c[1],3))
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv'); va = pd.read_csv('data/validation.csv')

chords=np.array([0.0254,0.0508,0.1016,0.1524,0.2286,0.3048])
vels=np.array([31.7,39.6,55.5,71.3])
def snap(x,arr):
    return arr[np.argmin(np.abs(arr-x))]
def cond(df):
    c=df['chord_length'].apply(lambda x: snap(x,chords))
    v=df['free_stream_velocity'].apply(lambda x: snap(x,vels))
    a=df['attack_angle'].round(1)
    return list(zip(a,c,v))
tr_c=set(cond(tr)); va_c=set(cond(va))
print('train conditions:', len(tr_c))
print('val conditions:', len(va_c))
print('overlap:', len(tr_c & va_c))
print('val-only:', len(va_c - tr_c))
# rows per condition
from collections import Counter
cc=Counter(cond(tr))
print('train rows/condition: min',min(cc.values()),'max',max(cc.values()),'mean',round(np.mean(list(cc.values())),1))
print('n train conditions', len(cc))
"

cd /app && cat > /tmp/evalharness.py << 'EOF'
import pandas as pd, numpy as np
chords=np.array([0.0254,0.0508,0.1016,0.1524,0.2286,0.3048])
vels=np.array([31.7,39.6,55.5,71.3])
def snap(x,arr): return arr[np.argmin(np.abs(arr-x))]
def cond_labels(df):
    c=df['chord_length'].apply(lambda x: snap(x,chords))
    v=df['free_stream_velocity'].apply(lambda x: snap(x,vels))
    a=df['attack_angle'].round(1)
    return pd.Series(list(zip(a,c,v)), index=df.index)
def metrics(y_true, y_pred, conds):
    y_true=np.asarray(y_true); y_pred=np.asarray(y_pred)
    err=y_pred-y_true
    rmse=np.sqrt(np.mean(err**2)); mae=np.mean(np.abs(err))
    df=pd.DataFrame({'e':err,'c':list(conds)})
    cr=df.groupby('c')['e'].apply(lambda s: np.sqrt(np.mean(s**2)))
    return dict(rmse=rmse, mae=mae,
               p90=np.percentile(cr,90), p95=np.percentile(cr,95),
               mx=cr.max(), ncond=len(cr))
TARGETS=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
def report(name,m):
    flags=''.join('OK ' if m[k]<=TARGETS[k] else 'XX ' for k in ['rmse','mae','p90','p95','mx'])
    print(f"{name:28s} rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f} | {flags}")
EOF
echo "harness written"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import pandas as pd, numpy as np
from evalharness import cond_labels, metrics, report, TARGETS
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
import train_and_predict as tp

tr=pd.read_csv('data/train.csv'); va=pd.read_csv('data/validation.csv')
model=make_pipeline(tp.AirfoilFeatures(), StandardScaler(), Ridge(alpha=10.0))
model.fit(tr[tp.FEATURES], tr[tp.TARGET])
pred=model.predict(va[tp.FEATURES])
m=metrics(va[tp.TARGET], pred, cond_labels(va))
report('starter Ridge (train->val)', m)
print('targets:', TARGETS)
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr=pd.read_csv('data/train.csv')
# correlation of target with raw and log features
import numpy as np
d=tr.copy()
d['logf']=np.log10(d['frequency'])
d['logt']=np.log10(d['suction_side_displacement_thickness'])
d['strouhal']=np.log10(d['frequency']*d['suction_side_displacement_thickness']/d['free_stream_velocity'])
d['logchord']=np.log10(d['chord_length'])
d['reynolds']=np.log10(d['free_stream_velocity']*d['chord_length'])
cols=['frequency','logf','attack_angle','chord_length','logchord','free_stream_velocity','suction_side_displacement_thickness','logt','strouhal','reynolds']
print(d[cols+['scaled_sound_pressure_level']].corr()['scaled_sound_pressure_level'].round(3))
"

cd /app && cat > /tmp/bench.py << 'EOF'
import sys; sys.path.insert(0,'/tmp'); sys.path.insert(0,'/app')
import pandas as pd, numpy as np
from evalharness import cond_labels, metrics, report
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import Ridge
from sklearn.ensemble import RandomForestRegressor, ExtraTreesRegressor, HistGradientBoostingRegressor, GradientBoostingRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.svm import SVR
from sklearn.neural_network import MLPRegressor
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel, Matern

RAW=["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]

def feats(df):
    f=df['frequency'].astype(float); a=df['attack_angle'].astype(float)
    c=df['chord_length'].astype(float); v=df['free_stream_velocity'].astype(float)
    t=df['suction_side_displacement_thickness'].astype(float)
    o=pd.DataFrame(index=df.index)
    o['logf']=np.log10(f); o['a']=a; o['c']=c; o['v']=v
    o['logt']=np.log10(t)
    o['strouhal']=np.log10(f*t/v)
    o['logchord']=np.log10(c)
    o['reynolds']=np.log10(v*c)
    o['logf2']=o['logf']**2
    o['strouhal2']=o['strouhal']**2
    o['a2']=a**2
    o['logf_v']=o['logf']*v
    o['logf_a']=o['logf']*a
    o['logf_t']=o['logf']*o['logt']
    o['logf_c']=o['logf']*o['logchord']
    return o.to_numpy(float)

tr=pd.read_csv('data/train.csv'); va=pd.read_csv('data/validation.csv')
Xtr=feats(tr); ytr=tr['scaled_sound_pressure_level'].values
Xva=feats(va); yva=va['scaled_sound_pressure_level'].values
cva=cond_labels(va)

def run(name, model):
    model.fit(Xtr,ytr)
    p=model.predict(Xva)
    report(name, metrics(yva,p,cva))

run('Ridge(a=1)+feats', make_pipeline(StandardScaler(), Ridge(alpha=1.0)))
run('RF(400)', RandomForestRegressor(n_estimators=400, n_jobs=-1, random_state=0))
run('ExtraTrees(400)', ExtraTreesRegressor(n_estimators=400, n_jobs=-1, random_state=0))
run('HGB', HistGradientBoostingRegressor(random_state=0))
run('HGB(lr.05,500)', HistGradientBoostingRegressor(learning_rate=0.05, max_iter=500, random_state=0))
run('GBR', GradientBoostingRegressor(random_state=0))
run('KNN(7)', make_pipeline(StandardScaler(), KNeighborsRegressor(n_neighbors=7, weights='distance')))
run('SVR(rbf)', make_pipeline(StandardScaler(), SVR(C=10, gamma='scale', epsilon=0.1)))
run('MLP', make_pipeline(StandardScaler(), MLPRegressor(hidden_layer_sizes=(128,128,64), max_iter=2000, random_state=0, early_stopping=True)))
k=ConstantKernel()*Matern(length_scale=np.ones(Xtr.shape[1]), nu=1.5)+WhiteKernel()
run('GPR(Matern)', make_pipeline(StandardScaler(), GaussianProcessRegressor(kernel=k, normalize_y=True, n_restarts_optimizer=0, random_state=0)))
EOF
python3 /tmp/bench.py

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
exec(open('/tmp/bench.py').read().split('def run')[0])  # load feats, data
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
k=ConstantKernel()*Matern(length_scale=np.ones(Xtr.shape[1]), nu=1.5)+WhiteKernel()
gp=make_pipeline(StandardScaler(), GaussianProcessRegressor(kernel=k, normalize_y=True, random_state=0))
gp.fit(Xtr,ytr); p=gp.predict(Xva)
err=p-yva
df=pd.DataFrame({'e':err,'c':list(cva)})
cr=df.groupby('c')['e'].apply(lambda s: np.sqrt(np.mean(s**2))).sort_values(ascending=False)
print('per-condition RMSE (GPR), worst first:')
print(cr.round(3))
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import pandas as pd, numpy as np
from evalharness import cond_labels
tr=pd.read_csv('data/train.csv')
tr['c']=cond_labels(tr)
# pick a few conditions, sort by frequency, look at SPL smoothness
for cond in list(tr['c'].unique())[:3]:
    sub=tr[tr['c']==cond].sort_values('frequency')
    print('cond',cond,'n=',len(sub))
    print('  freq:', np.round(sub['frequency'].values,0))
    print('  SPL :', np.round(sub['scaled_sound_pressure_level'].values,2))
    # fit smooth spline-ish: second difference magnitude
    lf=np.log10(sub['frequency'].values); y=sub['scaled_sound_pressure_level'].values
    if len(sub)>4:
        # residual from quadratic fit in logf
        co=np.polyfit(lf,y,3); res=y-np.polyval(co,lf)
        print('  cubic-fit resid std:', round(res.std(),3))
"

cd /app && cat > /tmp/cv.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
import pandas as pd, numpy as np
from evalharness import cond_labels, metrics
from sklearn.model_selection import GroupKFold

tr=pd.read_csv('/app/data/train.csv'); va=pd.read_csv('/app/data/validation.csv')
FULL=pd.concat([tr,va], ignore_index=True)
FULL['cond']=cond_labels(FULL)
groups=FULL['cond'].astype(str).values
y=FULL['scaled_sound_pressure_level'].values

def cv_eval(make_model, feat_fn, n_splits=10, seed=0):
    gkf=GroupKFold(n_splits=n_splits)
    oof=np.zeros(len(FULL))
    for tri,tei in gkf.split(FULL, y, groups):
        Xtr=feat_fn(FULL.iloc[tri]); Xte=feat_fn(FULL.iloc[tei])
        m=make_model(); m.fit(Xtr, y[tri])
        oof[tei]=m.predict(Xte)
    return metrics(y, oof, FULL['cond'])

def report(name,m):
    T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
    flags=''.join('OK ' if m[k]<=T[k] else 'XX ' for k in ['rmse','mae','p90','p95','mx'])
    print(f"{name:30s} rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f} nc={m['ncond']} | {flags}")
EOF
echo ok

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import cv_eval, report
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern, RBF

def feats_full(df):
    f=df['frequency'].astype(float).values; a=df['attack_angle'].astype(float).values
    c=df['chord_length'].astype(float).values; v=df['free_stream_velocity'].astype(float).values
    t=df['suction_side_displacement_thickness'].astype(float).values
    o=np.column_stack([np.log10(f), a, np.log10(c), v, np.log10(t), np.log10(f*t/v), np.log10(v*c)])
    return o
def feats_compact(df):
    f=df['frequency'].astype(float).values; a=df['attack_angle'].astype(float).values
    c=df['chord_length'].astype(float).values; v=df['free_stream_velocity'].astype(float).values
    t=df['suction_side_displacement_thickness'].astype(float).values
    o=np.column_stack([np.log10(f), a, np.log10(c), v, np.log10(t)])
    return o

def mk(nu):
    def f():
        d=5
        k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d), length_scale_bounds=(1e-2,1e3), nu=nu)+WhiteKernel(0.5,(1e-3,1e2))
        return make_pipeline(StandardScaler(), GaussianProcessRegressor(kernel=k, normalize_y=True, n_restarts_optimizer=2, random_state=0))
    return f
report('GPR Matern1.5 compact', cv_eval(mk(1.5), feats_compact, n_splits=10))
report('GPR Matern2.5 compact', cv_eval(mk(2.5), feats_compact, n_splits=10))
" 2>&1 | grep -v Warning

restart: true

cd /app && timeout 120 python3 -c "
import sys, time; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import cv_eval, report, FULL
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern

def feats_compact(df):
    f=df['frequency'].astype(float).values; a=df['attack_angle'].astype(float).values
    c=df['chord_length'].astype(float).values; v=df['free_stream_velocity'].astype(float).values
    t=df['suction_side_displacement_thickness'].astype(float).values
    return np.column_stack([np.log10(f), a, np.log10(c), v, np.log10(t)])

# Fit GPR once with optimizer on train->val to get good hyperparams, time it
t0=time.time()
def mk():
    d=5
    k=ConstantKernel(1.0)*Matern(length_scale=np.ones(d), nu=1.5)+WhiteKernel(0.3)
    return make_pipeline(StandardScaler(), GaussianProcessRegressor(kernel=k, normalize_y=True, n_restarts_optimizer=0, random_state=0))
report('GPR M1.5 compact 5fold', cv_eval(mk, feats_compact, n_splits=5))
print('time', round(time.time()-t0,1))
" 2>&1 | grep -vi warn

cd /app && cat > /tmp/ff.py << 'EOF'
import numpy as np, pandas as pd
def feats(df):
    f=df['frequency'].astype(float).values; a=df['attack_angle'].astype(float).values
    c=df['chord_length'].astype(float).values; v=df['free_stream_velocity'].astype(float).values
    t=df['suction_side_displacement_thickness'].astype(float).values
    logf=np.log10(f); logt=np.log10(t); logc=np.log10(c)
    st=np.log10(f*t/v); re=np.log10(v*c)
    o=np.column_stack([logf,a,logc,v,logt,st,re,logf*logf,st*st,a*a,logf*v,logf*a,logf*logt,a*logc,logf*logc])
    return o
EOF
cd /app && timeout 200 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from cv import cv_eval, report
from ff import feats
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor, ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.svm import SVR
report('ExtraTrees400', cv_eval(lambda: ExtraTreesRegressor(n_estimators=400,n_jobs=-1,random_state=0), feats, 8))
report('RF400', cv_eval(lambda: RandomForestRegressor(n_estimators=400,n_jobs=-1,random_state=0), feats, 8))
report('HGB', cv_eval(lambda: HistGradientBoostingRegressor(learning_rate=0.05,max_iter=600,max_leaf_nodes=31,l2_regularization=1.0,random_state=0), feats, 8))
report('KNNdist', cv_eval(lambda: make_pipeline(StandardScaler(), KNeighborsRegressor(n_neighbors=10,weights='distance')), feats, 8))
report('SVRrbf', cv_eval(lambda: make_pipeline(StandardScaler(), SVR(C=20,gamma='scale',epsilon=0.05)), feats, 8))
" 2>&1 | grep -vi warn

cd /app && timeout 60 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import pandas as pd, numpy as np
from evalharness import cond_labels
from cv import FULL
# within-condition variation of thickness relative to between-condition
g=FULL.groupby('cond')['suction_side_displacement_thickness']
within=g.std()/g.mean()
print('within-condition CV of thickness: median', round(within.median(),4), 'max', round(within.max(),4))
# check frequency coverage per condition
gf=FULL.groupby('cond')['frequency']
print('freqs per cond: min count', FULL.groupby('cond').size().min(), 'max', FULL.groupby('cond').size().max())
print('thickness is essentially constant within condition:', within.max()<0.02)
"

cd /app && cat > /tmp/oofgen.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd, time
from cv import FULL
from ff import feats
from evalharness import metrics
from sklearn.model_selection import GroupKFold
y=FULL['scaled_sound_pressure_level'].values
groups=FULL['cond'].astype(str).values
conds=FULL['cond']
NF=8
gkf=GroupKFold(n_splits=NF)
FOLDS=list(gkf.split(FULL,y,groups))
np.save('/tmp/y.npy', y)
pd.Series(conds.values).to_pickle('/tmp/conds.pkl')

def gen(name, make_model, feat_fn=feats):
    t0=time.time(); oof=np.zeros(len(FULL))
    for tri,tei in FOLDS:
        Xtr=feat_fn(FULL.iloc[tri]); Xte=feat_fn(FULL.iloc[tei])
        m=make_model(); m.fit(Xtr,y[tri]); oof[tei]=m.predict(Xte)
    np.save(f'/tmp/oof_{name}.npy', oof)
    mm=metrics(y,oof,conds)
    print(f"{name:16s} rmse={mm['rmse']:.3f} mae={mm['mae']:.3f} p90={mm['p90']:.3f} p95={mm['p95']:.3f} max={mm['mx']:.3f} t={time.time()-t0:.0f}s")
    return oof
EOF
cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from oofgen import gen
from ff import feats
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor, RandomForestRegressor
from sklearn.svm import SVR
gen('et', lambda: ExtraTreesRegressor(n_estimators=600,n_jobs=-1,random_state=0))
gen('hgb', lambda: HistGradientBoostingRegressor(learning_rate=0.05,max_iter=600,max_leaf_nodes=31,l2_regularization=1.0,random_state=0))
gen('rf', lambda: RandomForestRegressor(n_estimators=600,n_jobs=-1,random_state=0,min_samples_leaf=1))
gen('svr', lambda: make_pipeline(StandardScaler(), SVR(C=20,gamma='scale',epsilon=0.05)))
" 2>&1 | grep -vi warn

cd /app && cat >> /tmp/ff.py << 'EOF'

def feats_gpr(df):
    f=df['frequency'].astype(float).values; a=df['attack_angle'].astype(float).values
    c=df['chord_length'].astype(float).values; v=df['free_stream_velocity'].astype(float).values
    t=df['suction_side_displacement_thickness'].astype(float).values
    return np.column_stack([np.log10(f), a, np.log10(c), v, np.log10(t), np.log10(f*t/v)])
EOF
cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from oofgen import gen
from ff import feats_gpr
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
def mk():
    d=6
    k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d),length_scale_bounds=(1e-2,1e2),nu=1.5)+WhiteKernel(0.3,(1e-3,1e1))
    return make_pipeline(StandardScaler(), GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=1,random_state=0))
gen('gpr', mk, feats_gpr)
" 2>&1 | grep -vi warn

cd /app && timeout 120 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np
from cv import FULL
from ff import feats_gpr
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
y=FULL['scaled_sound_pressure_level'].values
X=feats_gpr(FULL); sc=StandardScaler().fit(X); Xs=sc.transform(X)
d=X.shape[1]
k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d),length_scale_bounds=(1e-2,1e2),nu=1.5)+WhiteKernel(0.3,(1e-3,1e1))
t0=time.time()
gp=GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=2,random_state=0).fit(Xs,y)
print('fit time', round(time.time()-t0,1))
print('learned kernel:', gp.kernel_)
" 2>&1 | grep -vi warn

cd /app && timeout 120 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np
from oofgen import gen
from ff import feats_gpr
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
def mk():
    d=6
    k=ConstantKernel(1.23**2)*Matern(length_scale=[1.86,2.22,1.26,4.6,1.97,1.51],nu=1.5)+WhiteKernel(0.0028)
    return make_pipeline(StandardScaler(), GaussianProcessRegressor(kernel=k,normalize_y=True,optimizer=None,random_state=0))
gen('gpr', mk, feats_gpr)
" 2>&1 | grep -vi warn

cd /app && cat > /tmp/ens.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from evalharness import metrics
y=np.load('/tmp/y.npy'); conds=pd.read_pickle('/tmp/conds.pkl')
names=['et','hgb','rf','svr','gpr']
OOF={n:np.load(f'/tmp/oof_{n}.npy') for n in names}
def ev(pred):
    return metrics(y,pred,conds)
def show(name,pred):
    m=ev(pred); T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
    flags=''.join('OK ' if m[k]<=T[k] else 'XX ' for k in ['rmse','mae','p90','p95','mx'])
    print(f"{name:26s} rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f} | {flags}")
    return m
EOF
cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
from ens import OOF, show
import numpy as np
show('gpr', OOF['gpr'])
show('et', OOF['et'])
show('gpr+et', 0.5*OOF['gpr']+0.5*OOF['et'])
show('gpr+et+hgb', (OOF['gpr']+OOF['et']+OOF['hgb'])/3)
show('2gpr+et', (2*OOF['gpr']+OOF['et'])/3)
show('gpr+et+hgb+svr', (OOF['gpr']+OOF['et']+OOF['hgb']+OOF['svr'])/4)
show('all5', sum(OOF.values())/5)
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import OOF, y, conds
from cv import FULL
pred=0.5*OOF['gpr']+0.5*OOF['et']
err=pred-y
df=pd.DataFrame({'e':err,'c':conds})
cr=df.groupby('c').apply(lambda s: pd.Series({'rmse':np.sqrt(np.mean(s['e']**2)),'bias':s['e'].mean(),'n':len(s)}))
cr=cr.sort_values('rmse',ascending=False)
print('WORST 12 conditions (cond = angle,chord,vel):')
print(cr.head(12).round(2))
print()
print('p90 threshold count: conditions with rmse>2.35:', (cr['rmse']>2.35).sum(), 'of', len(cr))
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd, itertools
from ens import OOF, y, conds
from evalharness import metrics
names=['et','hgb','rf','svr','gpr']
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
best=[]
grid=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,1.0]
from itertools import product
for w in product(grid, repeat=5):
    s=sum(w)
    if s<0.99: continue
    if (np.array(w)>0).sum()>3: continue  # at most 3 models for simplicity
    wn=np.array(w)/s
    pred=sum(wn[i]*OOF[names[i]] for i in range(5))
    m=metrics(y,pred,conds)
    # primary: minimize p90, tie-break rmse
    best.append((m['p90'], m['rmse'], m['mae'], m['p95'], m['mx'], w))
best.sort()
print('Top configs by p90 (et,hgb,rf,svr,gpr weights):')
for p90,rmse,mae,p95,mx,w in best[:12]:
    print(f'  p90={p90:.3f} rmse={rmse:.3f} mae={mae:.3f} p95={p95:.3f} max={mx:.3f}  w={w}')
"

cd /app && timeout 60 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
d=FULL.copy()
d['st']=np.log10(d['frequency']*d['suction_side_displacement_thickness']/d['free_stream_velocity'])
# For each condition compute peak SPL and peak-St; see if shape collapses
from evalharness import cond_labels
# correlation of SPL with st overall
print('corr SPL vs log-St:', np.corrcoef(d['st'], d['scaled_sound_pressure_level'])[0,1].round(3))
# Fit global cubic in st, check residual
for deg in [2,3,4,5]:
    co=np.polyfit(d['st'], d['scaled_sound_pressure_level'], deg)
    r=d['scaled_sound_pressure_level']-np.polyval(co,d['st'])
    print(f'global poly deg{deg} in st: resid std={r.std():.3f}')
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import OOF, y, conds
from evalharness import metrics
pred=0.5*OOF['gpr']+0.5*OOF['et']
df=pd.DataFrame({'p':pred,'y':y,'c':conds})
# oracle level correction: remove per-condition mean error
df['e']=df['p']-df['y']
bias=df.groupby('c')['e'].transform('mean')
pred_corr=df['p']-bias
m0=metrics(y,pred,conds); m1=metrics(y,pred_corr.values,conds)
print('original   p90=%.3f p95=%.3f max=%.3f rmse=%.3f'%(m0['p90'],m0['p95'],m0['mx'],m0['rmse']))
print('oracle-lvl p90=%.3f p95=%.3f max=%.3f rmse=%.3f'%(m1['p90'],m1['p95'],m1['mx'],m1['rmse']))
# distribution of per-condition bias magnitude
cb=df.groupby('c')['e'].mean()
print('per-condition |bias|: median %.2f p90 %.2f max %.2f'%(cb.abs().median(), cb.abs().quantile(.9), cb.abs().max()))
# shape-only RMSE per condition (after removing bias)
shape=df.assign(ec=df['e']-bias).groupby('c')['ec'].apply(lambda s: np.sqrt(np.mean(s**2)))
print('shape-only per-cond RMSE: median %.2f p90 %.2f max %.2f'%(shape.median(),shape.quantile(.9),shape.max()))
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import OOF, y, conds
from cv import FULL
pred=0.5*OOF['gpr']+0.5*OOF['et']
d=FULL.copy(); d['p']=pred; d['y']=y
for cond in [(12.6,0.1524,39.6),(8.9,0.1016,71.3),(19.7,0.0508,71.3)]:
    sub=d[d['cond']==cond].sort_values('frequency')
    print('cond',cond)
    print('  f   ', np.round(sub['frequency'].values,0).astype(int))
    print('  true', np.round(sub['y'].values,1))
    print('  pred', np.round(sub['p'].values,1))
"

cd /app && cat >> /tmp/ff.py << 'EOF'

def feats_nn(df):
    f=df['frequency'].astype(float).values; a=df['attack_angle'].astype(float).values
    c=df['chord_length'].astype(float).values; v=df['free_stream_velocity'].astype(float).values
    t=df['suction_side_displacement_thickness'].astype(float).values
    logf=np.log10(f); logt=np.log10(t); logc=np.log10(c); st=np.log10(f*t/v); re=np.log10(v*c)
    return np.column_stack([logf,a,logc,v,logt,st,re])
EOF
cd /app && timeout 250 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np
from oofgen import gen
from ff import feats_nn
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.neural_network import MLPRegressor
from sklearn.base import BaseEstimator, RegressorMixin
class BagMLP(BaseEstimator, RegressorMixin):
    def __init__(self, arch=(128,128,64), n=5, alpha=1e-3):
        self.arch=arch; self.n=n; self.alpha=alpha
    def fit(self,X,y):
        self.ms=[]
        for s in range(self.n):
            m=make_pipeline(StandardScaler(), MLPRegressor(hidden_layer_sizes=self.arch,alpha=self.alpha,
                max_iter=3000,random_state=s,early_stopping=True,n_iter_no_change=30,validation_fraction=0.12))
            m.fit(X,y); self.ms.append(m)
        return self
    def predict(self,X):
        return np.mean([m.predict(X) for m in self.ms],axis=0)
gen('mlp', lambda: BagMLP((128,128,64),5,1e-3), feats_nn)
" 2>&1 | grep -vi warn

cd /app && cat > /tmp/gpropt.py << 'EOF'
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
y=FULL['scaled_sound_pressure_level'].values
groups=FULL['cond'].astype(str).values
conds=FULL['cond']
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,groups))

def run(name, feat_fn, nu=1.5):
    X=feat_fn(FULL); d=X.shape[1]
    sc=StandardScaler().fit(X)
    k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d),length_scale_bounds=(5e-2,1e2),nu=nu)+WhiteKernel(0.1,(1e-3,1e1))
    t0=time.time()
    gp=GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=1,random_state=0).fit(sc.transform(X),y)
    theta=gp.kernel_
    # frozen CV
    oof=np.zeros(len(FULL))
    for tri,tei in FOLDS:
        sci=StandardScaler().fit(X[tri])
        g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
        oof[tei]=g.predict(sci.transform(X[tei]))
    m=metrics(y,oof,conds)
    print(f"{name:20s} rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f} t={time.time()-t0:.0f}s")
    np.save(f'/tmp/oof_{name}.npy', oof)
    return theta
EOF
echo ok

cd /app && timeout 280 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from gpropt import run
def f8(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(v*c),a*st])
def f9(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(v*c),a*st,a*logf])
t1=run('gpr8', f8, 1.5)
print('  kernel8:', t1)
t2=run('gpr9', f9, 1.5)
print('  kernel9:', t2)
" 2>&1 | grep -vi warn

cd /app && timeout 200 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
y=FULL['scaled_sound_pressure_level'].values
groups=FULL['cond'].astype(str).values; conds=FULL['cond']
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,groups))
def f8(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(v*c),a*st])
X=f8(FULL); d=X.shape[1]
# optimize on subsample
rng=np.random.RandomState(0); idx=rng.choice(len(X),600,replace=False)
sc=StandardScaler().fit(X[idx])
k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d),length_scale_bounds=(5e-2,1e2),nu=1.5)+WhiteKernel(0.1,(1e-3,1e1))
t0=time.time()
gp=GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=2,random_state=0).fit(sc.transform(X[idx]),y[idx])
print('opt time',round(time.time()-t0,1),'kernel',gp.kernel_)
theta=gp.kernel_
oof=np.zeros(len(FULL))
for tri,tei in FOLDS:
    sci=StandardScaler().fit(X[tri])
    g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
    oof[tei]=g.predict(sci.transform(X[tei]))
m=metrics(y,oof,conds)
print(f\"gpr8 rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f}\")
np.save('/tmp/oof_gpr8.npy', oof)
" 2>&1 | grep -vi warn

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from ens import y, conds, show
OOF={n:np.load(f'/tmp/oof_{n}.npy') for n in ['et','hgb','rf','svr','gpr','gpr8']}
show('gpr8', OOF['gpr8'])
show('gpr8+et', 0.5*OOF['gpr8']+0.5*OOF['et'])
show('gpr8+et+hgb', (OOF['gpr8']+OOF['et']+OOF['hgb'])/3)
show('gpr+gpr8+et', (OOF['gpr']+OOF['gpr8']+OOF['et'])/3)
show('gpr8+2et', (OOF['gpr8']+2*OOF['et'])/3)
show('.4gpr8+.4et+.2hgb', .4*OOF['gpr8']+.4*OOF['et']+.2*OOF['hgb'])
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds, show
from cv import FULL
OOF={n:np.load(f'/tmp/oof_{n}.npy') for n in ['et','gpr8']}
base=0.5*OOF['gpr8']+0.5*OOF['et']
logf=np.log10(FULL['frequency'].values)
def smooth(pred, deg, blend):
    out=pred.copy()
    for c in conds.unique():
        m=(conds==c).values
        if m.sum()<deg+2: continue
        x=logf[m]; p=pred[m]
        co=np.polyfit(x,p,deg); fit=np.polyval(co,x)
        out[m]=blend*fit+(1-blend)*p
    return out
show('base gpr8+et', base)
for deg in [3,4,5]:
    for bl in [0.5,1.0]:
        show(f'smooth deg{deg} bl{bl}', smooth(base,deg,bl))
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from evalharness import metrics
models={n:np.load(f'/tmp/oof_{n}.npy') for n in ['et','hgb','rf','svr','gpr','gpr8']}
models['gpr8+et']=0.5*models['gpr8']+0.5*models['et']
def shape_metrics(pred):
    df=pd.DataFrame({'e':pred-y,'c':conds})
    bias=df.groupby('c')['e'].transform('mean')
    sh=df['e']-bias
    cr=pd.DataFrame({'sh':sh,'c':conds}).groupby('c')['sh'].apply(lambda s:np.sqrt(np.mean(s**2)))
    return cr.median(), cr.quantile(.9), cr.max()
print('model            shape_med shape_p90 shape_max')
for n,p in models.items():
    med,p90,mx=shape_metrics(p)
    print(f'{n:12s}      {med:.3f}     {p90:.3f}     {mx:.3f}')
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from cv import FULL
pred=0.5*np.load('/tmp/oof_gpr8.npy')+0.5*np.load('/tmp/oof_et.npy')
df=pd.DataFrame({'e':pred-y,'c':conds,'a':FULL['attack_angle'].round(1),'ch':FULL['chord_length'].round(3),'v':FULL['free_stream_velocity'].round(0)})
cr=df.groupby('c').apply(lambda s: pd.Series({'rmse':np.sqrt(np.mean(s['e']**2)),'bias':s['e'].mean(),'a':s['a'].iloc[0],'ch':s['ch'].iloc[0],'v':s['v'].iloc[0],'n':len(s)}), include_groups=False).sort_values('rmse',ascending=False)
print('Worst 14 by full RMSE:')
print(cr.head(14)[['rmse','bias','a','ch','v','n']].round(2).to_string())
print()
print('angle distribution of worst 14:', sorted(cr.head(14)['a'].tolist()))
print('all conditions: n with a>=12:', (cr['a']>=12).sum(),'of',len(cr))
"

cd /app && timeout 120 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from oofgen import gen
from ff import feats
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import Ridge
from sklearn.kernel_ridge import KernelRidge
gen('pr3', lambda: make_pipeline(StandardScaler(), PolynomialFeatures(3,include_bias=False), Ridge(alpha=5.0)), feats)
gen('pr2', lambda: make_pipeline(StandardScaler(), PolynomialFeatures(2,include_bias=False), Ridge(alpha=2.0)), feats)
gen('krr', lambda: make_pipeline(StandardScaler(), KernelRidge(alpha=0.3, kernel='rbf', gamma=0.05)), feats)
" 2>&1 | grep -vi warn

cd /app && python3 -c "
import pandas as pd, numpy as np
tr=pd.read_csv('data/train.csv'); va=pd.read_csv('data/validation.csv')
full=pd.concat([tr,va],ignore_index=True)
# exact unique tuples
for cols in [['attack_angle','chord_length','free_stream_velocity']]:
    print('exact unique (a,c,v) tuples in full:', full.groupby(cols).ngroups)
print('exact unique a:', full['attack_angle'].nunique(),'c:', full['chord_length'].nunique(),'v:', full['free_stream_velocity'].nunique())
# so are a,c,v jittered (all unique) or repeated?
print()
print('sample raw chord values (first 6 unique):', sorted(full['chord_length'].unique())[:6])
print('counts of top chord value:', (full['chord_length']==full['chord_length'].mode()[0]).sum())
"

cd /app && cat > /tmp/gpr2.py << 'EOF'
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
y=FULL['scaled_sound_pressure_level'].values
groups=FULL['cond'].astype(str).values; conds=FULL['cond']
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,groups))
def f8(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(v*c),a*st])
def build(name, kernel, nsub=600):
    X=f8(FULL)
    rng=np.random.RandomState(0); idx=rng.choice(len(X),nsub,replace=False)
    sc=StandardScaler().fit(X[idx])
    t0=time.time()
    gp=GaussianProcessRegressor(kernel=kernel,normalize_y=True,n_restarts_optimizer=2,random_state=0).fit(sc.transform(X[idx]),y[idx])
    theta=gp.kernel_
    oof=np.zeros(len(FULL))
    for tri,tei in FOLDS:
        sci=StandardScaler().fit(X[tri])
        g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
        oof[tei]=g.predict(sci.transform(X[tei]))
    m=metrics(y,oof,conds)
    print(f"{name:10s} rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f} t={time.time()-t0:.0f}s")
    print('   kernel:', theta)
    np.save(f'/tmp/oof_{name}.npy', oof)
EOF
cd /app && timeout 200 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from gpr2 import build
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern, RBF
d=8
# two-scale sum kernel
k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d)*3,length_scale_bounds=(5e-1,1e2),nu=2.5) \
 + ConstantKernel(0.5,(1e-3,1e2))*Matern(length_scale=np.ones(d)*0.8,length_scale_bounds=(1e-1,5e0),nu=1.5) \
 + WhiteKernel(0.05,(1e-3,1e0))
build('gprsum', k)
" 2>&1 | grep -vi warn

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from ens import y, conds, show
O={n:np.load(f'/tmp/oof_{n}.npy') for n in ['et','hgb','gpr8','gprsum','svr']}
show('gprsum', O['gprsum'])
show('gprsum+et', 0.5*O['gprsum']+0.5*O['et'])
show('.6gprsum+.4et', .6*O['gprsum']+.4*O['et'])
show('.7gprsum+.3et', .7*O['gprsum']+.3*O['et'])
show('gprsum+et+hgb', (O['gprsum']+O['et']+O['hgb'])/3)
show('gprsum+gpr8+et', (O['gprsum']+O['gpr8']+O['et'])/3)
show('.5gprsum+.3et+.2gpr8', .5*O['gprsum']+.3*O['et']+.2*O['gpr8'])
"

cd /app && cat > /tmp/valproxy.py << 'EOF'
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from evalharness import cond_labels, metrics, report
from ff import feats
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
from sklearn.ensemble import ExtraTreesRegressor
tr=pd.read_csv('/app/data/train.csv'); va=pd.read_csv('/app/data/validation.csv')
ytr=tr['scaled_sound_pressure_level'].values; yva=va['scaled_sound_pressure_level'].values
cva=cond_labels(va)
def f8(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(v*c),a*st])
def make_gprsum():
    d=8
    return (ConstantKernel(1.44**2)*Matern(length_scale=[2.21,1.44,2.51,100,2.95,2.11,4.29,2.35],nu=2.5)
            + ConstantKernel(0.299**2)*Matern(length_scale=[0.434,1.08,0.282,5,5,5,0.534,2.5],nu=1.5)
            + WhiteKernel(0.001))
Xtr=f8(tr); Xva=f8(va); sc=StandardScaler().fit(Xtr)
gp=GaussianProcessRegressor(kernel=make_gprsum(),normalize_y=True,optimizer=None).fit(sc.transform(Xtr),ytr)
pg=gp.predict(sc.transform(Xva))
et=ExtraTreesRegressor(n_estimators=600,n_jobs=-1,random_state=0).fit(feats(tr),ytr)
pe=et.predict(feats(va))
report('val: gprsum', metrics(yva,pg,cva))
report('val: et', metrics(yva,pe,cva))
report('val: gprsum+et', metrics(yva,0.5*pg+0.5*pe,cva))
report('val: .6g+.4e', metrics(yva,0.6*pg+0.4*pe,cva))
EOF
cd /app && timeout 60 python3 /tmp/valproxy.py 2>&1 | grep -vi warn

cd /app && timeout 60 python3 -c "
import sys; sys.path.insert(0,'/tmp')
exec(open('/tmp/valproxy.py').read().split('report(')[0])
import numpy as np, pandas as pd
pred=0.5*pg+0.5*pe
df=pd.DataFrame({'e':pred-yva,'c':cva,'a':va['attack_angle'].round(1),'ch':va['chord_length'].round(3),'v':va['free_stream_velocity'].round(0)})
cr=df.groupby('c').apply(lambda s: pd.Series({'rmse':np.sqrt(np.mean(s['e']**2)),'bias':s['e'].mean(),'a':s['a'].iloc[0],'ch':s['ch'].iloc[0],'v':s['v'].iloc[0],'n':len(s)}),include_groups=False).sort_values('rmse',ascending=False)
print(cr[['rmse','bias','a','ch','v','n']].round(2).to_string())
" 2>&1 | grep -vi warn

cd /app && cat > /tmp/gpr3.py << 'EOF'
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
y=FULL['scaled_sound_pressure_level'].values
groups=FULL['cond'].astype(str).values; conds=FULL['cond']
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,groups))
def f9(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v); stc=np.log10(f*c/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,stc,np.log10(v*c),a*st])
def build(name, kernel, feat, nsub=800):
    X=feat(FULL)
    rng=np.random.RandomState(1); idx=rng.choice(len(X),nsub,replace=False)
    sc=StandardScaler().fit(X[idx])
    t0=time.time()
    gp=GaussianProcessRegressor(kernel=kernel,normalize_y=True,n_restarts_optimizer=1,random_state=0).fit(sc.transform(X[idx]),y[idx])
    theta=gp.kernel_
    oof=np.zeros(len(FULL))
    for tri,tei in FOLDS:
        sci=StandardScaler().fit(X[tri])
        g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
        oof[tei]=g.predict(sci.transform(X[tei]))
    m=metrics(y,oof,conds)
    print(f"{name:10s} rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f} t={time.time()-t0:.0f}s")
    print('   ',theta)
    np.save(f'/tmp/oof_{name}.npy', oof)
EOF
cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from gpr3 import build, f9
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
d=9
k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d)*3,length_scale_bounds=(5e-1,1e2),nu=2.5) \
 + ConstantKernel(0.5,(1e-3,1e2))*Matern(length_scale=np.ones(d)*0.8,length_scale_bounds=(8e-2,5e0),nu=1.5) \
 + WhiteKernel(0.05,(1e-3,1e0))
build('gpr9s', k, f9)
" 2>&1 | grep -vi warn

cd /app && timeout 250 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np
from gpr3 import f9, FULL, FOLDS, y, conds
from evalharness import metrics
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
d=9
k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d)*3,length_scale_bounds=(5e-1,1e2),nu=2.5) \
 + ConstantKernel(0.5,(1e-3,1e2))*Matern(length_scale=np.ones(d)*0.8,length_scale_bounds=(8e-2,5e0),nu=1.5) \
 + WhiteKernel(0.05,(1e-3,1e0))
X=f9(FULL)
rng=np.random.RandomState(1); idx=rng.choice(len(X),450,replace=False)
sc=StandardScaler().fit(X[idx])
t0=time.time()
gp=GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=0,random_state=0).fit(sc.transform(X[idx]),y[idx])
theta=gp.kernel_; print('opt',round(time.time()-t0),theta)
oof=np.zeros(len(FULL))
for tri,tei in FOLDS:
    sci=StandardScaler().fit(X[tri])
    g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
    oof[tei]=g.predict(sci.transform(X[tei]))
m=metrics(y,oof,conds)
print(f\"gpr9s rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f}\")
np.save('/tmp/oof_gpr9s.npy', oof)
" 2>&1 | grep -vi warn

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from ens import y, conds, show
O={n:np.load(f'/tmp/oof_{n}.npy') for n in ['et','hgb','gpr8','gprsum','gpr9s']}
show('gpr9s', O['gpr9s'])
show('gprsum+gpr9s', 0.5*O['gprsum']+0.5*O['gpr9s'])
show('gprsum+gpr9s+et', (O['gprsum']+O['gpr9s']+O['et'])/3)
show('.35sum+.35 9s+.3et', .35*O['gprsum']+.35*O['gpr9s']+.3*O['et'])
show('.4sum+.3 9s+.3et', .4*O['gprsum']+.3*O['gpr9s']+.3*O['et'])
show('.3sum+.3 9s+.25et+.15hgb', .3*O['gprsum']+.3*O['gpr9s']+.25*O['et']+.15*O['hgb'])
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from ens import y, conds
from evalharness import metrics
from itertools import product
O={n:np.load(f'/tmp/oof_{n}.npy') for n in ['et','hgb','gpr8','gprsum','gpr9s','rf']}
names=['gprsum','gpr9s','et','hgb','gpr8','rf']
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
grid=[0,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.5]
res=[]
for w in product(grid,repeat=len(names)):
    s=sum(w)
    if abs(s-1.0)>1e-6: continue
    if sum(1 for x in w if x>0)>4: continue
    pred=sum(w[i]*O[names[i]] for i in range(len(names)))
    m=metrics(y,pred,conds)
    if all(m[k]<=T[k] for k in T):
        res.append((m['p90'],m['mx'],m['rmse'],m['mae'],m['p95'],w))
res.sort()
print('Configs passing ALL metrics (full-CV). weights order:',names)
for r in res[:15]:
    print(f'  p90={r[0]:.3f} max={r[1]:.3f} rmse={r[2]:.3f} mae={r[3]:.3f} p95={r[4]:.3f} w={r[5]}')
print('total passing:',len(res))
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from evalharness import metrics
from itertools import product
O={n:np.load(f'/tmp/oof_{n}.npy') for n in ['et','hgb','gprsum','gpr9s']}
names=['gprsum','gpr9s','et','hgb']
grid=[0,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.6]
best_p90=[]; best_max=[]
for w in product(grid,repeat=4):
    if abs(sum(w)-1)>1e-6: continue
    pred=sum(w[i]*O[names[i]] for i in range(4))
    m=metrics(y,pred,conds)
    best_p90.append((m['p90'],m['mx'],m['rmse'],m['mae'],m['p95'],w))
# min max subject to p90<=2.45
cand=[r for r in best_p90 if r[0]<=2.45]
cand.sort(key=lambda r:r[1])
print('Min-max with p90<=2.45 (w order gprsum,gpr9s,et,hgb):')
for r in cand[:6]: print(f'  max={r[1]:.3f} p90={r[0]:.3f} rmse={r[2]:.3f} mae={r[3]:.3f} p95={r[4]:.3f} w={r[5]}')
# identify max-driving condition for equal 3-model
pred=(O['gprsum']+O['gpr9s']+O['et'])/3
df=pd.DataFrame({'e':pred-y,'c':conds})
cr=df.groupby('c')['e'].apply(lambda s:np.sqrt(np.mean(s**2))).sort_values(ascending=False)
print(); print('worst 5 conditions (3-model equal):'); print(cr.head(5).round(2))
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from evalharness import metrics
O={n:np.load(f'/tmp/oof_{n}.npy') for n in ['et','hgb','gprsum','gpr9s','gpr8']}
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
uc=pd.Series(conds.unique())
rng=np.random.RandomState(0)
def simulate(pred, k, n=3000):
    # precompute per-condition errors
    dfp=pd.DataFrame({'e':pred-y,'c':conds})
    groups={c:g['e'].values for c,g in dfp.groupby('c')}
    clist=list(groups)
    passes=0; p90s=[]; maxs=[]
    for _ in range(n):
        sel=rng.choice(len(clist), k, replace=False)
        errs=np.concatenate([groups[clist[i]] for i in sel])
        crs=np.array([np.sqrt(np.mean(groups[clist[i]]**2)) for i in sel])
        rmse=np.sqrt(np.mean(errs**2)); mae=np.mean(np.abs(errs))
        p90=np.percentile(crs,90); p95=np.percentile(crs,95); mx=crs.max()
        ok = rmse<=T['rmse'] and mae<=T['mae'] and p90<=T['p90'] and p95<=T['p95'] and mx<=T['mx']
        passes+=ok; p90s.append(p90); maxs.append(mx)
    return passes/n, np.median(p90s), np.percentile(p90s,75), np.median(maxs)
for name,w in [('sum+9s+et',(.34,.33,.33,0,0)),('.2/.4/.4',(.2,.4,.4,0,0)),('.1/.45/.45',(.1,.45,.45,0,0)),('.25/.4/.35',(.25,.4,.35,0,0))]:
    names=['gprsum','gpr9s','et','hgb','gpr8']
    pred=sum(w[i]*O[names[i]] for i in range(5))
    for k in [15,20]:
        pr,mp90,q75,mmax=simulate(pred,k)
        print(f'{name:12s} k={k}: pass={pr:.2f} median_p90={mp90:.3f} q75_p90={q75:.3f} median_max={mmax:.3f}')
"

cd /app && cat > /tmp/gprgen.py << 'EOF'
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
y=FULL['scaled_sound_pressure_level'].values
groups=FULL['cond'].astype(str).values; conds=FULL['cond']
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,groups))
def base(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    return f,a,c,v,t
def gen(name, feat, seed, nsub=500, nu_long=2.5):
    X=feat(FULL); d=X.shape[1]
    rng=np.random.RandomState(seed); idx=rng.choice(len(X),nsub,replace=False)
    sc=StandardScaler().fit(X[idx])
    k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d)*3,length_scale_bounds=(5e-1,1e2),nu=nu_long) \
     + ConstantKernel(0.4,(1e-3,1e2))*Matern(length_scale=np.ones(d)*0.8,length_scale_bounds=(8e-2,5e0),nu=1.5) \
     + WhiteKernel(0.05,(1e-3,1e0))
    t0=time.time()
    gp=GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=0,random_state=0).fit(sc.transform(X[idx]),y[idx])
    theta=gp.kernel_
    oof=np.zeros(len(FULL))
    for tri,tei in FOLDS:
        sci=StandardScaler().fit(X[tri])
        g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
        oof[tei]=g.predict(sci.transform(X[tei]))
    m=metrics(y,oof,conds)
    print(f"{name:8s} rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f} t={time.time()-t0:.0f}s")
    np.save(f'/tmp/oof_{name}.npy', oof)

def fA(df):
    f,a,c,v,t=base(df); logf=np.log10(f)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),np.log10(f*t/v),np.log10(f*c/v),np.log10(v*c),a*np.log10(f*t/v),a*a])
def fB(df):
    f,a,c,v,t=base(df); logf=np.log10(f)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),np.log10(f*t/v),np.log10(v*c),a*np.log10(f*t/v)])
def fC(df):
    f,a,c,v,t=base(df); logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(f*c/v),np.log10(v*c),a*st,logf*np.log10(c)])
EOF
cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from gprgen import gen, fA, fB, fC
gen('gprA', fA, seed=2, nsub=500)
gen('gprB', fB, seed=3, nsub=500)
gen('gprC', fC, seed=4, nsub=500)
" 2>&1 | grep -vi warn

cd /app && cat > /tmp/sim.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from evalharness import metrics
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
ALL=['et','hgb','rf','gpr8','gprsum','gpr9s','gprA','gprB','gprC']
O={n:np.load(f'/tmp/oof_{n}.npy') for n in ALL}
def fullcv(pred):
    return metrics(y,pred,conds)
dfc=pd.DataFrame({'c':conds})
_groupidx={c:np.where((conds==c).values)[0] for c in conds.unique()}
_clist=list(_groupidx)
def simulate(pred, k, n=4000, seed=0):
    rng=np.random.RandomState(seed)
    err=pred-y
    gerr={c:err[_groupidx[c]] for c in _clist}
    gcr={c:np.sqrt(np.mean(gerr[c]**2)) for c in _clist}
    passes=0; p90s=[]
    nc=len(_clist)
    for _ in range(n):
        sel=rng.choice(nc,k,replace=False)
        cs=[_clist[i] for i in sel]
        errs=np.concatenate([gerr[c] for c in cs])
        crs=np.array([gcr[c] for c in cs])
        rmse=np.sqrt(np.mean(errs**2)); mae=np.mean(np.abs(errs))
        p90=np.percentile(crs,90); p95=np.percentile(crs,95); mx=crs.max()
        p90s.append(p90)
        if rmse<=T['rmse'] and mae<=T['mae'] and p90<=T['p90'] and p95<=T['p95'] and mx<=T['mx']:
            passes+=1
    return passes/n, np.median(p90s)
def blend(weights):
    return sum(w*O[n] for n,w in weights.items())
EOF
cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
from sim import blend, simulate, fullcv, O
import numpy as np
gprs=['gprsum','gpr9s','gprA','gprB','gprC','gpr8']
allgpr=sum(O[n] for n in gprs)/len(gprs)
def show(name,pred):
    m=fullcv(pred); 
    pr15,mp15=simulate(pred,15); pr20,mp20=simulate(pred,20)
    print(f'{name:22s} cv:rmse={m[\"rmse\"]:.3f} mae={m[\"mae\"]:.3f} p90={m[\"p90\"]:.3f} p95={m[\"p95\"]:.3f} max={m[\"mx\"]:.3f} | pass15={pr15:.2f} mp90_15={mp15:.3f} pass20={pr20:.2f}')
show('allGPR', allgpr)
show('allGPR+et .7/.3', .7*allgpr+.3*O['et'])
show('allGPR+et .6/.4', .6*allgpr+.4*O['et'])
show('allGPR+et .5/.5', .5*allgpr+.5*O['et'])
show('allGPR+et+hgb', .6*allgpr+.25*O['et']+.15*O['hgb'])
"

cd /app && timeout 120 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel
# per-condition: fit SPL vs log-f with a flexible GP (in-sample), residual = irreducible noise
res=[]
for c,g in FULL.groupby('cond'):
    if len(g)<6: continue
    x=np.log10(g['frequency'].values).reshape(-1,1); y=g['scaled_sound_pressure_level'].values
    xs=(x-x.mean())/x.std()
    k=ConstantKernel(1.0)*RBF(0.5)+WhiteKernel(0.1)
    gp=GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=0).fit(xs,y)
    p=gp.predict(xs); r=y-p
    res.append((c,np.sqrt(np.mean(r**2)),len(g)))
rr=pd.DataFrame(res,columns=['c','rmse','n'])
print('In-sample per-condition curve-fit RMSE (irreducible-ish noise):')
print('  median %.3f  p75 %.3f  p90 %.3f  max %.3f'%(rr['rmse'].median(),rr['rmse'].quantile(.75),rr['rmse'].quantile(.9),rr['rmse'].max()))
print('  worst 5:'); print(rr.sort_values('rmse',ascending=False).head(5).to_string(index=False))
" 2>&1 | grep -vi warn

cd /app && timeout 280 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
y=FULL['scaled_sound_pressure_level'].values
groups=FULL['cond'].astype(str).values; conds=FULL['cond']
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,groups))
def f9(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(f*c/v),np.log10(v*c),a*st])
X=f9(FULL); d=X.shape[1]
# sample whole conditions to preserve curve structure
uc=conds.unique(); rng=np.random.RandomState(0)
selc=rng.choice(len(uc),48,replace=False)
mask=conds.isin(uc[selc]).values
sc=StandardScaler().fit(X[mask])
k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d)*3,length_scale_bounds=(3e-1,1e2),nu=2.5) \
 + ConstantKernel(0.4,(1e-3,1e2))*Matern(length_scale=np.ones(d)*0.7,length_scale_bounds=(8e-2,5e0),nu=1.5) \
 + WhiteKernel(0.03,(1e-3,1e0))
t0=time.time()
gp=GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=1,random_state=0).fit(sc.transform(X[mask]),y[mask])
theta=gp.kernel_; print('opt',round(time.time()-t0),'s'); print(theta)
oof=np.zeros(len(FULL))
for tri,tei in FOLDS:
    sci=StandardScaler().fit(X[tri])
    g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
    oof[tei]=g.predict(sci.transform(X[tei]))
m=metrics(y,oof,conds)
print(f\"gprW rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f}\")
np.save('/tmp/oof_gprW.npy', oof)
" 2>&1 | grep -vi warn

cd /app && timeout 200 python3 -c "
import sys,time,pickle; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
conds=FULL['cond']; y=FULL['scaled_sound_pressure_level'].values
def f9(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(f*c/v),np.log10(v*c),a*st])
X=f9(FULL); d=X.shape[1]
uc=conds.unique(); rng=np.random.RandomState(0)
selc=rng.choice(len(uc),40,replace=False)
mask=conds.isin(uc[selc]).values
sc=StandardScaler().fit(X[mask])
k=ConstantKernel(1.0,(1e-2,1e3))*Matern(length_scale=np.ones(d)*3,length_scale_bounds=(3e-1,1e2),nu=2.5) \
 + ConstantKernel(0.4,(1e-3,1e2))*Matern(length_scale=np.ones(d)*0.7,length_scale_bounds=(8e-2,5e0),nu=1.5) \
 + WhiteKernel(0.03,(1e-3,1e0))
t0=time.time()
gp=GaussianProcessRegressor(kernel=k,normalize_y=True,n_restarts_optimizer=0,random_state=0).fit(sc.transform(X[mask]),y[mask])
print('opt',round(time.time()-t0),'s'); print(gp.kernel_)
with open('/tmp/thetaW.pkl','wb') as fh: pickle.dump(gp.kernel_, fh)
" 2>&1 | grep -vi warn

cd /app && timeout 200 python3 -c "
import sys,pickle; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
conds=FULL['cond']; y=FULL['scaled_sound_pressure_level'].values
groups=FULL['cond'].astype(str).values
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,groups))
def f9(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(f*c/v),np.log10(v*c),a*st])
X=f9(FULL)
theta=pickle.load(open('/tmp/thetaW.pkl','rb'))
oof=np.zeros(len(FULL))
for tri,tei in FOLDS:
    sci=StandardScaler().fit(X[tri])
    g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
    oof[tei]=g.predict(sci.transform(X[tei]))
m=metrics(y,oof,conds)
print(f\"gprW rmse={m['rmse']:.3f} mae={m['mae']:.3f} p90={m['p90']:.3f} p95={m['p95']:.3f} max={m['mx']:.3f}\")
np.save('/tmp/oof_gprW.npy', oof)
" 2>&1 | grep -vi warn

cd /app && timeout 200 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np
from oofgen import gen
from ff import feats
from sklearn.ensemble import HistGradientBoostingRegressor, ExtraTreesRegressor
gen('hgb2', lambda: HistGradientBoostingRegressor(learning_rate=0.03,max_iter=1200,max_leaf_nodes=15,min_samples_leaf=15,l2_regularization=2.0,max_features=0.7,random_state=0), feats)
gen('et2', lambda: ExtraTreesRegressor(n_estimators=800,min_samples_leaf=1,max_features=0.6,n_jobs=-1,random_state=0), feats)
gen('et3', lambda: ExtraTreesRegressor(n_estimators=800,min_samples_leaf=2,max_features=0.8,n_jobs=-1,random_state=1), feats)
" 2>&1 | grep -vi warn

cd /app && timeout 200 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from evalharness import metrics
from itertools import product
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
gprs=['gprsum','gpr9s','gprA','gprB','gprC']
O={n:np.load(f'/tmp/oof_{n}.npy') for n in gprs+['et','hgb']}
gidx={c:np.where((conds==c).values)[0] for c in conds.unique()}
clist=list(gidx)
def sim(pred,k,n=2500,seed=1):
    rng=np.random.RandomState(seed); err=pred-y
    ge={c:err[gidx[c]] for c in clist}; gc={c:np.sqrt(np.mean(ge[c]**2)) for c in clist}
    nc=len(clist); p=0
    for _ in range(n):
        s=rng.choice(nc,k,replace=False); cs=[clist[i] for i in s]
        e=np.concatenate([ge[c] for c in cs]); crs=np.array([gc[c] for c in cs])
        if np.sqrt(np.mean(e**2))<=T['rmse'] and np.mean(np.abs(e))<=T['mae'] and np.percentile(crs,90)<=T['p90'] and np.percentile(crs,95)<=T['p95'] and crs.max()<=T['mx']: p+=1
    return p/n
# search: GPR mean weight wg split among gprs equally, et weight we, hgb wh
allg=sum(O[n] for n in gprs)/len(gprs)
best=[]
for we in [0,.1,.15,.2,.25,.3,.35,.4]:
  for wh in [0,.1,.15]:
    wg=1-we-wh
    if wg<=0: continue
    pred=wg*allg+we*O['et']+wh*O['hgb']
    m=metrics(y,pred,conds)
    p15=sim(pred,15); p20=sim(pred,20)
    best.append((p15+p20, p15,p20, m['p90'],m['mx'],m['rmse'],we,wh))
best.sort(reverse=True)
print('wg=1-we-wh split equally among 5 GPRs. top by pass15+pass20:')
for b in best[:10]:
    print(f'  pass15={b[1]:.2f} pass20={b[2]:.2f} | cv p90={b[3]:.3f} max={b[4]:.3f} rmse={b[5]:.3f} we={b[6]} wh={b[7]}')
"

cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
conds=FULL['cond']; y=FULL['scaled_sound_pressure_level'].values
groups=conds.astype(str).values
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,groups))
def f9(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float)
    c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float)
    t=df['suction_side_displacement_thickness'].values.astype(float)
    logf=np.log10(f); st=np.log10(f*t/v)
    return np.column_stack([logf,a,np.log10(c),v,np.log10(t),st,np.log10(f*c/v),np.log10(v*c),a*st])
X=f9(FULL)
# base length scales from gprW (whole-cond optimized), vary noise
long_ls=[4.97,4.34,100,14.6,4.88,3.85,3.04,6.14,100]
short_ls=[1.01,0.675,0.347,5,1.55,5,5,0.638,0.782]
for noise in [0.01,0.02,0.04,0.08]:
    theta=ConstantKernel(2.38**2)*Matern(length_scale=long_ls,nu=2.5)+ConstantKernel(0.405**2)*Matern(length_scale=short_ls,nu=1.5)+WhiteKernel(noise)
    oof=np.zeros(len(FULL))
    for tri,tei in FOLDS:
        sci=StandardScaler().fit(X[tri])
        g=GaussianProcessRegressor(kernel=theta,normalize_y=True,optimizer=None).fit(sci.transform(X[tri]),y[tri])
        oof[tei]=g.predict(sci.transform(X[tei]))
    m=metrics(y,oof,conds)
    print(f'noise={noise:.3f} rmse={m[\"rmse\"]:.3f} mae={m[\"mae\"]:.3f} p90={m[\"p90\"]:.3f} p95={m[\"p95\"]:.3f} max={m[\"mx\"]:.3f}')
    np.save(f'/tmp/oof_gprN{int(noise*1000)}.npy', oof)
" 2>&1 | grep -vi warn

cd /app && cat > /tmp/valproxy2.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from evalharness import cond_labels, metrics, report
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
from sklearn.ensemble import ExtraTreesRegressor
tr=pd.read_csv('/app/data/train.csv'); va=pd.read_csv('/app/data/validation.csv')
ytr=tr['scaled_sound_pressure_level'].values; yva=va['scaled_sound_pressure_level'].values
cva=cond_labels(va)
def f8(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float); c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float); t=df['suction_side_displacement_thickness'].values.astype(float)
    return np.column_stack([np.log10(f),a,np.log10(c),v,np.log10(t),np.log10(f*t/v),np.log10(v*c),a*np.log10(f*t/v)])
def f9(df):
    f=df['frequency'].values.astype(float); a=df['attack_angle'].values.astype(float); c=df['chord_length'].values.astype(float); v=df['free_stream_velocity'].values.astype(float); t=df['suction_side_displacement_thickness'].values.astype(float)
    return np.column_stack([np.log10(f),a,np.log10(c),v,np.log10(t),np.log10(f*t/v),np.log10(f*c/v),np.log10(v*c),a*np.log10(f*t/v)])
def feats(df):
    f=df['frequency'].astype(float).values;a=df['attack_angle'].astype(float).values;c=df['chord_length'].astype(float).values;v=df['free_stream_velocity'].astype(float).values;t=df['suction_side_displacement_thickness'].astype(float).values
    logf=np.log10(f);logt=np.log10(t);logc=np.log10(c);st=np.log10(f*t/v);re=np.log10(v*c)
    return np.column_stack([logf,a,logc,v,logt,st,re,logf*logf,st*st,a*a,logf*v,logf*a,logf*logt,a*logc,logf*logc])
Ksum=ConstantKernel(1.44**2)*Matern(length_scale=[2.21,1.44,2.51,100,2.95,2.11,4.29,2.35],nu=2.5)+ConstantKernel(0.299**2)*Matern(length_scale=[0.434,1.08,0.282,5,5,5,0.534,2.5],nu=1.5)+WhiteKernel(0.02)
K9s=ConstantKernel(1.66**2)*Matern(length_scale=[2.5,4.16,100,100,3.11,3.21,2.66,4.46,100],nu=2.5)+ConstantKernel(0.366**2)*Matern(length_scale=[5,0.314,0.434,5,5,5,5,0.552,0.569],nu=1.5)+WhiteKernel(0.02)
def gpr(K,featfn):
    Xt=featfn(tr);Xv=featfn(va);sc=StandardScaler().fit(Xt)
    g=GaussianProcessRegressor(kernel=K,normalize_y=True,optimizer=None).fit(sc.transform(Xt),ytr)
    return g.predict(sc.transform(Xv))
psum=gpr(Ksum,f8); p9s=gpr(K9s,f9)
et=ExtraTreesRegressor(n_estimators=800,n_jobs=-1,random_state=0).fit(feats(tr),ytr).predict(feats(va))
report('val psum',metrics(yva,psum,cva))
report('val p9s',metrics(yva,p9s,cva))
report('val gpr-avg',metrics(yva,(psum+p9s)/2,cva))
report('val .4sum+.4 9s+.2et',metrics(yva,.4*psum+.4*p9s+.2*et,cva))
report('val .35sum+.35 9s+.3et',metrics(yva,.35*psum+.35*p9s+.3*et,cva))
report('val .45sum+.35 9s+.2et',metrics(yva,.45*psum+.35*p9s+.2*et,cva))
EOF
cd /app && timeout 90 python3 /tmp/valproxy2.py 2>&1 | grep -vi warn

cd /app && timeout 120 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from cv import FULL
from evalharness import metrics
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, Matern
conds=FULL['cond']; y=FULL['scaled_sound_pressure_level'].values
FOLDS=list(GroupKFold(n_splits=8).split(FULL,y,conds.astype(str).values))
def f8(df):
    f=df['frequency'].values.astype(float);a=df['attack_angle'].values.astype(float);c=df['chord_length'].values.astype(float);v=df['free_stream_velocity'].values.astype(float);t=df['suction_side_displacement_thickness'].values.astype(float)
    return np.column_stack([np.log10(f),a,np.log10(c),v,np.log10(t),np.log10(f*t/v),np.log10(v*c),a*np.log10(f*t/v)])
def f9(df):
    f=df['frequency'].values.astype(float);a=df['attack_angle'].values.astype(float);c=df['chord_length'].values.astype(float);v=df['free_stream_velocity'].values.astype(float);t=df['suction_side_displacement_thickness'].values.astype(float)
    return np.column_stack([np.log10(f),a,np.log10(c),v,np.log10(t),np.log10(f*t/v),np.log10(f*c/v),np.log10(v*c),a*np.log10(f*t/v)])
Ksum=ConstantKernel(1.44**2)*Matern(length_scale=[2.21,1.44,2.51,100,2.95,2.11,4.29,2.35],nu=2.5)+ConstantKernel(0.299**2)*Matern(length_scale=[0.434,1.08,0.282,5,5,5,0.534,2.5],nu=1.5)+WhiteKernel(0.02)
K9s=ConstantKernel(1.66**2)*Matern(length_scale=[2.5,4.16,100,100,3.11,3.21,2.66,4.46,100],nu=2.5)+ConstantKernel(0.366**2)*Matern(length_scale=[5,0.314,0.434,5,5,5,5,0.552,0.569],nu=1.5)+WhiteKernel(0.02)
for nm,K,ff in [('psum',Ksum,f8),('p9s',K9s,f9)]:
    X=ff(FULL);oof=np.zeros(len(FULL))
    for tri,tei in FOLDS:
        sc=StandardScaler().fit(X[tri])
        g=GaussianProcessRegressor(kernel=K,normalize_y=True,optimizer=None).fit(sc.transform(X[tri]),y[tri])
        oof[tei]=g.predict(sc.transform(X[tei]))
    np.save(f'/tmp/oof_{nm}.npy',oof)
    m=metrics(y,oof,conds); print(nm, {k:round(m[k],3) for k in ['rmse','mae','p90','p95','mx']})
print('done')
" 2>&1 | grep -vi warn

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from evalharness import metrics
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
O={n:np.load(f'/tmp/oof_{n}.npy') for n in ['psum','p9s','et','hgb']}
gidx={c:np.where((conds==c).values)[0] for c in conds.unique()}; clist=list(gidx)
def sim(pred,k,n=4000,seed=7):
    rng=np.random.RandomState(seed); err=pred-y
    ge={c:err[gidx[c]] for c in clist}; gc={c:np.sqrt(np.mean(ge[c]**2)) for c in clist}
    nc=len(clist); p=0
    for _ in range(n):
        s=rng.choice(nc,k,replace=False); cs=[clist[i] for i in s]
        e=np.concatenate([ge[c] for c in cs]); crs=np.array([gc[c] for c in cs])
        if np.sqrt(np.mean(e**2))<=T['rmse'] and np.mean(np.abs(e))<=T['mae'] and np.percentile(crs,90)<=T['p90'] and np.percentile(crs,95)<=T['p95'] and crs.max()<=T['mx']: p+=1
    return p/n
for wsum,w9,we,wh in [(.5,.5,0,0),(.4,.4,.2,0),(.35,.35,.3,0),(.3,.3,.25,.15),(.45,.35,.2,0),(.4,.35,.25,0),(.33,.33,.34,0),(.35,.35,.2,.1)]:
    pred=wsum*O['psum']+w9*O['p9s']+we*O['et']+wh*O['hgb']
    m=metrics(y,pred,conds)
    print(f'w(sum,9s,et,hgb)=({wsum},{w9},{we},{wh}) cv p90={m[\"p90\"]:.3f} max={m[\"mx\"]:.3f} rmse={m[\"rmse\"]:.3f} mae={m[\"mae\"]:.3f} | pass15={sim(pred,15):.2f} pass20={sim(pred,20):.2f} pass12={sim(pred,12):.2f}')
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from evalharness import metrics
from itertools import product
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
names=['psum','p9s','gpr8','et']
O={n:np.load(f'/tmp/oof_{n}.npy') for n in names}
gidx={c:np.where((conds==c).values)[0] for c in conds.unique()}; clist=list(gidx)
def sim(pred,k,n=3000,seed=7):
    rng=np.random.RandomState(seed); err=pred-y
    ge={c:err[gidx[c]] for c in clist}; gc={c:np.sqrt(np.mean(ge[c]**2)) for c in clist}
    nc=len(clist); p=0
    for _ in range(n):
        s=rng.choice(nc,k,replace=False); cs=[clist[i] for i in s]
        e=np.concatenate([ge[c] for c in cs]); crs=np.array([gc[c] for c in cs])
        if np.sqrt(np.mean(e**2))<=T['rmse'] and np.mean(np.abs(e))<=T['mae'] and np.percentile(crs,90)<=T['p90'] and np.percentile(crs,95)<=T['p95'] and crs.max()<=T['mx']: p+=1
    return p/n
grid=[0,.1,.15,.2,.25,.3,.35,.4,.5]
res=[]
for w in product(grid,repeat=4):
    if abs(sum(w)-1)>1e-9: continue
    pred=sum(w[i]*O[names[i]] for i in range(4))
    m=metrics(y,pred,conds)
    res.append((min(sim(pred,15),sim(pred,20)), sim(pred,15),sim(pred,20),m['p90'],m['mx'],m['rmse'],w))
res.sort(reverse=True)
print('order',names,'| top by min(pass15,pass20):')
for r in res[:12]:
    print(f'  minpass={r[0]:.2f} p15={r[1]:.2f} p20={r[2]:.2f} cv_p90={r[3]:.3f} cv_max={r[4]:.3f} rmse={r[5]:.3f} w={r[6]}')
"

restart: true

cd /app && timeout 120 python3 -c "
import sys; sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
from ens import y, conds
from evalharness import metrics
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,mx=5.60)
names=['psum','p9s','gpr8','et']
O={n:np.load(f'/tmp/oof_{n}.npy') for n in names}
gidx={c:np.where((conds==c).values)[0] for c in conds.unique()}; clist=list(gidx); nc=len(clist)
def sim(pred,k,n=2000,seed=7):
    rng=np.random.RandomState(seed); err=pred-y
    ge={c:err[gidx[c]] for c in clist}; gc={c:np.sqrt(np.mean(ge[c]**2)) for c in clist}
    p=0
    for _ in range(n):
        s=rng.choice(nc,k,replace=False); cs=[clist[i] for i in s]
        e=np.concatenate([ge[c] for c in cs]); crs=np.array([gc[c] for c in cs])
        if np.sqrt(np.mean(e**2))<=T['rmse'] and np.mean(np.abs(e))<=T['mae'] and np.percentile(crs,90)<=T['p90'] and np.percentile(crs,95)<=T['p95'] and crs.max()<=T['mx']: p+=1
    return p/n
cands={'A .4/.4/0/.2':(.4,.4,0,.2),'B .3/.3/.2/.2':(.3,.3,.2,.2),'C .3/.3/.15/.25':(.3,.3,.15,.25),
       'D .25/.25/.25/.25':(.25,.25,.25,.25),'E .35/.3/.15/.2':(.35,.3,.15,.2),'F .3/.35/.15/.2':(.3,.35,.15,.2)}
for nm,w in cands.items():
    pred=sum(w[i]*O[names[i]] for i in range(4))
    m=metrics(y,pred,conds)
    print(f'{nm:20s} cv p90={m[\"p90\"]:.3f} max={m[\"mx\"]:.3f} rmse={m[\"rmse\"]:.3f} mae={m[\"mae\"]:.3f} p95={m[\"p95\"]:.3f} | p15={sim(pred,15):.2f} p20={sim(pred,20):.2f}')
"

command: create
path: /app/train_and_predict.py
file_text: #!/usr/bin/env python3
"""Airfoil self-noise surrogate.

Trains on the labeled public data (train + validation) and writes predictions
for the grader-staged evaluation features at /app/data/test_features.csv.

Modeling notes
--------------
The benchmark holds out complete aerodynamic *condition* groups
``(attack_angle, chord_length, free_stream_velocity)``, so the evaluation tests
generalization to unseen physical regimes.  Within a condition only frequency
varies, tracing a smooth sound-pressure spectrum in log-frequency.

The surrogate is a blend of two Gaussian-process regressors (two-scale Matern
kernels over physically motivated, log-scaled features including the Strouhal
number ``f*delta*/U``) and an ExtraTrees regressor.  The GPs capture the smooth
spectral level/shape and extrapolate across regimes; the tree model clamps the
GP's occasional wild extrapolation on extreme held-out corners, which controls
the worst per-condition error.  Kernel hyperparameters were tuned offline with
grouped (by-condition) cross-validation and are frozen here for a fast,
deterministic re-run.
"""

from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, Matern, WhiteKernel
from sklearn.preprocessing import StandardScaler

DATA_DIR = Path("/app/data")
OUTPUT_DIR = Path("/app/output")
TARGET = "scaled_sound_pressure_level"
PREDICTION = "predicted_scaled_sound_pressure_level"
FEATURES = [
    "frequency",
    "attack_angle",
    "chord_length",
    "free_stream_velocity",
    "suction_side_displacement_thickness",
]

# Blend weights (GP-heavy for low global/typical-condition error, with an
# ExtraTrees component to bound the worst held-out condition).
W_GP_SUM = 0.4
W_GP_9 = 0.4
W_ET = 0.2


def _raw(df: pd.DataFrame):
    f = df["frequency"].astype(float).to_numpy()
    a = df["attack_angle"].astype(float).to_numpy()
    c = df["chord_length"].astype(float).to_numpy()
    v = df["free_stream_velocity"].astype(float).to_numpy()
    t = df["suction_side_displacement_thickness"].astype(float).to_numpy()
    return f, a, c, v, t


def feats_gp_sum(df: pd.DataFrame) -> np.ndarray:
    """8 features used by the first Gaussian process."""
    f, a, c, v, t = _raw(df)
    st = np.log10(f * t / v)  # Strouhal number (log)
    return np.column_stack(
        [np.log10(f), a, np.log10(c), v, np.log10(t), st, np.log10(v * c), a * st]
    )


def feats_gp_9(df: pd.DataFrame) -> np.ndarray:
    """9 features used by the second Gaussian process (adds chord Strouhal)."""
    f, a, c, v, t = _raw(df)
    st = np.log10(f * t / v)
    stc = np.log10(f * c / v)  # chord-based Strouhal (log)
    return np.column_stack(
        [np.log10(f), a, np.log10(c), v, np.log10(t), st, stc, np.log10(v * c), a * st]
    )


def feats_tree(df: pd.DataFrame) -> np.ndarray:
    """Richer feature set (with interactions) for the ExtraTrees model."""
    f, a, c, v, t = _raw(df)
    logf = np.log10(f)
    logt = np.log10(t)
    logc = np.log10(c)
    st = np.log10(f * t / v)
    re = np.log10(v * c)
    return np.column_stack(
        [logf, a, logc, v, logt, st, re, logf * logf, st * st, a * a,
         logf * v, logf * a, logf * logt, a * logc, logf * logc]
    )


def kernel_gp_sum():
    return (
        ConstantKernel(1.44 ** 2)
        * Matern(length_scale=[2.21, 1.44, 2.51, 100.0, 2.95, 2.11, 4.29, 2.35], nu=2.5)
        + ConstantKernel(0.299 ** 2)
        * Matern(length_scale=[0.434, 1.08, 0.282, 5.0, 5.0, 5.0, 0.534, 2.5], nu=1.5)
        + WhiteKernel(0.02)
    )


def kernel_gp_9():
    return (
        ConstantKernel(1.66 ** 2)
        * Matern(length_scale=[2.5, 4.16, 100.0, 100.0, 3.11, 3.21, 2.66, 4.46, 100.0], nu=2.5)
        + ConstantKernel(0.366 ** 2)
        * Matern(length_scale=[5.0, 0.314, 0.434, 5.0, 5.0, 5.0, 5.0, 0.
…[truncated]

cd /app && rm -f output/predictions.csv && time python3 train_and_predict.py 2>&1 | grep -vi warn; echo "=== OUTPUT ==="; head -5 output/predictions.csv; echo "..."; wc -l output/predictions.csv

cd /app && python3 -c "
import pandas as pd, numpy as np
tf=pd.read_csv('data/test_features.csv'); pr=pd.read_csv('output/predictions.csv')
assert list(pr.columns)==['sample_id','predicted_scaled_sound_pressure_level'], pr.columns
assert len(pr)==len(tf), (len(pr),len(tf))
assert pr['sample_id'].tolist()==tf['sample_id'].astype(str).tolist(), 'id mismatch/order'
assert pr['sample_id'].is_unique
v=pr['predicted_scaled_sound_pressure_level'].to_numpy()
assert np.all(np.isfinite(v)), 'non-finite'
assert v.std()>1.0, 'looks constant'
print('FORMAT OK: rows',len(pr),'range [%.2f, %.2f] std %.2f'%(v.min(),v.max(),v.std()))
"

cd /app && timeout 90 python3 -c "
import sys; sys.path.insert(0,'/tmp'); sys.path.insert(0,'/app')
import numpy as np, pandas as pd
import train_and_predict as tp
from evalharness import cond_labels, metrics, report
tr=pd.read_csv('data/train.csv'); va=pd.read_csv('data/validation.csv')
y=tr[tp.TARGET].to_numpy()
sc_s,gp_s=tp._fit_gp(tp.kernel_gp_sum(),tp.feats_gp_sum,tr,y)
sc_9,gp_9=tp._fit_gp(tp.kernel_gp_9(),tp.feats_gp_9,tr,y)
from sklearn.ensemble import ExtraTreesRegressor
et=ExtraTreesRegressor(n_estimators=800,n_jobs=-1,random_state=0).fit(tp.feats_tree(tr),y)
ps=gp_s.predict(sc_s.transform(tp.feats_gp_sum(va)))
p9=gp_9.predict(sc_9.transform(tp.feats_gp_9(va)))
pe=et.predict(tp.feats_tree(va))
pred=tp.W_GP_SUM*ps+tp.W_GP_9*p9+tp.W_ET*pe
report('DEPLOYED train->val', metrics(va[tp.TARGET].to_numpy(), pred, cond_labels(va)))
print('targets: rmse<=2.25 mae<=1.55 p90<=2.35 p95<=4.75 max<=5.60')
" 2>&1 | grep -vi warn

command: str_replace
path: /app/train_and_predict.py
old_str: def main() -> None:
    labeled = load_labeled_data()
    evaluation = load_evaluation_features()
    y = labeled[TARGET].astype(float).to_numpy()
new_str: def main() -> None:
    labeled = load_labeled_data()
    evaluation = load_evaluation_features().copy()
    y = labeled[TARGET].astype(float).to_numpy()

    # Robustness: ensure required feature columns exist and are finite by
    # imputing any missing/non-finite evaluation values with training medians.
    for col in FEATURES:
        med = float(labeled[col].astype(float).median())
        if col not in evaluation.columns:
            evaluation[col] = med
        else:
            vals = pd.to_numeric(evaluation[col], errors="coerce")
            evaluation[col] = vals.where(np.isfinite(vals), med)

cd /app && rm -f output/predictions.csv && python3 train_and_predict.py 2>&1 | grep -vi warn; python3 -c "
import pandas as pd, numpy as np
tf=pd.read_csv('data/test_features.csv'); pr=pd.read_csv('output/predictions.csv')
assert list(pr.columns)==['sample_id','predicted_scaled_sound_pressure_level']
assert pr['sample_id'].tolist()==tf['sample_id'].astype(str).tolist()
assert len(pr)==len(tf) and pr['sample_id'].is_unique
v=pr['predicted_scaled_sound_pressure_level'].to_numpy()
assert np.all(np.isfinite(v)) and v.std()>1.0
print('OK after hardening: rows',len(pr),'std %.2f range [%.1f,%.1f]'%(v.std(),v.min(),v.max()))
"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/app')
import numpy as np, pandas as pd
import train_and_predict as tp
labeled=tp.load_labeled_data()
# craft a nasty evaluation frame: NaN, a non-numeric string, and a missing column
ev=pd.read_csv('data/test_features.csv').head(5).copy()
ev.loc[0,'frequency']=np.nan
ev.loc[1,'attack_angle']='bad'
ev=ev.drop(columns=['suction_side_displacement_thickness'])  # simulate missing col
# replicate main's imputation
for col in tp.FEATURES:
    med=float(labeled[col].astype(float).median())
    if col not in ev.columns: ev[col]=med
    else:
        vals=pd.to_numeric(ev[col],errors='coerce'); ev[col]=vals.where(np.isfinite(vals),med)
X=tp.feats_gp_sum(ev); X9=tp.feats_gp_9(ev); Xt=tp.feats_tree(ev)
print('all finite after impute:', np.all(np.isfinite(X)) and np.all(np.isfinite(X9)) and np.all(np.isfinite(Xt)))
print('shapes', X.shape, X9.shape, Xt.shape)
"

cd /app && sed -n '140,210p' train_and_predict.py

cd /app && md5sum data/*.csv > /tmp/before.md5
rm -f output/predictions.csv
python3 train_and_predict.py 2>&1 | grep -vi warn
cp output/predictions.csv /tmp/run1.csv
python3 train_and_predict.py 2>&1 | grep -vi warn
echo "--- determinism (diff of two runs; empty = identical) ---"
diff /tmp/run1.csv output/predictions.csv && echo "IDENTICAL across runs"
echo "--- input files unchanged? ---"
md5sum -c /tmp/before.md5
echo "--- final output head ---"
head -3 output/predictions.csv

cd /app && timeout 90 python3 -c "
import sys; sys.path.insert(0,'/tmp'); sys.path.insert(0,'/app')
import numpy as np, pandas as pd
import train_and_predict as tp
from evalharness import cond_labels, metrics
from sklearn.ensemble import ExtraTreesRegressor
tr=pd.read_csv('data/train.csv'); va=pd.read_csv('data/validation.csv')
y=tr[tp.TARGET].to_numpy(); yv=va[tp.TARGET].to_numpy(); cv=cond_labels(va)
sc_s,gp_s=tp._fit_gp(tp.kernel_gp_sum(),tp.feats_gp_sum,tr,y)
sc_9,gp_9=tp._fit_gp(tp.kernel_gp_9(),tp.feats_gp_9,tr,y)
et=ExtraTreesRegressor(n_estimators=800,n_jobs=-1,random_state=0).fit(tp.feats_tree(tr),y)
ps=gp_s.predict(sc_s.transform(tp.feats_gp_sum(va))); p9=gp_9.predict(sc_9.transform(tp.feats_gp_9(va))); pe=et.predict(tp.feats_tree(va))
for ws,w9,we in [(.5,.5,0),(.45,.45,.1),(.4,.4,.2),(.35,.35,.3),(.3,.3,.4)]:
    m=metrics(yv,ws*ps+w9*p9+we*pe,cv)
    print(f'({ws},{w9},{we}) val: p90={m[\"p90\"]:.3f} max={m[\"mx\"]:.3f} rmse={m[\"rmse\"]:.3f} mae={m[\"mae\"]:.3f} p95={m[\"p95\"]:.3f}')
" 2>&1 | grep -vi warn

cd /app && ls -la output/predictions.csv && echo "rows (incl header):" && wc -l < output/predictions.csv && echo "=== head ===" && head -3 output/predictions.csv