airfoil-self-noise

cd /app && ls -la && echo "---TREE---" && find . -maxdepth 3 -not -path '*/node_modules/*' | head -100

cd /app && cat train_and_predict.py

cd /app && 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 ===" && cat data/test_features.csv && wc -l data/test_features.csv

cd /app && python3 --version && python3 -c "import sklearn; print('sklearn', sklearn.__version__)" && python3 -c "import numpy; print('numpy', numpy.__version__)" && python3 -c "import pandas; print('pandas', pandas.__version__)" && python3 -c "import scipy; print('scipy', scipy.__version__)"

cd /app && for p in xgboost lightgbm catboost; do python3 -c "import $p; print('$p', $p.__version__)" 2>/dev/null || echo "$p NOT available"; done

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv')
va = pd.read_csv('data/validation.csv')
print('train shape', tr.shape, 'val shape', va.shape)
print()
print('=== TRAIN describe ===')
print(tr.describe().T)
print()
print('=== target summary ===')
print('train target', tr.scaled_sound_pressure_level.min(), tr.scaled_sound_pressure_level.max(), tr.scaled_sound_pressure_level.mean())
print('val target', va.scaled_sound_pressure_level.min(), va.scaled_sound_pressure_level.max(), va.scaled_sound_pressure_level.mean())
"

cd /app && python3 -c "
import pandas as pd, numpy as np
pd.set_option('display.width', 200)
tr = pd.read_csv('data/train.csv')
va = pd.read_csv('data/validation.csv')
for name, df in [('train', tr), ('val', va)]:
    print(f'=== {name} unique values ===')
    for c in ['frequency','attack_angle','chord_length','free_stream_velocity','suction_side_displacement_thickness']:
        print(f'{c}: n_unique={df[c].nunique()}, min={df[c].min():.5g}, max={df[c].max():.5g}')
    print()
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv')
va = pd.read_csv('data/validation.csv')
both = pd.concat([tr, va], ignore_index=True)
# nominal frequency: 1/3 octave bands
bands = [200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150,4000,5000,6300,8000,10000,12500,16000,20000]
def nearest(v, arr):
    arr = np.array(arr)
    return arr[np.argmin(np.abs(arr[None,:]-np.array(v)[:,None]), axis=1)]
both['nom_freq'] = nearest(both.frequency.values, bands)
# relative error
both['freq_relerr'] = (both.frequency - both.nom_freq)/both.nom_freq
print('freq rel err range:', both.freq_relerr.min(), both.freq_relerr.max())
# chord nominal
chords=[0.0254,0.0508,0.1016,0.1524,0.2286,0.3048]
both['nom_chord']=nearest(both.chord_length.values, chords)
print('chord rel err:', ((both.chord_length-both.nom_chord)/both.nom_chord).abs().max())
# velocity nominal
vels=[31.7,39.6,55.5,71.3]
both['nom_vel']=nearest(both.free_stream_velocity.values, vels)
print('vel abs err:', (both.free_stream_velocity-both.nom_vel).abs().max())
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv'); va = pd.read_csv('data/validation.csv')
both = pd.concat([tr, va], ignore_index=True)
# attack angle nominal candidates from UCI
angles=[0.0,1.5,2.0,3.0,4.0,5.3,6.7,7.3,8.4,8.9,9.9,11.2,12.3,12.6,15.4,15.6,17.4,19.7,22.2]
a=np.sort(both.attack_angle.unique())
# cluster by rounding to nearest 0.1 then group
from collections import Counter
r=np.round(both.attack_angle.values,1)
c=Counter(r)
print('distinct rounded-0.1 angles:', sorted(c.keys()))
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv'); va = pd.read_csv('data/validation.csv')
both = pd.concat([tr, va], ignore_index=True)
bands = [200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150,4000,5000,6300,8000,10000,12500,16000,20000]
chords=[0.0254,0.0508,0.1016,0.1524,0.2286,0.3048]; vels=[31.7,39.6,55.5,71.3]
def nearest(v, arr):
    arr=np.array(arr); return arr[np.argmin(np.abs(arr[None,:]-np.array(v)[:,None]),axis=1)]
both['nom_chord']=nearest(both.chord_length.values,chords)
both['nom_vel']=nearest(both.free_stream_velocity.values,vels)
# cluster attack angle: round to nearest known angle set
angset=np.array([0.0,1.5,3.0,4.0,5.3,7.3,8.9,9.9,12.3,12.6,15.4,15.6,17.4,19.7,22.2])
both['nom_ang']=nearest(both.attack_angle.values, angset)
# group and look at displacement thickness spread within group
g=both.groupby(['nom_ang','nom_chord','nom_vel'])
print('n groups:', g.ngroups)
spread=g['suction_side_displacement_thickness'].agg(lambda x:(x.max()-x.min())/x.mean())
print('within-group rel spread of thickness: max', spread.max(),'median',spread.median())
sz=g.size()
print('group sizes: min',sz.min(),'max',sz.max(),'median',sz.median())
print('n groups with size<=3:', (sz<=3).sum())
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv'); va = pd.read_csv('data/validation.csv')
chords=[0.0254,0.0508,0.1016,0.1524,0.2286,0.3048]; vels=[31.7,39.6,55.5,71.3]
angset=np.array([0.0,1.5,3.0,4.0,5.3,7.3,8.9,9.9,12.3,12.6,15.4,15.6,17.4,19.7,22.2])
def nearest(v, arr):
    arr=np.array(arr); return arr[np.argmin(np.abs(arr[None,:]-np.array(v)[:,None]),axis=1)]
def cond(df):
    return set(zip(nearest(df.attack_angle.values,angset), nearest(df.chord_length.values,chords), nearest(df.free_stream_velocity.values,vels)))
ctr=cond(tr); cva=cond(va)
print('train conds:', len(ctr), 'val conds:', len(cva))
print('overlap train&val:', len(ctr & cva))
print('val-only conds:', len(cva - ctr))
"

cd /app && python3 -c "
import pandas as pd, numpy as np
tr = pd.read_csv('data/train.csv'); va = pd.read_csv('data/validation.csv')
both=pd.concat([tr,va],ignore_index=True)
chords=[0.0254,0.0508,0.1016,0.1524,0.2286,0.3048]; vels=[31.7,39.6,55.5,71.3]
angset=np.array([0.0,1.5,3.0,4.0,5.3,7.3,8.9,9.9,12.3,12.6,15.4,15.6,17.4,19.7,22.2])
def nearest(v, arr):
    arr=np.array(arr); return arr[np.argmin(np.abs(arr[None,:]-np.array(v)[:,None]),axis=1)]
both['na']=nearest(both.attack_angle.values,angset)
# check within-group angle spread
both['nc']=nearest(both.chord_length.values,chords); both['nv']=nearest(both.free_stream_velocity.values,vels)
g=both.groupby(['na','nc','nv'])
aspread=g['attack_angle'].agg(lambda x:x.max()-x.min())
print('max within-group angle spread (deg):', aspread.max())
print('so angle jitter is small; nominal rounding is safe')
# Distribution of target within condition vs across
print('overall target std:', both.scaled_sound_pressure_level.std())
print('within-condition target std (mean):', g['scaled_sound_pressure_level'].std().mean())
"

cd /app && mkdir -p scratch && cat > scratch/evallib.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.model_selection import GroupKFold

FEATURES = ["frequency","attack_angle","chord_length","free_stream_velocity","suction_side_displacement_thickness"]
TARGET = "scaled_sound_pressure_level"
CHORDS=[0.0254,0.0508,0.1016,0.1524,0.2286,0.3048]; VELS=[31.7,39.6,55.5,71.3]
ANGSET=np.array([0.0,1.5,3.0,4.0,5.3,7.3,8.9,9.9,12.3,12.6,15.4,15.6,17.4,19.7,22.2])

def nearest(v, arr):
    arr=np.array(arr); return arr[np.argmin(np.abs(arr[None,:]-np.array(v)[:,None]),axis=1)]

def condition_ids(df):
    na=nearest(df.attack_angle.values,ANGSET); nc=nearest(df.chord_length.values,CHORDS); nv=nearest(df.free_stream_velocity.values,VELS)
    return pd.Series([f"{a:.2f}_{c:.4f}_{v:.1f}" for a,c,v in zip(na,nc,nv)], index=df.index)

def load():
    tr=pd.read_csv('/app/data/train.csv'); va=pd.read_csv('/app/data/validation.csv')
    return tr, va

def metrics_from(df_pred):
    # df_pred has columns: cond, err
    g=df_pred.groupby('cond')['err']
    cond_rmse=g.apply(lambda e: np.sqrt(np.mean(e**2)))
    allerr=df_pred['err'].values
    return {
        'rmse': float(np.sqrt(np.mean(allerr**2))),
        'mae': float(np.mean(np.abs(allerr))),
        'p90_cond': float(np.percentile(cond_rmse,90)),
        'p95_cond': float(np.percentile(cond_rmse,95)),
        'max_cond': float(cond_rmse.max()),
        'n_cond': len(cond_rmse),
    }

def grouped_cv(make_model, data, n_splits=10, seed=0, repeats=1):
    """Return pooled out-of-fold metrics across grouped CV."""
    rows=[]
    groups=condition_ids(data).values
    uniq=np.unique(groups)
    rng=np.random.RandomState(seed)
    for rep in range(repeats):
        # shuffle group assignment to folds
        gk=GroupKFold(n_splits=n_splits)
        # To get variety across repeats, remap groups via permutation
        perm=rng.permutation(len(uniq))
        gmap={u:i for i,u in enumerate(uniq[perm])}
        gidx=np.array([gmap[g] for g in groups])
        for tr_idx, te_idx in gk.split(data, data[TARGET], gidx):
            Xtr=data.iloc[tr_idx]; Xte=data.iloc[te_idx]
            m=make_model(); m.fit(Xtr[FEATURES], Xtr[TARGET].values)
            pred=m.predict(Xte[FEATURES])
            rows.append(pd.DataFrame({'cond':groups[te_idx],'err':pred-Xte[TARGET].values,'rep':rep}))
    allp=pd.concat(rows, ignore_index=True)
    # per-rep metrics then average for stability of percentile stats
    res=[]
    for rep in range(repeats):
        res.append(metrics_from(allp[allp.rep==rep]))
    out={k:np.mean([r[k] for r in res]) for k in res[0]}
    return out, allp

def val_holdout(make_model, tr, va):
    m=make_model(); m.fit(tr[FEATURES], tr[TARGET].values)
    pred=m.predict(va[FEATURES])
    dfp=pd.DataFrame({'cond':condition_ids(va).values,'err':pred-va[TARGET].values})
    return metrics_from(dfp)

def report(name, d):
    print(f"{name:28s} rmse={d['rmse']:.3f} mae={d['mae']:.3f} p90={d['p90_cond']:.3f} p95={d['p95_cond']:.3f} max={d['max_cond']:.3f} (nc={d.get('n_cond','?')})")
EOF
echo "written"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd
import evallib as E
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.ensemble import RandomForestRegressor, ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.neural_network import MLPRegressor

tr, va = E.load()
both = pd.concat([tr,va], ignore_index=True)

def ridge():
    return make_pipeline(StandardScaler(), Ridge(alpha=10.0))
def rf():
    return RandomForestRegressor(n_estimators=300, n_jobs=-1, random_state=0)
def et():
    return ExtraTreesRegressor(n_estimators=300, n_jobs=-1, random_state=0)
def hgb():
    return HistGradientBoostingRegressor(random_state=0)

for name, mk in [('ridge',ridge),('rf',rf),('et',et),('hgb',hgb)]:
    d,_=E.grouped_cv(mk, both, n_splits=10, seed=0, repeats=1)
    E.report('cv10 '+name, d)
" 2>&1 | grep -v Warning

