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 FEATURES---" && head -5 data/test_features.csv && wc -l data/test_features.csv

cd /app && python3 --version && echo "---" && pip list 2>/dev/null | grep -iE "scikit|numpy|pandas|scipy|xgboost|lightgbm|torch|tensorflow|catboost" && echo "---CPU---" && nproc && echo "---MEM---" && free -h 2>/dev/null | head -2

cd /app && ls -la && echo "---output---" && ls -la output/ 2>/dev/null && echo "---root---" && ls / && echo "---check write perms---" && touch /app/output/testwrite 2>&1 && echo "can write output" && rm -f /app/output/testwrite

cd /app && python3 -c "
import pandas as pd
import numpy as np
train = pd.read_csv('data/train.csv')
val = pd.read_csv('data/validation.csv')
print('TRAIN shape', train.shape)
print('VAL shape', val.shape)
print()
print('--- TRAIN describe ---')
print(train.describe().T)
print()
print('--- VAL describe ---')
print(val.describe().T)
"

cd /app && python3 -c "
import pandas as pd
import numpy as np
pd.set_option('display.width', 200)
pd.set_option('display.max_columns', 20)
train = pd.read_csv('data/train.csv')
val = pd.read_csv('data/validation.csv')
for name, df in [('TRAIN', train), ('VAL', val)]:
    print('===', name, '===')
    for c in df.columns:
        if c=='sample_id': continue
        print(f'{c:40s} min={df[c].min():.5g} max={df[c].max():.5g} n_unique={df[c].nunique()}')
    print()
"

cd /app && python3 -c "
import pandas as pd
import numpy as np
train = pd.read_csv('data/train.csv')
val = pd.read_csv('data/validation.csv')
alld = pd.concat([train,val])
# Original UCI canonical values
print('chord_length sorted unique-ish (rounded 3):')
print(np.sort(alld['chord_length'].round(3).unique()))
print()
print('free_stream_velocity rounded 0:')
print(np.sort(alld['free_stream_velocity'].round(0).unique()))
print()
print('Check clustering of chord_length - histogram of rounded to 3 decimals:')
print(alld['chord_length'].round(3).value_counts().sort_index())
print()
print('free_stream_velocity rounded 1:')
print(alld['free_stream_velocity'].round(1).value_counts().sort_index())
"

cd /app && python3 -c "
import pandas as pd
import numpy as np
train = pd.read_csv('data/train.csv')
val = pd.read_csv('data/validation.csv')
alld = pd.concat([train,val])
print('attack_angle rounded 1 value counts:')
print(alld['attack_angle'].round(1).value_counts().sort_index())
print()
print('n unique rounded:', alld['attack_angle'].round(1).nunique())
"

cd /app && python3 -c "
import pandas as pd
import numpy as np
train = pd.read_csv('data/train.csv')
val = pd.read_csv('data/validation.csv')
train['src']='train'; val['src']='val'
alld = pd.concat([train,val], ignore_index=True)

# Snap chord and velocity to canonical grids
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, grid):
    return grid[np.argmin(np.abs(grid[None,:]-np.asarray(x)[:,None]),axis=1)]
alld['chord_s']=snap(alld['chord_length'].values, chords)
alld['vel_s']=snap(alld['free_stream_velocity'].values, vels)
# For angle, cluster: sort unique and group within tolerance
print('Check residual of chord snap (max abs rel):', np.max(np.abs(alld['chord_length']-alld['chord_s'])/alld['chord_s']))
print('Check residual of vel snap (max abs):', np.max(np.abs(alld['free_stream_velocity']-alld['vel_s'])))
# thickness per (chord_s, vel_s) - is it determined also by angle? check variance
print()
print('Number of (chord_s, vel_s) combos:', alld.groupby(['chord_s','vel_s']).ngroups)
"

cd /app && python3 -c "
import pandas as pd
import numpy as np
train = pd.read_csv('data/train.csv')
val = pd.read_csv('data/validation.csv')
alld = pd.concat([train,val], ignore_index=True)
a = np.sort(alld['attack_angle'].values)
# find gaps
gaps = np.diff(a)
# print values where gap > 0.2 to identify clusters
clusters=[[a[0]]]
for i in range(1,len(a)):
    if a[i]-a[i-1] > 0.3:
        clusters.append([])
    clusters[-1].append(a[i])
print('num clusters with gap>0.3:', len(clusters))
for c in clusters:
    print(f'  center~{np.mean(c):7.3f}  n={len(c):3d}  range=[{min(c):.3f},{max(c):.3f}] width={max(c)-min(c):.3f}')
"

cd /app && python3 -c "
import pandas as pd
import numpy as np
train = pd.read_csv('data/train.csv')
val = pd.read_csv('data/validation.csv')
train['src']='train'; val['src']='val'
alld = pd.concat([train,val], ignore_index=True)
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, grid):
    return grid[np.argmin(np.abs(grid[None,:]-np.asarray(x)[:,None]),axis=1)]
alld['chord_s']=snap(alld['chord_length'].values, chords)
alld['vel_s']=snap(alld['free_stream_velocity'].values, vels)
# snap angle by rounding to nearest 0.1 then cluster centers
a = alld['attack_angle'].values
order = np.argsort(a)
asort = a[order]
grp = np.zeros(len(a),dtype=int)
g=0
for i in range(1,len(asort)):
    if asort[i]-asort[i-1] > 0.09:
        g+=1
    grp[order[i]]=g
alld['angle_g']=grp
print('num angle groups (gap>0.09):', alld['angle_g'].nunique())
# condition = (chord_s, vel_s, angle_g)
alld['cond']=alld['chord_s'].astype(str)+'_'+alld['vel_s'].astype(str)+'_'+alld['angle_g'].astype(str)
print('num conditions total:', alld['cond'].nunique())
# within-condition thickness spread
g = alld.groupby('cond')['suction_side_displacement_thickness']
rel_spread = ((g.max()-g.min())/g.mean())
print('within-cond thickness rel spread: median={:.4f} p95={:.4f} max={:.4f}'.format(rel_spread.median(), rel_spread.quantile(.95), rel_spread.max()))
# overlap train/val
tc = set(alld[alld.src=='train']['cond']); vc=set(alld[alld.src=='val']['cond'])
print('train conds:', len(tc), 'val conds:', len(vc), 'overlap:', len(tc&vc))
print('val-only conds:', len(vc-tc))
"

cd /app && cat > /tmp/harness.py << 'EOF'
import pandas as pd, numpy as np
from sklearn.model_selection import GroupKFold

DATA='/app/data'
TARGET='scaled_sound_pressure_level'
FEATURES=["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]

def load():
    tr=pd.read_csv(f'{DATA}/train.csv'); va=pd.read_csv(f'{DATA}/validation.csv')
    return tr, va

def recover_conditions(df):
    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,grid): return grid[np.argmin(np.abs(grid[None,:]-np.asarray(x)[:,None]),axis=1)]
    cs=snap(df['chord_length'].values,chords)
    vs=snap(df['free_stream_velocity'].values,vels)
    a=df['attack_angle'].values
    order=np.argsort(a); asort=a[order]; grp=np.zeros(len(a),dtype=int); g=0
    for i in range(1,len(asort)):
        if asort[i]-asort[i-1]>0.09: g+=1
        grp[order[i]]=g
    cond=[f'{c:.4f}_{v:.1f}_{gg}' for c,v,gg in zip(cs,vs,grp)]
    return np.array(cond)

def per_condition_rmse(y_true,y_pred,cond):
    df=pd.DataFrame({'y':y_true,'p':y_pred,'c':cond})
    out=df.groupby('c').apply(lambda d: np.sqrt(np.mean((d['y']-d['p'])**2)))
    return out

def summarize(y_true,y_pred,cond,label=''):
    rmse=np.sqrt(np.mean((y_true-y_pred)**2))
    mae=np.mean(np.abs(y_true-y_pred))
    pc=per_condition_rmse(y_true,y_pred,cond)
    p90=np.percentile(pc,90); p95=np.percentile(pc,95); mx=pc.max()
    print(f'{label:22s} RMSE={rmse:.3f} MAE={mae:.3f} | condRMSE p90={p90:.3f} p95={p95:.3f} max={mx:.3f} (nconds={len(pc)})')
    return dict(rmse=rmse,mae=mae,p90=p90,p95=p95,mx=mx)

def cv_eval(make_model, X, y, groups, nsplits=8, label=''):
    gkf=GroupKFold(n_splits=nsplits)
    oof=np.zeros(len(y))
    for tr,te in gkf.split(X,y,groups):
        m=make_model()
        m.fit(X.iloc[tr],y.iloc[tr])
        oof[te]=m.predict(X.iloc[te])
    return summarize(y.values,oof,groups,label=label)
EOF
echo "harness written"

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
import numpy as np, pandas as pd
tr,va=load()
alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld)
print('total conds', len(set(cond)), 'rows', len(alld))
X=alld[FEATURES]; y=alld[TARGET]

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

# Starter features
class AF(BaseEstimator,TransformerMixin):
    def fit(self,X,y=None): return self
    def transform(self,X):
        f=pd.DataFrame(X,columns=FEATURES).astype(float)
        lf=np.log10(f['frequency']); lt=np.log10(f['suction_side_displacement_thickness'])
        o=pd.DataFrame(index=f.index)
        o['frequency']=f['frequency']; o['attack_angle']=f['attack_angle']; o['chord_length']=f['chord_length']
        o['free_stream_velocity']=f['free_stream_velocity']; o['thk']=f['suction_side_displacement_thickness']
        o['lf']=lf; o['lt']=lt; o['fv']=lf*f['free_stream_velocity']; o['fa']=lf*f['attack_angle']
        o['av']=f['attack_angle']*f['free_stream_velocity']; o['ta']=lt*f['attack_angle']; o['cv']=f['chord_length']*f['free_stream_velocity']
        return o.to_numpy(float)

cv_eval(lambda: make_pipeline(AF(),StandardScaler(),Ridge(alpha=10.0)), X,y,cond,label='Starter Ridge')

from sklearn.ensemble import HistGradientBoostingRegressor, RandomForestRegressor, ExtraTreesRegressor
cv_eval(lambda: HistGradientBoostingRegressor(random_state=0), X,y,cond,label='HGB default')
cv_eval(lambda: RandomForestRegressor(n_estimators=300,random_state=0,n_jobs=-1), X,y,cond,label='RF default')
cv_eval(lambda: ExtraTreesRegressor(n_estimators=300,random_state=0,n_jobs=-1), X,y,cond,label='ET default')
from sklearn.neighbors import KNeighborsRegressor
" 2>&1 | grep -v Warning

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
import numpy as np, pandas as pd
tr,va=load()
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,grid): return grid[np.argmin(np.abs(grid[None,:]-np.asarray(x)[:,None]),axis=1)]
for col,grid in [('chord_length',chords),('free_stream_velocity',vels)]:
    t=set(np.round(snap(tr[col].values,grid),4)); v=set(np.round(snap(va[col].values,grid),4))
    print(col,'val-only snapped:', sorted(v-t))
# angle: round to 0.1 after clustering is messy; check raw nearest
ta=np.sort(tr['attack_angle'].round(1).unique()); vaa=np.sort(va['attack_angle'].round(1).unique())
print('train angles(round1):',ta)
print('val angles(round1):',vaa)
# for each val angle, min dist to a train angle
for x in vaa:
    d=np.min(np.abs(ta-x))
    if d>0.15: print(f'  val angle {x} min dist to train {d:.2f} (possibly unseen)')
print('frequency range train:',tr.frequency.min(),tr.frequency.max(),'val:',va.frequency.min(),va.frequency.max())
"

