{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.14"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":70367,"databundleVersionId":9188054}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":5293.975235,"end_time":"2024-10-08T07:51:20.091348","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-10-08T06:23:06.116113","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom tqdm import tqdm\nimport joblib\n\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\nimport itertools\n\nfrom scipy.optimize import minimize\nfrom scipy import optimize\n\nfrom astropy.stats import sigma_clip","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":3.72139,"end_time":"2024-10-08T06:23:13.351214","exception":false,"start_time":"2024-10-08T06:23:09.629824","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-11-10T07:52:54.189971Z","iopub.execute_input":"2024-11-10T07:52:54.19054Z","iopub.status.idle":"2024-11-10T07:52:54.460798Z","shell.execute_reply.started":"2024-11-10T07:52:54.190478Z","shell.execute_reply":"2024-11-10T07:52:54.459681Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"dataset = 'test'\nadc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/'+f'{dataset}_adc_info.csv',index_col='planet_id')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')","metadata":{"papermill":{"duration":0.235059,"end_time":"2024-10-08T06:23:13.593638","exception":false,"start_time":"2024-10-08T06:23:13.358579","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-11-10T07:52:54.462104Z","iopub.execute_input":"2024-11-10T07:52:54.462775Z","iopub.status.idle":"2024-11-10T07:52:54.719362Z","shell.execute_reply.started":"2024-11-10T07:52:54.462732Z","shell.execute_reply":"2024-11-10T07:52:54.717869Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if dataset == \"train\":\n    adc_info = adc_info[:10]","metadata":{"execution":{"iopub.status.busy":"2024-11-10T07:52:54.721846Z","iopub.execute_input":"2024-11-10T07:52:54.722269Z","iopub.status.idle":"2024-11-10T07:52:54.728192Z","shell.execute_reply.started":"2024-11-10T07:52:54.72219Z","shell.execute_reply":"2024-11-10T07:52:54.72658Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def apply_linear_corr(linear_corr,clean_signal):\n    linear_corr = np.flip(linear_corr, axis=0)\n    for x, y in itertools.product(\n                range(clean_signal.shape[1]), range(clean_signal.shape[2])\n            ):\n        poli = np.poly1d(linear_corr[:, x, y])\n        clean_signal[:, x, y] = poli(clean_signal[:, x, y])\n    return clean_signal\n\ndef clean_dark(signal, dark, dt):\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n    signal -= dark* dt[:, np.newaxis, np.newaxis]\n    return signal\n\ndef preproc(dataset, adc_info, sensor, binning = 15):\n    cut_inf, cut_sup = 39, 321\n    sensor_sizes_dict = {\"AIRS-CH0\":[[11250, 32, 356], [1, 32, cut_sup-cut_inf]], \"FGS1\":[[135000, 32, 32], [1, 32, 32]]}\n    binned_dict = {\"AIRS-CH0\":[11250 // binning // 2, 282], \"FGS1\":[135000 // binning // 2]}\n    linear_corr_dict = {\"AIRS-CH0\":(6, 32, 356), \"FGS1\":(6, 32, 32)}\n    planet_ids = adc_info.index\n    \n    feats = []\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/{sensor}_signal.parquet').to_numpy()\n        dark_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/dark.parquet', engine='pyarrow').to_numpy()\n        dead_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/dead.parquet', engine='pyarrow').to_numpy()\n        flat_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/flat.parquet', engine='pyarrow').to_numpy()\n        linear_corr = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/linear_corr.parquet').values.astype(np.float64).reshape(linear_corr_dict[sensor])\n\n        signal = signal.reshape(sensor_sizes_dict[sensor][0]) \n        gain = adc_info[f'{sensor}_adc_gain'].values[i]\n        offset = adc_info[f'{sensor}_adc_offset'].values[i]\n        signal = signal / gain + offset\n        \n        hot = sigma_clip(\n            dark_frame, sigma=5, maxiters=5\n        ).mask\n        \n        if sensor != \"FGS1\":\n            signal = signal[:, :, cut_inf:cut_sup] \n            dt = np.ones(len(signal))*0.1 \n            dt[1::2] += 4.5 #@bilzard idea\n            linear_corr = linear_corr[:, :, cut_inf:cut_sup]\n            dark_frame = dark_frame[:, cut_inf:cut_sup]\n            dead_frame = dead_frame[:, cut_inf:cut_sup]\n            flat_frame = flat_frame[:, cut_inf:cut_sup]\n            hot = hot[:, cut_inf:cut_sup]\n        else:\n            dt = np.ones(len(signal))*0.1\n            dt[1::2] += 0.1\n            \n        # signal = signal.clip(0) #@graySnow idea\n        linear_corr_signal = apply_linear_corr(linear_corr, signal)\n        signal = clean_dark(linear_corr_signal, dark_frame, dt)\n        \n        flat = flat_frame.reshape(sensor_sizes_dict[sensor][1])\n        flat[dead_frame.reshape(sensor_sizes_dict[sensor][1])] = np.nan\n        #flat[hot.reshape(sensor_sizes_dict[sensor][1])] = np.nan\n        signal = signal / flat\n        \n        \n        if sensor == \"FGS1\":\n            signal = signal[:,10:22,10:22] # **** updates ****\n            signal = signal.reshape(sensor_sizes_dict[sensor][0][0],144) # # **** updates ****\n\n        if sensor != \"FGS1\":\n            # backgrounds are [0:8] and [24:32]\n            signal_bg = np.nanmean(\n                np.concatenate([signal[:, 0:8, :], signal[:, 24:32, :]], axis=1), axis=1\n            )\n            signal_bg[np.isnan(signal_bg)] = 0\n            signal = signal[:, 8:24, :]  # **** updates ****\n\n        mean_signal = np.nanmean(signal, axis=1)\n        cds_signal = mean_signal[1::2] - mean_signal[0::2]\n        cds_signal_bg = signal_bg[1::2] - signal_bg[0::2]\n\n        cds_signal_bg = np.nanmean(cds_signal_bg, axis=0, keepdims=True)\n\n        cds_signal -= cds_signal_bg\n\n        binned = np.zeros((binned_dict[sensor]))\n        for j in range(cds_signal.shape[0] // binning):\n            binned[j] = cds_signal[j * binning : j * binning + binning].mean(axis=0)\n\n        if sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0], 1))\n\n        feats.append(binned)\n        \n    return np.stack(feats)\n\nsignal_list = preproc(f'{dataset}', adc_info, \"AIRS-CH0\", 1)","metadata":{"papermill":{"duration":5151.643632,"end_time":"2024-10-08T07:49:05.388291","exception":false,"start_time":"2024-10-08T06:23:13.744659","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-11-10T07:52:54.729987Z","iopub.execute_input":"2024-11-10T07:52:54.730489Z","iopub.status.idle":"2024-11-10T07:53:02.138082Z","shell.execute_reply.started":"2024-11-10T07:52:54.730434Z","shell.execute_reply":"2024-11-10T07:53:02.136772Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Modelization","metadata":{"papermill":{"duration":0.072544,"end_time":"2024-10-08T07:49:05.536937","exception":false,"start_time":"2024-10-08T07:49:05.464393","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import numpy as np\nfrom scipy.optimize import minimize\nimport matplotlib.pyplot as plt\n\n# コスト関数\ndef cost_function(params, data):\n    t1, t2, a, b = params\n    t1 = int(t1)\n    t2 = int(t2)\n    \n    # 制約違反\n    if t1 > t2:\n        return(t1 - t2) * 1e9\n    \n    y = np.full((data.shape[0], ), a)\n    y[t1:t2] = np.linspace(a, b, t2-t1)\n    y[t2:] = b\n    \n    cost = np.sum((data - y) ** 2)\n    \n    return cost\n\ndef phase_sec_detector(data):\n    \"\"\"\n    Return\n    t1:入り始め\n    t2:入り終わり\n    t3:出始め\n    t4:出終わり\n    \"\"\"\n    \n    window_size = 100\n    kernel = np.ones(window_size) / window_size\n    # 移動平均で平す\n    moving_average = np.convolve(data, kernel, mode='valid')\n    \n    # 微分を取って平す\n    diff = np.diff(moving_average)\n    moving_average_diff = np.convolve(diff, kernel, mode='valid')\n    \n    min_idx = np.argmin(moving_average_diff)\n    max_idx = np.argmax(moving_average_diff)\n    \n    # min_idx と max_idxを元のデータのインデックスに対応させる\n    min_idx_on_data = min_idx + window_size\n    max_idx_on_data = max_idx + window_size\n    \n    signal = data.copy()\n    \n    # データを300点前後にトリミング\n    width = 300\n    data = signal[min_idx_on_data-width:min_idx_on_data+width]\n\n    # 初期パラメータ (t1, t2, a, b)\n    initial_params = [width//2, width*3//2, np.mean(data[:width//2]), np.mean(data[width*3//2:])]\n\n    # 境界条件\n    bounds = [(1, len(data) - 2),  # t1 の範囲\n              (2, len(data) - 1),  # t2 の範囲\n              (min(data), max(data)),  # a の範囲\n              (min(data), max(data))]  # b の範囲\n\n    # 最適化\n    result = minimize(cost_function, initial_params, args=(data,), bounds=bounds, method='Nelder-Mead')\n\n    # 最適化結果\n    t1_opt, t2_opt, a_opt, b_opt = result.x\n    t1_opt = int(t1_opt)\n    t2_opt = int(t2_opt)\n    \n    # データを300点前後にトリミング\n    width = 300\n    data = signal[max_idx_on_data-width:max_idx_on_data+width]\n\n    # 初期パラメータ (t1, t2, a, b)\n    initial_params = [width//2, width*3//2, np.mean(data[:width//2]), np.mean(data[width*3//2:])]\n\n    # 境界条件\n    bounds = [(1, len(data) - 2),  # t1 の範囲\n              (2, len(data) - 1),  # t2 の範囲\n              (min(data), max(data)),  # a の範囲\n              (min(data), max(data))]  # b の範囲\n\n    # 最適化\n    result = minimize(cost_function, initial_params, args=(data,), bounds=bounds, method='Nelder-Mead')\n\n    # 最適化結果\n    t3_opt, t4_opt, a_opt, b_opt = result.x\n    t3_opt = int(t3_opt)\n    t4_opt = int(t4_opt)\n\n    return min_idx_on_data-width+t1_opt, min_idx_on_data-width+t2_opt, max_idx_on_data-width+t3_opt, max_idx_on_data-width+t4_opt","metadata":{"papermill":{"duration":0.102073,"end_time":"2024-10-08T07:49:05.709817","exception":false,"start_time":"2024-10-08T07:49:05.607744","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-11-10T07:53:10.230026Z","iopub.execute_input":"2024-11-10T07:53:10.231064Z","iopub.status.idle":"2024-11-10T07:53:10.251399Z","shell.execute_reply.started":"2024-11-10T07:53:10.231012Z","shell.execute_reply":"2024-11-10T07:53:10.250041Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from scipy.optimize import least_squares\n\n\ndef objective(s, signal, t1, t2, t3, t4):\n    y = signal[:t1].tolist() + (signal[t2:t3] / (1 - s)).tolist() + signal[t4:].tolist()\n    x = list(range(0, t1)) + list(range(t2, t3)) + list(range(t4, signal.shape[0]))\n\n    z = np.polyfit(x, y, deg=4)\n    p = np.poly1d(z)\n    q = ((p(x) - y) ** 2).mean()\n\n    return q\n\n\n# 目的関数の定義\ndef model(params, t, w, d_t, d_w, spectrum):\n    # f の係数抽出\n    coeffs_t = params[: d_t + 1]\n    f = np.polyval(coeffs_t, t)\n    f_sum = np.sum(f)\n    coeffs_t[-1] -= f_sum / t.shape[0]\n    f = np.polyval(coeffs_t, t)\n\n    # g の係数抽出\n    coeffs_w = params[d_t + 1 : d_t + d_w + 2]\n    g = np.polyval(coeffs_w, w)\n\n    # 恒星強度\n    star = params[-1]\n\n    # モデル式\n    return star * spectrum * (1 + np.outer(f, g))\n\n\n# 誤差関数の定義\ndef loss_function(params, t, w, t1, t2, t3, t4, d_t, d_w, spectrum, y_true):\n    y_pred = model(params, t, w, d_t, d_w, spectrum)\n    opt_ratio = np.sum(y_pred[t2:t3] * y_true[t2:t3], axis=0) / np.sum(\n        y_pred[t2:t3] ** 2, axis=0\n    )\n    y_pred[t2:t3] *= opt_ratio\n    y_p = np.concatenate([y_pred[:t1], y_pred[t2:t3], y_pred[t4:]])\n    y_t = np.concatenate([y_true[:t1], y_true[t2:t3], y_true[t4:]])\n\n    return ((y_p - y_t) / scale).ravel()\n\n\n# 誤差関数の定義\ndef loss_function2(params, t, w, t1, t2, t3, t4, d_t, d_w, spectrum, y_true):\n    y_pred = model(params, t, w, d_t, d_w, spectrum)\n\n    # :t1 and :t4\n    y_pred2 = np.concatenate([y_pred[:t1], y_pred[t4:]], axis=0)\n    y_true2 = np.concatenate([y_true[:t1], y_true[t4:]], axis=0)\n    opt_ratio2 = np.sum(y_pred2 * y_true2, axis=0) / np.sum(y_pred2**2, axis=0)\n    y_pred[:] *= opt_ratio2\n\n    opt_ratio = np.sum(y_pred[t2:t3] * y_true[t2:t3], axis=0) / np.sum(\n        y_pred[t2:t3] ** 2, axis=0\n    )\n    y_pred[t2:t3] *= opt_ratio\n    y_p = np.concatenate([y_pred[:t1], y_pred[t2:t3], y_pred[t4:]])\n    y_t = np.concatenate([y_true[:t1], y_true[t2:t3], y_true[t4:]])\n\n    return ((y_p - y_t) / scale).ravel()\n\n\ndef get_dip(params, t, w, t1, t2, t3, t4, d_t, d_w, spectrum, y_true):\n    y_pred = model(params, t, w, d_t, d_w, spectrum)\n\n    # :t1 and :t4\n    y_pred2 = np.concatenate([y_pred[:t1], y_pred[t4:]], axis=0)\n    y_true2 = np.concatenate([y_true[:t1], y_true[t4:]], axis=0)\n    opt_ratio2 = np.sum(y_pred2 * y_true2, axis=0) / np.sum(y_pred2**2, axis=0)\n    y_pred[:] *= opt_ratio2\n\n    opt_ratio = np.sum(y_pred[t2:t3] * y_true[t2:t3], axis=0) / np.sum(\n        y_pred[t2:t3] ** 2, axis=0\n    )\n    all_opt_ratio = np.sum(y_pred[t2:t3] * y_true[t2:t3]) / np.sum(y_pred[t2:t3] ** 2)\n\n    return 1 - opt_ratio, 1 - all_opt_ratio\n\n\ndef get_error(params, t, w, t1, t2, t3, t4, d_t, d_w, spectrum, y_true):\n    y_pred = model(params, t, w, d_t, d_w, spectrum)\n\n    # :t1 and :t4\n    y_pred2 = np.concatenate([y_pred[:t1], y_pred[t4:]], axis=0)\n    y_true2 = np.concatenate([y_true[:t1], y_true[t4:]], axis=0)\n    opt_ratio2 = np.sum(y_pred2 * y_true2, axis=0) / np.sum(y_pred2**2, axis=0)\n    y_pred[:] *= opt_ratio2\n\n    y_p = y_pred[t2:t3, :]\n    y_t = y_true[t2:t3, :]\n    num_boost = 1000\n    sample_index = np.random.randint(0, y_p.shape[0], size=(num_boost, y_p.shape[0]))\n    y_p_sample = np.take(y_p, sample_index, axis=0)\n    y_t_sample = np.take(y_t, sample_index, axis=0)\n    opt_ratio_list = np.sum(y_p_sample * y_t_sample, axis=1) / np.sum(\n        y_p_sample**2, axis=1\n    )\n    error = np.std(opt_ratio_list, axis=0)\n\n    return error\n\n\ndef fit(signal, d_t=4, d_w=4):\n    t1, t2, t3, t4 = phase_sec_detector(signal.mean(axis=1))\n\n    t = np.arange(signal.shape[0])\n    w = np.arange(signal.shape[1])\n\n    # 初期パラメータの設定\n    r = minimize(\n        objective,\n        [\n            1\n            - signal[t2:t3, :].mean()\n            / np.concatenate([signal[:t1, :], signal[:t4, :]]).mean()\n        ],\n        args=(signal.mean(axis=1), t1, t2, t3, t4),\n        method=\"Nelder-Mead\",\n    )\n    dip_mean = r.x[0]\n    y = (\n        (signal.mean(axis=1))[:t1].tolist()\n        + ((signal.mean(axis=1))[t2:t3] / (1 - dip_mean)).tolist()\n        + (signal.mean(axis=1))[t4:].tolist()\n    )\n    x = list(range(0, t1)) + list(range(t2, t3)) + list(range(t4, signal.shape[0]))\n    z_init = np.polyfit(x, y, deg=d_t) / signal.mean()\n    z_init[-1] -= 1\n    spectrum = signal.mean(axis=0)\n\n    # binでまとめて高速化\n    time_bin = 30\n    wave_bin = 1\n    signal_bin = signal[\n        : signal.shape[0] // time_bin * time_bin,\n        : signal.shape[1] // wave_bin * wave_bin,\n    ]\n    signal_bin = signal_bin.reshape(\n        signal.shape[0] // time_bin, time_bin, signal.shape[1] // wave_bin, wave_bin\n    ).mean(axis=(1, 3))\n    t_bin = (\n        t[: t.shape[0] // time_bin * time_bin]\n        .reshape(t.shape[0] // time_bin, time_bin)\n        .mean(axis=1)\n    )\n    w_bin = (\n        w[: w.shape[0] // wave_bin * wave_bin]\n        .reshape(w.shape[0] // wave_bin, wave_bin)\n        .mean(axis=1)\n    )\n    t1_bin = t1 // time_bin\n    t2_bin = t2 // time_bin\n    t3_bin = t3 // time_bin\n    t4_bin = t4 // time_bin\n    t1_bin -= 2\n    t2_bin += 2\n    t3_bin -= 2\n    t4_bin += 2\n    spectrum_bin = (\n        spectrum[: spectrum.shape[0] // wave_bin * wave_bin]\n        .reshape(spectrum.shape[0] // wave_bin, wave_bin)\n        .mean(axis=1)\n    )\n\n    initial_params = np.concatenate([z_init, np.repeat(0, d_w), np.array([1, 1])])\n\n    # 最適化の実行\n    result = least_squares(\n        loss_function,\n        initial_params,\n        args=(\n            t_bin,\n            w_bin,\n            t1_bin,\n            t2_bin,\n            t3_bin,\n            t4_bin,\n            d_t,\n            d_w,\n            spectrum_bin,\n            signal_bin,\n        ),\n        x_scale=\"jac\",\n    )\n\n    result = least_squares(\n        loss_function2,\n        result.x,\n        args=(\n            t_bin,\n            w_bin,\n            t1_bin,\n            t2_bin,\n            t3_bin,\n            t4_bin,\n            d_t,\n            d_w,\n            spectrum_bin,\n            signal_bin,\n        ),\n        x_scale=\"jac\",\n    )\n\n    # 最適化されたパラメータ\n    optimized_params = result.x\n    dip, dip_mean = get_dip(\n        optimized_params, t, w, t1, t2, t3, t4, d_t, d_w, spectrum, signal\n    )\n    err = get_error(optimized_params, t, w, t1, t2, t3, t4, d_t, d_w, spectrum, signal)\n\n    return optimized_params, dip, dip_mean, err","metadata":{"execution":{"iopub.status.busy":"2024-11-10T07:53:15.951306Z","iopub.execute_input":"2024-11-10T07:53:15.95173Z","iopub.status.idle":"2024-11-10T07:53:16.002491Z","shell.execute_reply.started":"2024-11-10T07:53:15.95169Z","shell.execute_reply":"2024-11-10T07:53:16.001264Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from sklearn.gaussian_process import GaussianProcessRegressor\nfrom sklearn.gaussian_process.kernels import RBF, ConstantKernel as C\nfrom sklearn.gaussian_process.kernels import Matern, ConstantKernel as C\nfrom sklearn.gaussian_process.kernels import RationalQuadratic, ConstantKernel as C\nimport warnings\nfrom sklearn.gaussian_process.kernels import ConvergenceWarning\n\n# ConvergenceWarning を無視する\nwarnings.filterwarnings(\"ignore\", category=ConvergenceWarning)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:55:38.302193Z","iopub.execute_input":"2024-10-31T17:55:38.303191Z","iopub.status.idle":"2024-10-31T17:55:38.316019Z","shell.execute_reply.started":"2024-10-31T17:55:38.303142Z","shell.execute_reply":"2024-10-31T17:55:38.314697Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"all_s = []\nall_sigma = []\ndip_list = []\nfor idx in tqdm(range(len(adc_info))):\n    d_t = 4\n    d_w = 4\n    scale = (signal_list[idx, 1:, :] - signal_list[idx, :-1, :]).std(axis=0)\n    optimized_params, dip, dip_mean, err = fit(signal_list[idx].copy(), d_t=d_t, d_w=d_w)\n    dip_list.append(dip[::-1])\n    \n    dip_ave = np.mean(dip)\n    \n    # ガウス過程回帰\n    data = signal_list[idx].copy()\n    x = np.arange(282)\n    y = dip.copy() - dip_ave\n    s = err * 1.6  # NEW\n    X = x.reshape(-1, 1)\n\n    kernel1 = C(y.max() - y.min(), (1e-9, 1e3)) * RBF(10, (1, 1e5))\n    kernel2 = C(y.max() - y.min(), (1e-9, 1e3)) * Matern(\n        length_scale=10, length_scale_bounds=(1, 1e5), nu=1.5\n    )\n    kernel = kernel1 + kernel2\n    \n    # ガウス過程回帰\n    gp = GaussianProcessRegressor(kernel=kernel, alpha=s**2, n_restarts_optimizer=10)\n    gp.fit(X, y)\n\n    # 新しい波長データに対する予測\n    x_pred = np.arange(283).reshape(-1, 1)  # 予測したい波長範囲\n    y_pred, y_std = gp.predict(x_pred, return_std=True)\n    y_pred = y_pred[::-1] + dip_ave\n    y_std = y_std[::-1]\n    all_s.append(y_pred)\n    \n    # 標準偏差推定\n    window_size = 20\n    \n    kernel = np.ones(window_size) / window_size\n    moving_average = np.convolve(dip, kernel, mode='valid')\n    # all_sigma.append(np.repeat(np.std(moving_average), 283))\n    \n    sigma = np.std(moving_average)\n    coef = 0.6\n    p = 1\n    sigma = (sigma**p + coef * y_std**p) ** (1 / p)\n    all_sigma.append(sigma)\n\nall_s = np.array(all_s)\n# all_sigma = np.array(all_sigma)*0.59\n\n# 追加\noffset_sigma = 7e-6\nall_sigma = np.array(all_sigma)*0.35 + offset_sigma","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:55:38.317544Z","iopub.execute_input":"2024-10-31T17:55:38.317945Z","iopub.status.idle":"2024-10-31T17:57:49.03109Z","shell.execute_reply.started":"2024-10-31T17:55:38.317905Z","shell.execute_reply":"2024-10-31T17:57:49.029132Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from keras.layers import Input, Dense\nfrom keras.models import Model\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint\nfrom sklearn.model_selection import KFold\nimport numpy as np\nimport pandas as pd\nfrom keras import regularizers\nfrom sklearn.decomposition import NMF\n\nif len(adc_info) > 11:\n    # データの前処理\n    new_list = []\n    for dip in dip_list:\n        window_size = 5\n        new_dip = pd.Series(dip.clip(0)).rolling(window=window_size, min_periods=1).median().to_numpy()\n        new_list.append(new_dip)\n\n    data = np.array(new_list)\n\n    # autoencoder\n    # 平均値を引いてデータを中心化\n    data_mean = np.mean(data, axis=1, keepdims=True)\n    data_centered = data - data_mean\n    data_std = np.std(data, axis=1, keepdims=True)\n    data_centered /= data_std\n\n    data_std2 = np.std(data[:,1:] - data[:,:-1], axis=0)\n    data_std2 = np.concatenate([data_std2, data_std2[-1:]])\n    data_centered /= data_std2\n\n    n_components = 4\n    input_dim = data_centered.shape[1]\n    encoding_dim = n_components  # 次元数を設定\n    \n    # Autoencoderの設計\n    input_data = Input(shape=(input_dim,))\n    encoded = Dense(encoding_dim, activation='relu')(input_data)\n    decoded = Dense(input_dim, activation='linear')(encoded)\n\n    autoencoder = Model(input_data, decoded)\n    autoencoder.compile(optimizer='adamw', loss='mse')\n\n    # EarlyStoppingとModelCheckpointコールバックを設定\n    early_stopping = EarlyStopping(monitor='val_loss', patience=5, restore_best_weights=True)\n    model_checkpoint = ModelCheckpoint('best_autoencoder.keras', save_best_only=True, monitor='val_loss')\n\n    # 学習\n    autoencoder.fit(data_centered, data_centered, epochs=500, batch_size=128, shuffle=True, verbose=0, \n                    validation_data=(data_centered, data_centered),\n                    callbacks=[early_stopping, model_checkpoint])\n\n    # 検証データの再構成\n    reconstructed_data_autoencoder = autoencoder.predict(data_centered) * data_std2 * data_std + data_mean\n    \n    # NMF\n    data = np.array(new_list)\n    # スケールを合わせる\n    data_std2 = np.std(data[:, 1:] - data[:, :-1], axis=0)\n    data_std2 = np.concatenate([data_std2, data_std2[-1:]])\n    data_centered = data/data_std2\n\n    # NMFのインスタンス作成\n    n_components = 5\n    nmf = NMF(n_components=n_components, init='random', random_state=0, max_iter=10000)\n\n    # NMFの適用\n    W = nmf.fit_transform(data_centered)\n    H = nmf.components_\n\n    # 再構成されたスペクトル\n    reconstructed_data_nmf = np.dot(W, H)*data_std2\n\n    all_s[:,0] = np.min(all_s[:,:283], axis=1)*0.6 + np.min(reconstructed_data_autoencoder[:,:282], axis=1)*0.2 + np.min(reconstructed_data_nmf[:,:282], axis=1)*0.2\n    all_s[:,1:283] = all_s[:,1:283]*0.6 + reconstructed_data_autoencoder[:,:282]*0.2 + reconstructed_data_nmf[:,:282]*0.2","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:57:49.033237Z","iopub.execute_input":"2024-10-31T17:57:49.033818Z","iopub.status.idle":"2024-10-31T17:57:49.075724Z","shell.execute_reply.started":"2024-10-31T17:57:49.033758Z","shell.execute_reply":"2024-10-31T17:57:49.073878Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import scipy\ndef score(\n        solution,\n        submission,\n        naive_mean = 0.0025517145902829823,\n        naive_sigma = 0.0017261973793536417,\n        sigma_true = 1e-5\n    ) -> float:\n    '''\n    This is a Gaussian Log Likelihood based metric. For a submission, which contains the predicted mean (x_hat) and variance (x_hat_std),\n    we calculate the Gaussian Log-likelihood (GLL) value to the provided ground truth (x). We treat each pair of x_hat,\n    x_hat_std as a 1D gaussian, meaning there will be 283 1D gaussian distributions, hence 283 values for each test spectrum,\n    the GLL value for one spectrum is the sum of all of them.\n\n    Inputs:\n        - solution: Ground Truth spectra (from test set)\n            - shape: (nsamples, n_wavelengths)\n        - submission: Predicted spectra and errors (from participants)\n            - shape: (nsamples, n_wavelengths*2)\n        naive_mean: (float) mean from the train set.\n        naive_sigma: (float) standard deviation from the train set.\n        sigma_true: (float) essentially sets the scale of the outputs.\n    '''\n\n    n_wavelengths = 283\n\n    y_pred = submission[:, :n_wavelengths]\n    # Set a non-zero minimum sigma pred to prevent division by zero errors.\n    sigma_pred = np.clip(submission[:, n_wavelengths:], a_min=10**-15, a_max=None)\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(scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true * np.ones_like(y_true)))\n    GLL_mean = np.sum(scipy.stats.norm.logpdf(y_true, loc=naive_mean * np.ones_like(y_true), scale=naive_sigma * np.ones_like(y_true)))\n\n    submit_score = (GLL_pred - GLL_mean)/(GLL_true - GLL_mean)\n    return float(np.clip(submit_score, 0.0, 1.0))","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:57:49.077786Z","iopub.execute_input":"2024-10-31T17:57:49.078269Z","iopub.status.idle":"2024-10-31T17:57:49.101727Z","shell.execute_reply.started":"2024-10-31T17:57:49.078213Z","shell.execute_reply":"2024-10-31T17:57:49.100294Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if dataset == \"train\":\n    label = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2024/train_labels.csv\")\n    label = label.to_numpy()[:,1:]\n    sub = np.concatenate([all_s, all_sigma], axis=1)\n    score_sigma_estimate = score(label[:len(adc_info)], sub)\n    print(score_sigma_estimate)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:57:49.108719Z","iopub.execute_input":"2024-10-31T17:57:49.110804Z","iopub.status.idle":"2024-10-31T17:57:49.195244Z","shell.execute_reply.started":"2024-10-31T17:57:49.110734Z","shell.execute_reply":"2024-10-31T17:57:49.19383Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\n\nif dataset == \"train\" and len(adc_info) <= 100:\n    # if dataset == \"train\":\n    # サブプロットの行数と列数を設定 (例: 4行5列)\n    n_rows = (len(adc_info) - 1) // 5 + 1\n    fig, axes = plt.subplots(\n        n_rows, 5, figsize=(20, n_rows * 4)\n    )  # figsizeで全体のサイズを調整\n\n    label = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2024/train_labels.csv\")\n    # label = pd.read_csv(f\"../../input/ariel-data-challenge-2024/train_labels.csv\")\n    label = label.to_numpy()[:, 1:]\n\n    # 20個のプロットを作成\n    for idx in range(len(adc_info)):\n        if len(adc_info) < 5:\n            ax = axes[idx]\n        else:\n            ax = axes[idx // 5, idx % 5]  # 4行5列のgrid内で位置を指定\n        dip = dip_list[idx]\n        pred = all_s[idx, 1:]\n        l = label[idx, 1:]\n        ax.scatter(np.arange(len(dip)), dip[::-1], s=3)\n        ax.plot(np.arange(len(pred)), pred, label=f\"taurex {idx}\", color=\"green\")\n        ax.plot(np.arange(len(l)), l, label=f\"label {idx}\", color=\"red\")\n        ax.legend()\n        score_idx = score(label[idx : idx + 1], sub[idx : idx + 1])\n        ax.set_title(f\"Score {score_idx}\")  # 各プロットにタイトルを設定\n\n    # 余白を調整して見やすく\n    plt.tight_layout()\n\n    # ファイルに保存\n    plt.savefig(\"combined_plots.png\")\n\n    # グラフを表示\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:58:46.393829Z","iopub.execute_input":"2024-10-31T17:58:46.394339Z","iopub.status.idle":"2024-10-31T17:58:52.342295Z","shell.execute_reply.started":"2024-10-31T17:58:46.394295Z","shell.execute_reply":"2024-10-31T17:58:52.340836Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"ss = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/sample_submission.csv')\nsigma = all_sigma\npred = all_s.clip(0) \nsubmission = pd.DataFrame(np.concatenate([pred,sigma], axis=1), columns=ss.columns[1:])\nsubmission.index = adc_info.index\nsubmission.to_csv('submission.csv')\nsubmission","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:57:53.99008Z","iopub.execute_input":"2024-10-31T17:57:53.99051Z","iopub.status.idle":"2024-10-31T17:57:54.07476Z","shell.execute_reply.started":"2024-10-31T17:57:53.990467Z","shell.execute_reply":"2024-10-31T17:57:54.073441Z"},"trusted":true},"outputs":[],"execution_count":null}]}