cd /app/scratch && cat > feats.py << 'EOF'
import numpy as np, pandas as pd
C_SOUND=340.0
def phys_features(X):
    f=X['frequency'].astype(float).values
    a=X['attack_angle'].astype(float).values
    c=X['chord_length'].astype(float).values
    U=X['free_stream_velocity'].astype(float).values
    d=X['suction_side_displacement_thickness'].astype(float).values
    logf=np.log10(f); logd=np.log10(d); logc=np.log10(c); logU=np.log10(U)
    St=f*d/U                     # Strouhal
    logSt=np.log10(St)
    M=U/C_SOUND
    logM=np.log10(M)
    Re_c=U*c/1.5e-5              # chord Reynolds
    logRe=np.log10(Re_c)
    out=pd.DataFrame(index=X.index)
    out['logf']=logf; out['a']=a; out['logc']=logc; out['logU']=logU; out['logd']=logd
    out['logSt']=logSt; out['logSt2']=logSt**2; out['logSt3']=logSt**3
    out['logM']=logM; out['logRe']=logRe
    out['f_d_U']=St
    out['a2']=a**2
    out['logf2']=logf**2
    out['logf_a']=logf*a
    out['logf_logc']=logf*logc
    out['logf_logU']=logf*logU
    out['logd_a']=logd*a
    out['logSt_a']=logSt*a
    out['logSt_logc']=logSt*logc
    out['logSt_logU']=logSt*logU
    out['a_logU']=a*logU
    return out
EOF
python3 -c "
import numpy as np, pandas as pd, evallib as E, feats
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.preprocessing import FunctionTransformer
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True)
ft=FunctionTransformer(feats.phys_features)
def ridge(alpha=1.0):
    return make_pipeline(ft, StandardScaler(), Ridge(alpha=alpha))
for al in [0.3,1.0,3.0,10.0]:
    d,_=E.grouped_cv(lambda al=al: ridge(al), both, n_splits=10, seed=0)
    E.report(f'ridge-phys a={al}', d)
" 2>&1 | grep -v Warning

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.neural_network import MLPRegressor
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True)
ft=FunctionTransformer(feats.phys_features)
def mlp(h=(128,128),a=1e-3):
    return make_pipeline(ft, StandardScaler(), MLPRegressor(hidden_layer_sizes=h, alpha=a, max_iter=2000, random_state=0, early_stopping=False))
def et():
    return make_pipeline(ft, ExtraTreesRegressor(n_estimators=400,n_jobs=-1,random_state=0))
def hgb():
    return make_pipeline(ft, HistGradientBoostingRegressor(random_state=0, max_iter=500, learning_rate=0.05))
for name,mk in [('et-phys',et),('hgb-phys',hgb),('mlp-128x128',lambda:mlp((128,128)))]:
    d,_=E.grouped_cv(mk, both, n_splits=10, seed=0)
    E.report(name, d)
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
tr,va=E.load()
ft=FunctionTransformer(feats.phys_features)
def et():
    return make_pipeline(ft, ExtraTreesRegressor(n_estimators=600,n_jobs=-1,random_state=0))
def hgb():
    return make_pipeline(ft, HistGradientBoostingRegressor(random_state=0, max_iter=500, learning_rate=0.05))
for name,mk in [('et-phys',et),('hgb-phys',hgb)]:
    d=E.val_holdout(mk, tr, va)
    E.report('VAL '+name, d)
print('val n conditions:', E.condition_ids(va).nunique())
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True)
both['cond']=E.condition_ids(both)
# pick a condition with many rows
for cond, g in list(both.groupby('cond'))[:3]:
    g=g.sort_values('frequency')
    print('COND', cond, 'n=',len(g))
    for _,r in g.iterrows():
        print(f'  f={r.frequency:8.1f}  St={r.frequency*r.suction_side_displacement_thickness/r.free_stream_velocity:8.4f}  spl={r.scaled_sound_pressure_level:8.3f}')
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True)
both['cond']=E.condition_ids(both)
both['logSt']=np.log10(both.frequency*both.suction_side_displacement_thickness/both.free_stream_velocity)
res=[]
for cond,g in both.groupby('cond'):
    if len(g)<6: continue
    x=g.logSt.values; y=g.scaled_sound_pressure_level.values
    # fit cubic in logSt
    c=np.polyfit(x,y,min(4,len(g)-1)); yp=np.polyval(c,x)
    res.append(y-yp)
res=np.concatenate(res)
print('within-condition residual std after poly4(logSt):', np.std(res), 'max abs', np.max(np.abs(res)))
"

cd /app/scratch && cat > feats.py << 'EOF'
import numpy as np, pandas as pd
C_SOUND=343.0; NU=1.5e-5
def phys_features(X):
    f=X['frequency'].astype(float).values
    a=X['attack_angle'].astype(float).values
    c=X['chord_length'].astype(float).values
    U=X['free_stream_velocity'].astype(float).values
    d=X['suction_side_displacement_thickness'].astype(float).values
    logf=np.log10(f); logd=np.log10(d); logc=np.log10(c); logU=np.log10(U)
    M=U/C_SOUND; logM=np.log10(M)
    St=f*d/U; logSt=np.log10(St)
    logRe=np.log10(U*c/NU)
    logRed=np.log10(U*d/NU)
    out=pd.DataFrame(index=X.index)
    # base
    out['logf']=logf; out['a']=a; out['logc']=logc; out['logU']=logU; out['logd']=logd
    out['M']=M; out['logM']=logM
    # strouhal family
    out['logSt']=logSt; out['logSt2']=logSt**2; out['logSt3']=logSt**3; out['logSt4']=logSt**4
    out['logRe']=logRe; out['logRed']=logRed
    # angle family
    out['a2']=a**2; out['a3']=a**3; out['sqrt_a']=np.sqrt(np.abs(a))
    # interactions (peak-shift & level)
    out['logSt_logM']=logSt*logM
    out['logSt_a']=logSt*a
    out['logSt_logc']=logSt*logc
    out['logSt_logU']=logSt*logU
    out['logSt_logd']=logSt*logd
    out['a_logU']=a*logU
    out['a_logc']=a*logc
    out['logd_logU']=logd*logU
    out['logc_logU']=logc*logU
    out['a_logd']=a*logd
    return out
EOF
python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True)
ft=FunctionTransformer(feats.phys_features)
for ne in [400,800]:
  for mf in [0.5,0.7,1.0]:
    for ms in [1,2,3]:
        def mk(ne=ne,mf=mf,ms=ms):
            return make_pipeline(ft, ExtraTreesRegressor(n_estimators=ne,max_features=mf,min_samples_leaf=ms,n_jobs=-1,random_state=0))
        d,_=E.grouped_cv(mk, both, n_splits=10, seed=0)
        dv=E.val_holdout(mk, tr, va)
        print(f'ET ne={ne} mf={mf} ms={ms} | CV rmse={d[\"rmse\"]:.3f} p90={d[\"p90_cond\"]:.3f} max={d[\"max_cond\"]:.3f} | VAL rmse={dv[\"rmse\"]:.3f} mae={dv[\"mae\"]:.3f} p90={dv[\"p90_cond\"]:.3f} max={dv[\"max_cond\"]:.3f}')
"

restart: true

cd /app/scratch && timeout 250 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.kernel_ridge import KernelRidge
from sklearn.svm import SVR
tr,va=E.load()
ft=FunctionTransformer(feats.phys_features)
def kr(alpha,gamma):
    return make_pipeline(ft, StandardScaler(), KernelRidge(kernel='rbf', alpha=alpha, gamma=gamma))
def svr(C,g,e):
    return make_pipeline(ft, StandardScaler(), SVR(kernel='rbf', C=C, gamma=g, epsilon=e))
print('--- KernelRidge RBF (val holdout) ---')
for alpha in [0.01,0.03,0.1]:
  for gamma in [0.02,0.05,0.1]:
    dv=E.val_holdout(lambda: kr(alpha,gamma), tr, va)
    print(f'KR a={alpha} g={gamma} | VAL rmse={dv[\"rmse\"]:.3f} mae={dv[\"mae\"]:.3f} p90={dv[\"p90_cond\"]:.3f} p95={dv[\"p95_cond\"]:.3f} max={dv[\"max_cond\"]:.3f}')
"

cd /app/scratch && cat > mlpens.py << 'EOF'
import numpy as np
from sklearn.base import BaseEstimator, RegressorMixin, clone
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.neural_network import MLPRegressor
import feats

class MLPEnsemble(BaseEstimator, RegressorMixin):
    def __init__(self, hidden=(96,96), alpha=3e-3, n=6, max_iter=1500):
        self.hidden=hidden; self.alpha=alpha; self.n=n; self.max_iter=max_iter
    def fit(self, X, y):
        self.models_=[]
        for s in range(self.n):
            m=make_pipeline(FunctionTransformer(feats.phys_features), StandardScaler(),
                MLPRegressor(hidden_layer_sizes=self.hidden, alpha=self.alpha, max_iter=self.max_iter,
                    random_state=s, activation='tanh', solver='lbfgs'))
            m.fit(X, y); self.models_.append(m)
        return self
    def predict(self, X):
        return np.mean([m.predict(X) for m in self.models_], axis=0)
EOF
timeout 250 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings
warnings.filterwarnings('ignore')
from mlpens import MLPEnsemble
tr,va=E.load()
for hid in [(64,64),(96,96),(128,64)]:
  for al in [1e-3,3e-3]:
    dv=E.val_holdout(lambda hid=hid,al=al: MLPEnsemble(hidden=hid, alpha=al, n=6), tr, va)
    print(f'MLP-lbfgs {hid} a={al} | VAL rmse={dv[\"rmse\"]:.3f} mae={dv[\"mae\"]:.3f} p90={dv[\"p90_cond\"]:.3f} p95={dv[\"p95_cond\"]:.3f} max={dv[\"max_cond\"]:.3f}')
"

cd /app/scratch && timeout 200 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor
tr,va=E.load()
ft=FunctionTransformer(feats.phys_features)
m=make_pipeline(ft, ExtraTreesRegressor(n_estimators=600,max_features=1.0,min_samples_leaf=1,n_jobs=-1,random_state=0))
m.fit(tr[E.FEATURES], tr[E.TARGET].values)
va=va.copy(); va['pred']=m.predict(va[E.FEATURES]); va['err']=va.pred-va[E.TARGET]; va['cond']=E.condition_ids(va)
cr=va.groupby('cond').apply(lambda g: pd.Series({'rmse':np.sqrt(np.mean(g.err**2)),'n':len(g),'ang':g.attack_angle.mean(),'chord':g.chord_length.mean(),'vel':g.free_stream_velocity.mean()}))
cr=cr.sort_values('rmse',ascending=False)
pd.set_option('display.width',200)
print(cr)
print('worst rows:')
print(va.reindex(va.err.abs().sort_values(ascending=False).index)[['frequency','attack_angle','chord_length','free_stream_velocity','scaled_sound_pressure_level','pred','err']].head(10).to_string())
"

cd /app/scratch && cat > gam.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.preprocessing import SplineTransformer, StandardScaler
from sklearn.linear_model import Ridge
C_SOUND=343.0; NU=1.5e-5
def base_vars(X):
    f=X['frequency'].astype(float).values; a=X['attack_angle'].astype(float).values
    c=X['chord_length'].astype(float).values; U=X['free_stream_velocity'].astype(float).values
    d=X['suction_side_displacement_thickness'].astype(float).values
    logf=np.log10(f); logd=np.log10(d); logc=np.log10(c); logU=np.log10(U)
    M=U/C_SOUND; logM=np.log10(M); logSt=np.log10(f*d/U); logRe=np.log10(U*c/NU)
    return dict(f=f,a=a,c=c,U=U,d=d,logf=logf,logd=logd,logc=logc,logU=logU,M=M,logM=logM,logSt=logSt,logRe=logRe)

class GAM(BaseEstimator, RegressorMixin):
    def __init__(self, n_knots=10, degree=3, alpha=1.0, modulate=True):
        self.n_knots=n_knots; self.degree=degree; self.alpha=alpha; self.modulate=modulate
    def _design(self, X, fit=False):
        v=base_vars(X)
        logSt=v['logSt'].reshape(-1,1)
        if fit:
            self.spl_=SplineTransformer(n_knots=self.n_knots, degree=self.degree, extrapolation='continue')
            S=self.spl_.fit_transform(logSt)
        else:
            S=self.spl_.transform(logSt)
        a=v['a']; logU=v['logU']; logd=v['logd']; logc=v['logc']; logRe=v['logRe']; M=v['M']; logM=v['logM']
        # level block (smooth overall offset per condition)
        level=np.column_stack([a, a**2, a**3, logc, logU, logd, logRe, M,
                               a*logU, a*logc, a*logd, logc*logU, logd*logU, logc*logd, a*M])
        blocks=[S, level]
        if self.modulate:
            # let spectral shape shift/scale with Mach, angle, Re
            blocks.append(S*M.reshape(-1,1))
            blocks.append(S*a.reshape(-1,1))
            blocks.append(S*logRe.reshape(-1,1))
        D=np.hstack(blocks)
        return D
    def fit(self, X, y):
        D=self._design(X, fit=True)
        self.sc_=StandardScaler().fit(D)
        self.rid_=Ridge(alpha=self.alpha).fit(self.sc_.transform(D), y)
        return self
    def predict(self, X):
        D=self._design(X, fit=False)
        return self.rid_.predict(self.sc_.transform(D))