cd /app && cat > /tmp/feats.py << 'EOF'
import numpy as np, pandas as pd
FEATURES=["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]
NU=1.5e-5; C=340.0
def phys(X):
    f=pd.DataFrame(X,columns=FEATURES).astype(float).reset_index(drop=True)
    freq=f['frequency']; ang=f['attack_angle']; chord=f['chord_length']
    vel=f['free_stream_velocity']; thk=f['suction_side_displacement_thickness']
    o=pd.DataFrame(index=f.index)
    lf=np.log10(freq); lt=np.log10(thk); lv=np.log10(vel); lc=np.log10(chord)
    o['lf']=lf; o['lt']=lt; o['lv']=lv; o['lc']=lc
    o['ang']=ang
    o['logSt']=lf+lt-lv                 # Strouhal number
    o['logRe_c']=lv+lc                  # Reynolds (chord)
    o['logRe_d']=lv+lt                  # Reynolds (thickness)
    o['mach']=vel/C
    o['ang_sin']=np.sin(np.radians(ang))
    o['freq']=freq; o['chord']=chord; o['vel']=vel; o['thk']=thk
    o['logSt2']=o['logSt']**2
    o['lf_ang']=lf*ang
    o['lf_vel']=lf*vel
    return o
PHYS_COLS=None
EOF
python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld)
X=alld[FEATURES]; y=alld[TARGET]
from sklearn.ensemble import HistGradientBoostingRegressor, ExtraTreesRegressor, RandomForestRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer

def mk(est): return make_pipeline(FunctionTransformer(phys), est)
cv_eval(lambda: mk(HistGradientBoostingRegressor(random_state=0)), X,y,cond,label='HGB phys')
cv_eval(lambda: mk(ExtraTreesRegressor(n_estimators=400,random_state=0,n_jobs=-1)), X,y,cond,label='ET phys')
cv_eval(lambda: mk(RandomForestRegressor(n_estimators=400,random_state=0,n_jobs=-1)), X,y,cond,label='RF phys')
" 2>&1 | grep -v Warning

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld)
X=alld[FEATURES]; y=alld[TARGET]
gkf=GroupKFold(n_splits=8); oof=np.zeros(len(y))
for t,e in gkf.split(X,y,cond):
    m=make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=400,random_state=0,n_jobs=-1))
    m.fit(X.iloc[t],y.iloc[t]); oof[e]=m.predict(X.iloc[e])
pc=per_condition_rmse(y.values,oof,cond).sort_values(ascending=False)
print('Worst 10 conditions:'); print(pc.head(10))
# characterize worst conditions
worst=pc.head(6).index
for w in worst:
    sub=alld[cond==w]
    print(f'cond {w}: n={len(sub)} ang={sub.attack_angle.mean():.1f} chord={sub.chord_length.mean():.4f} vel={sub.free_stream_velocity.mean():.1f} thk={sub.suction_side_displacement_thickness.mean():.4f} SPLrange=[{sub[TARGET].min():.1f},{sub[TARGET].max():.1f}]')
" 2>&1 | grep -v Warning

cd /app && timeout 280 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel, Matern
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]

# GP
def mk_gp():
    k=ConstantKernel(50.0)*Matern(length_scale=np.ones(18),nu=1.5)+WhiteKernel(1.0)
    return make_pipeline(FunctionTransformer(phys),StandardScaler(),GaussianProcessRegressor(kernel=k,normalize_y=True,alpha=1e-6,n_restarts_optimizer=0))
cv_eval(mk_gp, X,y,cond,label='GP Matern1.5',nsplits=8)

# tuned ET
for msl in [1,2,3]:
    cv_eval(lambda msl=msl: make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=600,min_samples_leaf=msl,random_state=0,n_jobs=-1)), X,y,cond,label=f'ET msl={msl}')
# tuned HGB
cv_eval(lambda: make_pipeline(FunctionTransformer(phys),HistGradientBoostingRegressor(max_iter=500,learning_rate=0.05,max_leaf_nodes=31,l2_regularization=1.0,random_state=0)), X,y,cond,label='HGB tuned')
" 2>&1 | grep -v Warning

cd /app && timeout 290 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel, Matern
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
nf=phys(X).shape[1]; print('n phys feats',nf)

def mk_gp(nu):
    k=ConstantKernel(50.0)*Matern(length_scale=np.ones(nf),nu=nu)+WhiteKernel(1.0)
    return make_pipeline(FunctionTransformer(phys),StandardScaler(),GaussianProcessRegressor(kernel=k,normalize_y=True,alpha=1e-6,n_restarts_optimizer=0))
cv_eval(lambda: mk_gp(1.5), X,y,cond,label='GP Matern1.5',nsplits=8)
cv_eval(lambda: mk_gp(2.5), X,y,cond,label='GP Matern2.5',nsplits=8)
for msl in [1,2,3]:
    cv_eval(lambda msl=msl: make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=600,min_samples_leaf=msl,random_state=0,n_jobs=-1)), X,y,cond,label=f'ET msl={msl}')
cv_eval(lambda: make_pipeline(FunctionTransformer(phys),HistGradientBoostingRegressor(max_iter=500,learning_rate=0.05,max_leaf_nodes=31,l2_regularization=1.0,random_state=0)), X,y,cond,label='HGB tuned')
" 2>&1 | grep -v Warning

cd /app && timeout 200 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, WhiteKernel, Matern
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
nf=phys(X).shape[1]
t0=time.time()
def mk_gp(nu):
    k=ConstantKernel(50.0)*Matern(length_scale=np.ones(nf),nu=nu)+WhiteKernel(1.0)
    return make_pipeline(FunctionTransformer(phys),StandardScaler(),GaussianProcessRegressor(kernel=k,normalize_y=True,alpha=1e-6,n_restarts_optimizer=0))
cv_eval(lambda: mk_gp(1.5), X,y,cond,label='GP Matern1.5 (4fold)',nsplits=4)
print('time', time.time()-t0)
" 2>&1 | grep -v Warning

cd /app && timeout 250 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
for msl in [1,2,3,5]:
    cv_eval(lambda msl=msl: make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=600,min_samples_leaf=msl,random_state=0,n_jobs=-1)), X,y,cond,label=f'ET msl={msl}')
for lr,mi in [(0.05,400),(0.03,800),(0.08,300)]:
    cv_eval(lambda lr=lr,mi=mi: make_pipeline(FunctionTransformer(phys),HistGradientBoostingRegressor(max_iter=mi,learning_rate=lr,max_leaf_nodes=31,l2_regularization=1.0,random_state=0)), X,y,cond,label=f'HGB lr={lr} mi={mi}')
" 2>&1 | grep -v Warning

cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.neural_network import MLPRegressor
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]

def mk_mlp(hls,alpha):
    return make_pipeline(FunctionTransformer(phys),StandardScaler(),
        MLPRegressor(hidden_layer_sizes=hls,activation='tanh',alpha=alpha,max_iter=2000,random_state=0,early_stopping=False))
cv_eval(lambda: mk_mlp((64,64),1e-3), X,y,cond,label='MLP 64x64 a=1e-3')
cv_eval(lambda: mk_mlp((128,64),1e-2), X,y,cond,label='MLP 128x64 a=1e-2')
cv_eval(lambda: mk_mlp((100,100,50),1e-3), X,y,cond,label='MLP 100x100x50')
" 2>&1 | grep -v Warning

cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor, RandomForestRegressor
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]

def models():
    return {
      'et': make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=600,min_samples_leaf=1,random_state=0,n_jobs=-1)),
      'hgb': make_pipeline(FunctionTransformer(phys),HistGradientBoostingRegressor(max_iter=400,learning_rate=0.05,max_leaf_nodes=31,l2_regularization=1.0,random_state=0)),
      'rf': make_pipeline(FunctionTransformer(phys),RandomForestRegressor(n_estimators=500,min_samples_leaf=1,random_state=0,n_jobs=-1)),
    }
gkf=GroupKFold(n_splits=8)
preds={k:np.zeros(len(y)) for k in models()}
for t,e in gkf.split(X,y,cond):
    ms=models()
    for k,m in ms.items():
        m.fit(X.iloc[t],y.iloc[t]); preds[k][e]=m.predict(X.iloc[e])
for k in preds: summarize(y.values,preds[k],cond,label=k)
import itertools
for combo in [('et','hgb'),('et','rf'),('et','hgb','rf'),('hgb','rf')]:
    avg=np.mean([preds[k] for k in combo],axis=0)
    summarize(y.values,avg,cond,label='avg:'+'+'.join(combo))
# weighted et heavy
for w in [0.6,0.7]:
    avg=w*preds['et']+(1-w)*preds['hgb']
    summarize(y.values,avg,cond,label=f'et{w}+hgb{1-w:.1f}')
" 2>&1 | grep -v Warning

cd /app && python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
import numpy as np, pandas as pd
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); alld['cond']=cond
alld['lf']=np.log10(alld['frequency'])
resids=[]
for c,d in alld.groupby('cond'):
    if len(d)<6: continue
    d=d.sort_values('lf')
    # fit poly deg 3 in lf
    for deg in [3]:
        coef=np.polyfit(d['lf'],d[TARGET],deg)
        pred=np.polyval(coef,d['lf'])
        resids.append((d[TARGET]-pred).values)
r=np.concatenate(resids)
print(f'Within-condition poly3 residual: std={r.std():.3f} mae={np.mean(np.abs(r)):.3f} max|r|={np.max(np.abs(r)):.3f}')
# deg4
resids=[]
for c,d in alld.groupby('cond'):
    if len(d)<7: continue
    d=d.sort_values('lf')
    coef=np.polyfit(d['lf'],d[TARGET],4)
    pred=np.polyval(coef,d['lf'])
    resids.append((d[TARGET]-pred).values)
r=np.concatenate(resids)
print(f'Within-condition poly4 residual: std={r.std():.3f} mae={np.mean(np.abs(r)):.3f}')
print('n rows per condition: min/med/max', alld.groupby('cond').size().min(), int(alld.groupby('cond').size().median()), alld.groupby('cond').size().max())
"

cd /app && timeout 200 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
tr,va=load()
condv=recover_conditions(va)
Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
print('val conditions:', len(set(condv)))
et=make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=600,min_samples_leaf=1,random_state=0,n_jobs=-1))
hgb=make_pipeline(FunctionTransformer(phys),HistGradientBoostingRegressor(max_iter=400,learning_rate=0.05,max_leaf_nodes=31,l2_regularization=1.0,random_state=0))
et.fit(Xtr,ytr); hgb.fit(Xtr,ytr)
pe=et.predict(Xva); ph=hgb.predict(Xva)
summarize(yva.values,pe,condv,label='ET val')
summarize(yva.values,ph,condv,label='HGB val')
summarize(yva.values,0.7*pe+0.3*ph,condv,label='0.7ET+0.3HGB val')
summarize(yva.values,0.5*pe+0.5*ph,condv,label='0.5ET+0.5HGB val')
" 2>&1 | grep -v Warning

cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.kernel_ridge import KernelRidge
from sklearn.neighbors import KNeighborsRegressor
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]

for g in [0.05,0.1,0.2,0.3]:
  for a in [0.1,0.5,1.0]:
    cv_eval(lambda g=g,a=a: make_pipeline(FunctionTransformer(phys),StandardScaler(),KernelRidge(kernel='rbf',gamma=g,alpha=a)), X,y,cond,label=f'KRR g={g} a={a}')
" 2>&1 | grep -v Warning

cd /app && timeout 250 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
gkf=GroupKFold(n_splits=8); oof=np.zeros(len(y))
for t,e in gkf.split(X,y,cond):
    m=make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=400,random_state=0,n_jobs=-1))
    m.fit(X.iloc[t],y.iloc[t]); oof[e]=m.predict(X.iloc[e])
res=y.values-oof
d=alld.copy(); d['res']=res; d['lf']=np.log10(d['frequency'])
# residual vs features (binned)
for col in ['lf','attack_angle','chord_length','free_stream_velocity','suction_side_displacement_thickness']:
    q=pd.qcut(d[col],6,duplicates='drop')
    g=d.groupby(q,observed=True)['res'].agg(['mean','std','count'])
    print('---',col); print(g.round(3))
print('overall res mean',res.mean().round(3),'std',res.std().round(3))
# feature importance
m=make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=400,random_state=0,n_jobs=-1)); m.fit(X,y)
imp=m.named_steps['extratreesregressor'].feature_importances_
names=list(phys(X).columns)
for n,i in sorted(zip(names,imp),key=lambda t:-t[1]): print(f'  {n:10s} {i:.3f}')
" 2>&1 | grep -v Warning

