{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"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":"nvidiaTeslaT4","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":9833462,"sourceType":"datasetVersion","datasetId":6031236}],"dockerImageVersionId":30787,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# NeurIPS ARIEL2024 5th place solution [TRAIN kernel]\n\nThis kernel requires ~64GB of RAM, so it wont run on full data in the kernel. Set `LOW_MEM = False` when running on a larger machine.\n\nPlease find our solution description here: https://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/543760\n\nInference: https://www.kaggle.com/code/ilu000/neurips-ariel24-5th-place-solution-inference","metadata":{}},{"cell_type":"code","source":"LOW_MEM = True","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:10:17.668125Z","iopub.execute_input":"2024-11-07T14:10:17.668957Z","iopub.status.idle":"2024-11-07T14:10:17.677973Z","shell.execute_reply.started":"2024-11-07T14:10:17.668906Z","shell.execute_reply":"2024-11-07T14:10:17.677095Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import gc\nimport numpy as np\nimport pandas as pd\nfrom time import time\nfrom sklearn.model_selection import StratifiedKFold\nfrom scipy.signal import savgol_filter\nfrom sklearn.linear_model import LinearRegression\nfrom scipy.optimize import minimize\nimport scipy.stats\n\nt0 = time()\n\n\ndef score(solution, submission, naive_mean, naive_sigma):\n    # naive_mean = 0.002551714\n    # naive_sigma = 0.00172619\n    sigma_true = 0.00001\n\n    n_wavelengths = 283\n    y_pred = submission[:, :n_wavelengths]\n    sigma_pred = submission[:, n_wavelengths:]\n    y_true = solution\n\n    GLL_pred = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred))\n    GLL_true = np.sum(\n        scipy.stats.norm.logpdf(\n            y_true, loc=y_true, scale=sigma_true * np.ones_like(y_true)\n        )\n    )\n    GLL_mean = np.sum(\n        scipy.stats.norm.logpdf(\n            y_true,\n            loc=naive_mean * np.ones_like(y_true),\n            scale=naive_sigma * np.ones_like(y_true),\n        )\n    )\n\n    submit_score = (GLL_pred - GLL_mean) / (GLL_true - GLL_mean)\n    return submit_score\n\n\ndir1 = \"/kaggle/input/ariel-data-challenge-2024/\"\ndir2 = \"/kaggle/input/neurips-ariel24-5th-place-solution-data/\"\nlog = True\nVERSION = 10","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:10:18.347335Z","iopub.execute_input":"2024-11-07T14:10:18.348085Z","iopub.status.idle":"2024-11-07T14:10:19.427182Z","shell.execute_reply.started":"2024-11-07T14:10:18.348049Z","shell.execute_reply":"2024-11-07T14:10:19.426353Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"gc.enable()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:10:19.428601Z","iopub.execute_input":"2024-11-07T14:10:19.429023Z","iopub.status.idle":"2024-11-07T14:10:19.433016Z","shell.execute_reply.started":"2024-11-07T14:10:19.42899Z","shell.execute_reply":"2024-11-07T14:10:19.432071Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# read data\nx = np.load(dir2 + \"train_signal_v26.npy\")  # 673, 5625, 357\nstar = pd.read_csv(dir1 + \"train_adc_info.csv\", index_col=\"planet_id\")[\"star\"].values\ny = pd.read_csv(dir1 + \"train_labels.csv\", index_col=\"planet_id\").values\nif LOW_MEM:\n    x = x[:100].copy()\n    star = star[:100].copy()\n    y = y[:100].copy()\n    gc.collect()\nif log:\n    print(\"read data\", x.shape, int(time() - t0), \"sec\")\n\n# augmented data\nif LOW_MEM:\n    x_aug = np.load(dir2 + \"train_signal_v26_aug7.npy\")[:100].copy()  # 673, 5625, 357\n    gc.collect()\n    x_aug2 = np.load(dir2 + \"train_signal_v26_aug19.npy\")[:100].copy()  # 673, 5625, 357\n    gc.collect()\n    x_aug3 = np.load(dir2 + \"train_signal_v26_aug20.npy\")[:100].copy()  # 673, 5625, 357\n    gc.collect()\n    x = np.concatenate((x, x_aug, x_aug2, x_aug3), axis=0)\n    gc.collect()\n    y_aug = pd.read_csv(dir2 + \"train_labels_aug7.csv\", index_col=\"planet_id\").values[:100].copy()\n    gc.collect()\n    y_aug2 = pd.read_csv(dir2 + \"train_labels_aug19.csv\", index_col=\"planet_id\").values[:100].copy()\n    gc.collect()\n    y_aug3 = pd.read_csv(dir2 + \"train_labels_aug20.csv\", index_col=\"planet_id\").values[:100].copy()\n    gc.collect()\n    y = np.concatenate((y, y_aug, y_aug2, y_aug3), axis=0)\n    gc.collect()\nelse:\n    x_aug = np.load(dir2 + \"train_signal_v26_aug7.npy\")  # 673, 5625, 357\n    x_aug2 = np.load(dir2 + \"train_signal_v26_aug19.npy\")  # 673, 5625, 357\n    x_aug3 = np.load(dir2 + \"train_signal_v26_aug20.npy\")  # 673, 5625, 357\n    x = np.concatenate((x, x_aug, x_aug2, x_aug3), axis=0)\n    y_aug = pd.read_csv(dir2 + \"train_labels_aug7.csv\", index_col=\"planet_id\").values\n    y_aug2 = pd.read_csv(dir2 + \"train_labels_aug19.csv\", index_col=\"planet_id\").values\n    y_aug3 = pd.read_csv(dir2 + \"train_labels_aug20.csv\", index_col=\"planet_id\").values\n    y = np.concatenate((y, y_aug, y_aug2, y_aug3), axis=0)\n\nNUM_TRAINS = 4\nstar = np.concatenate([star] * NUM_TRAINS, axis=0)\n\nlabels = y.copy()\nif log:\n    print(\"read data_aug\", x.shape, int(time() - t0), \"sec\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:10:19.434211Z","iopub.execute_input":"2024-11-07T14:10:19.434537Z","iopub.status.idle":"2024-11-07T14:13:28.260541Z","shell.execute_reply.started":"2024-11-07T14:10:19.434506Z","shell.execute_reply":"2024-11-07T14:13:28.259515Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# group data\nsh = (x.shape[0], x.shape[1], 1)\nx1 = np.concatenate(\n    (x[:, :, 0].reshape(sh), x[:, :, 1:].mean(-1).reshape(sh)), axis=-1\n)  # 1=F, 2=mean(A)\nx1 = np.concatenate(\n    (x1, x[:, :, 1 + 37 : -37].mean(-1).reshape(sh)), axis=-1\n)  # 3=mean2(A) - true spectrum\n\n# scale all stars to the same (uniform) spectrum.\n# After totals - unscaled data is a lot more stable, use it for bp/all_s\nstar_list = pd.Series(star).unique()\nfor s in star_list:\n    x[star == s, :, :] /= x[star == s, :, :].mean(0).mean(0).reshape(1, 1, -1)\n\n\n# add 3 components to capture the slope: level, and 2 linear ones\n# 4=level - Am3\nx1 = np.concatenate((x1, x[:, :, 1 + 37 : -37].mean(-1).reshape(sh)), axis=-1).astype(\n    np.float32\n)\n\n# construct spectral components.\nx = np.concatenate(\n    (x[:, :, 0].reshape(sh), np.flip(x[:, :, 1 + 37 : -37], -1)), -1\n)  # flip and drop tails\nwl = (\n    pd.read_csv(dir1 + \"wavelengths.csv\").values.astype(np.float32).ravel()\n)  # WL: 283, 0.71 + 1.95 to 3.90\n\n# Spectral lines of gasses from RB table 7-1: co2=2.03, co=2.35, h2o=2.69, 4.25(main), nh3=3.0(main), ch4=3.3(main)\ncc = [2.00, 2.38, 2.58, 2.76, 3.00, 3.32, 3.50, 3.98]\nstds = [0.08, 0.12, 0.13, 0.11, 0.16, 0.16, 0.19, 0.19]\nfor i in range(len(cc)):  # loop over components\n    c = cc[i]\n    std = stds[i]\n\n    # construct gaussian component with mean and std\n    comp = np.exp(-((wl - c) ** 2) / 2 / std**2)\n    comp /= comp.max()  # norm it, just in case\n\n    # blend\n    x2 = (x * comp.reshape(1, 1, -1)).sum(-1)\n    x1 = np.concatenate((x1, x2.reshape(sh)), axis=-1).astype(np.float32)\nx = x1  # final grouped data\nif log:\n    print(\"end grouping\", x.shape, int(time() - t0), \"sec\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:13:28.262552Z","iopub.execute_input":"2024-11-07T14:13:28.262835Z","iopub.status.idle":"2024-11-07T14:14:00.01941Z","shell.execute_reply.started":"2024-11-07T14:13:28.262804Z","shell.execute_reply":"2024-11-07T14:14:00.018272Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# find transit zones: use the derivative\ndef phase_detector(signal):\n    MIN = signal.shape[0] // 2\n    signal1 = signal[:MIN]\n    signal2 = signal[MIN:]\n    first_derivative1 = np.gradient(signal1)\n    first_derivative2 = np.gradient(signal2)\n    phase1 = np.argmin(first_derivative1)\n    phase2 = np.argmax(first_derivative2) + MIN\n    return phase1, phase2\n\n\nbp = np.zeros(x.shape[0], dtype=np.int32)\nbp2 = np.zeros(x.shape[0], dtype=np.int32)\n\nfor i in range(x.shape[0]):\n    signal = (\n        x[i, :, 2] - 0.2 * x[i, :, 1]\n    )\n    signal = signal.reshape(-1, 15).mean(-1)  # collapse by 15; 375\n    signal = savgol_filter(signal, 11, 2)  # smooth\n    p1, p2 = phase_detector(signal)\n    bp[i] = 15 * p1\n    bp2[i] = 15 * p2\nif log:\n    print(\n        \"end bp\",\n        x.shape,\n        int(0.49 + np.mean(np.minimum(1842, np.maximum(1800, bp)))),\n        int(time() - t0),\n        \"sec\",\n    )\n\n# # save bps\n# np.save(\"bp_\" + str(VERSION) + \".npy\", bp)\n# np.save(\"bp2_\" + str(VERSION) + \".npy\", bp2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:14:00.020713Z","iopub.execute_input":"2024-11-07T14:14:00.021037Z","iopub.status.idle":"2024-11-07T14:14:00.288817Z","shell.execute_reply.started":"2024-11-07T14:14:00.021003Z","shell.execute_reply":"2024-11-07T14:14:00.287893Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def try_s(s):\n    y0 = signal * (1 + s * mult)\n    z = np.polyfit(xx0, y0, deg)\n    p = np.poly1d(z)\n    d = p(xx0) - y0\n    q = np.abs(d).mean()  # MAE - minimize it. Seems better than MSE\n    return q * 1000  # this increases solver precision\n\n\ndeg = 4\nNS2 = 70  # number of periods used for smoothing. param.\nw = 70  # width of transition period\nall_s = []\nfor i in range(x.shape[0]):\n    xx0 = np.arange(-x.shape[1] // 2, x.shape[1] // 2)\n    mult = np.zeros(x.shape[1], dtype=np.float32)\n    mult[bp[i] : bp2[i]] = 1\n\n    signal = x[i, :, 2] - 0.2 * x[i, :, 1]\n    signal = savgol_filter(signal, NS2, 2)\n\n    idx = [\n        j\n        for j in range(x.shape[1])\n        if j < bp[i] - w or j > bp2[i] + w or (j > bp[i] + w and j < bp2[i] - w)\n    ]\n    xx0 = xx0[idx]\n    signal = signal[idx]\n    mult = mult[idx]\n\n    r = minimize(try_s, [0.0003], method=\"Nelder-Mead\")\n    all_s.append(r.x[0])\nall_s = np.array(all_s).astype(np.float32)\nall_s2 = []\nfor i in range(x.shape[0]):\n    xx0 = np.arange(-x.shape[1] // 2, x.shape[1] // 2)\n    mult = np.zeros(x.shape[1], dtype=np.float32)\n    mult[bp[i] : bp2[i]] = 1\n\n    signal = x[i, :, 3]\n    signal = savgol_filter(signal, NS2, 2)\n\n    idx = [\n        j\n        for j in range(x.shape[1])\n        if j < bp[i] - w or j > bp2[i] + w or (j > bp[i] + w and j < bp2[i] - w)\n    ]\n    xx0 = xx0[idx]\n    signal = signal[idx]\n    mult = mult[idx]\n\n    r = minimize(try_s, [0.0003], method=\"Nelder-Mead\")\n    all_s2.append(r.x[0])\nall_s2 = np.array(all_s2).astype(np.float32)\nif log:\n    print(\"end all_s\", int(time() - t0), \"sec\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:14:00.28981Z","iopub.execute_input":"2024-11-07T14:14:00.29011Z","iopub.status.idle":"2024-11-07T14:14:59.093871Z","shell.execute_reply.started":"2024-11-07T14:14:00.290079Z","shell.execute_reply":"2024-11-07T14:14:59.092884Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"w = 10\nunobscured = []\nobscured = []\nfor i in range(x.shape[0]):\n    obscuredi = x[i, bp[i] + w : bp2[i] - w, :].mean(axis=0)\n    unobscuredi = (\n        x[i, : bp[i] - w, :].mean(axis=0) + x[i, bp2[i] + w :, :].mean(axis=0)\n    ) / 2\n    obscured.append(obscuredi)\n    unobscured.append(unobscuredi)\nunobscured = np.array(unobscured)\ndd = (unobscured - np.array(obscured)) / unobscured  # ratio\nif log:\n    print(\"end diff of obscured\", int(time() - t0), \"sec\")\n\n\ndeg2 = 2  # degree of polyfit for dd. param\nfor i in range(x.shape[0]):\n    xx0 = np.arange(-x.shape[1] // 2, x.shape[1] // 2)\n    mult = np.zeros(x.shape[1], dtype=np.float32)\n    mult[bp[i] : bp2[i]] = 1\n\n    signal = x[i, :, :]  # all components\n    signal = savgol_filter(signal, NS2, 2, axis=0)\n    y0 = signal * (1 + all_s[i] * mult.reshape(-1, 1))\n\n    idx = [\n        j\n        for j in range(x.shape[1])\n        if j < bp[i] - w or j > bp2[i] + w or (j > bp[i] + w and j < bp2[i] - w)\n    ]\n    xx0 = xx0[idx]\n    y0 = y0[idx]\n\n    z = np.polyfit(xx0, y0, deg2)  # use deg2 here\n    pred = 0\n    xx0 = np.arange(-x.shape[1] // 2, x.shape[1] // 2)\n    for j in range(deg2):  # excl const\n        pred = xx0.reshape(-1, 1) * (z[j, :] + pred)\n    x[i, :, :] -= pred\n\n\nN = 250\nx = savgol_filter(x, N, 2, axis=1)\nN2 = int(\n    ((N * 2.5) // 2) * 2\n)  # fully capture the transition\nx = x[:, N2:, :] - x[:, :-N2, :]  # diff: 673, 5625-2*N, 17 (cuts N from each end)\nd = []\nfor i in range(x.shape[0]):\n    R = int(N * 0.75)  # window for finding peaks\n    d1i = x[i, bp2[i] - N2 // 2 - R : bp2[i] - N2 // 2 + R, :].mean(\n        axis=0\n    )\n    d2i = -x[i, bp[i] - N2 // 2 - R : bp[i] - N2 // 2 + R, :].mean(axis=0)\n    di = (d1i + d2i) / 2 / unobscured[i]  # ratio\n    d.append(di)\ndel x\ngc.collect()\nd = np.array(d)\n\n\n# define new all_s as best blend of all vars incl max\n# use N best features from all_s/d/dd: e = 27.0\nall_s = (\n    np.maximum(all_s, d[:, 3]) * 0.333\n    + np.minimum(all_s, d[:, 2]) * 0.204\n    + np.minimum(all_s2, dd[:, 3]) * 0.035\n    + np.maximum(all_s2, d[:, 1]) * 0.152\n    + np.maximum(d[:, 1], d[:, 3]) * 0.28\n)\n\nif log:\n    print(\n        \"all_s_blend\",\n        np.round(np.sqrt(((y.mean(-1) - all_s) ** 2).mean()) * 1e6, 1),\n        int(time() - t0),\n        \"sec\",\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:14:59.09505Z","iopub.execute_input":"2024-11-07T14:14:59.095354Z","iopub.status.idle":"2024-11-07T14:15:40.60321Z","shell.execute_reply.started":"2024-11-07T14:14:59.09532Z","shell.execute_reply":"2024-11-07T14:15:40.6023Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# # define all_s as best blend of all vars incl max - solver\n# x = np.zeros([d.shape[0], 1000], dtype=np.float32)\n# x[:, 0] = all_s\n# x[:, 1] = all_s2\n# x[:, 2:4] = d[:, :2]\n# x[:, 4:6] = dd[:, :2]\n# llf = [\"all_s0\", \"all_s1\", \"d0\", \"d1\", \"dd0\", \"dd1\"]\n# k0, k = 6, 6\n# for i in range(k0 - 1):  # add max of all combinations\n#     for j in range(i + 1, k0):\n#         print(k, i, j, llf[i], llf[j])\n#         x[:, k] = np.maximum(x[:, i], x[:, j])\n#         k += 1\n# x = x[:, :k]\n# x0 = x.copy()\n# # select larger and larger set of features - best approach\n# ll0 = [j for j in range(x0.shape[1])]\n# llb_prior = set()\n# import itertools\n\n# for i1 in range(1, 6):  # use so many features\n#     s0, n = 10000, 0\n#     for ll in itertools.combinations(ll0, i1):\n#         n += 1\n#         x = x0[:, ll].copy()\n#         reg = LinearRegression(fit_intercept=False, positive=True).fit(\n#             x, y.mean(-1)\n#         )  # positive?\n#         pred = reg.predict(x)\n#         s = np.round(np.sqrt(((y.mean(-1) - pred) ** 2).mean()) * 1e6, 1)\n#         if s < s0:\n#             s0, llb, ccb = s, list(ll).copy(), np.round(reg.coef_, 3)\n#     print(i1, s0, n, llb, ccb, int(time() - t0), \"sec\")\n#     llb_prior = set(llb.copy())\n# stop","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:15:40.606257Z","iopub.execute_input":"2024-11-07T14:15:40.606644Z","iopub.status.idle":"2024-11-07T14:15:40.611871Z","shell.execute_reply.started":"2024-11-07T14:15:40.606589Z","shell.execute_reply":"2024-11-07T14:15:40.610972Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# redefine d/dd (0=Am1, 1=Am2, 2-6=c1-5) in terms of all_s/Am2\nind = [4, 5, 6, 7, 8, 9, 10, 11]  # center c1-c5 on Am3\nd[:, ind] = d[:, ind] - d[:, 3].reshape(-1, 1)\ndd[:, ind] = dd[:, ind] - dd[:, 3].reshape(-1, 1)\nind = [0, 1, 2, 3]  # center F/Am on all_s\nd[:, ind] = d[:, ind] - all_s.reshape(-1, 1)\ndd[:, ind] = dd[:, ind] - all_s.reshape(-1, 1)\n\n# combine d and dd\nx = np.concatenate((d, dd), axis=1)\n\n# drop some features to improve CV\nind = [0, 1, 4, 5, 7, 11, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23]\nx = x[:, ind]\n\nif log:\n    print(\"end filter\", x.shape, int(time() - t0), \"sec\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:15:40.613159Z","iopub.execute_input":"2024-11-07T14:15:40.613451Z","iopub.status.idle":"2024-11-07T14:15:40.623547Z","shell.execute_reply.started":"2024-11-07T14:15:40.61342Z","shell.execute_reply":"2024-11-07T14:15:40.622726Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# model training\ny0 = y.copy()\ny = y - all_s.reshape(-1, 1)  # subtract\npred1 = np.zeros([x.shape[0], y.shape[1]], dtype=np.float32)\nc12 = []\n\nfor fold in range(1 + star.max()):\n    reg = LinearRegression(fit_intercept=False).fit(\n        x[star != fold, :], y[star != fold, :]\n    )\n    pred1[star == fold, :] = reg.predict(x[star == fold, :])\n    c12.append(reg.coef_)\nnp.save(\"c12_\" + str(VERSION) + \".npy\", np.array(c12).mean(axis=0))\n\n# assign sigma\npred1 = pred1 + all_s.reshape(-1, 1)  # add it back\npred = np.concatenate((pred1, pred1), axis=-1)\npred[:, 283:] = np.sqrt(((y0 - pred[:, :283]) ** 2).mean())\nprint(\n    \"best score0\",\n    np.round(\n        score(\n            labels[: len(labels) // NUM_TRAINS],\n            pred[: len(labels) // NUM_TRAINS],\n            labels[: len(labels) // NUM_TRAINS].mean(),\n            labels[: len(labels) // NUM_TRAINS].std(),\n        ),\n        4,\n    ),\n    np.round(np.sqrt(((y0 - pred[:, :283]) ** 2).mean()) * 1e6, 1),\n    int(time() - t0),\n    \"sec*****************************************\",\n)\nprint(\n    \"best score0\",\n    np.round(\n        score(\n            labels[len(labels) // NUM_TRAINS :],\n            pred[len(labels) // NUM_TRAINS :],\n            labels[len(labels) // NUM_TRAINS :].mean(),\n            labels[len(labels) // NUM_TRAINS :].std(),\n        ),\n        4,\n    ),\n    np.round(np.sqrt(((y0 - pred[:, :283]) ** 2).mean()) * 1e6, 1),\n    int(time() - t0),\n    \"sec*****************************************\",\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:15:40.624889Z","iopub.execute_input":"2024-11-07T14:15:40.625225Z","iopub.status.idle":"2024-11-07T14:15:40.722Z","shell.execute_reply.started":"2024-11-07T14:15:40.625193Z","shell.execute_reply":"2024-11-07T14:15:40.721234Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# # drop features 1 at a time ****************************************************************************\n# # instead, select larger and larger set of features - best approach\n# x0 = x.copy()\n# ll0 = [j for j in range(x0.shape[1])]\n# llb_prior = set()\n# import itertools\n\n# for i1 in range(1, x0.shape[1]):  # use so many features\n#     s0, n = -1000, 0\n#     for ll in itertools.combinations(ll0, i1):\n#         if (\n#             len(ll) > 1 and len(set(ll).intersection(llb_prior)) < len(llb_prior) - 1\n#         ):  # allow drop of no more than 1 var from prior best - quickest possible approach\n#             continue\n#         n += 1\n#         x = x0[:, ll].copy()\n#         for fold in [0, 1]:\n#             reg = LinearRegression(fit_intercept=False).fit(\n#                 x[star == fold, :], y[star == fold, :]\n#             )\n#             pred1[star != fold, :] = reg.predict(x[star != fold, :])\n#         pred1 = pred1 + all_s.reshape(-1, 1)\n#         pred = np.concatenate((pred1, pred1), axis=-1)\n#         pred[:, 283:] = np.sqrt(((y0 - pred1) ** 2).mean()).reshape(-1, 1)\n#         s = np.round(\n#             score(\n#                 labels[:len(labels) // 2],\n#                 pred[:len(labels) // 2],\n#                 labels[:len(labels) // 2].mean(),\n#                 labels[:len(labels) // 2].std(),\n#             ),\n#             4,\n#         )\n#         if s > s0:\n#             s0, llb = s, list(ll).copy()\n#     print(i1, s0, n, llb, int(time() - t0), \"sec\")\n#     llb_prior = set(llb.copy())\n# stop","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:15:40.723093Z","iopub.execute_input":"2024-11-07T14:15:40.724616Z","iopub.status.idle":"2024-11-07T14:15:40.73318Z","shell.execute_reply.started":"2024-11-07T14:15:40.724569Z","shell.execute_reply":"2024-11-07T14:15:40.732346Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# predict 3 models - stars 0, 1, 0+1 - LB seems better this way\nfolds = 5\npred1 = np.zeros([x.shape[0], y.shape[1]], dtype=np.float32)\nkf = StratifiedKFold(n_splits=folds, shuffle=True, random_state=13)\nyy = y0.mean(axis=1) * 1000  # 0.4 to 7.4\nyy = 1 * (yy > 1.23) + 1 * (yy > 1.92) + 1 * (yy > 3.6)  # split into 4 equal groups\nc = []\nfor fold, (tr_ind, va_ind) in enumerate(kf.split(x, yy)):\n    reg = LinearRegression(fit_intercept=False).fit(x[tr_ind, :], y[tr_ind, :])\n    c.append(reg.coef_)\nnp.save(\"c_\" + str(VERSION) + \".npy\", np.array(c).mean(axis=0))\n# star 0 only\nxs0 = x[star == 0, :]\nys0 = y[star == 0, :]\npreds0 = np.zeros([xs0.shape[0], y.shape[1]], dtype=np.float32)\nyy = y0[star == 0, :].mean(axis=1) * 1000  # 1.2 to 7.4\nyy = 1 * (yy > 2.4) + 1 * (yy > 3.5) + 1 * (yy > 5.1)  # split into 4 equal groups\nc0 = []\nfor fold, (tr_ind, va_ind) in enumerate(kf.split(xs0, yy)):\n    reg = LinearRegression(fit_intercept=False).fit(xs0[tr_ind, :], ys0[tr_ind, :])\n    c0.append(reg.coef_)\nnp.save(\"c0_\" + str(VERSION) + \".npy\", np.array(c0).mean(axis=0))\n# star 1 only\nxs1 = x[star == 1, :]\nys1 = y[star == 1, :]\npreds1 = np.zeros([xs1.shape[0], y.shape[1]], dtype=np.float32)\nyy = y0[star == 1, :].mean(axis=1) * 1000  # 0.4 to 2.6\nyy = 1 * (yy > 0.85) + 1 * (yy > 1.21) + 1 * (yy > 1.66)  # split into 4 equal groups\nc1 = []\nfor fold, (tr_ind, va_ind) in enumerate(kf.split(xs1, yy)):\n    reg = LinearRegression(fit_intercept=False).fit(xs1[tr_ind, :], ys1[tr_ind, :])\n    c1.append(reg.coef_)\nnp.save(\"c1_\" + str(VERSION) + \".npy\", np.array(c1).mean(axis=0))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:15:40.734826Z","iopub.execute_input":"2024-11-07T14:15:40.736119Z","iopub.status.idle":"2024-11-07T14:15:40.880533Z","shell.execute_reply.started":"2024-11-07T14:15:40.736076Z","shell.execute_reply":"2024-11-07T14:15:40.879788Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# predict sigma: A * std + B\nstd_mult0 = np.load(dir2 + \"std_mult0.npy\")  # .77 to 2.1\n\n\ndef obj3(m):\n    global pred\n    pred[:, 283:] = m[0] * pred[:, :283].std(-1).reshape(-1, 1) + m[1] * 1e-6\n    return -score(\n        labels[: len(labels) // NUM_TRAINS],\n        pred[: len(labels) // NUM_TRAINS],\n        labels[: len(labels) // NUM_TRAINS].mean(),\n        labels[: len(labels) // NUM_TRAINS].std(),\n    )\n\n\nr = minimize(obj3, [0.3, 23], method=\"Nelder-Mead\")\nprint(\n    \"best score1\",\n    np.round(-r.fun, 4),\n    np.round(r.x, 3),\n    int(time() - t0),\n    \"sec**********************************************\",\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:22:26.6771Z","iopub.execute_input":"2024-11-07T14:22:26.677981Z","iopub.status.idle":"2024-11-07T14:22:26.957199Z","shell.execute_reply.started":"2024-11-07T14:22:26.677925Z","shell.execute_reply":"2024-11-07T14:22:26.956121Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# NN prep\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras.callbacks import Callback, LearningRateScheduler, EarlyStopping\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.optimizers.schedules import CosineDecay\nfrom tensorflow.keras.layers import Input, Dense, Concatenate, Dropout\n\n# norm x\nx = np.concatenate((x, 0.1 * all_s.reshape(-1, 1)), axis=1)  # incl all_s, JIC\nx *= 6000\n\n# NN predict sigma only\np2 = np.zeros([x.shape[0], y.shape[1]], dtype=np.float32)\nmse = (y0 - pred[:, :283]) ** 2  # target\n\nNL = 7\nD = 64\nact = \"swish\"\nBS = 64  # out of 673\nfor fold in [0, 1]:\n    x_t, x_v = x[star == fold, :], x[star != fold, :]\n    y_t, y_v = mse[star == fold, :], mse[star != fold, :]\n    with tf.device(\"/GPU:0\"):\n\n        def customloss1(mse, p):\n            var = p * p\n            d2 = (tf.math.log(var) + mse / var + 11.724) / 11.3017\n            loss = tf.math.reduce_mean(d2, axis=-1)\n            return loss  # =-1 * score, so minimize this\n\n        i1 = Input(shape=(x.shape[1],))\n        d1 = Dense(D, activation=act)(i1)\n        for i in range(NL - 1):\n            d1 = Dense(D, activation=act)(d1)\n        o = 7e-6 + 12e-6 * Dense(y_t.shape[1], activation=\"softplus\")(\n            d1\n        )  # low mult here avoids a blow-up\n        model = tf.keras.Model(inputs=i1, outputs=o)\n        sc = CosineDecay(initial_learning_rate=5e-4, decay_steps=520)\n        model.compile(optimizer=Adam(learning_rate=sc), loss=customloss1)\n        history = model.fit(\n            x=x_t,\n            y=y_t,\n            validation_data=(x_v, y_v),\n            epochs=52,\n            callbacks=[],\n            verbose=0,\n            batch_size=BS,\n        )\n        model.save(\n            \"model_\" + str(fold) + \"_\" + str(VERSION) + \".keras\"\n        )\n        p2[star != fold, :] = model.predict(x_v, batch_size=BS, verbose=0)\npred[:, 283:] = p2\nprint(\n    \"best score3\",\n    np.round(\n        score(\n            labels[: len(labels) // NUM_TRAINS],\n            pred[: len(labels) // NUM_TRAINS],\n            labels[: len(labels) // NUM_TRAINS].mean(),\n            labels[: len(labels) // NUM_TRAINS].std(),\n        ),\n        3,\n    ),\n    int(time() - t0),\n    \"sec\",\n)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2024-11-07T14:22:30.410657Z","iopub.execute_input":"2024-11-07T14:22:30.41107Z","iopub.status.idle":"2024-11-07T14:22:57.327389Z","shell.execute_reply.started":"2024-11-07T14:22:30.411032Z","shell.execute_reply":"2024-11-07T14:22:57.326437Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}