EOF
timeout 200 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings
warnings.filterwarnings('ignore')
from gam import GAM
tr,va=E.load()
for nk in [8,12,16]:
  for al in [0.3,1.0,3.0]:
    dv=E.val_holdout(lambda nk=nk,al=al: GAM(n_knots=nk, alpha=al), tr, va)
    print(f'GAM nk={nk} a={al} | VAL rmse={dv[\"rmse\"]:.3f} mae={dv[\"mae\"]:.3f} p90={dv[\"p90_cond\"]:.3f} p95={dv[\"p95_cond\"]:.3f} max={dv[\"max_cond\"]:.3f}')
"

cd /app/scratch && timeout 260 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor, RandomForestRegressor
from mlpens import MLPEnsemble
tr,va=E.load()
ft=FunctionTransformer(feats.phys_features)
def fit_pred(m):
    m.fit(tr[E.FEATURES], tr[E.TARGET].values); return m.predict(va[E.FEATURES])
et=make_pipeline(ft, ExtraTreesRegressor(n_estimators=600,max_features=1.0,min_samples_leaf=1,n_jobs=-1,random_state=0))
hgb=make_pipeline(ft, HistGradientBoostingRegressor(random_state=0,max_iter=600,learning_rate=0.03,max_leaf_nodes=31,l2_regularization=1.0))
rf=make_pipeline(ft, RandomForestRegressor(n_estimators=600,max_features=0.5,n_jobs=-1,random_state=0))
mlp=MLPEnsemble(hidden=(128,64),alpha=3e-3,n=6)
P={}
for nm,m in [('et',et),('hgb',hgb),('rf',rf),('mlp',mlp)]:
    P[nm]=fit_pred(m)
y=va[E.TARGET].values; cond=E.condition_ids(va).values
def rep(nm,pred):
    dfp=pd.DataFrame({'cond':cond,'err':pred-y}); d=E.metrics_from(dfp)
    print(f'{nm:22s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
for nm in P: rep(nm,P[nm])
print('--- blends ---')
rep('et+hgb', 0.5*P['et']+0.5*P['hgb'])
rep('et+rf', 0.5*P['et']+0.5*P['rf'])
rep('et+hgb+rf', (P['et']+P['hgb']+P['rf'])/3)
rep('et+hgb+mlp', (P['et']+P['hgb']+P['mlp'])/3)
rep('et.5 hgb.3 mlp.2', 0.5*P['et']+0.3*P['hgb']+0.2*P['mlp'])
rep('all4', (P['et']+P['hgb']+P['rf']+P['mlp'])/4)
"

cd /app/scratch && cat >> feats.py << 'EOF'

def phys_plus_raw(X):
    out=phys_features(X)
    out['r_frequency']=X['frequency'].astype(float).values
    out['r_chord']=X['chord_length'].astype(float).values
    out['r_U']=X['free_stream_velocity'].astype(float).values
    out['r_thick']=X['suction_side_displacement_thickness'].astype(float).values
    return out
EOF
cat > logo.py << 'EOF'
import numpy as np, pandas as pd, evallib as E
from sklearn.model_selection import LeaveOneGroupOut
def logo_eval(make_model, data, groups):
    logo=LeaveOneGroupOut(); rows=[]
    for tri,tei in logo.split(data, data[E.TARGET], groups):
        Xtr=data.iloc[tri]; Xte=data.iloc[tei]
        m=make_model(); m.fit(Xtr[E.FEATURES], Xtr[E.TARGET].values)
        rows.append(pd.DataFrame({'cond':groups[tei],'err':m.predict(Xte[E.FEATURES])-Xte[E.TARGET].values}))
    return E.metrics_from(pd.concat(rows,ignore_index=True))
EOF
timeout 280 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings, logo
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True)
groups=E.condition_ids(both).values
ftp=FunctionTransformer(feats.phys_features); ftr=FunctionTransformer(feats.phys_plus_raw)
for nm,ft,ms in [('phys ms1',ftp,1),('phys ms2',ftp,2),('phys+raw ms1',ftr,1)]:
    def mk(ft=ft,ms=ms): return make_pipeline(ft, ExtraTreesRegressor(n_estimators=400,max_features=1.0,min_samples_leaf=ms,n_jobs=-1,random_state=0))
    d=logo.logo_eval(mk, both, groups)
    print(f'LOGO ET {nm:14s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f} nc={d[\"n_cond\"]}')
"

cd /app/scratch && timeout 280 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.model_selection import LeaveOneGroupOut
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
groups=E.condition_ids(both).values
ft=FunctionTransformer(feats.phys_features)
logo=LeaveOneGroupOut(); recs=[]
both['logSt']=np.log10(both.frequency*both.suction_side_displacement_thickness/both.free_stream_velocity)
allerr=np.zeros(len(both))
for tri,tei in logo.split(both, both[E.TARGET], groups):
    m=make_pipeline(ft, ExtraTreesRegressor(n_estimators=250,max_features=1.0,n_jobs=-1,random_state=0))
    m.fit(both.iloc[tri][E.FEATURES], both.iloc[tri][E.TARGET].values)
    allerr[tei]=m.predict(both.iloc[tei][E.FEATURES])-both.iloc[tei][E.TARGET].values
both['err']=allerr
cr=both.groupby('cond' if 'cond' in both else E.condition_ids(both)).apply(lambda g: pd.Series({
    'rmse':np.sqrt(np.mean(g.err**2)),'bias':g.err.mean(),'shapestd':g.err.std(),'n':len(g),
    'ang':g.attack_angle.mean(),'chord':g.chord_length.mean(),'vel':g.free_stream_velocity.mean()}))
cr=cr.sort_values('rmse',ascending=False)
pd.set_option('display.width',220)
print('Top 12 worst conditions (LOGO):')
print(cr.head(12).to_string())
print()
print('fraction of per-cond RMSE^2 explained by bias: %.2f'%((cr.bias**2).sum()/(cr.rmse**2).sum()))
print('overall mae=%.3f rmse=%.3f'%(both.err.abs().mean(), np.sqrt((both.err**2).mean())))
" 

cd /app/scratch && timeout 280 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.model_selection import LeaveOneGroupOut
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both)
groups=both['cond'].values
ft=FunctionTransformer(feats.phys_features)
logo=LeaveOneGroupOut(); allpred=np.zeros(len(both))
for tri,tei in logo.split(both, both[E.TARGET], groups):
    m=make_pipeline(ft, ExtraTreesRegressor(n_estimators=300,max_features=1.0,n_jobs=-1,random_state=0))
    m.fit(both.iloc[tri][E.FEATURES], both.iloc[tri][E.TARGET].values)
    allpred[tei]=m.predict(both.iloc[tei][E.FEATURES])
both['pred']=allpred; both['err']=both.pred-both[E.TARGET]
both['logSt']=np.log10(both.frequency*both.suction_side_displacement_thickness/both.free_stream_velocity)
both.to_csv('logo_err.csv',index=False)
print('saved. overall mae=%.3f rmse=%.3f'%(both.err.abs().mean(), np.sqrt((both.err**2).mean())))
# residual vs logSt bins
both['stbin']=pd.qcut(both.logSt,12)
print(both.groupby('stbin').err.agg(['mean','std','count']).to_string())
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd
b=pd.read_csv('logo_err.csv')
print('=== residual vs attack angle bins ===')
b['abin']=pd.cut(b.attack_angle,[-1,1,3,5,8,11,14,18,23])
print(b.groupby('abin').err.agg(['mean','std','count']).to_string())
print()
print('=== residual vs velocity ===')
print(b.groupby(np.round(b.free_stream_velocity,0)).err.agg(['mean','std','count']).to_string())
print()
print('=== rmse by chord ===')
print(b.groupby(np.round(b.chord_length,4)).err.agg(lambda e:np.sqrt((e**2).mean())).to_string())
"

cd /app/scratch && cat > fastcv.py << 'EOF'
import numpy as np, pandas as pd, evallib as E
from sklearn.model_selection import GroupKFold
def fast_eval(make_model, data, groups, n_splits=20, seed=0):
    uniq=np.unique(groups); rng=np.random.RandomState(seed); perm=rng.permutation(len(uniq))
    gmap={u:i for i,u in enumerate(uniq[perm])}; gidx=np.array([gmap[g] for g in groups])
    gk=GroupKFold(n_splits=n_splits); rows=[]
    for tri,tei in gk.split(data, data[E.TARGET], gidx):
        m=make_model(); m.fit(data.iloc[tri][E.FEATURES], data.iloc[tri][E.TARGET].values)
        rows.append(pd.DataFrame({'cond':groups[tei],'err':m.predict(data.iloc[tei][E.FEATURES])-data.iloc[tei][E.TARGET].values}))
    return E.metrics_from(pd.concat(rows,ignore_index=True))
EOF
timeout 280 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, warnings, fastcv
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
ft=FunctionTransformer(feats.phys_features)
def mk(ms=1,ne=400,mf=1.0,boot=False):
    return make_pipeline(ft, ExtraTreesRegressor(n_estimators=ne,max_features=mf,min_samples_leaf=ms,bootstrap=boot,n_jobs=-1,random_state=0))