cd /app && timeout 260 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.kernel_ridge import KernelRidge
from sklearn.compose import TransformedTargetRegressor
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
ymean=y.mean()
def ctr(est): return TransformedTargetRegressor(regressor=est, func=lambda y:y-ymean, inverse_func=lambda y:y+ymean)
for g in [0.02,0.05,0.1]:
  for a in [0.1,0.3,1.0]:
    cv_eval(lambda g=g,a=a: ctr(make_pipeline(FunctionTransformer(phys),StandardScaler(),KernelRidge(kernel='rbf',gamma=g,alpha=a))), X,y,cond,label=f'KRRc g={g} a={a}')
" 2>&1 | grep -v Warning

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

def fit_predict_val(make_model, Xtr,ytr,Xva):
    m=make_model(); m.fit(Xtr,ytr); return m.predict(Xva)
EOF
timeout 280 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]

for mf in [0.5,0.7,1.0,'sqrt']:
    lab=f'ET mf={mf}'
    r=cv_eval(lambda mf=mf: make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=500,max_features=mf,random_state=0,n_jobs=-1)), X,y,cond,label='CV '+lab)
    m=make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=500,max_features=mf,random_state=0,n_jobs=-1)); m.fit(Xtr,ytr)
    summarize(yva.values,m.predict(Xva),condv,label='VAL '+lab)
    print()
" 2>&1 | grep -v Warning

cd /app && timeout 280 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.neighbors import KNeighborsRegressor
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
# weighted phys: scale some features up
class WPhys(FunctionTransformer):
    pass
def wphys(X):
    o=phys(X)
    return o
for k in [5,8,12,20]:
  for w in ['distance']:
    cv_eval(lambda k=k,w=w: make_pipeline(FunctionTransformer(phys),StandardScaler(),KNeighborsRegressor(n_neighbors=k,weights=w)), X,y,cond,label=f'KNN k={k} {w}')
" 2>&1 | grep -v Warning

cd /app && timeout 290 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler, SplineTransformer
from sklearn.ensemble import HistGradientBoostingRegressor, ExtraTreesRegressor
from sklearn.linear_model import Ridge
from sklearn.base import BaseEstimator, RegressorMixin, clone
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]

# HGB absolute loss
for mi,lr in [(600,0.03),(800,0.02)]:
    cv_eval(lambda mi=mi,lr=lr: make_pipeline(FunctionTransformer(phys),HistGradientBoostingRegressor(loss='absolute_error',max_iter=mi,learning_rate=lr,max_leaf_nodes=31,l2_regularization=1.0,random_state=0)), X,y,cond,label=f'HGB-abs mi={mi} lr={lr}')

# Hybrid: smooth spline ridge + ET residual
class Hybrid(BaseEstimator,RegressorMixin):
    def __init__(self, stage1, stage2): self.stage1=stage1; self.stage2=stage2
    def fit(self,X,y):
        self.s1=clone(self.stage1).fit(X,y); r=y-self.s1.predict(X)
        self.s2=clone(self.stage2).fit(X,r); return self
    def predict(self,X): return self.s1.predict(X)+self.s2.predict(X)
def smooth():
    return make_pipeline(FunctionTransformer(phys),SplineTransformer(n_knots=6,degree=3),Ridge(alpha=1.0))
cv_eval(smooth, X,y,cond,label='Smooth spline ridge')
et=lambda: make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=500,random_state=0,n_jobs=-1))
cv_eval(lambda: Hybrid(smooth(), et()), X,y,cond,label='Hybrid smooth+ET')
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/hybrid.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.base import BaseEstimator, RegressorMixin, clone
class Hybrid(BaseEstimator,RegressorMixin):
    def __init__(self, stage1, stage2): self.stage1=stage1; self.stage2=stage2
    def fit(self,X,y):
        self.s1=clone(self.stage1).fit(X,y); r=y-self.s1.predict(X)
        self.s2=clone(self.stage2).fit(X,r); return self
    def predict(self,X): return self.s1.predict(X)+self.s2.predict(X)
EOF
timeout 290 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from hybrid import Hybrid
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.linear_model import Ridge
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def smooth(nk,al):
    return make_pipeline(FunctionTransformer(phys),SplineTransformer(n_knots=nk,degree=3),Ridge(alpha=al))
et=lambda: make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=500,random_state=0,n_jobs=-1))
for nk,al in [(5,1.0),(6,1.0),(8,1.0),(6,0.3),(6,3.0),(10,1.0)]:
    mk=lambda nk=nk,al=al: Hybrid(smooth(nk,al), et())
    cv_eval(mk, X,y,cond,label=f'CV Hyb nk={nk} al={al}')
    m=mk(); m.fit(Xtr,ytr); summarize(yva.values,m.predict(Xva),condv,label=f'VAL Hyb nk={nk} al={al}')
    print()
" 2>&1 | grep -v Warning

cd /app && timeout 200 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from hybrid import Hybrid
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
def smooth(): return make_pipeline(FunctionTransformer(phys),SplineTransformer(n_knots=6,degree=3),Ridge(alpha=0.3))
et=lambda: make_pipeline(FunctionTransformer(phys),ExtraTreesRegressor(n_estimators=500,random_state=0,n_jobs=-1))
gkf=GroupKFold(n_splits=8); oof=np.zeros(len(y))
for t,e in gkf.split(X,y,cond):
    m=Hybrid(smooth(),et()); m.fit(X.iloc[t],y.iloc[t]); oof[e]=m.predict(X.iloc[e])
pc=per_condition_rmse(y.values,oof,cond).sort_values()
print('percentiles of per-cond RMSE: p50=%.2f p75=%.2f p90=%.2f p95=%.2f max=%.2f'%(pc.median(),pc.quantile(.75),pc.quantile(.9),pc.quantile(.95),pc.max()))
print('frac conds <=2.35:', (pc<=2.35).mean())
print('Worst 8:'); 
for c,v in pc.tail(8).items():
    sub=alld[cond==c]; print(f'  {c}: RMSE={v:.2f} ang={sub.attack_angle.mean():.1f} chord={sub.chord_length.mean():.4f} vel={sub.free_stream_velocity.mean():.0f}')
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/feats2.py << 'EOF'
import numpy as np, pandas as pd
FEATURES=["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]
NU=1.5e-5; C=340.0
def phys2(X):
    f=pd.DataFrame(X,columns=FEATURES).astype(float).reset_index(drop=True)
    freq=f['frequency']; ang=f['attack_angle']; chord=f['chord_length']
    vel=f['free_stream_velocity']; thk=f['suction_side_displacement_thickness']
    o=pd.DataFrame(index=f.index)
    lf=np.log10(freq); lt=np.log10(thk); lv=np.log10(vel); lc=np.log10(chord)
    logSt=lf+lt-lv
    o['lf']=lf; o['lt']=lt; o['lv']=lv; o['lc']=lc; o['ang']=ang
    o['logSt']=logSt; o['logRe_c']=lv+lc; o['logRe_d']=lv+lt; o['mach']=vel/C
    o['ang_sin']=np.sin(np.radians(ang))
    o['logSt_ang']=logSt*ang
    o['logSt_lv']=logSt*lv
    o['logSt_lc']=logSt*lc
    o['lf_ang']=lf*ang
    o['ang2']=ang**2
    return o
EOF
timeout 295 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats2 import phys2
from hybrid import Hybrid
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.linear_model import Ridge
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def smooth(al=0.3,nk=6): return make_pipeline(FunctionTransformer(phys2),SplineTransformer(n_knots=nk,degree=3),Ridge(alpha=al))
et=lambda: make_pipeline(FunctionTransformer(phys2),ExtraTreesRegressor(n_estimators=500,random_state=0,n_jobs=-1))
# Hybrid with enriched feats
for al in [0.3,1.0]:
    mk=lambda al=al: Hybrid(smooth(al),et())
    cv_eval(mk,X,y,cond,label=f'CV Hyb2 al={al}')
    m=mk(); m.fit(Xtr,ytr); summarize(yva.values,m.predict(Xva),condv,label=f'VAL Hyb2 al={al}'); print()
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/blend.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from feats2 import phys2
from hybrid import Hybrid
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.linear_model import Ridge
from sklearn.model_selection import GroupKFold

tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]

def sm(feat,al,nk=6): return make_pipeline(FunctionTransformer(feat),SplineTransformer(n_knots=nk,degree=3),Ridge(alpha=al))
def et(feat,seed=0,msl=1): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=500,min_samples_leaf=msl,random_state=seed,n_jobs=-1))
def hgb(feat): return make_pipeline(FunctionTransformer(feat),HistGradientBoostingRegressor(loss='absolute_error',max_iter=600,learning_rate=0.03,max_leaf_nodes=31,l2_regularization=1.0,random_state=0))

MODELS={
 'A': lambda: Hybrid(sm(phys2,0.3), et(phys2,0)),
 'B': lambda: Hybrid(sm(phys,0.3), et(phys,0)),
 'C': lambda: et(phys,0),
 'D': lambda: Hybrid(sm(phys2,0.3), et(phys,0)),
 'E': lambda: Hybrid(sm(phys2,0.3), hgb(phys2)),
 'F': lambda: Hybrid(sm(phys,1.0), et(phys2,0)),
}

def compute():
    gkf=GroupKFold(n_splits=8)
    oof={k:np.zeros(len(y)) for k in MODELS}
    valp={k:np.zeros(len(yva)) for k in MODELS}
    for k,mk in MODELS.items():
        for t,e in gkf.split(X,y,cond):
            m=mk(); m.fit(X.iloc[t],y.iloc[t]); oof[k][e]=m.predict(X.iloc[e])
        m=mk(); m.fit(Xtr,ytr); valp[k]=m.predict(Xva)
    return oof,valp

if __name__=='__main__':
    oof,valp=compute()
    np.savez('/tmp/preds.npz', y=y.values, cond=cond, yva=yva.values, condv=condv,
             **{f'oof_{k}':oof[k] for k in MODELS}, **{f'val_{k}':valp[k] for k in MODELS})
    for k in MODELS:
        summarize(y.values,oof[k],cond,label=f'CV {k}')
        summarize(yva.values,valp[k],condv,label=f'VAL {k}')
        print()
EOF
timeout 295 python3 /tmp/blend.py 2>&1 | grep -v Warning

cd /app && timeout 250 python3 -c "
import numpy as np, itertools, sys; sys.path.insert(0,'/tmp')
from harness import per_condition_rmse
d=np.load('/tmp/preds.npz',allow_pickle=True)
keys=['A','B','C','D','E','F']
y=d['y']; cond=d['cond']; yva=d['yva']; condv=d['condv']
oof=np.vstack([d[f'oof_{k}'] for k in keys])
val=np.vstack([d[f'val_{k}'] for k in keys])
def metr(pred,yy,cc):
    rmse=np.sqrt(np.mean((yy-pred)**2)); mae=np.mean(np.abs(yy-pred))
    pc=per_condition_rmse(yy,pred,cc)
    return rmse,mae,np.percentile(pc,90),np.percentile(pc,95),pc.max()
# search simplex weights in steps
best=[]
grid=np.arange(0,1.01,0.2)
import itertools
for w in itertools.product(grid,repeat=len(keys)):
    s=sum(w)
    if abs(s-1)>1e-6: continue
    w=np.array(w)
    po=w@oof
    r,m,p90,p95,mx=metr(po,y,cond)
    # objective: must-pass margins, combine p90 and mae
    score=max(r/2.25,m/1.55,p90/2.35,p95/4.75,mx/5.60)
    best.append((score,w,r,m,p90,p95,mx))
best.sort(key=lambda t:t[0])
print('Top CV blends (by worst-normalized-metric):')
for score,w,r,m,p90,p95,mx in best[:8]:
    pv=w@val
    rv,mv,p90v,p95v,mxv=metr(pv,yva,condv)
    ws=','.join(f'{k}:{ww:.1f}' for k,ww in zip(keys,w) if ww>0)
    print(f'  [{ws}] CV r={r:.2f} m={m:.2f} p90={p90:.2f} p95={p95:.2f} mx={mx:.2f} score={score:.3f} || VAL r={rv:.2f} m={mv:.2f} p90={p90v:.2f} mx={mxv:.2f}')
" 2>&1 | grep -v Warning

