{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":202299006,"sourceType":"kernelVersion"},{"sourceId":202375230,"sourceType":"kernelVersion"}],"dockerImageVersionId":30786,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# great ideas borrowed from 'ariel_only_correlation', https://www.kaggle.com/code/sergeifironov/ariel-only-correlation","metadata":{"execution":{"iopub.status.busy":"2024-10-17T02:46:11.005038Z","iopub.execute_input":"2024-10-17T02:46:11.005512Z","iopub.status.idle":"2024-10-17T02:46:11.061488Z","shell.execute_reply.started":"2024-10-17T02:46:11.005467Z","shell.execute_reply":"2024-10-17T02:46:11.060187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport numpy as np, pandas as pd, matplotlib.pyplot as plt\nimport random\n\nfrom ariel_calib import train_planets, test_planets, read_signal, CalibConfig \n\ninpdir = '../input/ariel-data-challenge-2024'\ndatadir = './ariel-2024-calibrated'\nif not os.path.isdir(datadir):\n    !unzip -q ../input/ariel-2024-calibrated/ariel-2024-calibrated.zip -d {datadir}\n\n!ls -Llah {datadir}\n\ntime_binning = 30","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:13:33.874704Z","iopub.execute_input":"2024-10-24T05:13:33.875809Z","iopub.status.idle":"2024-10-24T05:13:46.963216Z","shell.execute_reply.started":"2024-10-24T05:13:33.875756Z","shell.execute_reply":"2024-10-24T05:13:46.961842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_labels = pd.read_csv(f'{datadir}/train_labels.csv', index_col='planet_id')","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:13:53.242996Z","iopub.execute_input":"2024-10-24T05:13:53.243508Z","iopub.status.idle":"2024-10-24T05:13:53.319468Z","shell.execute_reply.started":"2024-10-24T05:13:53.243457Z","shell.execute_reply":"2024-10-24T05:13:53.31828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import itertools\nfrom scipy.optimize import minimize\nfrom functools import partial\nimport random, os\n\n# phase_detector, try_s, calibrate_signal, calibrate_train, feature_extract\n# from https://www.kaggle.com/code/sergeifironov/ariel-only-correlation\ndef phase_detector(signal):\n    phase1, phase2 = None, None\n    best_drop = 0\n    for i in range(50//2,150//2):        \n        t1 = signal[i:i+20//2].max() - signal[i:i+20//2].min()\n        if t1 > best_drop:\n            phase1 = i+(20+5)//2\n            best_drop = t1\n    \n    best_drop = 0\n    for i in range(200//2,250//2):\n        t1 = signal[i:i+20//2].max() - signal[i:i+20//2].min()\n        if t1 > best_drop:\n            phase2 = i-5//2\n            best_drop = t1\n    \n    return phase1, phase2\n\ndef try_s(signal, p1, p2, deg, s):\n    out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n    x, y = out, signal[out].tolist()\n    x = x + list(range(p1,p2))\n\n    y = y + (signal[p1:p2] * (1 + s[0])).tolist()\n    z = np.polyfit(x, y, deg)\n    p = np.poly1d(z)\n    q = np.abs(p(x) - y).mean()\n\n    if s < 1e-4:\n        return q + 1e3\n\n    return q\n    \ndef calibrate_signal(signal):\n    p1,p2 = phase_detector(signal)\n\n    best_deg, best_score = 1, 1e12\n    for deg in range(1, 4):\n        f = partial(try_s, signal, p1, p2, deg)\n        r = minimize(f, [0.001], method = 'Nelder-Mead')\n        s = r.x[0]\n\n        out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n        x, y = out, signal[out].tolist()\n        x = x + list(range(p1,p2))\n        y = y + (signal[p1:p2] * (1 + s)).tolist()\n    \n        z = np.polyfit(x, y, deg)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        \n        if q < best_score:\n            best_score = q\n            best_deg = deg\n        \n        print(deg, q)\n            \n    z = np.polyfit(x, y, best_deg)\n    p = np.poly1d(z)\n\n    return s, x, y, p(x)\n\ndef calibrate_train(signal):\n    p1,p2 = phase_detector(signal)\n\n    best_deg, best_score = 1, 1e12\n    for deg in range(1, 4):\n        f = partial(try_s, signal, p1, p2, deg)\n        r = minimize(f, [0.0001], method = 'Nelder-Mead')\n        s = r.x[0]\n\n        out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n        x, y = out, signal[out].tolist()\n        x = x + list(range(p1,p2))\n        y = y + (signal[p1:p2] * (1 + s)).tolist()\n    \n        z = np.polyfit(x, y, deg)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        \n        if q < best_score:\n            best_score = q\n            best_deg = deg\n            \n    z = np.polyfit(x, y, best_deg)\n    p = np.poly1d(z)\n    \n    return s, p(np.arange(signal.shape[0])), p1, p2\n\ndef feature_extract_1(signal, p1, p2):\n    # get best_deg\n    best_deg, best_score, best_s = 1, 1e12, None\n    for deg in range(1, 4):\n        f = partial(try_s, signal, p1, p2, deg)\n        r = minimize(f, [0.0001], method = 'Nelder-Mead')\n        s = r.x[0]\n    \n        out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n        x, y = out, signal[out].tolist()\n        x = x + list(range(p1,p2))\n        y = y + (signal[p1:p2] * (1 + s)).tolist()\n    \n        z = np.polyfit(x, y, deg)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        \n        if q < best_score:\n            best_score = q\n            best_deg = deg\n            best_s = s\n\n    best_s = np.clip(best_s, 1e-4, 0.1)   \n    return best_s\n\ndef feature_extract(airs_signal, lambda_labels=None):\n    mean_signal = airs_signal.mean(1)\n    p1, p2 = phase_detector(mean_signal)\n    \n    if lambda_labels is None:\n        s1 = feature_extract_1(mean_signal, p1, p2)\n        s = [s1]\n    else:\n        if not isinstance(lambda_labels, np.ndarray):\n            lambda_labels = np.array(lambda_labels)\n\n        s = []\n        for k in range(0, lambda_labels.max()+1):\n            y = airs_signal[:, lambda_labels==k].mean(1)\n            s1 = feature_extract_1(y, p1, p2)\n            s.append(s1)\n\n    return s","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:23:46.945821Z","iopub.execute_input":"2024-10-24T05:23:46.946287Z","iopub.status.idle":"2024-10-24T05:23:46.983343Z","shell.execute_reply.started":"2024-10-24T05:23:46.946243Z","shell.execute_reply":"2024-10-24T05:23:46.982132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# partition lambdas\nfrom sklearn.cluster import SpectralClustering\nY = train_labels.values\nY1 = Y / Y.mean(1, keepdims=True) - 1\nclst = SpectralClustering(n_clusters=5).fit(Y1.T[1:])\ndmap = {c:k for k, c in enumerate(np.flip( np.argsort(np.unique(clst.labels_, return_counts=True)[1]) ))}\nlambda_labels = list(map(dmap.get, clst.labels_))","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:23:32.666277Z","iopub.execute_input":"2024-10-24T05:23:32.66745Z","iopub.status.idle":"2024-10-24T05:23:32.983097Z","shell.execute_reply.started":"2024-10-24T05:23:32.667379Z","shell.execute_reply":"2024-10-24T05:23:32.982126Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.linear_model import QuantileRegressor, Ridge\nfrom copy import deepcopy\n\nfit_cloned = []\n\nclass SigmaRegressor:\n    def __init__(self, alpha=0):\n        self.alpha = alpha\n        self.l_regrs = None\n\n    def get_params(self, deep=True):\n        return {'alpha': self.alpha}\n\n    def fit(self, X, Y, verbose=True):\n        self.l_regrs = []\n        s = X[:, :1]\n        X1 = X[:,1:] / s\n        Y = Y / s\n        \n        for k in tqdm.tqdm(range(Y.shape[1]), disable=not verbose):\n            y = Y[:, k]\n            regr0 = QuantileRegressor(quantile=1-0.953, alpha=self.alpha, solver='highs').fit(X1, y) # -2 * sigma\n            regr1 = QuantileRegressor(quantile=0.5, alpha=self.alpha, solver='highs').fit(X1, y)\n            regr2 = QuantileRegressor(quantile=0.953, alpha=self.alpha, solver='highs').fit(X1, y) # +2 * sigma\n            self.l_regrs.append((regr0, regr1, regr2))  \n\n        fit_cloned.append( deepcopy(self) )\n        \n        return self\n\n    def predict(self, X):\n\n        s = X[:, :1]\n        X1 = X[:,1:] / s\n\n        yp0, yp1, yp2 = [], [], []\n        for k, (regr0, regr1, regr2) in enumerate(self.l_regrs):\n            yp0.append( regr0.predict(X1) )\n            yp1.append( regr1.predict(X1) )        \n            yp2.append( regr2.predict(X1) )        \n\n        yp0 = np.stack(yp0, 1) * s \n        yp1 = np.stack(yp1, 1) * s\n        yp2 = np.stack(yp2, 1) * s\n\n        return np.concatenate([yp1, (yp2 - yp0)/4], 1)  # mean, sigma","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:23:33.565606Z","iopub.execute_input":"2024-10-24T05:23:33.566105Z","iopub.status.idle":"2024-10-24T05:23:33.58201Z","shell.execute_reply.started":"2024-10-24T05:23:33.56606Z","shell.execute_reply":"2024-10-24T05:23:33.580887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tqdm\n\nX = []\nfor planet_id in tqdm.tqdm(train_labels.index.values):\n    airs_signal = np.load(f'{datadir}/train_{time_binning}/{planet_id}.npz')['airs']\n    X.append( feature_extract(airs_signal, lambda_labels) )\n\nX = np.stack(X, 0)\nY = train_labels.values","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:23:53.842687Z","iopub.execute_input":"2024-10-24T05:23:53.843205Z","iopub.status.idle":"2024-10-24T05:25:33.982718Z","shell.execute_reply.started":"2024-10-24T05:23:53.843161Z","shell.execute_reply":"2024-10-24T05:25:33.981429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import KFold, cross_val_predict\n\nfit_cloned = []  # fit_cloned after cross_val_predict will contain fitted estimators for evaluating test dataset.\n\nntrial = 5\nYP, YVar = None, None\nfor k in range(ntrial):\n    print(f\"k={k}\")\n    cv_split = KFold(n_splits=5, random_state=k+1234, shuffle=True).split(X,Y)\n    YP1 = cross_val_predict(SigmaRegressor(), X, Y, cv=cv_split)\n    YSigma1 = YP1[:, 283:]\n    YP1 = YP1[:, :283]\n    if YP is None:\n        YP = YP1.copy()\n        YVar = YSigma1**2\n    else:\n        YP += YP1\n        YVar += YSigma1**2\n        \nYP /= ntrial\nYSigma = np.sqrt( YVar / ntrial )","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:26:26.751699Z","iopub.execute_input":"2024-10-24T05:26:26.752184Z","iopub.status.idle":"2024-10-24T05:39:25.259861Z","shell.execute_reply.started":"2024-10-24T05:26:26.752138Z","shell.execute_reply":"2024-10-24T05:39:25.258717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(fit_cloned)","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:40:11.639855Z","iopub.execute_input":"2024-10-24T05:40:11.640355Z","iopub.status.idle":"2024-10-24T05:40:11.648394Z","shell.execute_reply.started":"2024-10-24T05:40:11.640313Z","shell.execute_reply":"2024-10-24T05:40:11.647084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sample plots\nrandom.seed(23871)\n_, axs = plt.subplots(5, 3, figsize=(3*4, 5*3))\nfor ax in axs.flatten():\n    k = random.randint(0, len(Y)-1)\n    planet_id = train_labels.index[k]\n    ax.plot(Y[k], label='ytrue')\n    ax.plot(YP[k] + 2.5*YSigma[k], linewidth=0.5, label='yp+2.5s')\n    ax.plot(YP[k], label='yp', alpha=0.8)\n    ax.plot(YP[k] - 2.5*YSigma[k], linewidth=0.5, label='yp-2.5s')\n    ax.legend()\n    ax.set_title(f'planet={planet_id}')\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:40:12.634452Z","iopub.execute_input":"2024-10-24T05:40:12.635486Z","iopub.status.idle":"2024-10-24T05:40:18.049017Z","shell.execute_reply.started":"2024-10-24T05:40:12.635424Z","shell.execute_reply":"2024-10-24T05:40:18.047676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# find sigma multiplier\nimport scipy \n\ndef gll(yp, ysigma, yt):\n    return scipy.stats.norm.logpdf(yt, loc=yp, scale=ysigma).sum(1)\n\nmin_sigma = 1e-5\nL_ideal = gll(Y, min_sigma, Y).mean()  # 2998.0983016896525\nL_ref = gll(Y.flatten().mean(), Y.flatten().std(), Y).mean() # 1398.8395122870256\n\n# find sigma multiplier\nscores = []\nmultipliers = np.linspace(0.2, 2.0, 100)\nfor m in multipliers:\n    score = (gll(YP, np.clip(YSigma * m, min_sigma * 5, None), Y) - L_ref) / (L_ideal - L_ref)\n    score = np.clip(score, 0, None).mean()\n    scores.append(score)\n\nbest_score, multiplier = np.max(scores), multipliers[np.argmax(scores)]\nprint(f'best_score={best_score:.4f}, sigma multiplier={multiplier:.4f}')","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:41:06.935511Z","iopub.execute_input":"2024-10-24T05:41:06.935937Z","iopub.status.idle":"2024-10-24T05:41:07.713009Z","shell.execute_reply.started":"2024-10-24T05:41:06.935899Z","shell.execute_reply":"2024-10-24T05:41:07.71177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(multipliers, scores)","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:41:19.214652Z","iopub.execute_input":"2024-10-24T05:41:19.215762Z","iopub.status.idle":"2024-10-24T05:41:19.527392Z","shell.execute_reply.started":"2024-10-24T05:41:19.215712Z","shell.execute_reply":"2024-10-24T05:41:19.52607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import json\nwith open(f'{datadir}/train_{time_binning}/calib-config.json', 'rt') as f:\n    conf = CalibConfig(**json.load(f))\n\nXtest = []\nfor planet_id in test_planets:\n    airs_signal = read_signal(planet_id, 'airs', conf)\n    Xtest.append( feature_extract(airs_signal, lambda_labels) )\n    \nXtest = np.stack(Xtest, 0)","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:42:37.660303Z","iopub.execute_input":"2024-10-24T05:42:37.661371Z","iopub.status.idle":"2024-10-24T05:42:44.100686Z","shell.execute_reply.started":"2024-10-24T05:42:37.661304Z","shell.execute_reply":"2024-10-24T05:42:44.099382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"YPtest, YVtest = None, None\nfor regr in fit_cloned:\n    YP1 = regr.predict(Xtest)\n    YS1 = YP1[:, 283:]\n    YP1 = YP1[:, :283]\n    if YPtest is None:\n        YPtest = YP1.copy()\n        YVtest = YS1**2\n    else:\n        YPtest += YP1\n        YVtest += YS1**2\n\nYPtest /= len(fit_cloned)\nYPtest = np.clip(YPtest, 0, None)\nYStest = np.sqrt( YVtest / len(fit_cloned) )\nYStest = np.clip(YStest * multiplier, min_sigma * 5, None) ","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:43:04.056263Z","iopub.execute_input":"2024-10-24T05:43:04.056757Z","iopub.status.idle":"2024-10-24T05:43:05.68491Z","shell.execute_reply.started":"2024-10-24T05:43:04.056711Z","shell.execute_reply":"2024-10-24T05:43:05.683649Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submit = pd.read_csv(f'{inpdir}/sample_submission.csv')\ndf = pd.concat([\n    pd.DataFrame({'planet_id':list(test_planets)}), \n    pd.DataFrame(np.concatenate([YPtest, YStest], 1), columns=submit.columns[1:])], axis=1)\ndf.to_csv('submission.csv', index=None)","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:43:06.490927Z","iopub.execute_input":"2024-10-24T05:43:06.491499Z","iopub.status.idle":"2024-10-24T05:43:06.531993Z","shell.execute_reply.started":"2024-10-24T05:43:06.491451Z","shell.execute_reply":"2024-10-24T05:43:06.53078Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# remove training data cache\n!rm -rf {datadir}","metadata":{"execution":{"iopub.status.busy":"2024-10-24T05:43:09.667828Z","iopub.execute_input":"2024-10-24T05:43:09.668282Z","iopub.status.idle":"2024-10-24T05:43:11.065035Z","shell.execute_reply.started":"2024-10-24T05:43:09.668239Z","shell.execute_reply":"2024-10-24T05:43:11.063192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}