for cfg in [dict(),dict(ne=800),dict(mf=0.7),dict(boot=True),dict(boot=True,mf=0.7)]:
    d=fastcv.fast_eval(lambda cfg=cfg: mk(**cfg), both, groups, n_splits=20)
    print(f'ET {str(cfg):30s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
"

cd /app/scratch && cat > feats2.py << 'EOF'
import numpy as np, pandas as pd
C=343.0; NU=1.5e-5
def make(X, use_angle=True, extra=True):
    f=X['frequency'].astype(float).values; a=X['attack_angle'].astype(float).values
    c=X['chord_length'].astype(float).values; U=X['free_stream_velocity'].astype(float).values
    d=X['suction_side_displacement_thickness'].astype(float).values
    logf=np.log10(f); logd=np.log10(d); logc=np.log10(c); logU=np.log10(U)
    M=U/C; logSt=np.log10(f*d/U); logRe=np.log10(U*c/NU); logRed=np.log10(U*d/NU)
    dc=d/c  # thickness-to-chord
    o=pd.DataFrame(index=X.index)
    o['logf']=logf; o['logc']=logc; o['logU']=logU; o['logd']=logd; o['M']=M
    o['logSt']=logSt; o['logSt2']=logSt**2; o['logSt3']=logSt**3
    o['logRe']=logRe; o['logRed']=logRed; o['logdc']=np.log10(dc)
    o['logSt_logM']=logSt*np.log10(M); o['logSt_logc']=logSt*logc
    o['logSt_logU']=logSt*logU; o['logSt_logd']=logSt*logd; o['logSt_logRe']=logSt*logRe
    o['logc_logU']=logc*logU; o['logd_logU']=logd*logU; o['logc_logd']=logc*logd
    if use_angle:
        o['a']=a; o['a2']=a**2; o['a3']=a**3
        o['a_logU']=a*logU; o['a_logc']=a*logc; o['a_logd']=a*logd
        o['logSt_a']=logSt*a
    if extra:
        o['logSt_logdc']=logSt*np.log10(dc)
        o['logRed2']=logRed**2
    return o
def with_angle(X): return make(X, True, True)
def no_angle(X): return make(X, False, True)
EOF
timeout 280 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats2, warnings, fastcv
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
for nm,fn in [('with_angle',feats2.with_angle),('no_angle',feats2.no_angle)]:
    ft=FunctionTransformer(fn)
    def mk(ft=ft): return make_pipeline(ft, ExtraTreesRegressor(n_estimators=500,max_features=0.7,n_jobs=-1,random_state=0))
    d=fastcv.fast_eval(mk, both, groups, n_splits=20)
    print(f'{nm:12s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings
warnings.filterwarnings('ignore')
b=pd.read_csv('logo_err.csv')  # has cond, pred (LOGO ET), logSt, target
y=b[E.TARGET].values
def metrics(pred):
    dfp=pd.DataFrame({'cond':b['cond'].values,'err':pred-y}); return E.metrics_from(dfp)
def rep(nm,pred):
    d=metrics(pred); print(f'{nm:22s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
rep('ET raw', b.pred.values)
# within-condition polynomial smoothing of the ET predictions
for deg in [2,3,4,5]:
    sm=b.pred.values.copy()
    for c,g in b.groupby('cond'):
        idx=g.index.values; x=g.logSt.values; p=g.pred.values
        dd=min(deg,len(g)-1)
        if dd>=1:
            coef=np.polyfit(x,p,dd); sm[idx]=np.polyval(coef,x)
    rep(f'ET smooth poly{deg}', sm)
"

cd /app/scratch && cat > gen_oof.py << 'EOF'
import numpy as np, pandas as pd, evallib as E, feats, feats2, warnings
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.model_selection import GroupKFold
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
uniq=np.unique(groups); rng=np.random.RandomState(1); perm=rng.permutation(len(uniq))
gmap={u:i for i,u in enumerate(uniq[perm])}; gidx=np.array([gmap[g] for g in groups])
gk=GroupKFold(n_splits=15)
def ET(fn,mf=0.7,ne=400,ms=1): return make_pipeline(FunctionTransformer(fn), ExtraTreesRegressor(n_estimators=ne,max_features=mf,min_samples_leaf=ms,n_jobs=-1,random_state=0))
def HGB(fn): return make_pipeline(FunctionTransformer(fn), HistGradientBoostingRegressor(random_state=0,max_iter=600,learning_rate=0.03,max_leaf_nodes=31,l2_regularization=1.0,min_samples_leaf=15))
bases={
 'et_phys': lambda: ET(feats.phys_features,1.0),
 'et_phys07': lambda: ET(feats.phys_features,0.7),
 'et_f2': lambda: ET(feats2.with_angle,0.7),
 'et_noang': lambda: ET(feats2.no_angle,0.7),
 'hgb_phys': lambda: HGB(feats.phys_features),
 'hgb_f2': lambda: HGB(feats2.with_angle),
}
oof={k:np.zeros(len(both)) for k in bases}
for tri,tei in gk.split(both, both[E.TARGET], gidx):
    for k,mkf in bases.items():
        m=mkf(); m.fit(both.iloc[tri][E.FEATURES], both.iloc[tri][E.TARGET].values)
        oof[k][tei]=m.predict(both.iloc[tei][E.FEATURES])
out=pd.DataFrame(oof); out['cond']=groups; out['y']=both[E.TARGET].values
out.to_csv('oof.csv',index=False)
print('done', out.shape)
EOF
nohup timeout 290 python3 gen_oof.py > gen_oof.log 2>&1 &
echo "started"; sleep 5; cat gen_oof.log

cd /app/scratch && sleep 120 && cat gen_oof.log && ls -la oof.csv 2>/dev/null

cd /app/scratch && sleep 90 && cat gen_oof.log && echo "---" && ls -la oof.csv 2>/dev/null

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E, itertools
o=pd.read_csv('oof.csv'); y=o.y.values; cond=o['cond'].values
cols=['et_phys','et_phys07','et_f2','et_noang','hgb_phys','hgb_f2']
def metrics(pred):
    return E.metrics_from(pd.DataFrame({'cond':cond,'err':pred-y}))
def rep(nm,pred):
    d=metrics(pred); print(f'{nm:30s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
for c in cols: rep(c, o[c].values)
print('--- simple blends ---')
rep('all ET avg', o[['et_phys','et_phys07','et_f2']].mean(1).values)
rep('all ET+noang', o[['et_phys','et_phys07','et_f2','et_noang']].mean(1).values)
rep('ET3+HGB2', (o[['et_phys','et_phys07','et_f2']].mean(1)*0.7+o[['hgb_phys','hgb_f2']].mean(1)*0.3).values)
rep('all6 avg', o[cols].mean(1).values)
rep('ET-heavy', (o[['et_phys','et_phys07','et_f2','et_noang']].mean(1)*0.8+o[['hgb_phys','hgb_f2']].mean(1)*0.2).values)
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
o=pd.read_csv('oof.csv'); y=o.y.values; cond=o['cond'].values
cols=['et_phys','et_phys07','et_f2','et_noang','hgb_phys','hgb_f2']
pred=o[cols].mean(1).values
o['err']=pred-y
cr=o.groupby('cond').err.apply(lambda e: np.sqrt(np.mean(e**2))).sort_values(ascending=False)
print('n conditions:', len(cr))
print('percentiles of per-cond rmse: p50=%.3f p75=%.3f p90=%.3f p95=%.3f max=%.3f'%(
    np.percentile(cr,50),np.percentile(cr,75),np.percentile(cr,90),np.percentile(cr,95),cr.max()))
print('# conditions > 2.35:', (cr>2.35).sum(), 'out of', len(cr))
print('# conditions > 3.0:', (cr>3.0).sum())
print()
print('worst 10:'); print(cr.head(10).to_string())
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
o=pd.read_csv('oof.csv')
# parse cond string a_c_v
parts=o['cond'].str.split('_',expand=True).astype(float); parts.columns=['a','c','v']
o=pd.concat([o,parts],axis=1)
condtab=o.drop_duplicates('cond')[['cond','a','c','v']]
print('conditions per nominal angle:')
print(condtab.groupby('a').size().to_string())
print()
print('conditions per (angle) with which chords exist:')
for a in sorted(condtab.a.unique()):
    sub=condtab[condtab.a==a]
    print(f'  ang {a:5.1f}: {len(sub):2d} conds, chords={sorted(set(sub.c))}, vels={sorted(set(sub.v))}')
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
o=pd.read_csv('oof.csv'); y=o.y.values
cols=['et_phys','et_phys07','et_f2','et_noang','hgb_phys','hgb_f2']
o['pred']=o[cols].mean(1).values; o['err']=o.pred-y
# per-condition stats
g=o.groupby('cond')
cr=g.err.apply(lambda e: np.sqrt(np.mean(e**2)))
rows_by_cond={c:sub.err.values for c,sub in g}
conds=list(cr.index)
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,max=5.60)
rng=np.random.RandomState(0)
for size in [15,20,25,30]:
    cnt={k:0 for k in T}; passall=0; N=4000
    for _ in range(N):
        S=rng.choice(conds,size=size,replace=False)
        errs=np.concatenate([rows_by_cond[c] for c in S])
        crr=np.array([cr[c] for c in S])
        m=dict(rmse=np.sqrt(np.mean(errs**2)),mae=np.mean(np.abs(errs)),
               p90=np.percentile(crr,90),p95=np.percentile(crr,95),max=crr.max())
        ok=True
        for k in T:
            if m[k]<=T[k]: cnt[k]+=1
            else: ok=False
        passall+=ok
    print(f'subset={size}: '+' '.join(f'{k}={cnt[k]/N:.2f}' for k in T)+f' | ALL={passall/N:.3f}')
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
from scipy.interpolate import UnivariateSpline
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True)
both['cond']=E.condition_ids(both)
both['logSt']=np.log10(both.frequency*both.suction_side_displacement_thickness/both.free_stream_velocity)
res=[]
res_sm=[]
for c,g in both.groupby('cond'):
    if len(g)<8: continue
    g=g.sort_values('logSt'); x=g.logSt.values; y=g.scaled_sound_pressure_level.values
    # smoothing spline with small s -> flexible
    try:
        sp=UnivariateSpline(x,y,k=4,s=len(x)*0.15)  # moderate smoothing
        res_sm.append(y-sp(x))
    except Exception as ex:
        pass
    # local second difference based noise estimate
    d2=y[2:]-2*y[1:-1]+y[:-2]
    res.append(d2/np.sqrt(6))  # scaled 2nd diff ~ noise if trend smooth
rs=np.concatenate(res); rsm=np.concatenate(res_sm)
print('2nd-diff noise estimate std: %.3f'%np.std(rs))
print('flexible spline residual std: %.3f  max abs %.3f'%(np.std(rsm), np.max(np.abs(rsm))))
print('fraction |spline resid|>1.5: %.3f'%np.mean(np.abs(rsm)>1.5))
"

cd /app/scratch && cat > gam2.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.preprocessing import SplineTransformer, StandardScaler
from sklearn.linear_model import Ridge
C=343.0; NU=1.5e-5
def bv(X):
    f=X['frequency'].astype(float).values; a=X['attack_angle'].astype(float).values
    c=X['chord_length'].astype(float).values; U=X['free_stream_velocity'].astype(float).values
    d=X['suction_side_displacement_thickness'].astype(float).values
    return dict(a=a,logc=np.log10(c),logU=np.log10(U),logd=np.log10(d),M=U/C,
                logSt=np.log10(f*d/U),logRe=np.log10(U*c/NU),logRed=np.log10(U*d/NU))
class GAM2(BaseEstimator, RegressorMixin):
    def __init__(self, n_knots=12, degree=3, alpha=3.0):
        self.n_knots=n_knots; self.degree=degree; self.alpha=alpha
    def _design(self, X, fit=False):
        v=bv(X); St=v['logSt'].reshape(-1,1)
        if fit:
            self.spl_=SplineTransformer(n_knots=self.n_knots, degree=self.degree, extrapolation='linear')
            S=self.spl_.fit_transform(St)
        else: S=self.spl_.transform(St)
        a=v['a']; logU=v['logU']; logd=v['logd']; logc=v['logc']; logRe=v['logRe']; M=v['M']; logRed=v['logRed']
        level=np.column_stack([a,a**2,a**3,logc,logU,logd,logRe,logRed,M,
            a*logU,a*logc,a*logd,logc*logU,logd*logU,logc*logd,a*M,logRe*a])
        mod=np.column_stack([a,M,logRe,logc,logd,logU])  # shape modulators
        blocks=[S, level]
        for j in range(mod.shape[1]):
            blocks.append(S*mod[:,j:j+1])
        return np.hstack(blocks)
    def fit(self, X, y):
        D=self._design(X, True); self.sc_=StandardScaler().fit(D)
        self.rid_=Ridge(alpha=self.alpha).fit(self.sc_.transform(D), y); return self
    def predict(self, X):
        return self.rid_.predict(self.sc_.transform(self._design(X, False)))
EOF
timeout 120 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings, fastcv
warnings.filterwarnings('ignore')
from gam2 import GAM2
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
for nk in [8,12,16]:
  for al in [1.0,3.0,10.0]:
    d=fastcv.fast_eval(lambda nk=nk,al=al: GAM2(n_knots=nk,alpha=al), both, groups, n_splits=15)
    print(f'GAM2 nk={nk} a={al:4.1f} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
"

cd /app/scratch && timeout 200 python3 -c "
import numpy as np, pandas as pd, evallib as E, feats, feats2, warnings, fastcv
warnings.filterwarnings('ignore')
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.neighbors import KNeighborsRegressor
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
for k in [5,8,12]:
  for w in ['distance']:
    def mk(k=k,w=w): return make_pipeline(FunctionTransformer(feats2.with_angle), StandardScaler(), KNeighborsRegressor(n_neighbors=k,weights=w))
    d=fastcv.fast_eval(mk, both, groups, n_splits=15)
    print(f'KNN k={k} {w:8s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
"

cd /app/scratch && cat > twostage.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.base import BaseEstimator, RegressorMixin, clone
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.linear_model import Ridge
from sklearn.ensemble import ExtraTreesRegressor
import feats, feats2
class TwoStage(BaseEstimator, RegressorMixin):
    def __init__(self, feat_fn=feats2.with_angle, alpha=1.0, ne=500, mf=0.7, ms=1):
        self.feat_fn=feat_fn; self.alpha=alpha; self.ne=ne; self.mf=mf; self.ms=ms
    def fit(self, X, y):
        self.lin_=make_pipeline(FunctionTransformer(self.feat_fn), StandardScaler(), Ridge(alpha=self.alpha))
        self.lin_.fit(X,y)
        r=y-self.lin_.predict(X)
        self.tree_=make_pipeline(FunctionTransformer(self.feat_fn), ExtraTreesRegressor(
            n_estimators=self.ne,max_features=self.mf,min_samples_leaf=self.ms,n_jobs=-1,random_state=0))
        self.tree_.fit(X,r); return self
    def predict(self, X):
        return self.lin_.predict(X)+self.tree_.predict(X)
EOF
timeout 200 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings, fastcv, feats, feats2
warnings.filterwarnings('ignore')
from twostage import TwoStage
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
for al in [0.3,1.0,3.0,10.0]:
    d=fastcv.fast_eval(lambda al=al: TwoStage(feats2.with_angle, alpha=al), both, groups, n_splits=15)
    print(f'TwoStage f2 a={al:5.1f} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
"

cd /app/scratch && cat >> twostage.py << 'EOF'

C=343.0; NU=1.5e-5
def level_feats(X):
    c=X['chord_length'].astype(float).values; U=X['free_stream_velocity'].astype(float).values
    d=X['suction_side_displacement_thickness'].astype(float).values
    a=X['attack_angle'].astype(float).values
    o=pd.DataFrame(index=X.index)
    o['logc']=np.log10(c); o['logU']=np.log10(U); o['logd']=np.log10(d)
    o['logRe']=np.log10(U*c/NU); o['logRed']=np.log10(U*d/NU); o['M']=U/C
    o['a']=a
    return o
class TwoStageSafe(BaseEstimator, RegressorMixin):
    def __init__(self, tree_feat=feats2.with_angle, alpha=5.0, ne=500, mf=0.7, ms=1):
        self.tree_feat=tree_feat; self.alpha=alpha; self.ne=ne; self.mf=mf; self.ms=ms
    def fit(self, X, y):
        self.lin_=make_pipeline(FunctionTransformer(level_feats), StandardScaler(), Ridge(alpha=self.alpha))
        self.lin_.fit(X,y); r=y-self.lin_.predict(X)
        self.tree_=make_pipeline(FunctionTransformer(self.tree_feat), ExtraTreesRegressor(
            n_estimators=self.ne,max_features=self.mf,min_samples_leaf=self.ms,n_jobs=-1,random_state=0))
        self.tree_.fit(X,r); return self
    def predict(self, X):
        return self.lin_.predict(X)+self.tree_.predict(X)
EOF
timeout 200 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings, fastcv, feats2
warnings.filterwarnings('ignore')
from twostage import TwoStageSafe
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
for al in [1.0,3.0,10.0,30.0]:
    d=fastcv.fast_eval(lambda al=al: TwoStageSafe(feats2.with_angle, alpha=al), both, groups, n_splits=15)
    print(f'TwoStageSafe a={al:5.1f} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
"

cd /app/scratch && cat >> twostage.py << 'EOF'

from sklearn.ensemble import HistGradientBoostingRegressor
class TwoStageGen(BaseEstimator, RegressorMixin):
    """Safe linear level + arbitrary tree stage on residual."""
    def __init__(self, tree_feat=feats2.with_angle, alpha=20.0, kind='et', ne=500, mf=0.7, ms=1):
        self.tree_feat=tree_feat; self.alpha=alpha; self.kind=kind; self.ne=ne; self.mf=mf; self.ms=ms
    def _tree(self):
        if self.kind=='et':
            return ExtraTreesRegressor(n_estimators=self.ne,max_features=self.mf,min_samples_leaf=self.ms,n_jobs=-1,random_state=0)
        if self.kind=='hgb':
            return HistGradientBoostingRegressor(random_state=0,max_iter=700,learning_rate=0.03,max_leaf_nodes=31,l2_regularization=1.0,min_samples_leaf=20)
    def fit(self, X, y):
        self.lin_=make_pipeline(FunctionTransformer(level_feats), StandardScaler(), Ridge(alpha=self.alpha))
        self.lin_.fit(X,y); r=y-self.lin_.predict(X)
        self.tree_=make_pipeline(FunctionTransformer(self.tree_feat), self._tree()); self.tree_.fit(X,r); return self
    def predict(self, X):
        return self.lin_.predict(X)+self.tree_.predict(X)
EOF
cat > gen_oof2.py << 'EOF'
import numpy as np, pandas as pd, evallib as E, feats, feats2, warnings
warnings.filterwarnings('ignore')
from sklearn.model_selection import GroupKFold
from twostage import TwoStageGen
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
uniq=np.unique(groups); rng=np.random.RandomState(1); perm=rng.permutation(len(uniq))
gmap={u:i for i,u in enumerate(uniq[perm])}; gidx=np.array([gmap[g] for g in groups])
gk=GroupKFold(n_splits=15)
bases={
 'ts_et_f2':   lambda: TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=0.7),
 'ts_et_f2b':  lambda: TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=1.0),
 'ts_et_phys': lambda: TwoStageGen(feats.phys_features, alpha=30, kind='et', mf=0.7),
 'ts_et_noang':lambda: TwoStageGen(feats2.no_angle, alpha=30, kind='et', mf=0.7),
 'ts_hgb_f2':  lambda: TwoStageGen(feats2.with_angle, alpha=30, kind='hgb'),
 'ts_hgb_phys':lambda: TwoStageGen(feats.phys_features, alpha=30, kind='hgb'),
}
oof={k:np.zeros(len(both)) for k in bases}
for tri,tei in gk.split(both, both[E.TARGET], gidx):
    for k,mkf in bases.items():
        m=mkf(); m.fit(both.iloc[tri][E.FEATURES], both.iloc[tri][E.TARGET].values)
        oof[k][tei]=m.predict(both.iloc[tei][E.FEATURES])
out=pd.DataFrame(oof); out['cond']=groups; out['y']=both[E.TARGET].values
out.to_csv('oof2.csv',index=False); print('done', out.shape)
EOF
nohup timeout 290 python3 gen_oof2.py > gen_oof2.log 2>&1 &
echo started; sleep 4; cat gen_oof2.log

cd /app/scratch && sleep 180 && cat gen_oof2.log && ls -la oof2.csv 2>/dev/null

cd /app/scratch && sleep 60 && cat gen_oof2.log && echo "---" && ls -la oof2.csv 2>/dev/null

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
o=pd.read_csv('oof2.csv'); y=o.y.values; cond=o['cond'].values
cols=['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']
def metrics(p): return E.metrics_from(pd.DataFrame({'cond':cond,'err':p-y}))
def rep(nm,p):
    d=metrics(p); print(f'{nm:26s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
for c in cols: rep(c,o[c].values)
print('--- blends ---')
rep('all6', o[cols].mean(1).values)
rep('et4', o[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang']].mean(1).values)
rep('et4+hgb2', (o[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang']].mean(1)*0.7 + o[['ts_hgb_f2','ts_hgb_phys']].mean(1)*0.3).values)
rep('et-f2+hgb', (o[['ts_et_f2','ts_et_f2b']].mean(1)*0.6 + o[['ts_hgb_f2','ts_hgb_phys']].mean(1)*0.4).values)
rep('4way', o[['ts_et_f2','ts_et_phys','ts_hgb_f2','ts_hgb_phys']].mean(1).values)
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
o2=pd.read_csv('oof2.csv')
y=o2.y.values; cond=o2['cond'].values
# candidate blends (from two-stage models)
cand={
 'all6': o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1).values,
 'et4hgb2': (o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang']].mean(1)*0.7+o2[['ts_hgb_f2','ts_hgb_phys']].mean(1)*0.3).values,
 'etf2_hgb': (o2[['ts_et_f2','ts_et_f2b']].mean(1)*0.6+o2[['ts_hgb_f2','ts_hgb_phys']].mean(1)*0.4).values,
 'et_f2_only': o2[['ts_et_f2','ts_et_f2b']].mean(1).values,
}
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,max=5.60)
for nm,pred in cand.items():
    o2['err']=pred-y; g=o2.groupby('cond'); cr=g.err.apply(lambda e:np.sqrt(np.mean(e**2)))
    rows_by={c:s.err.values for c,s in g}; conds=list(cr.index)
    rng=np.random.RandomState(0)
    for size in [20,25]:
        cnt={k:0 for k in T}; pa=0; N=3000
        for _ in range(N):
            S=rng.choice(conds,size=size,replace=False)
            errs=np.concatenate([rows_by[c] for c in S]); crr=np.array([cr[c] for c in S])
            m=dict(rmse=np.sqrt(np.mean(errs**2)),mae=np.mean(np.abs(errs)),p90=np.percentile(crr,90),p95=np.percentile(crr,95),max=crr.max())
            ok=True
            for k in T:
                if m[k]<=T[k]: cnt[k]+=1
                else: ok=False
            pa+=ok
        print(f'{nm:12s} sz={size}: '+' '.join(f'{k}={cnt[k]/N:.2f}' for k in T)+f' ALL={pa/N:.3f}')
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
o1=pd.read_csv('oof.csv'); o2=pd.read_csv('oof2.csv')
assert (o1['cond'].values==o2['cond'].values).all() and (o1.y.values==o2.y.values).all()
y=o2.y.values; cond=o2['cond'].values
pred=o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1).values
df=pd.DataFrame({'cond':cond,'err':pred-y,
  'a':pd.Series(cond).str.split('_').str[0].astype(float).values,
  'c':pd.Series(cond).str.split('_').str[1].astype(float).values,
  'v':pd.Series(cond).str.split('_').str[2].astype(float).values})
cr=df.groupby('cond').apply(lambda g: pd.Series({'rmse':np.sqrt(np.mean(g.err**2)),'bias':g.err.mean(),'n':len(g),'a':g.a.iloc[0],'c':g.c.iloc[0],'v':g.v.iloc[0]}))
cr=cr.sort_values('rmse',ascending=False)
pd.set_option('display.width',200)
print('worst 14 (two-stage all6 blend):'); print(cr.head(14).to_string())
print('count > 2.35:', (cr.rmse>2.35).sum(), ' >2.6:', (cr.rmse>2.6).sum())
"

cd /app/scratch && cat >> twostage.py << 'EOF'

from sklearn.preprocessing import SplineTransformer
def _sbv(X):
    f=X['frequency'].astype(float).values; U=X['free_stream_velocity'].astype(float).values
    d=X['suction_side_displacement_thickness'].astype(float).values; c=X['chord_length'].astype(float).values
    a=X['attack_angle'].astype(float).values
    return (np.log10(f*d/U), U/343.0, np.log10(c), np.log10(U), np.log10(d), np.log10(U*c/1.5e-5), a)
class BackboneTree(BaseEstimator, RegressorMixin):
    """Stage1: smooth extrapolating backbone (level + Strouhal spline + Mach/Re mod). Stage2: ET residual."""
    def __init__(self, tree_feat=feats2.with_angle, alpha=20.0, n_knots=8, mf=0.7, ne=500):
        self.tree_feat=tree_feat; self.alpha=alpha; self.n_knots=n_knots; self.mf=mf; self.ne=ne
    def _design(self, X, fit=False):
        logSt,M,logc,logU,logd,logRe,a=_sbv(X); St=logSt.reshape(-1,1)
        if fit:
            self.spl_=SplineTransformer(n_knots=self.n_knots,degree=3,extrapolation='linear')
            S=self.spl_.fit_transform(St)
        else: S=self.spl_.transform(St)
        level=np.column_stack([logc,logU,logd,logRe,M])
        blocks=[S, S*M.reshape(-1,1), S*logU.reshape(-1,1), level]
        return np.hstack(blocks)
    def fit(self, X, y):
        D=self._design(X,True); self.sc_=StandardScaler().fit(D)
        self.lin_=Ridge(alpha=self.alpha).fit(self.sc_.transform(D), y)
        r=y-self.lin_.predict(self.sc_.transform(D))
        self.tree_=make_pipeline(FunctionTransformer(self.tree_feat),
            ExtraTreesRegressor(n_estimators=self.ne,max_features=self.mf,n_jobs=-1,random_state=0))
        self.tree_.fit(X,r); return self
    def predict(self, X):
        D=self._design(X,False); return self.lin_.predict(self.sc_.transform(D))+self.tree_.predict(X)
EOF
timeout 220 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings, fastcv, feats2
warnings.filterwarnings('ignore')
from twostage import BackboneTree
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
for nk in [6,8]:
  for al in [10.0,30.0,100.0]:
    d=fastcv.fast_eval(lambda nk=nk,al=al: BackboneTree(feats2.with_angle,alpha=al,n_knots=nk), both, groups, n_splits=15)
    print(f'Backbone nk={nk} a={al:5.1f} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
"

cd /app/scratch && cat > gen_oof3.py << 'EOF'
import numpy as np, pandas as pd, evallib as E, feats, feats2, warnings
warnings.filterwarnings('ignore')
from sklearn.model_selection import GroupKFold
from twostage import BackboneTree
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
uniq=np.unique(groups); rng=np.random.RandomState(1); perm=rng.permutation(len(uniq))
gmap={u:i for i,u in enumerate(uniq[perm])}; gidx=np.array([gmap[g] for g in groups])
gk=GroupKFold(n_splits=15)
bases={
 'bb_f2':   lambda: BackboneTree(feats2.with_angle, alpha=100, n_knots=8, mf=0.7),
 'bb_f2b':  lambda: BackboneTree(feats2.with_angle, alpha=30, n_knots=8, mf=1.0),
 'bb_phys': lambda: BackboneTree(feats.phys_features, alpha=100, n_knots=8, mf=0.7),
}
oof={k:np.zeros(len(both)) for k in bases}
for tri,tei in gk.split(both, both[E.TARGET], gidx):
    for k,mkf in bases.items():
        m=mkf(); m.fit(both.iloc[tri][E.FEATURES], both.iloc[tri][E.TARGET].values)
        oof[k][tei]=m.predict(both.iloc[tei][E.FEATURES])
out=pd.DataFrame(oof); out.to_csv('oof3.csv',index=False); print('done', out.shape)
EOF
nohup timeout 290 python3 gen_oof3.py > gen_oof3.log 2>&1 &
echo started; sleep 150; cat gen_oof3.log; ls -la oof3.csv 2>/dev/null

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
o2=pd.read_csv('oof2.csv'); o3=pd.read_csv('oof3.csv')
y=o2.y.values; cond=o2['cond'].values
M=pd.concat([o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']], o3], axis=1)
cols=list(M.columns)
def metrics(p): return E.metrics_from(pd.DataFrame({'cond':cond,'err':p-y}))
def rep(nm,p):
    d=metrics(p); print(f'{nm:30s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
# random search weights to minimize p90 s.t. constraints
rng=np.random.RandomState(0); X=M.values
best=None
for it in range(60000):
    w=rng.dirichlet(np.ones(len(cols)))
    p=X@w; o=pd.DataFrame({'cond':cond,'err':p-y}); g=o.groupby('cond').err
    cr=g.apply(lambda e:np.sqrt(np.mean(e**2)))
    rmse=np.sqrt(np.mean((p-y)**2)); mae=np.mean(np.abs(p-y))
    p90=np.percentile(cr,90); p95=np.percentile(cr,95); mx=cr.max()
    if rmse<=2.15 and mae<=1.55 and p95<=4.0 and mx<=5.0:
        score=p90
        if best is None or score<best[0]: best=(score,w.copy(),rmse,mae,p90,p95,mx)
print('best p90 blend:', {c:round(wv,3) for c,wv in zip(cols,best[1]) if wv>0.02})
print(f'  rmse={best[2]:.3f} mae={best[3]:.3f} p90={best[4]:.3f} p95={best[5]:.3f} max={best[6]:.3f}')
rep('equal all9', M.mean(1).values)
rep('ts6 + bb3 (0.6/0.4)', (o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1)*0.6 + o3.mean(1)*0.4).values)
"

restart: true

cd /app/scratch && timeout 200 python3 -c "
import numpy as np, pandas as pd, evallib as E
o2=pd.read_csv('oof2.csv'); o3=pd.read_csv('oof3.csv')
y=o2.y.values; cond=o2['cond'].values
M=pd.concat([o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']], o3], axis=1)
cols=list(M.columns); X=M.values
# precompute condition index groups
uc, inv = np.unique(cond, return_inverse=True)
nC=len(uc); counts=np.bincount(inv)
def fast_metrics(p):
    e=p-y; rmse=np.sqrt(np.mean(e**2)); mae=np.mean(np.abs(e))
    sse=np.bincount(inv, weights=e*e); cr=np.sqrt(sse/counts)
    return rmse,mae,np.percentile(cr,90),np.percentile(cr,95),cr.max()
def rep(nm,p):
    r=fast_metrics(p); print(f'{nm:24s} rmse={r[0]:.3f} mae={r[1]:.3f} p90={r[2]:.3f} p95={r[3]:.3f} max={r[4]:.3f}')
rng=np.random.RandomState(0); best=None
for it in range(8000):
    w=rng.dirichlet(np.ones(len(cols))); r=fast_metrics(X@w)
    if r[0]<=2.15 and r[1]<=1.55 and r[3]<=4.2 and r[4]<=5.1:
        if best is None or r[2]<best[0]: best=(r[2],w.copy(),r)
print('best p90 blend weights:', {c:round(wv,2) for c,wv in zip(cols,best[1]) if wv>0.03})
print(f'  -> rmse={best[2][0]:.3f} mae={best[2][1]:.3f} p90={best[2][2]:.3f} p95={best[2][3]:.3f} max={best[2][4]:.3f}')
rep('equal all9', M.mean(1).values)
rep('ts6*.6+bb3*.4', o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1).values*0.6 + o3.mean(1).values*0.4)
rep('ts6*.7+bb3*.3', o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1).values*0.7 + o3.mean(1).values*0.3)
"

cd /app/scratch && timeout 120 python3 -c "
import numpy as np, pandas as pd
o2=pd.read_csv('oof2.csv'); o3=pd.read_csv('oof3.csv')
y=o2.y.values; cond=o2['cond'].values
ts=o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1).values
bb=o3.mean(1).values
pred=0.6*ts+0.4*bb
uc,inv=np.unique(cond,return_inverse=True); counts=np.bincount(inv)
e=pred-y; sse=np.bincount(inv,weights=e*e); cr=np.sqrt(sse/counts)
rows_by=[np.where(inv==i)[0] for i in range(len(uc))]
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,max=5.60)
rng=np.random.RandomState(1)
for size in [15,20,25]:
    cnt={k:0 for k in T}; pa=0; N=5000
    for _ in range(N):
        S=rng.choice(len(uc),size=size,replace=False)
        idx=np.concatenate([rows_by[i] for i in S]); ee=pred[idx]-y[idx]; crr=cr[S]
        m=dict(rmse=np.sqrt(np.mean(ee**2)),mae=np.mean(np.abs(ee)),p90=np.percentile(crr,90),p95=np.percentile(crr,95),max=crr.max())
        ok=True
        for k in T:
            if m[k]<=T[k]: cnt[k]+=1
            else: ok=False
        pa+=ok
    print(f'sz={size}: '+' '.join(f'{k}={cnt[k]/N:.2f}' for k in T)+f' ALL={pa/N:.3f}')
"

cd /app/scratch && cat > blend.py << 'EOF'
import numpy as np, pandas as pd
from sklearn.base import BaseEstimator, RegressorMixin
import feats, feats2
from twostage import TwoStageGen, BackboneTree
def build_models():
    return {
     'ts_et_f2':   (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=0.7), 0.6/6),
     'ts_et_f2b':  (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=1.0), 0.6/6),
     'ts_et_phys': (TwoStageGen(feats.phys_features, alpha=30, kind='et', mf=0.7), 0.6/6),
     'ts_et_noang':(TwoStageGen(feats2.no_angle, alpha=30, kind='et', mf=0.7), 0.6/6),
     'ts_hgb_f2':  (TwoStageGen(feats2.with_angle, alpha=30, kind='hgb'), 0.6/6),
     'ts_hgb_phys':(TwoStageGen(feats.phys_features, alpha=30, kind='hgb'), 0.6/6),
     'bb_f2':      (BackboneTree(feats2.with_angle, alpha=100, n_knots=8, mf=0.7), 0.4/3),
     'bb_f2b':     (BackboneTree(feats2.with_angle, alpha=30, n_knots=8, mf=1.0), 0.4/3),
     'bb_phys':    (BackboneTree(feats.phys_features, alpha=100, n_knots=8, mf=0.7), 0.4/3),
    }
class Blend(BaseEstimator, RegressorMixin):
    def fit(self, X, y):
        self.models_=build_models()
        for k,(m,w) in self.models_.items(): m.fit(X,y)
        return self
    def predict(self, X):
        tot=0.0
        for k,(m,w) in self.models_.items(): tot=tot+w*m.predict(X)
        return tot
EOF
timeout 200 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings
warnings.filterwarnings('ignore')
from blend import Blend
tr,va=E.load()
m=Blend().fit(tr[E.FEATURES], tr[E.TARGET].values)
pred=m.predict(va[E.FEATURES])
d=E.metrics_from(pd.DataFrame({'cond':E.condition_ids(va).values,'err':pred-va[E.TARGET].values}))
print('TRAIN->VAL holdout (15 disjoint conditions):')
print(f\"  rmse={d['rmse']:.3f} mae={d['mae']:.3f} p90={d['p90_cond']:.3f} p95={d['p95_cond']:.3f} max={d['max_cond']:.3f}\")
print('Targets: rmse<=2.25 mae<=1.55 p90<=2.35 p95<=4.75 max<=5.60')
"

cd /app/scratch && timeout 200 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings
warnings.filterwarnings('ignore')
from scipy.interpolate import UnivariateSpline
from blend import Blend
tr,va=E.load()
m=Blend().fit(tr[E.FEATURES], tr[E.TARGET].values)
va=va.copy(); va['pred']=m.predict(va[E.FEATURES]); va['cond']=E.condition_ids(va)
va['logSt']=np.log10(va.frequency*va.suction_side_displacement_thickness/va.free_stream_velocity)
def metrics(col):
    d=E.metrics_from(pd.DataFrame({'cond':va['cond'].values,'err':va[col].values-va[E.TARGET].values}))
    return d
def rep(nm,col):
    d=metrics(col); print(f'{nm:20s} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
rep('blend raw','pred')
for sfac in [0.3,0.6,1.0]:
    sm=va.pred.values.copy()
    for c,g in va.groupby('cond'):
        if len(g)<6: continue
        gg=g.sort_values('logSt'); x=gg.logSt.values; p=gg.pred.values
        try:
            sp=UnivariateSpline(x,p,k=3,s=len(x)*sfac); sm[gg.index]=sp(x)
        except Exception: pass
    va['sm']=sm; rep(f'smooth s={sfac}','sm')
"

cd /app/scratch && timeout 250 python3 -c "
import time, numpy as np, pandas as pd, evallib as E, warnings
warnings.filterwarnings('ignore')
from blend import Blend
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True)
t0=time.time()
m=Blend().fit(both[E.FEATURES], both[E.TARGET].values)
t1=time.time()
p=m.predict(va[E.FEATURES])
t2=time.time()
print(f'fit {t1-t0:.1f}s predict {t2-t1:.1f}s  pred range [{p.min():.1f},{p.max():.1f}] finite={np.isfinite(p).all()}')
"

cd /app/scratch && timeout 150 python3 -c "
import numpy as np, pandas as pd
o2=pd.read_csv('oof2.csv'); o3=pd.read_csv('oof3.csv')
y=o2.y.values; cond=o2['cond'].values
M=pd.concat([o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']], o3], axis=1)
cols=list(M.columns); X=M.values
uc,inv=np.unique(cond,return_inverse=True); counts=np.bincount(inv)
def fm(p):
    e=p-y; sse=np.bincount(inv,weights=e*e); cr=np.sqrt(sse/counts)
    return np.sqrt(np.mean(e**2)),np.mean(np.abs(e)),np.percentile(cr,90),np.percentile(cr,95),cr.max()
rng=np.random.RandomState(3); best=None
for it in range(20000):
    w=rng.dirichlet(np.ones(len(cols))*0.5)
    r=fm(X@w)
    if r[0]<=2.20 and r[3]<=4.3 and r[4]<=5.3:
        obj=max(r[2]/2.35, r[1]/1.55)  # binding metrics normalized
        if best is None or obj<best[0]: best=(obj,w.copy(),r)
w=best[1]
print('weights:', {c:round(v,3) for c,v in zip(cols,w) if v>0.02})
r=best[1]; rr=best[2]
print(f'obj={best[0]:.3f} rmse={rr[0]:.3f} mae={rr[1]:.3f} p90={rr[2]:.3f} p95={rr[3]:.3f} max={rr[4]:.3f}')
np.save('best_w.npy', w)
import json; json.dump(cols, open('cols.json','w'))
"

command: str_replace
path: /app/scratch/blend.py
old_str: def build_models():
    return {
     'ts_et_f2':   (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=0.7), 0.6/6),
     'ts_et_f2b':  (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=1.0), 0.6/6),
     'ts_et_phys': (TwoStageGen(feats.phys_features, alpha=30, kind='et', mf=0.7), 0.6/6),
     'ts_et_noang':(TwoStageGen(feats2.no_angle, alpha=30, kind='et', mf=0.7), 0.6/6),
     'ts_hgb_f2':  (TwoStageGen(feats2.with_angle, alpha=30, kind='hgb'), 0.6/6),
     'ts_hgb_phys':(TwoStageGen(feats.phys_features, alpha=30, kind='hgb'), 0.6/6),
     'bb_f2':      (BackboneTree(feats2.with_angle, alpha=100, n_knots=8, mf=0.7), 0.4/3),
     'bb_f2b':     (BackboneTree(feats2.with_angle, alpha=30, n_knots=8, mf=1.0), 0.4/3),
     'bb_phys':    (BackboneTree(feats.phys_features, alpha=100, n_knots=8, mf=0.7), 0.4/3),
    }
new_str: def build_models():
    return {
     'ts_et_f2':   (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=0.7), 0.22),
     'ts_et_f2b':  (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=1.0), 0.30),
     'ts_et_phys': (TwoStageGen(feats.phys_features, alpha=30, kind='et', mf=0.7), 0.14),
     'ts_et_noang':(TwoStageGen(feats2.no_angle, alpha=30, kind='et', mf=0.7), 0.06),
     'ts_hgb_f2':  (TwoStageGen(feats2.with_angle, alpha=30, kind='hgb'), 0.10),
     'ts_hgb_phys':(TwoStageGen(feats.phys_features, alpha=30, kind='hgb'), 0.05),
     'bb_f2':      (BackboneTree(feats2.with_angle, alpha=100, n_knots=8, mf=0.7), 0.08),
     'bb_f2b':     (BackboneTree(feats2.with_angle, alpha=30, n_knots=8, mf=1.0), 0.05),
    }

cd /app/scratch && timeout 120 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings
warnings.filterwarnings('ignore')
from blend import Blend
tr,va=E.load()
w=sum(wt for _,(_,wt) in Blend().__class__.__dict__.items() if False)  # noop
m=Blend().fit(tr[E.FEATURES], tr[E.TARGET].values)
pred=m.predict(va[E.FEATURES])
d=E.metrics_from(pd.DataFrame({'cond':E.condition_ids(va).values,'err':pred-va[E.TARGET].values}))
print('TRAIN->VAL (tuned weights):')
print(f\"  rmse={d['rmse']:.3f} mae={d['mae']:.3f} p90={d['p90_cond']:.3f} p95={d['p95_cond']:.3f} max={d['max_cond']:.3f}\")
import blend
print('weight sum =', sum(wt for (_,wt) in blend.build_models().values()))
"

cd /app/scratch && timeout 120 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings
warnings.filterwarnings('ignore')
import blend
from blend import Blend
tr,va=E.load()
m=Blend().fit(tr[E.FEATURES], tr[E.TARGET].values)
pred=m.predict(va[E.FEATURES])
d=E.metrics_from(pd.DataFrame({'cond':E.condition_ids(va).values,'err':pred-va[E.TARGET].values}))
print('TRAIN->VAL (tuned weights):')
print(f\"  rmse={d['rmse']:.3f} mae={d['mae']:.3f} p90={d['p90_cond']:.3f} p95={d['p95_cond']:.3f} max={d['max_cond']:.3f}\")
print('weight sum =', round(sum(wt for (_,wt) in blend.build_models().values()),4))
"

command: str_replace
path: /app/scratch/blend.py
old_str: def build_models():
    return {
     'ts_et_f2':   (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=0.7), 0.22),
     'ts_et_f2b':  (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=1.0), 0.30),
     'ts_et_phys': (TwoStageGen(feats.phys_features, alpha=30, kind='et', mf=0.7), 0.14),
     'ts_et_noang':(TwoStageGen(feats2.no_angle, alpha=30, kind='et', mf=0.7), 0.06),
     'ts_hgb_f2':  (TwoStageGen(feats2.with_angle, alpha=30, kind='hgb'), 0.10),
     'ts_hgb_phys':(TwoStageGen(feats.phys_features, alpha=30, kind='hgb'), 0.05),
     'bb_f2':      (BackboneTree(feats2.with_angle, alpha=100, n_knots=8, mf=0.7), 0.08),
     'bb_f2b':     (BackboneTree(feats2.with_angle, alpha=30, n_knots=8, mf=1.0), 0.05),
    }
new_str: def build_models():
    ts = 0.6 / 6.0
    bb = 0.4 / 3.0
    return {
     'ts_et_f2':   (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=0.7), ts),
     'ts_et_f2b':  (TwoStageGen(feats2.with_angle, alpha=30, kind='et', mf=1.0), ts),
     'ts_et_phys': (TwoStageGen(feats.phys_features, alpha=30, kind='et', mf=0.7), ts),
     'ts_et_noang':(TwoStageGen(feats2.no_angle, alpha=30, kind='et', mf=0.7), ts),
     'ts_hgb_f2':  (TwoStageGen(feats2.with_angle, alpha=30, kind='hgb'), ts),
     'ts_hgb_phys':(TwoStageGen(feats.phys_features, alpha=30, kind='hgb'), ts),
     'bb_f2':      (BackboneTree(feats2.with_angle, alpha=100, n_knots=8, mf=0.7), bb),
     'bb_f2b':     (BackboneTree(feats2.with_angle, alpha=30, n_knots=8, mf=1.0), bb),
     'bb_phys':    (BackboneTree(feats.phys_features, alpha=100, n_knots=8, mf=0.7), bb),
    }

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

Trains a physics-informed blended regressor on the labeled public data
(train.csv + validation.csv) and writes predictions for the grader-staged
feature-only rows at /app/data/test_features.csv.

Design notes
------------
The split holds out whole aerodynamic condition groups
``(attack_angle, chord_length, free_stream_velocity)``, so the model must
generalize to unseen physical regimes.  Key ideas that make this work:

* Physics features derived from the Brooks-Pope-Marcolini (BPM) trailing-edge
  noise model: the Strouhal number ``St = f * delta* / U`` collapses the
  spectral shape, while log scalings of velocity / chord / thickness capture
  the overall level.
* A two-stage "safe" structure: a linear model on *frequency-independent*
  log-level features predicts the overall condition level (this extrapolates
  gracefully to grid-boundary conditions such as the highest velocity),
  and tree ensembles fit the remaining smooth spectral residual.
* A smooth spline "backbone" variant that extrapolates the Strouhal spectral
  curve with linear tails, taming the single worst held-out condition.
* A diversified blend of these models to reduce per-condition tail errors
  (the p90 / p95 / max condition-RMSE metrics).
"""

from pathlib import Path
import warnings

import numpy as np
import pandas as pd

from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler, SplineTransformer
from sklearn.linear_model import Ridge
from sklearn.ensemble import ExtraTreesRegressor, HistGradientBoostingRegressor

warnings.filterwarnings("ignore")

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",
]

C_SOUND = 343.0          # speed of sound (m/s)
NU = 1.5e-5              # kinematic viscosity of air (m^2/s)
# Safety clip: training targets span ~105-140 dB; guard against any wild
# extrapolation while never touching physically plausible predictions.
CLIP_LO, CLIP_HI = 95.0, 148.0


# --------------------------------------------------------------------------- #
# Feature engineering                                                         #
# --------------------------------------------------------------------------- #
def _cols(X):
    f = np.maximum(X["frequency"].astype(float).values, 1.0)
    a = X["attack_angle"].astype(float).values
    c = np.maximum(X["chord_length"].astype(float).values, 1e-6)
    U = np.maximum(X["free_stream_velocity"].astype(float).values, 1e-3)
    d = np.maximum(X["suction_side_displacement_thickness"].astype(float).values, 1e-9)
    return f, a, c, U, d


def phys_features(X):
    """Rich physics feature set including raw attack angle."""
    f, a, c, U, d = _cols(X)
    logf, logd, logc, logU = np.log10(f), np.log10(d), np.log10(c), np.log10(U)
    M = U / C_SOUND
    logSt = np.log10(f * d / U)
    logRe = np.log10(U * c / NU)
    logRed = np.log10(U * d / NU)
    o = pd.DataFrame(index=X.index)
    o["logf"] = logf; o["logc"] = logc; o["logU"] = logU; o["logd"] = logd; o["M"] = M
    o["logSt"] = logSt; o["logSt2"] = logSt ** 2; o["logSt3"] = logSt ** 3
    o["logRe"] = logRe; o["logRed"] = logRed; o["logdc"] = np.log10(d / c)
    o["a"] = a; o["a2"] = a ** 2; o["a3"] = a ** 3
    o["logSt_logM"] = logSt * np.log10(M); o["logSt_logc"] = logSt * logc
    o["logSt_logU"] = logSt * logU; o["logSt_logd"] = logSt * logd
    o["logSt_logRe"] = logSt * logRe; o["logSt_a"] = logSt * a
    o["a_logU"] = a * logU; o["a_logc"] = a * logc; o["a_logd"] = a * logd
    o["logc_logU"] = logc * logU; o["logd_logU"] = logd * logU; o["logc_logd"] = logc * logd
    return o


d
…[truncated]

cd /app && rm -f output/predictions.csv && time python3 train_and_predict.py && echo "---HEAD---" && head -5 output/predictions.csv && echo "---VALIDATE---" && python3 -c "
import pandas as pd, numpy as np
pred=pd.read_csv('output/predictions.csv'); test=pd.read_csv('data/test_features.csv')
assert list(pred.columns)==['sample_id','predicted_scaled_sound_pressure_level'], pred.columns.tolist()
assert len(pred)==len(test), (len(pred),len(test))
assert set(pred.sample_id.astype(str))==set(test.sample_id.astype(str)), 'id mismatch'
assert np.isfinite(pred.predicted_scaled_sound_pressure_level).all()
assert pred.predicted_scaled_sound_pressure_level.nunique()>5, 'not constant'
print('rows', len(pred), 'unique preds', pred.predicted_scaled_sound_pressure_level.nunique())
print('OK: columns, alignment, finiteness, non-constant all pass')
"

cd /app && timeout 150 python3 -c "
import numpy as np, pandas as pd, warnings, importlib.util
warnings.filterwarnings('ignore')
spec=importlib.util.spec_from_file_location('tap','/app/train_and_predict.py')
tap=importlib.util.module_from_spec(spec); spec.loader.exec_module(tap)
tr=pd.read_csv('data/train.csv'); va=pd.read_csv('data/validation.csv')
# nominal condition ids for val (for per-condition metrics)
CH=[0.0254,0.0508,0.1016,0.1524,0.2286,0.3048]; VE=[31.7,39.6,55.5,71.3]
AN=np.array([0.0,1.5,3.0,4.0,5.3,7.3,8.9,9.9,12.3,12.6,15.4,15.6,17.4,19.7,22.2])
def nr(v,a): a=np.array(a); return a[np.argmin(np.abs(a[None,:]-np.array(v)[:,None]),axis=1)]
cond=[f'{x:.2f}_{y:.4f}_{z:.1f}' for x,y,z in zip(nr(va.attack_angle,AN),nr(va.chord_length,CH),nr(va.free_stream_velocity,VE))]
X=tr[tap.FEATURES]; y=tr[tap.TARGET].astype(float).values
preds=np.zeros(len(va)); ws=0.0
for m,w in tap.build_blend():
    m.fit(X,y); preds+=w*m.predict(va[tap.FEATURES]); ws+=w
preds=np.clip(preds/ws, tap.CLIP_LO, tap.CLIP_HI)
err=preds-va[tap.TARGET].values
d=pd.DataFrame({'cond':cond,'err':err})
cr=d.groupby('cond').err.apply(lambda e:np.sqrt(np.mean(e**2)))
print('CONSOLIDATED script, TRAIN->VAL (honest, disjoint conditions):')
print(f'  rmse={np.sqrt(np.mean(err**2)):.3f}  mae={np.mean(np.abs(err)):.3f}  p90={np.percentile(cr,90):.3f}  p95={np.percentile(cr,95):.3f}  max={cr.max():.3f}')
print('  targets: rmse<=2.25 mae<=1.55 p90<=2.35 p95<=4.75 max<=5.60')
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd, evallib as E
o2=pd.read_csv('oof2.csv'); o3=pd.read_csv('oof3.csv')
y=o2.y.values; cond=o2['cond'].values
ts=o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1).values
bb=o3.mean(1).values
pred=0.6*ts+0.4*bb
df=pd.DataFrame({'cond':cond,'err':pred-y})
df['ang']=pd.Series(cond).str.split('_').str[0].astype(float).values
# count conditions per nominal angle
ang_counts=df.drop_duplicates('cond').groupby('ang').size()
singleton_angles=set(ang_counts[ang_counts<=1].index)
print('singleton angles (in only 1 condition):', sorted(singleton_angles))
def metrics(sub):
    g=sub.groupby('cond').err; cr=g.apply(lambda e:np.sqrt(np.mean(e**2)))
    return (np.sqrt(np.mean(sub.err**2)),np.mean(np.abs(sub.err)),np.percentile(cr,90),np.percentile(cr,95),cr.max(),len(cr))
for nm,sub in [('ALL 82 conditions',df),('excl singleton-angle conds',df[~df.ang.isin(singleton_angles)])]:
    r=metrics(sub)
    print(f'{nm:32s} rmse={r[0]:.3f} mae={r[1]:.3f} p90={r[2]:.3f} p95={r[3]:.3f} max={r[4]:.3f} (nc={r[5]})')
"

cd /app/scratch && timeout 120 python3 -c "
import numpy as np, pandas as pd
o2=pd.read_csv('oof2.csv'); o3=pd.read_csv('oof3.csv')
y=o2.y.values; cond=o2['cond'].values
pred=0.6*o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1).values+0.4*o3.mean(1).values
df=pd.DataFrame({'cond':cond,'err':pred-y}); df['ang']=pd.Series(cond).str.split('_').str[0].astype(float).values
ac=df.drop_duplicates('cond').groupby('ang').size(); singles=set(ac[ac<=1].index)
sub=df[~df.ang.isin(singles)].reset_index(drop=True)
uc,inv=np.unique(sub['cond'].values,return_inverse=True); counts=np.bincount(inv)
e=sub.err.values; cr=np.sqrt(np.bincount(inv,weights=e*e)/counts)
rows_by=[np.where(inv==i)[0] for i in range(len(uc))]
T=dict(rmse=2.25,mae=1.55,p90=2.35,p95=4.75,max=5.60); rng=np.random.RandomState(7)
for size in [15,20,25]:
    cnt={k:0 for k in T}; pa=0; N=5000
    for _ in range(N):
        S=rng.choice(len(uc),size=min(size,len(uc)),replace=False)
        idx=np.concatenate([rows_by[i] for i in S]); ee=sub.err.values[idx]; crr=cr[S]
        m=dict(rmse=np.sqrt(np.mean(ee**2)),mae=np.mean(np.abs(ee)),p90=np.percentile(crr,90),p95=np.percentile(crr,95),max=crr.max())
        ok=True
        for k in T:
            if m[k]<=T[k]: cnt[k]+=1
            else: ok=False
        pa+=ok
    print(f'sz={size}: '+' '.join(f'{k}={cnt[k]/N:.2f}' for k in T)+f' ALL={pa/N:.3f}')
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd
o2=pd.read_csv('oof2.csv'); o3=pd.read_csv('oof3.csv')
y=o2.y.values; cond=o2['cond'].values
cols2=['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']
A=np.column_stack([o2[cols2].values, o3.values])  # 9 models
w=np.array([0.1]*6+[0.4/3]*3)
uc,inv=np.unique(cond,return_inverse=True); counts=np.bincount(inv)
def m(pred):
    e=pred-y; cr=np.sqrt(np.bincount(inv,weights=e*e)/counts)
    return np.sqrt(np.mean(e**2)),np.mean(np.abs(e)),np.percentile(cr,90),np.percentile(cr,95),cr.max()
def rep(nm,pred):
    r=m(pred); print(f'{nm:22s} rmse={r[0]:.3f} mae={r[1]:.3f} p90={r[2]:.3f} p95={r[3]:.3f} max={r[4]:.3f}')
rep('weighted mean', A@w)
rep('plain mean', A.mean(1))
rep('median', np.median(A,1))
# trimmed mean (drop min & max per row)
As=np.sort(A,1); rep('trimmed mean(1)', As[:,1:-1].mean(1))
# mean of the 4 ET-f2/phys two-stage (best individuals) + median
rep('ET-ts mean', o2[['ts_et_f2','ts_et_f2b','ts_et_phys']].mean(1).values)
"

cd /app && echo "=== /app contents ===" && ls -la && echo "=== *.py in /app (top level) ===" && ls *.py && echo "=== confirm script imports only std/sklearn ===" && grep -n "^import\|^from" train_and_predict.py | grep -iv "sklearn\|numpy\|pandas\|pathlib\|warnings" || echo "(no non-standard imports)"

cd /app/scratch && timeout 260 python3 -c "
import numpy as np, pandas as pd, evallib as E, warnings, fastcv, feats, feats2
warnings.filterwarnings('ignore')
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer, StandardScaler
from sklearn.linear_model import Ridge
from sklearn.neural_network import MLPRegressor
from twostage import level_feats
class TSMLP(BaseEstimator, RegressorMixin):
    def __init__(self, feat_fn=feats2.with_angle, alpha=30.0, hidden=(64,64), nn_alpha=1e-2, n=4):
        self.feat_fn=feat_fn; self.alpha=alpha; self.hidden=hidden; self.nn_alpha=nn_alpha; self.n=n
    def fit(self,X,y):
        y=np.asarray(y,float)
        self.lin_=make_pipeline(FunctionTransformer(level_feats),StandardScaler(),Ridge(alpha=self.alpha)); self.lin_.fit(X,y)
        r=y-self.lin_.predict(X)
        self.ms_=[]
        for s in range(self.n):
            m=make_pipeline(FunctionTransformer(self.feat_fn),StandardScaler(),
                MLPRegressor(hidden_layer_sizes=self.hidden,alpha=self.nn_alpha,activation='tanh',solver='lbfgs',max_iter=1200,random_state=s))
            m.fit(X,r); self.ms_.append(m)
        return self
    def predict(self,X):
        return self.lin_.predict(X)+np.mean([m.predict(X) for m in self.ms_],axis=0)
tr,va=E.load(); both=pd.concat([tr,va],ignore_index=True).reset_index(drop=True)
both['cond']=E.condition_ids(both); groups=both['cond'].values
for hid in [(48,48),(64,64)]:
    d=fastcv.fast_eval(lambda hid=hid: TSMLP(feats2.with_angle,hidden=hid,nn_alpha=1e-2,n=4), both, groups, n_splits=12)
    print(f'TS-MLP {hid} rmse={d[\"rmse\"]:.3f} mae={d[\"mae\"]:.3f} p90={d[\"p90_cond\"]:.3f} p95={d[\"p95_cond\"]:.3f} max={d[\"max_cond\"]:.3f}')
"

cd /app && rm -rf __pycache__ && rm -f output/predictions.csv && echo "=== FINAL RUN (as grader invokes) ===" && python3 train_and_predict.py && echo && echo "=== OUTPUT CHECK ===" && python3 -c "
import pandas as pd, numpy as np
p=pd.read_csv('output/predictions.csv'); t=pd.read_csv('data/test_features.csv')
assert list(p.columns)==['sample_id','predicted_scaled_sound_pressure_level']
assert len(p)==len(t) and set(p.sample_id.astype(str))==set(t.sample_id.astype(str))
v=p['predicted_scaled_sound_pressure_level']
assert np.isfinite(v).all() and v.nunique()>5
assert v.round(6).duplicated().mean()<0.5  # not implausibly exact/constant
print('rows:',len(p),'| range: [%.2f, %.2f]'%(v.min(),v.max()),'| unique:',v.nunique())
print('ALL OUTPUT CHECKS PASSED')
"

cd /app/scratch && python3 -c "
import numpy as np, pandas as pd
o2=pd.read_csv('oof2.csv'); o3=pd.read_csv('oof3.csv')
y=o2.y.values; cond=o2['cond'].values
ts=o2[['ts_et_f2','ts_et_f2b','ts_et_phys','ts_et_noang','ts_hgb_f2','ts_hgb_phys']].mean(1).values
bb=o3.mean(1).values
uc,inv=np.unique(cond,return_inverse=True); counts=np.bincount(inv)
def m(pred):
    e=pred-y; cr=np.sqrt(np.bincount(inv,weights=e*e)/counts)
    return np.sqrt(np.mean(e**2)),np.mean(np.abs(e)),np.percentile(cr,90),np.percentile(cr,95),cr.max()
for wts in [0.5,0.6,0.7,0.8,1.0]:
    r=m(wts*ts+(1-wts)*bb)
    print(f'ts={wts:.1f}/bb={1-wts:.1f}: rmse={r[0]:.3f} mae={r[1]:.3f} p90={r[2]:.3f} p95={r[3]:.3f} max={r[4]:.3f}')
"

command: str_replace
path: /app/train_and_predict.py
old_str:     ts = 0.6 / 6.0
    bb = 0.4 / 3.0
new_str:     ts = 0.5 / 6.0
    bb = 0.5 / 3.0

cd /app && rm -rf __pycache__ && rm -f output/predictions.csv && python3 train_and_predict.py && echo "---" && head -3 output/predictions.csv && python3 -c "
import pandas as pd, numpy as np
p=pd.read_csv('output/predictions.csv'); t=pd.read_csv('data/test_features.csv')
assert list(p.columns)==['sample_id','predicted_scaled_sound_pressure_level']
assert len(p)==len(t)==p.sample_id.nunique()
assert np.isfinite(p.iloc[:,1]).all() and p.iloc[:,1].nunique()>5
print('FINAL OK: rows',len(p))
"

cd /app && rm -rf __pycache__ && timeout 150 python3 -c "
import numpy as np, pandas as pd, warnings, importlib.util
warnings.filterwarnings('ignore')
spec=importlib.util.spec_from_file_location('tap','/app/train_and_predict.py'); tap=importlib.util.module_from_spec(spec); spec.loader.exec_module(tap)
tr=pd.read_csv('data/train.csv'); va=pd.read_csv('data/validation.csv')
CH=[0.0254,0.0508,0.1016,0.1524,0.2286,0.3048]; VE=[31.7,39.6,55.5,71.3]
AN=np.array([0.0,1.5,3.0,4.0,5.3,7.3,8.9,9.9,12.3,12.6,15.4,15.6,17.4,19.7,22.2])
def nr(v,a): a=np.array(a); return a[np.argmin(np.abs(a[None,:]-np.array(v)[:,None]),axis=1)]
cond=[f'{x:.2f}_{y:.4f}_{z:.1f}' for x,y,z in zip(nr(va.attack_angle,AN),nr(va.chord_length,CH),nr(va.free_stream_velocity,VE))]
X=tr[tap.FEATURES]; y=tr[tap.TARGET].astype(float).values
preds=np.zeros(len(va)); ws=0.0
for m,w in tap.build_blend(): m.fit(X,y); preds+=w*m.predict(va[tap.FEATURES]); ws+=w
err=np.clip(preds/ws,tap.CLIP_LO,tap.CLIP_HI)-va[tap.TARGET].values
cr=pd.DataFrame({'c':cond,'e':err}).groupby('c').e.apply(lambda e:np.sqrt(np.mean(e**2)))
print('FINAL train->val proxy (conservative, 67 train conds):')
print(f'  rmse={np.sqrt(np.mean(err**2)):.3f}(<=2.25)  mae={np.mean(np.abs(err)):.3f}(<=1.55)  p90={np.percentile(cr,90):.3f}(<=2.35)  p95={np.percentile(cr,95):.3f}(<=4.75)  max={cr.max():.3f}(<=5.60)')
print('Note: OOF estimate (76-cond train, closer to production 82): rmse~2.01 mae~1.49 p90~2.63 p95~3.63 max~4.42')
"

cd /app && rm -rf __pycache__ && echo "=== deliverables ===" && ls -la train_and_predict.py output/predictions.csv && echo && sed -n '/^def build_blend/,/^    ]/p' train_and_predict.py