cd /app && timeout 200 python3 -c "
import numpy as np, sys; sys.path.insert(0,'/tmp')
from harness import per_condition_rmse, load, recover_conditions, FEATURES, TARGET
d=np.load('/tmp/preds.npz',allow_pickle=True)
keys=['A','B','C','D','E','F']
y=d['y']; cond=d['cond']
oof=np.vstack([d[f'oof_{k}'] for k in keys])
# best blend weights from search
w=np.array([0,0.4,0.2,0,0.2,0.2]) # B,C,E,F
po=w@oof
tr,va=load(); alld=__import__('pandas').concat([tr,va],ignore_index=True)
lf=np.log10(alld['frequency'].values)
def metr(pred):
    rmse=np.sqrt(np.mean((y-pred)**2)); mae=np.mean(np.abs(y-pred))
    pc=per_condition_rmse(y,pred,cond); 
    return f'RMSE={rmse:.3f} MAE={mae:.3f} p90={np.percentile(pc,90):.3f} p95={np.percentile(pc,95):.3f} max={pc.max():.3f} frac<=2.35={np.mean(pc<=2.35):.2f}'
print('raw blend :', metr(po))
import pandas as pd
for deg in [2,3,4]:
    sm=po.copy()
    df=pd.DataFrame({'lf':lf,'p':po,'c':cond})
    for c,idx in df.groupby('c').groups.items():
        idx=np.array(idx); 
        if len(idx)<deg+2: continue
        x=lf[idx]; coef=np.polyfit(x,po[idx],deg); sm[idx]=np.polyval(coef,x)
    print(f'smooth d{deg}:', metr(sm))
" 2>&1 | grep -v Warning

cd /app && timeout 120 python3 -c "
import numpy as np, pandas as pd, sys; sys.path.insert(0,'/tmp')
from harness import load, recover_conditions, TARGET
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); alld['cond']=cond
alld['lf']=np.log10(alld['frequency'])
# model-free noise: second difference along sorted freq; Var(2nd diff)=6*sigma^2 if smooth
sig=[]; 
perc=[]
for c,d in alld.groupby('cond'):
    if len(d)<5: continue
    d=d.sort_values('lf'); v=d[TARGET].values
    dd=v[2:]-2*v[1:-1]+v[:-2]
    s=np.sqrt(np.mean(dd**2)/6.0)
    sig.append(s)
sig=np.array(sig)
print(f'Per-condition noise sigma (2nd-diff): median={np.median(sig):.3f} mean={sig.mean():.3f} p90={np.percentile(sig,90):.3f} max={sig.max():.3f}')
print(f'Implied best-case per-cond RMSE ~ sigma. median noise->RMSE~{np.median(sig):.2f}')
# correlate noise with angle
info=[]
for c,d in alld.groupby('cond'):
    if len(d)<5: continue
    dd=d.sort_values('lf')[TARGET].values; dd=dd[2:]-2*dd[1:-1]+dd[:-2]
    info.append((d['attack_angle'].mean(), d['chord_length'].mean(), np.sqrt(np.mean(dd**2)/6.0)))
info=pd.DataFrame(info,columns=['ang','chord','sig'])
print(info.sort_values('sig',ascending=False).head(8).round(3).to_string(index=False))
"

cd /app && cat > /tmp/hybrid2.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.base import BaseEstimator, RegressorMixin, clone
class HybridBag(BaseEstimator,RegressorMixin):
    '''stage1 smooth; stage2 = average of several ET with different seeds.'''
    def __init__(self, stage1, stage2_factory, n_seed=5):
        self.stage1=stage1; self.stage2_factory=stage2_factory; self.n_seed=n_seed
    def fit(self,X,y):
        self.s1=clone(self.stage1).fit(X,y); r=y-self.s1.predict(X)
        self.s2=[self.stage2_factory(s).fit(X,r) for s in range(self.n_seed)]
        return self
    def predict(self,X):
        return self.s1.predict(X)+np.mean([m.predict(X) for m in self.s2],axis=0)
EOF
timeout 295 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from feats2 import phys2
from hybrid2 import HybridBag
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def sm(feat,al): return make_pipeline(FunctionTransformer(feat),SplineTransformer(n_knots=6,degree=3),Ridge(alpha=al))
def etf(feat): return lambda s: make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=300,random_state=s,n_jobs=-1))
for feat,fn,al in [(phys2,'phys2',0.3),(phys,'phys',0.3)]:
    mk=lambda feat=feat,al=al: HybridBag(sm(feat,al), etf(feat), n_seed=6)
    cv_eval(mk,X,y,cond,label=f'CV HybBag {fn}')
    m=mk(); m.fit(Xtr,ytr); summarize(yva.values,m.predict(Xva),condv,label=f'VAL HybBag {fn}'); print()
" 2>&1 | grep -v Warning

cd /app && timeout 120 python3 -c "
import numpy as np, pandas as pd, sys; sys.path.insert(0,'/tmp')
from harness import per_condition_rmse, load, recover_conditions, TARGET
d=np.load('/tmp/preds.npz',allow_pickle=True)
keys=['A','B','C','D','E','F']; y=d['y']; cond=d['cond']
oof=np.vstack([d[f'oof_{k}'] for k in keys])
w=np.array([0,0.4,0.2,0,0.2,0.2]); po=w@oof
res=y-po
df=pd.DataFrame({'res':res,'c':cond})
g=df.groupby('c')['res']
bias=g.mean(); spread=g.std(); rmse=np.sqrt(g.apply(lambda v:np.mean(v**2)))
summ=pd.DataFrame({'rmse':rmse,'bias':bias,'spread':spread})
summ['bias_frac']=summ['bias']**2/summ['rmse']**2
hard=summ.sort_values('rmse',ascending=False)
print('For worst 12 conditions: fraction of MSE from bias(level offset):')
print(hard.head(12).round(2).to_string())
print()
print('Avg bias_frac over all conds: %.2f'%summ['bias_frac'].mean())
print('Avg bias_frac over worst 20: %.2f'%hard.head(20)['bias_frac'].mean())
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/feats_gp.py << 'EOF'
import numpy as np, pandas as pd
FEATURES=["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]
def gpfeats(X):
    f=pd.DataFrame(X,columns=FEATURES).astype(float).reset_index(drop=True)
    lf=np.log10(f['frequency']); lt=np.log10(f['suction_side_displacement_thickness'])
    lv=np.log10(f['free_stream_velocity']); lc=np.log10(f['chord_length'])
    o=pd.DataFrame(index=f.index)
    o['logSt']=lf+lt-lv; o['lf']=lf; o['ang']=f['attack_angle']; o['lv']=lv; o['lc']=lc; o['lt']=lt
    return o.to_numpy(float)
EOF
timeout 200 python3 -c "
import sys,time; sys.path.insert(0,'/tmp')
from harness import *
from feats_gp import gpfeats
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel as C, WhiteKernel, Matern, DotProduct
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
nf=gpfeats(X).shape[1]
def mkgp():
    k=C(100.0)*Matern(length_scale=np.ones(nf),nu=2.5)+DotProduct(sigma_0=1.0)+WhiteKernel(1.0,(1e-2,10))
    return make_pipeline(FunctionTransformer(gpfeats),StandardScaler(),GaussianProcessRegressor(kernel=k,normalize_y=True,alpha=1e-8,n_restarts_optimizer=0))
t0=time.time()
cv_eval(mkgp, X,y,cond,label='GP6feat',nsplits=6)
print('time',round(time.time()-t0,1))
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/feats3.py << 'EOF'
import numpy as np, pandas as pd
FEATURES=["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]
C=340.0
def phys3(X):
    f=pd.DataFrame(X,columns=FEATURES).astype(float).reset_index(drop=True)
    freq=f['frequency']; ang=f['attack_angle']; chord=f['chord_length']
    vel=f['free_stream_velocity']; thk=f['suction_side_displacement_thickness']
    lf=np.log10(freq); lt=np.log10(thk); lv=np.log10(vel); lc=np.log10(chord)
    logSt=lf+lt-lv
    o=pd.DataFrame(index=f.index)
    o['lf']=lf; o['lt']=lt; o['lv']=lv; o['lc']=lc; o['ang']=ang
    o['logSt']=logSt; o['logRe_c']=lv+lc; o['logRe_d']=lv+lt; o['mach']=vel/C
    o['ang_sin']=np.sin(np.radians(ang))
    o['logSt_ang']=logSt*ang
    o['logSt_lc']=logSt*lc
    o['lf_ang']=lf*ang
    o['ang2']=ang**2
    o['lt_ang']=lt*ang
    o['thk_chord']=thk/chord          # bl thickness ratio
    o['logStc']=lf+lc-lv              # chord-based Strouhal
    return o
EOF
timeout 295 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats3 import phys3
from feats2 import phys2
from hybrid import Hybrid
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def sm(feat,al,nk=6): return make_pipeline(FunctionTransformer(feat),SplineTransformer(n_knots=nk,degree=3),Ridge(alpha=al))
def et(feat,msl=1): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=500,min_samples_leaf=msl,random_state=0,n_jobs=-1))
for al in [0.3,0.6]:
  for msl in [1,2]:
    mk=lambda al=al,msl=msl: Hybrid(sm(phys3,al),et(phys3,msl))
    cv_eval(mk,X,y,cond,label=f'CV Hyb3 al={al} msl={msl}')
    m=mk(); m.fit(Xtr,ytr); summarize(yva.values,m.predict(Xva),condv,label=f'VAL Hyb3 al={al} msl={msl}'); print()
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/blend2.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from feats2 import phys2
from feats3 import phys3
from hybrid import Hybrid
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
from sklearn.model_selection import GroupKFold

tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
NE=300
def sm(feat,al,nk=6): return make_pipeline(FunctionTransformer(feat),SplineTransformer(n_knots=nk,degree=3),Ridge(alpha=al))
def et(feat,msl=1): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=NE,min_samples_leaf=msl,random_state=0,n_jobs=-1))
MODELS={
 'M1': lambda: Hybrid(sm(phys3,0.3), et(phys3)),
 'M2': lambda: Hybrid(sm(phys2,0.3), et(phys2)),
 'M3': lambda: Hybrid(sm(phys,0.3), et(phys)),
 'M4': lambda: et(phys3),
 'M5': lambda: Hybrid(sm(phys3,0.6), et(phys3,2)),
 'M6': lambda: Hybrid(sm(phys2,0.3), et(phys3)),
}
gkf=GroupKFold(n_splits=8)
oof={k:np.zeros(len(y)) for k in MODELS}; valp={k:np.zeros(len(yva)) for k in MODELS}
for k,mk in MODELS.items():
    for t,e in gkf.split(X,y,cond):
        m=mk(); m.fit(X.iloc[t],y.iloc[t]); oof[k][e]=m.predict(X.iloc[e])
    m=mk(); m.fit(Xtr,ytr); valp[k]=m.predict(Xva)
np.savez('/tmp/preds2.npz', y=y.values, cond=cond, yva=yva.values, condv=condv,
         keys=list(MODELS.keys()),
         **{f'oof_{k}':oof[k] for k in MODELS}, **{f'val_{k}':valp[k] for k in MODELS})
for k in MODELS:
    r=summarize(y.values,oof[k],cond,label=f'CV {k}')
print('saved')
EOF
timeout 298 python3 /tmp/blend2.py 2>&1 | grep -v Warning

cd /app && timeout 250 python3 -c "
import numpy as np, itertools, sys; sys.path.insert(0,'/tmp')
from harness import per_condition_rmse
d=np.load('/tmp/preds2.npz',allow_pickle=True)
keys=list(d['keys'])
y=d['y']; cond=d['cond']; yva=d['yva']; condv=d['condv']
oof=np.vstack([d[f'oof_{k}'] for k in keys]); val=np.vstack([d[f'val_{k}'] for k in keys])
T=dict(r=2.25,m=1.55,p90=2.35,p95=4.75,mx=5.60)
def metr(pred,yy,cc):
    rmse=np.sqrt(np.mean((yy-pred)**2)); mae=np.mean(np.abs(yy-pred)); pc=per_condition_rmse(yy,pred,cc)
    return dict(r=rmse,m=mae,p90=np.percentile(pc,90),p95=np.percentile(pc,95),mx=pc.max())
def score(md): return max(md['r']/T['r'],md['m']/T['m'],md['p90']/T['p90'],md['p95']/T['p95'],md['mx']/T['mx'])
res=[]
grid=np.arange(0,1.001,0.1)
for w in itertools.product(grid,repeat=len(keys)):
    if abs(sum(w)-1)>1e-6: continue
    w=np.array(w)
    mc=metr(w@oof,y,cond); mv=metr(w@val,yva,condv)
    # combined objective: emphasize CV (more conditions) but require val ok too
    comb=0.6*score(mc)+0.4*score(mv)
    res.append((comb,score(mc),score(mv),w,mc,mv))
res.sort(key=lambda t:t[0])
for comb,sc,sv,w,mc,mv in res[:6]:
    ws=','.join(f'{k}:{ww:.1f}' for k,ww in zip(keys,w) if ww>0)
    print(f'[{ws}]')
    print(f'   CV  r={mc[\"r\"]:.2f} m={mc[\"m\"]:.2f} p90={mc[\"p90\"]:.2f} p95={mc[\"p95\"]:.2f} mx={mc[\"mx\"]:.2f} (sc={sc:.2f})')
    print(f'   VAL r={mv[\"r\"]:.2f} m={mv[\"m\"]:.2f} p90={mv[\"p90\"]:.2f} p95={mv[\"p95\"]:.2f} mx={mv[\"mx\"]:.2f} (sc={sv:.2f})')
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/blend3.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from feats2 import phys2
from feats3 import phys3
from hybrid import Hybrid
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor, RandomForestRegressor, GradientBoostingRegressor
from sklearn.linear_model import Ridge
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def sm(feat,al,nk=6): return make_pipeline(FunctionTransformer(feat),SplineTransformer(n_knots=nk,degree=3),Ridge(alpha=al))
def rf(feat): return make_pipeline(FunctionTransformer(feat),RandomForestRegressor(n_estimators=400,min_samples_leaf=1,random_state=0,n_jobs=-1))
def gb(feat): return make_pipeline(FunctionTransformer(feat),GradientBoostingRegressor(n_estimators=400,learning_rate=0.03,max_depth=3,subsample=0.8,random_state=0))
def et(feat): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=300,random_state=1,n_jobs=-1))
MODELS={
 'M7': lambda: Hybrid(sm(phys3,0.3), rf(phys3)),
 'M8': lambda: Hybrid(sm(phys3,0.3), gb(phys3)),
 'M9': lambda: et(phys2),
}
gkf=GroupKFold(n_splits=8)
oof={k:np.zeros(len(y)) for k in MODELS}; valp={k:np.zeros(len(yva)) for k in MODELS}
for k,mk in MODELS.items():
    for t,e in gkf.split(X,y,cond):
        m=mk(); m.fit(X.iloc[t],y.iloc[t]); oof[k][e]=m.predict(X.iloc[e])
    m=mk(); m.fit(Xtr,ytr); valp[k]=m.predict(Xva)
    summarize(y.values,oof[k],cond,label=f'CV {k}')
np.savez('/tmp/preds3.npz', keys=list(MODELS.keys()),
         **{f'oof_{k}':oof[k] for k in MODELS}, **{f'val_{k}':valp[k] for k in MODELS})
print('saved')
EOF
timeout 298 python3 /tmp/blend3.py 2>&1 | grep -v Warning

cd /app && timeout 260 python3 -c "
import numpy as np, itertools, sys; sys.path.insert(0,'/tmp')
from harness import per_condition_rmse
d2=np.load('/tmp/preds2.npz',allow_pickle=True); d3=np.load('/tmp/preds3.npz',allow_pickle=True)
y=d2['y']; cond=d2['cond']; yva=d2['yva']; condv=d2['condv']
keys=list(d2['keys'])+list(d3['keys'])
def get(k,pfx): return d2[f'{pfx}_{k}'] if f'{pfx}_{k}' in d2 else d3[f'{pfx}_{k}']
oof=np.vstack([get(k,'oof') for k in keys]); val=np.vstack([get(k,'val') for k in keys])
T=dict(r=2.25,m=1.55,p90=2.35,p95=4.75,mx=5.60)
def metr(pred,yy,cc):
    pc=per_condition_rmse(yy,pred,cc)
    return dict(r=np.sqrt(np.mean((yy-pred)**2)),m=np.mean(np.abs(yy-pred)),p90=np.percentile(pc,90),p95=np.percentile(pc,95),mx=pc.max())
# coordinate-descent weight optimization minimizing CV p90 + 0.5*val p90 with penalties for other metrics
def obj(mc,mv):
    pen=0
    for kk in T:
        pen+=max(0,mc[kk]/T[kk]-1)*3 + max(0,mv[kk]/T[kk]-1)*2
    return (mc['p90']+mv['p90']) + 10*pen
w=np.ones(len(keys))/len(keys)
import random; random.seed(0)
best=(1e9,None)
for it in range(4000):
    ww=np.random.dirichlet(np.ones(len(keys))*0.5)
    mc=metr(ww@oof,y,cond); mv=metr(ww@val,yva,condv)
    o=obj(mc,mv)
    if o<best[0]: best=(o,ww,mc,mv)
o,ww,mc,mv=best
ws=','.join(f'{k}:{w:.2f}' for k,w in zip(keys,ww) if w>0.02)
print('Best for p90:',ws)
print(f'  CV  r={mc[\"r\"]:.3f} m={mc[\"m\"]:.3f} p90={mc[\"p90\"]:.3f} p95={mc[\"p95\"]:.3f} mx={mc[\"mx\"]:.3f}')
print(f'  VAL r={mv[\"r\"]:.3f} m={mv[\"m\"]:.3f} p90={mv[\"p90\"]:.3f} p95={mv[\"p95\"]:.3f} mx={mv[\"mx\"]:.3f}')
# also equal-weight all
ww=np.ones(len(keys))/len(keys); mc=metr(ww@oof,y,cond); mv=metr(ww@val,yva,condv)
print('Equal weight:')
print(f'  CV  r={mc[\"r\"]:.3f} m={mc[\"m\"]:.3f} p90={mc[\"p90\"]:.3f} p95={mc[\"p95\"]:.3f} mx={mc[\"mx\"]:.3f}')
print(f'  VAL r={mv[\"r\"]:.3f} m={mv[\"m\"]:.3f} p90={mv[\"p90\"]:.3f} p95={mv[\"p95\"]:.3f} mx={mv[\"mx\"]:.3f}')
" 2>&1 | grep -v Warning

cd /app && timeout 200 python3 -c "
import numpy as np, pandas as pd, sys; sys.path.insert(0,'/tmp')
from harness import *
from feats3 import phys3
import numpy as np
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); alld['cond']=cond; alld['lf']=np.log10(alld['frequency'])
# cleaner noise: fit local quadratic (Savitzky-Golay-like) in window, residual
from numpy.polynomial import polynomial as P
res=[]
for c,d in alld.groupby('cond'):
    d=d.sort_values('lf'); v=d[TARGET].values; x=d['lf'].values; n=len(v)
    if n<7: continue
    for i in range(n):
        lo=max(0,i-3); hi=min(n,i+4)
        xx=x[lo:hi]-x[i]; yy=v[lo:hi]
        deg=min(3,len(xx)-1)
        coef=np.polyfit(xx,yy,deg)
        res.append(v[i]-np.polyval(coef,0.0))
res=np.array(res)
print(f'Local-cubic residual noise: std={res.std():.3f} mae={np.mean(np.abs(res)):.3f}')
# per-condition noise via local fit
pcn=[]
for c,d in alld.groupby('cond'):
    d=d.sort_values('lf'); v=d[TARGET].values; x=d['lf'].values; n=len(v)
    if n<7: continue
    r=[]
    for i in range(n):
        lo=max(0,i-3); hi=min(n,i+4); xx=x[lo:hi]-x[i]; yy=v[lo:hi]
        coef=np.polyfit(xx,yy,min(3,len(xx)-1)); r.append(v[i]-np.polyval(coef,0.0))
    pcn.append(np.std(r))
pcn=np.array(pcn)
print(f'Per-cond noise floor: median={np.median(pcn):.2f} p90={np.percentile(pcn,90):.2f} max={pcn.max():.2f}')
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/vc.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.base import BaseEstimator, TransformerMixin
FEATURES=["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]
class VCFeatures(BaseEstimator, TransformerMixin):
    '''Varying-coefficient design: spline basis in logSt and in lf, each interacted
    with condition descriptors [1, ang, ang2, lv, lc, lt], plus main effects.'''
    def __init__(self, n_knots=7):
        self.n_knots=n_knots
    def _base(self, X):
        f=pd.DataFrame(X,columns=FEATURES).astype(float).reset_index(drop=True)
        lf=np.log10(f['frequency']); lt=np.log10(f['suction_side_displacement_thickness'])
        lv=np.log10(f['free_stream_velocity']); lc=np.log10(f['chord_length']); ang=f['attack_angle']
        logSt=lf+lt-lv
        return lf,lt,lv,lc,ang,logSt
    def fit(self, X, y=None):
        lf,lt,lv,lc,ang,logSt=self._base(X)
        self.st_knots_=np.quantile(logSt, np.linspace(0,1,self.n_knots))
        self.lf_knots_=np.quantile(lf, np.linspace(0,1,self.n_knots))
        return self
    def _spl(self, v, knots):
        # natural-ish cubic via truncated power basis
        cols=[v, v**2, v**3]
        k=knots[1:-1]
        for kk in k:
            cols.append(np.clip(v-kk,0,None)**3)
        return np.vstack(cols).T
    def transform(self, X):
        lf,lt,lv,lc,ang,logSt=self._base(X)
        ang2=ang**2
        desc=np.vstack([np.ones(len(lf)), ang, ang2, lv, lc, lt]).T  # condition descriptors
        Sst=self._spl(logSt.values, self.st_knots_)
        Slf=self._spl(lf.values, self.lf_knots_)
        blocks=[desc]  # main effects (level as function of condition)
        # interactions: each spline basis column * each descriptor
        for S in [Sst, Slf]:
            for j in range(S.shape[1]):
                blocks.append(S[:,[j]]*desc)
        return np.hstack(blocks)
EOF
timeout 290 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats3 import phys3
from hybrid import Hybrid
from vc import VCFeatures
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def vcsm(al,nk=7): return make_pipeline(VCFeatures(nk),StandardScaler(),Ridge(alpha=al))
def et(feat): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=400,random_state=0,n_jobs=-1))
# VC smooth alone
for al in [1.0,5.0,20.0]:
    cv_eval(lambda al=al: vcsm(al), X,y,cond,label=f'CV VCsmooth al={al}')
print()
for al in [5.0,20.0]:
    mk=lambda al=al: Hybrid(vcsm(al), et(phys3))
    cv_eval(mk, X,y,cond,label=f'CV VChybrid al={al}')
    m=mk(); m.fit(Xtr,ytr); summarize(yva.values,m.predict(Xva),condv,label=f'VAL VChybrid al={al}'); print()
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/vc2.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.preprocessing import SplineTransformer
FEATURES=["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]
class VCFeatures2(BaseEstimator, TransformerMixin):
    def __init__(self, n_knots=7, degree=3):
        self.n_knots=n_knots; self.degree=degree
    def _base(self, X):
        f=pd.DataFrame(X,columns=FEATURES).astype(float).reset_index(drop=True)
        lf=np.log10(f['frequency']); lt=np.log10(f['suction_side_displacement_thickness'])
        lv=np.log10(f['free_stream_velocity']); lc=np.log10(f['chord_length']); ang=f['attack_angle']
        logSt=(lf+lt-lv)
        return lf.values,lt.values,lv.values,lc.values,ang.values,logSt.values
    def fit(self, X, y=None):
        lf,lt,lv,lc,ang,logSt=self._base(X)
        self.st_=SplineTransformer(n_knots=self.n_knots,degree=self.degree,extrapolation='linear').fit(logSt.reshape(-1,1))
        self.lf_=SplineTransformer(n_knots=self.n_knots,degree=self.degree,extrapolation='linear').fit(lf.reshape(-1,1))
        return self
    def _desc(self, ang, lv, lc, lt):
        return np.vstack([np.ones_like(ang), ang, ang**2, lv, lc, lt]).T
    def transform(self, X):
        lf,lt,lv,lc,ang,logSt=self._base(X)
        desc=self._desc(ang,lv,lc,lt)
        Sst=self.st_.transform(logSt.reshape(-1,1))
        Slf=self.lf_.transform(lf.reshape(-1,1))
        blocks=[desc]
        for S in [Sst, Slf]:
            for j in range(S.shape[1]):
                blocks.append(S[:,[j]]*desc)
        return np.hstack(blocks)
EOF
timeout 295 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats3 import phys3
from hybrid import Hybrid
from vc2 import VCFeatures2
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def vcsm(al,nk=7): return make_pipeline(VCFeatures2(nk),StandardScaler(with_mean=False),Ridge(alpha=al))
def et(feat): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=400,random_state=0,n_jobs=-1))
for al in [2.0,5.0,10.0]:
  for nk in [6,7]:
    mk=lambda al=al,nk=nk: Hybrid(vcsm(al,nk), et(phys3))
    cv_eval(mk, X,y,cond,label=f'CV VC2hyb al={al} nk={nk}')
    m=mk(); m.fit(Xtr,ytr); summarize(yva.values,m.predict(Xva),condv,label=f'VAL VC2hyb al={al} nk={nk}'); print()
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/blend4.py << 'EOF'
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats3 import phys3
from feats2 import phys2
from hybrid import Hybrid
from vc import VCFeatures
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def vcsm(al,nk=7): return make_pipeline(VCFeatures(nk),StandardScaler(),Ridge(alpha=al))
def et(feat,seed=0): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=400,random_state=seed,n_jobs=-1))
MODELS={
 'V1': lambda: Hybrid(vcsm(5.0,7), et(phys3)),
 'V2': lambda: Hybrid(vcsm(20.0,7), et(phys2)),
}
gkf=GroupKFold(n_splits=8)
oof={k:np.zeros(len(y)) for k in MODELS}; valp={k:np.zeros(len(yva)) for k in MODELS}
for k,mk in MODELS.items():
    for t,e in gkf.split(X,y,cond):
        m=mk(); m.fit(X.iloc[t],y.iloc[t]); oof[k][e]=m.predict(X.iloc[e])
    m=mk(); m.fit(Xtr,ytr); valp[k]=m.predict(Xva)
    summarize(y.values,oof[k],cond,label=f'CV {k}')
np.savez('/tmp/preds4.npz', keys=list(MODELS.keys()),
         **{f'oof_{k}':oof[k] for k in MODELS}, **{f'val_{k}':valp[k] for k in MODELS})
print('saved')
EOF
timeout 298 python3 /tmp/blend4.py 2>&1 | grep -v Warning

cd /app && timeout 280 python3 -c "
import numpy as np, sys; sys.path.insert(0,'/tmp')
from harness import per_condition_rmse
ds=[np.load(f'/tmp/preds{i}.npz',allow_pickle=True) for i in [2,3,4]]
y=ds[0]['y']; cond=ds[0]['cond']; yva=ds[0]['yva']; condv=ds[0]['condv']
keys=[]; 
for d in ds: keys+=list(d['keys'])
def get(k,pfx):
    for d in ds:
        if f'{pfx}_{k}' in d: return d[f'{pfx}_{k}']
oof=np.vstack([get(k,'oof') for k in keys]); val=np.vstack([get(k,'val') for k in keys])
T=dict(r=2.25,m=1.55,p90=2.35,p95=4.75,mx=5.60)
def metr(pred,yy,cc):
    pc=per_condition_rmse(yy,pred,cc)
    return dict(r=np.sqrt(np.mean((yy-pred)**2)),m=np.mean(np.abs(yy-pred)),p90=np.percentile(pc,90),p95=np.percentile(pc,95),mx=pc.max())
np.random.seed(0)
best=(1e9,None)
for it in range(60000):
    ww=np.random.dirichlet(np.ones(len(keys))*0.4)
    mc=metr(ww@oof,y,cond)
    # hard constraints on CV tail + RMSE/MAE, minimize p90 (CV + val avg)
    if mc['mx']>5.55 or mc['p95']>4.7 or mc['r']>2.24 or mc['m']>1.54: continue
    mv=metr(ww@val,yva,condv)
    if mv['mx']>5.55 or mv['r']>2.24 or mv['m']>1.60: continue
    sc=0.5*mc['p90']+0.5*mv['p90']
    if sc<best[0]: best=(sc,ww,mc,mv)
if best[1] is None:
    print('no feasible blend under constraints')
else:
    sc,ww,mc,mv=best
    ws=','.join(f'{k}:{w:.2f}' for k,w in zip(keys,ww) if w>0.02)
    print('Blend:',ws)
    print(f'  CV  r={mc[\"r\"]:.3f} m={mc[\"m\"]:.3f} p90={mc[\"p90\"]:.3f} p95={mc[\"p95\"]:.3f} mx={mc[\"mx\"]:.3f}')
    print(f'  VAL r={mv[\"r\"]:.3f} m={mv[\"m\"]:.3f} p90={mv[\"p90\"]:.3f} p95={mv[\"p95\"]:.3f} mx={mv[\"mx\"]:.3f}')
" 2>&1 | grep -v Warning

cd /app && timeout 200 python3 -c "
import numpy as np, sys; sys.path.insert(0,'/tmp')
ds=[np.load(f'/tmp/preds{i}.npz',allow_pickle=True) for i in [2,3,4]]
y=ds[0]['y']; cond=ds[0]['cond']; yva=ds[0]['yva']; condv=ds[0]['condv']
keys=[]
for d in ds: keys+=list(d['keys'])
def get(k,pfx):
    for d in ds:
        if f'{pfx}_{k}' in d: return d[f'{pfx}_{k}']
oof=np.vstack([get(k,'oof') for k in keys]); val=np.vstack([get(k,'val') for k in keys])
# factorize conditions for fast per-cond rmse
def factor(c):
    u,inv=np.unique(c,return_inverse=True); return inv,len(u)
ci,nc=factor(cond); cvi,ncv=factor(condv)
def pcrmse(err2, idx, n):
    s=np.bincount(idx,weights=err2,minlength=n); cnt=np.bincount(idx,minlength=n)
    return np.sqrt(s/cnt)
def metr(pred,yy,idx,n):
    e=yy-pred; pc=pcrmse(e*e,idx,n)
    return np.sqrt(np.mean(e*e)),np.mean(np.abs(e)),np.percentile(pc,90),np.percentile(pc,95),pc.max()
np.random.seed(0)
best=(1e9,None)
for it in range(20000):
    ww=np.random.dirichlet(np.ones(len(keys))*0.4)
    r,m,p90,p95,mx=metr(ww@oof,y,ci,nc)
    if mx>5.55 or p95>4.7 or r>2.24 or m>1.54: continue
    rv,mv,p90v,p95v,mxv=metr(ww@val,yva,cvi,ncv)
    if mxv>5.55 or rv>2.24 or mv>1.60: continue
    sc=0.5*p90+0.5*p90v
    if sc<best[0]: best=(sc,ww.copy(),(r,m,p90,p95,mx),(rv,mv,p90v,p95v,mxv))
if best[1] is None:
    print('no feasible blend')
else:
    sc,ww,mc,mv=best
    print('Blend:',', '.join(f'{k}:{w:.2f}' for k,w in zip(keys,ww) if w>0.02))
    print('  CV  r=%.3f m=%.3f p90=%.3f p95=%.3f mx=%.3f'%mc)
    print('  VAL r=%.3f m=%.3f p90=%.3f p95=%.3f mx=%.3f'%mv)
    np.save('/tmp/best_w.npy', ww); np.save('/tmp/keys.npy', np.array(keys))
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/hybrid_clip.py << 'EOF'
import numpy as np
from sklearn.base import BaseEstimator, RegressorMixin, clone
class HybridClip(BaseEstimator,RegressorMixin):
    def __init__(self, stage1, stage2, pad=6.0):
        self.stage1=stage1; self.stage2=stage2; self.pad=pad
    def fit(self,X,y):
        self.s1=clone(self.stage1).fit(X,y)
        self.lo=float(np.min(y))-self.pad; self.hi=float(np.max(y))+self.pad
        p1=np.clip(self.s1.predict(X),self.lo,self.hi)
        r=y-p1; self.s2=clone(self.stage2).fit(X,r); return self
    def predict(self,X):
        p1=np.clip(self.s1.predict(X),self.lo,self.hi)
        return p1+self.s2.predict(X)
EOF
timeout 295 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats3 import phys3
from feats2 import phys2
from hybrid_clip import HybridClip
from vc import VCFeatures
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
condv=recover_conditions(va); Xtr,ytr=tr[FEATURES],tr[TARGET]; Xva,yva=va[FEATURES],va[TARGET]
def vcsm(al,nk=7): return make_pipeline(VCFeatures(nk),StandardScaler(),Ridge(alpha=al))
def et(feat,seed=0): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=400,random_state=seed,n_jobs=-1))
for al,pad in [(5.0,4.0),(10.0,4.0),(20.0,3.0),(10.0,2.0)]:
    mk=lambda al=al,pad=pad: HybridClip(vcsm(al), et(phys3), pad=pad)
    cv_eval(mk, X,y,cond,label=f'CV VCclip al={al} pad={pad}')
    m=mk(); m.fit(Xtr,ytr); summarize(yva.values,m.predict(Xva),condv,label=f'VAL VCclip al={al} pad={pad}'); print()
" 2>&1 | grep -v Warning

cd /app && timeout 295 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from feats3 import phys3
from hybrid import Hybrid
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
def sm(feat,al): return make_pipeline(FunctionTransformer(feat),SplineTransformer(n_knots=6,degree=3),Ridge(alpha=al))
def et(feat): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=400,random_state=0,n_jobs=-1))
MODELS={'M1':lambda:Hybrid(sm(phys3,0.3),et(phys3)),'M3':lambda:Hybrid(sm(phys,0.3),et(phys)),'M4':lambda:et(phys3)}
for nsp in [8,15,25]:
    gkf=GroupKFold(n_splits=nsp); oof={k:np.zeros(len(y)) for k in MODELS}
    for k,mk in MODELS.items():
        for t,e in gkf.split(X,y,cond):
            m=mk(); m.fit(X.iloc[t],y.iloc[t]); oof[k][e]=m.predict(X.iloc[e])
    blend=0.2*oof['M1']+0.55*oof['M3']+0.25*oof['M4']
    summarize(y.values,blend,cond,label=f'blend nsplits={nsp}')
" 2>&1 | grep -v Warning

cd /app && timeout 298 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from feats3 import phys3
from hybrid import Hybrid
from vc import VCFeatures
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
def sm(feat,al): return make_pipeline(FunctionTransformer(feat),SplineTransformer(n_knots=6,degree=3),Ridge(alpha=al))
def vcsm(al,nk=7): return make_pipeline(VCFeatures(nk),StandardScaler(),Ridge(alpha=al))
def et(feat): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=400,random_state=0,n_jobs=-1))
MODELS={'M1':lambda:Hybrid(sm(phys3,0.3),et(phys3)),'M3':lambda:Hybrid(sm(phys,0.3),et(phys)),
        'M4':lambda:et(phys3),'V1':lambda:Hybrid(vcsm(10.0),et(phys3)),'V2':lambda:Hybrid(vcsm(20.0),et(phys))}
for nsp in [15,25]:
    gkf=GroupKFold(n_splits=nsp); oof={k:np.zeros(len(y)) for k in MODELS}
    for k,mk in MODELS.items():
        for t,e in gkf.split(X,y,cond):
            m=mk(); m.fit(X.iloc[t],y.iloc[t]); oof[k][e]=m.predict(X.iloc[e])
    for k in ['V1','V2']: summarize(y.values,oof[k],cond,label=f'{k} nsp={nsp}')
    b1=0.2*oof['M1']+0.55*oof['M3']+0.25*oof['M4']
    summarize(y.values,b1,cond,label=f'blend(no VC) nsp={nsp}')
    for wv in [0.2,0.35,0.5]:
        b=(1-wv)*b1+wv*(0.5*oof['V1']+0.5*oof['V2'])
        summarize(y.values,b,cond,label=f'blend+VC w={wv} nsp={nsp}')
    print()
" 2>&1 | grep -v Warning

cd /app && timeout 298 python3 -c "
import sys; sys.path.insert(0,'/tmp')
from harness import *
from feats import phys
from feats3 import phys3
from hybrid import Hybrid
from vc import VCFeatures
import numpy as np, pandas as pd
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
from sklearn.model_selection import GroupKFold
tr,va=load(); alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[FEATURES]; y=alld[TARGET]
def sm(feat,al): return make_pipeline(FunctionTransformer(feat),SplineTransformer(n_knots=6,degree=3),Ridge(alpha=al))
def vcsm(al,nk=7): return make_pipeline(VCFeatures(nk),StandardScaler(),Ridge(alpha=al))
def et(feat): return make_pipeline(FunctionTransformer(feat),ExtraTreesRegressor(n_estimators=300,random_state=0,n_jobs=-1))
MODELS={'M1':lambda:Hybrid(sm(phys3,0.3),et(phys3)),'M3':lambda:Hybrid(sm(phys,0.3),et(phys)),
        'M4':lambda:et(phys3),'V1':lambda:Hybrid(vcsm(10.0),et(phys3)),'V2':lambda:Hybrid(vcsm(20.0),et(phys))}
gkf=GroupKFold(n_splits=15); oof={k:np.zeros(len(y)) for k in MODELS}
for k,mk in MODELS.items():
    for t,e in gkf.split(X,y,cond):
        m=mk(); m.fit(X.iloc[t],y.iloc[t]); oof[k][e]=m.predict(X.iloc[e])
np.savez('/tmp/oof15.npz',y=y.values,cond=cond,**{k:oof[k] for k in MODELS})
for k in ['V1','V2']: summarize(y.values,oof[k],cond,label=f'{k}')
b1=0.2*oof['M1']+0.55*oof['M3']+0.25*oof['M4']
summarize(y.values,b1,cond,label='blend(no VC)')
for wv in [0.25,0.4,0.55]:
    b=(1-wv)*b1+wv*(0.5*oof['V1']+0.5*oof['V2'])
    summarize(y.values,b,cond,label=f'blend+VC w={wv}')
" 2>&1 | grep -v Warning

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.csv + validation.csv) and writes
predictions for the grader-staged evaluation features in test_features.csv.

Modeling strategy
-----------------
The hidden split holds out whole aerodynamic condition groups
``(attack_angle, chord_length, free_stream_velocity)``, so the surrogate must
generalize to unseen physical regimes rather than interpolate nearby rows.

Two ideas drive the design:

1. Physics-informed features.  Airfoil self-noise collapses onto a master
   spectrum when expressed through the boundary-layer Strouhal number
   ``St = f * delta* / U`` together with Reynolds / Mach style amplitude terms.
   Building features in log space around these quantities lets the model
   extrapolate the spectral *shape* across conditions.

2. Semi-parametric hybrids.  Each base model is a smooth stage-1 regressor
   (spline / varying-coefficient ridge) that captures the global, smoothly
   varying physics and therefore extrapolates gracefully to unseen condition
   groups, followed by an ExtraTrees stage-2 fit on the residual that cleans up
   local structure.  Several such hybrids (plus a plain ExtraTrees) are blended.
   The blend averages away the occasional catastrophic per-condition miss that
   any single smooth model can make on an extreme held-out regime, which keeps
   the worst / tail per-condition errors bounded.

The whole workflow is deterministic and self-contained (no internet, no
external data, and the raw source labels are never used).
"""

from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.base import BaseEstimator, RegressorMixin, TransformerMixin, clone
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, SplineTransformer, 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",
]

SPEED_OF_SOUND = 340.0
N_TREES = 500
RANDOM_STATE = 0


# --------------------------------------------------------------------------- #
# Physics-informed feature builders
# --------------------------------------------------------------------------- #
def _raw(X):
    f = pd.DataFrame(X, columns=FEATURES).astype(float).reset_index(drop=True)
    freq = f["frequency"].clip(lower=1e-6)
    ang = f["attack_angle"]
    chord = f["chord_length"].clip(lower=1e-9)
    vel = f["free_stream_velocity"].clip(lower=1e-6)
    thk = f["suction_side_displacement_thickness"].clip(lower=1e-12)
    lf = np.log10(freq)
    lt = np.log10(thk)
    lv = np.log10(vel)
    lc = np.log10(chord)
    return freq, ang, chord, vel, thk, lf, lt, lv, lc


def phys(X):
    """Compact physics feature set."""
    freq, ang, chord, vel, thk, lf, lt, lv, lc = _raw(X)
    logSt = lf + lt - lv
    o = pd.DataFrame(index=range(len(lf)))
    o["lf"] = lf
    o["lt"] = lt
    o["lv"] = lv
    o["lc"] = lc
    o["ang"] = ang
    o["logSt"] = logSt
    o["logRe_c"] = lv + lc
    o["logRe_d"] = lv + lt
    o["mach"] = vel / SPEED_OF_SOUND
    o["ang_sin"] = np.sin(np.radians(ang))
    o["logSt2"] = logSt ** 2
    o["lf_ang"] = lf * ang
    o["lf_vel"] = lf * vel
    return o.to_numpy(float)


def phys2(X):
    """Physics features with extra Strouhal / angle interactions."""
    freq, ang, chord, vel, thk, lf, lt, lv, lc = _raw(X)
    logSt = lf + lt - lv
    o = pd.DataFrame(index=range(len(lf)))
    o["lf"] = lf
    o["lt"] = lt
    o["lv"] = lv
    o["lc"] = lc
    o["ang"] = ang
    o["logSt"] = logSt
    o["logRe_c"] = lv + lc
    o["logRe_d"] = lv + lt
    o["mach"] = vel / SPEED_OF_SOUND
    o["ang_sin"] = np.si
…[truncated]

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

cd /app && timeout 200 python3 -c "
import sys; sys.path.insert(0,'/app'); sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
import train_and_predict as T
from harness import recover_conditions, summarize
tr=pd.read_csv('/app/data/train.csv'); va=pd.read_csv('/app/data/validation.csv')
Xtr,ytr=tr[T.FEATURES],tr[T.TARGET].astype(float)
Xva,yva=va[T.FEATURES],va[T.TARGET].astype(float)
condv=recover_conditions(va)
tot=0.0; pred=np.zeros(len(Xva))
for name,fac,w in T.build_models():
    m=fac(); m.fit(Xtr,ytr); p=m.predict(Xva); pred+=w*p; tot+=w
    summarize(yva.values,p,condv,label=name)
pred/=tot
print('='*60)
summarize(yva.values,pred,condv,label='FINAL BLEND (train->val)')
" 2>&1 | grep -v Warning

cd /app && cat > /tmp/finalcv.py << 'EOF'
import sys; sys.path.insert(0,'/app'); sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
import train_and_predict as T
from harness import recover_conditions
from sklearn.model_selection import GroupKFold
tr=pd.read_csv('/app/data/train.csv'); va=pd.read_csv('/app/data/validation.csv')
alld=pd.concat([tr,va],ignore_index=True)
cond=recover_conditions(alld); X=alld[T.FEATURES]; y=alld[T.TARGET].astype(float)
condv=recover_conditions(va); Xtr,ytr=tr[T.FEATURES],tr[T.TARGET].astype(float)
Xva,yva=va[T.FEATURES],va[T.TARGET].astype(float)
names=[n for n,_,_ in T.build_models()]; facs={n:f for n,f,_ in T.build_models()}
gkf=GroupKFold(n_splits=8)
oof={n:np.zeros(len(y)) for n in names}; valp={n:np.zeros(len(yva)) for n in names}
for n in names:
    for t,e in gkf.split(X,y,cond):
        m=facs[n](); m.fit(X.iloc[t],y.iloc[t]); oof[n][e]=m.predict(X.iloc[e])
    m=facs[n](); m.fit(Xtr,ytr); valp[n]=m.predict(Xva)
np.savez('/tmp/final_oof.npz', y=y.values, cond=cond, yva=yva.values, condv=condv,
         names=names, **{f'o_{n}':oof[n] for n in names}, **{f'v_{n}':valp[n] for n in names})
print('done', names)
EOF
timeout 298 python3 /tmp/finalcv.py 2>&1 | grep -v Warning

cd /app && timeout 200 python3 -c "
import numpy as np
d=np.load('/tmp/final_oof.npz',allow_pickle=True)
names=list(d['names']); y=d['y']; cond=d['cond']; yva=d['yva']; condv=d['condv']
O=np.vstack([d[f'o_{n}'] for n in names]); V=np.vstack([d[f'v_{n}'] for n in names])
def fac(c): u,inv=np.unique(c,return_inverse=True); return inv,len(u)
ci,nc=fac(cond); cvi,ncv=fac(condv)
def pc(e2,idx,n): 
    s=np.bincount(idx,weights=e2,minlength=n); c=np.bincount(idx,minlength=n); return np.sqrt(s/c)
def metr(pred,yy,idx,n):
    e=yy-pred; p=pc(e*e,idx,n)
    return dict(r=np.sqrt(np.mean(e*e)),m=np.mean(np.abs(e)),p90=np.percentile(p,90),p95=np.percentile(p,95),mx=p.max())
T=dict(r=2.25,m=1.55,p90=2.35,p95=4.75,mx=5.60)
def score(md): return max(v/T[k] for k,v in md.items())
np.random.seed(0); best=(1e9,None)
for _ in range(40000):
    w=np.random.dirichlet(np.ones(len(names))*0.5)
    mc=metr(w@O,y,ci,nc)
    if mc['mx']>5.55 or mc['p95']>4.7: continue
    mv=metr(w@V,yva,cvi,ncv)
    if mv['mx']>5.55: continue
    # robust objective: emphasize p90 & mae on both, with CV weighted more
    obj=0.55*(0.5*mc['p90']/T['p90']+0.5*mc['m']/T['m'])+0.45*(0.5*mv['p90']/T['p90']+0.5*mv['m']/T['m'])
    if obj<best[0]: best=(obj,w.copy(),mc,mv)
obj,w,mc,mv=best
print('weights:',{n:round(float(x),3) for n,x in zip(names,w)})
print('CV :',{k:round(v,3) for k,v in mc.items()})
print('VAL:',{k:round(v,3) for k,v in mv.items()})
np.save('/tmp/fw.npy',w)
" 2>&1 | grep -v Warning

cd /app && timeout 200 python3 -c "
import numpy as np
d=np.load('/tmp/final_oof.npz',allow_pickle=True)
names=list(d['names']); y=d['y']; cond=d['cond']; yva=d['yva']; condv=d['condv']
O=np.vstack([d[f'o_{n}'] for n in names]); V=np.vstack([d[f'v_{n}'] for n in names])
def fac(c): u,inv=np.unique(c,return_inverse=True); return inv,len(u)
ci,nc=fac(cond); cvi,ncv=fac(condv)
def pc(e2,idx,n): s=np.bincount(idx,weights=e2,minlength=n); c=np.bincount(idx,minlength=n); return np.sqrt(s/c)
def metr(pred,yy,idx,n):
    e=yy-pred; p=pc(e*e,idx,n)
    return dict(r=np.sqrt(np.mean(e*e)),m=np.mean(np.abs(e)),p90=np.percentile(p,90),p95=np.percentile(p,95),mx=p.max())
np.random.seed(1); best=(1e9,None)
for _ in range(80000):
    w=np.random.dirichlet(np.ones(len(names))*0.5)
    mc=metr(w@O,y,ci,nc)
    if mc['m']>1.53 or mc['mx']>5.55 or mc['p95']>4.6 or mc['r']>2.24: continue
    mv=metr(w@V,yva,cvi,ncv)
    if mv['m']>1.58 or mv['mx']>5.0 or mv['r']>2.2: continue
    obj=0.6*mc['p90']+0.4*mv['p90']
    if obj<best[0]: best=(obj,w.copy(),mc,mv)
obj,w,mc,mv=best
print('weights:',{n:round(float(x),3) for n,x in zip(names,w)})
print('CV :',{k:round(v,3) for k,v in mc.items()})
print('VAL:',{k:round(v,3) for k,v in mv.items()})
np.save('/tmp/fw2.npy',w)
" 2>&1 | grep -v Warning

cd /app && timeout 298 python3 -c "
import sys; sys.path.insert(0,'/app'); sys.path.insert(0,'/tmp')
import numpy as np, pandas as pd
import train_and_predict as T
from harness import recover_conditions
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.model_selection import GroupKFold
tr=pd.read_csv('/app/data/train.csv'); va=pd.read_csv('/app/data/validation.csv')
alld=pd.concat([tr,va],ignore_index=True); cond=recover_conditions(alld)
X=alld[T.FEATURES]; y=alld[T.TARGET].astype(float); condv=recover_conditions(va)
Xtr,ytr=tr[T.FEATURES],tr[T.TARGET].astype(float); Xva,yva=va[T.FEATURES],va[T.TARGET].astype(float)
def rf(): return make_pipeline(FunctionTransformer(T.phys3),RandomForestRegressor(n_estimators=400,random_state=0,n_jobs=-1))
def gb(): return make_pipeline(FunctionTransformer(T.phys3),GradientBoostingRegressor(n_estimators=400,learning_rate=0.03,max_depth=3,subsample=0.8,random_state=0))
M={'rf_h':lambda:T.Hybrid(T._spline_ridge(T.phys3,0.3),rf()),'gb_h':lambda:T.Hybrid(T._spline_ridge(T.phys3,0.3),gb())}
gkf=GroupKFold(n_splits=8); oof={k:np.zeros(len(y)) for k in M}; vp={k:np.zeros(len(yva)) for k in M}
for k,mk in M.items():
    for t,e in gkf.split(X,y,cond):
        m=mk(); m.fit(X.iloc[t],y.iloc[t]); oof[k][e]=m.predict(X.iloc[e])
    m=mk(); m.fit(Xtr,ytr); vp[k]=m.predict(Xva)
np.savez('/tmp/extra_oof.npz',**{f'o_{k}':oof[k] for k in M},**{f'v_{k}':vp[k] for k in M},names=list(M.keys()))
print('done')
" 2>&1 | grep -v Warning

cd /app && timeout 220 python3 -c "
import numpy as np
d=np.load('/tmp/final_oof.npz',allow_pickle=True); e=np.load('/tmp/extra_oof.npz',allow_pickle=True)
names=list(d['names'])+list(e['names'])
y=d['y']; cond=d['cond']; yva=d['yva']; condv=d['condv']
def oget(n): return d[f'o_{n}'] if f'o_{n}' in d else e[f'o_{n}']
def vget(n): return d[f'v_{n}'] if f'v_{n}' in d else e[f'v_{n}']
O=np.vstack([oget(n) for n in names]); V=np.vstack([vget(n) for n in names])
def fac(c): u,inv=np.unique(c,return_inverse=True); return inv,len(u)
ci,nc=fac(cond); cvi,ncv=fac(condv)
def pc(e2,idx,n): s=np.bincount(idx,weights=e2,minlength=n); c=np.bincount(idx,minlength=n); return np.sqrt(s/c)
def metr(pred,yy,idx,n):
    er=yy-pred; p=pc(er*er,idx,n)
    return dict(r=np.sqrt(np.mean(er*er)),m=np.mean(np.abs(er)),p90=np.percentile(p,90),p95=np.percentile(p,95),mx=p.max())
np.random.seed(3); best=(1e9,None)
for _ in range(120000):
    w=np.random.dirichlet(np.ones(len(names))*0.45)
    mc=metr(w@O,y,ci,nc)
    if mc['m']>1.52 or mc['mx']>5.5 or mc['p95']>4.5 or mc['r']>2.2: continue
    mv=metr(w@V,yva,cvi,ncv)
    if mv['m']>1.57 or mv['mx']>5.0 or mv['r']>2.15: continue
    obj=0.6*mc['p90']+0.4*mv['p90']
    if obj<best[0]: best=(obj,w.copy(),mc,mv)
obj,w,mc,mv=best
print('names:',names)
print('weights:',[round(float(x),3) for x in w])
print('CV :',{k:round(v,3) for k,v in mc.items()})
print('VAL:',{k:round(v,3) for k,v in mv.items()})
np.save('/tmp/fw3.npy',w); 
import json; open('/tmp/fw3_names.json','w').write(json.dumps(names))
" 2>&1 | grep -v Warning

cd /app && timeout 120 python3 -c "
import numpy as np
d=np.load('/tmp/final_oof.npz',allow_pickle=True)
names=list(d['names']); y=d['y']; cond=d['cond']; yva=d['yva']; condv=d['condv']
W={'hyb_phys3':0.25,'hyb_phys':0.14,'et_phys3':0.22,'hyb_phys2':0.12,'vc_phys3':0.27,'vc_phys':0.0}
def blend(pfx,n):
    s=sum(W[nm]*d[f'{pfx}_{nm}'] for nm in names); return s/sum(W.values())
def fac(c): u,inv=np.unique(c,return_inverse=True); return inv,len(u)
def show(pred,yy,c,lab):
    idx,n=fac(c); e=yy-pred
    s=np.bincount(idx,weights=e*e,minlength=n); cnt=np.bincount(idx,minlength=n); p=np.sqrt(s/cnt)
    print(f'{lab}: RMSE={np.sqrt(np.mean(e*e)):.3f} MAE={np.mean(np.abs(e)):.3f} p90={np.percentile(p,90):.3f} p95={np.percentile(p,95):.3f} max={p.max():.3f}')
show(blend('o',names),y,cond,'CV  (8-fold, train~72 conds)')
show(blend('v',names),yva,condv,'VAL (train 67 conds)')
print('Targets:            RMSE<=2.25 MAE<=1.55 p90<=2.35 p95<=4.75 max<=5.60')
" 2>&1 | grep -v Warning

command: str_replace
path: /app/train_and_predict.py
old_str:     return [
        ("hyb_phys3", lambda: Hybrid(_spline_ridge(phys3, 0.3), _et(phys3)), 0.20),
        ("hyb_phys", lambda: Hybrid(_spline_ridge(phys, 0.3), _et(phys)), 0.34),
        ("et_phys3", lambda: _et(phys3), 0.18),
        ("hyb_phys2", lambda: Hybrid(_spline_ridge(phys2, 0.3), _et(phys3)), 0.12),
        ("vc_phys3", lambda: Hybrid(_vc_ridge(10.0), _et(phys3)), 0.08),
        ("vc_phys", lambda: Hybrid(_vc_ridge(20.0), _et(phys)), 0.08),
    ]
new_str:     return [
        ("hyb_phys3", lambda: Hybrid(_spline_ridge(phys3, 0.3), _et(phys3)), 0.25),
        ("hyb_phys", lambda: Hybrid(_spline_ridge(phys, 0.3), _et(phys)), 0.14),
        ("et_phys3", lambda: _et(phys3), 0.22),
        ("hyb_phys2", lambda: Hybrid(_spline_ridge(phys2, 0.3), _et(phys3)), 0.12),
        ("vc_phys3", lambda: Hybrid(_vc_ridge(10.0), _et(phys3)), 0.27),
    ]

cd /app && rm -f output/predictions.csv && time python3 train_and_predict.py 2>&1 | grep -viE "warn" | tail -5
echo "=== VALIDATION ==="
python3 -c "
import pandas as pd, numpy as np
feat=pd.read_csv('/app/data/test_features.csv')
pred=pd.read_csv('/app/output/predictions.csv')
print('columns:', list(pred.columns))
assert list(pred.columns)==['sample_id','predicted_scaled_sound_pressure_level'], 'bad columns'
print('n_pred:', len(pred), 'n_feat:', len(feat))
assert len(pred)==len(feat), 'row count mismatch'
assert list(pred['sample_id'].astype(str))==list(feat['sample_id'].astype(str)), 'sample_id misalignment'
assert pred['predicted_scaled_sound_pressure_level'].notna().all(), 'NaNs'
assert np.isfinite(pred['predicted_scaled_sound_pressure_level']).all(), 'non-finite'
assert pred['predicted_scaled_sound_pressure_level'].nunique()>1, 'constant'
print('pred range: [%.2f, %.2f]'%(pred['predicted_scaled_sound_pressure_level'].min(),pred['predicted_scaled_sound_pressure_level'].max()))
print('ALL FORMAT CHECKS PASSED')
"

cd /app && cp data/test_features.csv /tmp/tf_backup.csv && md5sum data/test_features.csv
echo "--- create synthetic hidden test (7 novel rows, new ids, shuffled cols order) ---"
python3 -c "
import pandas as pd
tr=pd.read_csv('/app/data/train.csv')
# take some rows, perturb slightly, assign NEW sample_ids, drop target
s=tr.sample(7,random_state=42).reset_index(drop=True)
s['sample_id']=['hidden_%04d'%i for i in range(7)]
cols=['sample_id','frequency','attack_angle','chord_length','free_stream_velocity','suction_side_displacement_thickness']
s[cols].to_csv('/app/data/test_features.csv',index=False)
print('wrote synthetic test_features with', len(s), 'rows')
"
rm -f output/predictions.csv && python3 train_and_predict.py 2>&1 | grep -viE "warn" | tail -3
python3 -c "
import pandas as pd, numpy as np
f=pd.read_csv('/app/data/test_features.csv'); p=pd.read_csv('/app/output/predictions.csv')
assert list(p.columns)==['sample_id','predicted_scaled_sound_pressure_level']
assert list(p.sample_id.astype(str))==list(f.sample_id.astype(str)), 'misaligned'
assert len(p)==7 and np.isfinite(p.iloc[:,1]).all()
print('Synthetic hidden run OK. ids:', list(p.sample_id))
print(p.to_string(index=False))
"
echo "--- restore original ---"
cp /tmp/tf_backup.csv data/test_features.csv && md5sum data/test_features.csv

cd /app && rm -f output/predictions.csv && python3 train_and_predict.py 2>/dev/null && cp output/predictions.csv /tmp/run1.csv && python3 train_and_predict.py 2>/dev/null && cp output/predictions.csv /tmp/run2.csv
if diff -q /tmp/run1.csv /tmp/run2.csv >/dev/null; then echo "DETERMINISTIC: two runs identical"; else echo "NON-DETERMINISTIC"; diff /tmp/run1.csv /tmp/run2.csv | head; fi
echo "--- confirm no /tmp imports in script ---"
grep -nE "sys.path|/tmp|import harness|import feats|import vc|import hybrid" train_and_predict.py || echo "clean: no external/dev imports"
echo "--- final head of predictions ---"; head -3 output/predictions.csv

cd /app && python3 -c "import ast; ast.parse(open('train_and_predict.py').read()); print('syntax OK')" && sed -n '250,320p' train_and_predict.py

cd /app && ls -la output/predictions.csv && echo && timeout 120 python3 -c "
import numpy as np
d=np.load('/tmp/final_oof.npz',allow_pickle=True)
names=['hyb_phys3','hyb_phys','et_phys3','hyb_phys2','vc_phys3']
W={'hyb_phys3':0.25,'hyb_phys':0.14,'et_phys3':0.22,'hyb_phys2':0.12,'vc_phys3':0.27}
tw=sum(W.values())
def fac(c): u,inv=np.unique(c,return_inverse=True); return inv,len(u)
def m(pred,yy,c):
    idx,n=fac(c); e=yy-pred
    s=np.bincount(idx,weights=e*e,minlength=n); cnt=np.bincount(idx,minlength=n); p=np.sqrt(s/cnt)
    return np.sqrt(np.mean(e*e)),np.mean(np.abs(e)),np.percentile(p,90),np.percentile(p,95),p.max()
ocv=sum(W[n]*d[f'o_{n}'] for n in names)/tw
oval=sum(W[n]*d[f'v_{n}'] for n in names)/tw
T=[2.25,1.55,2.35,4.75,5.60]
labels=['RMSE','MAE','p90_cond','p95_cond','max_cond']
print('%-26s %8s %8s %8s %8s %8s'%('Scenario','RMSE','MAE','p90','p95','max'))
print('%-26s %8.3f %8.3f %8.3f %8.3f %8.3f'%(('CV 8-fold (train~72 conds)',)+m(ocv,d['y'],d['cond'])))
print('%-26s %8.3f %8.3f %8.3f %8.3f %8.3f'%(('TRAIN->VAL (67 conds)',)+m(oval,d['yva'],d['condv'])))
print('%-26s %8.2f %8.2f %8.2f %8.2f %8.2f'%(('TARGET (<=)',)+tuple(T)))
print()
print('Final model trains on all 82 conds -> expect further improvement')
print('(fold-count study: 72->79 conds moved p90 2.73->2.61, MAE 1.50->1.42)')
"