{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","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"}],"dockerImageVersionId":30747,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom pandas.plotting import autocorrelation_plot\nfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf\n\nfrom scipy.interpolate import griddata\nfrom tqdm import tqdm\nimport os\nimport shutil\nimport time\nimport warnings\nwarnings.filterwarnings('ignore', module='statsmodels')\n\nfrom statsmodels.tsa.stattools import adfuller, kpss\nfrom scipy.signal import butter, filtfilt\nfrom scipy.fft import fft, fftfreq\nimport pywt\n\nimport pandas.api.types\nimport scipy.stats\n\nimport glob\nimport joblib\nimport gc ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-08-06T17:26:58.337229Z","iopub.execute_input":"2024-08-06T17:26:58.337744Z","iopub.status.idle":"2024-08-06T17:27:16.137202Z","shell.execute_reply.started":"2024-08-06T17:26:58.337689Z","shell.execute_reply":"2024-08-06T17:27:16.135811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Task Understanding\nThe task of this competition is to extract the atmospheric spectra for each observations with an **estimate of their level of uncertainty**. (Performing this detrending process to extract ***atmospheric spectra*** and their ***associated errorbars*** from raw observational data is a crucial and common prerequisite step for any modern astronomical instrument before the data can undergo scientific analysis.)\n\n# Data Understanding\n\n## train_labels.csv\nGround truth spectra. (Ground truth spectroscopy refers to known spectral data for a particular object, such as a planet or a star. These data are usually derived from observations from high-precision instruments or validated model calculations.)","metadata":{}},{"cell_type":"code","source":"train_labels = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2024/train_labels.csv\").set_index(\"planet_id\")\nprint(train_labels.shape)\ntrain_labels.head(5)","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:27:16.13905Z","iopub.execute_input":"2024-08-06T17:27:16.140191Z","iopub.status.idle":"2024-08-06T17:27:16.291233Z","shell.execute_reply.started":"2024-08-06T17:27:16.140147Z","shell.execute_reply":"2024-08-06T17:27:16.289978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10, 4))\n\nplt.plot(train_labels.columns, train_labels.iloc[0], label = f'id_{train_labels.index[0]}')\nplt.plot(train_labels.columns, train_labels.iloc[1], label = f'id_{train_labels.index[1]}')\nplt.plot(train_labels.columns, train_labels.iloc[2], label = f'id_{train_labels.index[2]}')\nplt.plot(train_labels.columns, train_labels.iloc[3], label = f'id_{train_labels.index[3]}')\nplt.plot(train_labels.columns, train_labels.iloc[4], label = f'id_{train_labels.index[4]}')\n\nplt.title('Spectral Data of Samples')\nplt.xlabel('Wavelength')\nplt.ylabel('Value')\nplt.legend()\nplt.grid(True)\n\ntick_positions = range(0, len(train_labels.columns), 40)\ntick_labels = [train_labels.columns[i] for i in tick_positions]\nplt.xticks(ticks=tick_positions, labels=tick_labels, rotation=45)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-06T15:41:16.057583Z","iopub.execute_input":"2024-08-06T15:41:16.058254Z","iopub.status.idle":"2024-08-06T15:41:16.589091Z","shell.execute_reply.started":"2024-08-06T15:41:16.058188Z","shell.execute_reply":"2024-08-06T15:41:16.587453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that the spectrum of each planet fluctuates in a small range, which means that the spectrum can be used to distinguish each planet. Fluctuations are due to some atmospheric phenomenon, and planets with close spectral curves may be at similar distances from their stars. \n\n*So I think the mean of all the spectra can also be used as a simple marker to distinguish planets.*","metadata":{}},{"cell_type":"code","source":"correlation_matrix = train_labels.corr()\n\nplt.figure(figsize=(5, 4))\nsns.heatmap(correlation_matrix, annot=False, cmap='coolwarm', fmt='.2f')\nplt.title('Correlation Heatmap in Samples')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It appears that there are many planets whose ground truth spectra is similar.","metadata":{}},{"cell_type":"code","source":"# for i in range(train_labels.shape[1]):\n#     plt.figure(figsize=(10, 6))\n#     autocorrelation_plot(train_labels.iloc[:, i])\n#     plt.title(f'Autocorrelation Plot for Feature {i+1}')\n#     plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean_values = train_labels.mean(axis=1)\n\nprint(\"max_wl_mean\",np.max(mean_values), \"min_wl_mean\", np.min(mean_values))\n\n# mean_df = pd.DataFrame({\n#     'planet_id': data_subset.index,\n#     'mean_gts': mean_values\n# }).set_index('planet_id')\n\n# mean_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting the MSE distribution as a histogram\nplt.figure(figsize=(12, 2))\nplt.hist(mean_values, bins=20, color='skyblue', edgecolor='black')\nplt.title('Distribution of Mean')\nplt.xlabel('Mean')\nplt.ylabel('Frequency')\nplt.grid(True)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## wavelength.csv\nThe wavelength grid for each ground truth spectrum in the dataset.","metadata":{}},{"cell_type":"code","source":"wavelengths = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2024/wavelengths.csv\")\nwavelengths","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:27:16.292946Z","iopub.execute_input":"2024-08-06T17:27:16.293471Z","iopub.status.idle":"2024-08-06T17:27:16.335232Z","shell.execute_reply.started":"2024-08-06T17:27:16.293428Z","shell.execute_reply":"2024-08-06T17:27:16.333914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Comparing the ground truth spectrum with the wavelength grid does not make sense per se. What really matters is:\n* The ground truth spectra of different planets are compared under a unified wavelength grid to understand their spectral characteristics.\n* The observed spectra and ground truth spectra of the same planet were compared under a unified wavelength grid to verify the accuracy and validity of the observed data.\n\n## [train/test]_adc_info.csv\nContains analog-to-digital (ADC) conversion parameters (gain and offset) for restoring the original dynamic range of the data. Also includes a star column identifying which star was used for that planet's simulation. \n\n*To convert analog data into digital signals, ADCs apply specific gain and offset parameters. These parameters can restore the data to its original dynamic range, which is very important to ensure the accuracy and consistency of the data:*\n* **Gain**: Adjusts the amplitude of a signal, usually by a multiplicative factor.\n* **Offset**: The base point at which the signal is adjusted, usually by an additive factor.\n\n*The original observations are corrected using gains and offsets to recover their original dynamic range:*\n$Corrected Value=(Raw Value−Offset)\\times Gain$","metadata":{}},{"cell_type":"code","source":"train_adc_info = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv\").set_index(\"planet_id\")\nprint(train_adc_info.shape)\ntrain_adc_info.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10, 4))\nsns.countplot(data=train_adc_info, x='star', order=train_adc_info['star'].value_counts().index)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* **star**: The column is a binary distribution of 0 and 1, meaning that only two stars were used in the simulation to generate the data. 0 and 1 are the identifiers of two different stars.\n\nOther details:\n1. **FGS1** refers to data from the First Fine Guidance Sensor (FGS). FGS are commonly used to precisely locate and track stars or celestial bodies to ensure a stable and accurate pointing of the telescope. FGS1 is one of the channels responsible for collecting a specific band or type of optical data.\n1. **AIRS-CH0** refers to Channel 0 of the Atmospheric Infrared Sounder (AIRS). AIRS is a hyperspectral resolution infrared sounder commonly used to measure atmospheric temperature, humidity, cloud cover, and other meteorological parameters. AIRS-CH0 is one of its multiple channels, collecting a specific band or type of infrared data.","metadata":{}},{"cell_type":"code","source":"# data_subset = train_adc_info.head(5).drop([\"star\"],axis=1)\n# data_subset = pd.concat([data_subset, mean_df], axis=1)\n\n# fig, axes = plt.subplots(nrows=1, ncols=data_subset.shape[1], figsize=(14, 4))\n# fig.subplots_adjust(wspace=0.5)\n\n# for idx, col in enumerate(data_subset.columns):\n#     sns.barplot(x=data_subset.index, y=data_subset[col], ax=axes[idx])\n#     axes[idx].set_title({col})\n#     axes[idx].set_xlabel('planet_id')\n#     axes[idx].set_ylabel(col)\n#     axes[idx].grid(True, linestyle='--', linewidth=0.5)\n#     axes[idx].set_xticklabels(axes[idx].get_xticklabels(), rotation=90)\n\n# plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sorted_data = pd.DataFrame()\n\n# for col in data_subset.columns:\n#     sorted_series = data_subset[[col]].copy()\n#     sorted_series['abs_value'] = sorted_series[col].abs()\n#     sorted_series = sorted_series.sort_values(by='abs_value', ascending=False)\n#     sorted_series = sorted_series.drop(columns='abs_value')\n#     sorted_data[col] = sorted_series.index\n\n# sorted_data","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## axis_info.parquet\nAxis information for both instruments.\n\n* **AIRS -CH0-AXis0-H:** Represents the value of axis0 of the AIRS-CH0 instrument at a point in time in units of h (usually a physical quantity related to the instrument, such as height or other measurement).\n* **AIRS -CH0-AXis2-UM**: Represents the value of axis2 of the AIRS-CH0 instrument at a point in time, in μm (micrometers, usually representing the wavelength of the spectrum).\n* **AIRS -CH0-Integration_time**: Indicates the integration time of the AIRS-CH0 instrument at a certain point in time, that is, the length of time for each measurement.\n* **Fgs1-axis0-h**: Represents the value of axis0 of the FGS1 instrument at a point in time in units of h.\n\n*--Idea--*\n* Data registration: This axis information can be used to register image data with actual physical measurements, such as correlating spectral data with wavelength, time, and other information.\n* Data analysis: Knowing the measurement conditions at different points in time, such as integration time, can help analyze data quality and instrument performance.\n* Correction: This information may be used to correct image data, especially when NaN values appear in some columns of the data, which can be completed or adjusted using information from other columns.","metadata":{}},{"cell_type":"code","source":"axis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\nprint(axis_info.shape)\naxis_info.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(10, 6))\n\ncolumns = ['AIRS-CH0-axis0-h', 'AIRS-CH0-axis2-um', 'AIRS-CH0-integration_time', 'FGS1-axis0-h']\nfor ax, col in zip(axs.flat, columns):\n    ax.scatter(axis_info.index, axis_info[col], s=1)  \n    ax.set_title(col)\n    #ax.set_xlabel('Index')\n    #ax.set_ylabel(col)\n    ax.grid(True)\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## [train/test]/[planet_id]/FGS1_signal.parquet\nSignal data from the FGS1 instrument. Each file contains 135,000 rows of images at 0.1 second time steps. Each 32x32 image has been flattened into 1024 columns. You can un-flatten the data with numpy.reshape(135000, 32, 32). Similar to AIR-CH0, the data is generated in uint16. To restore its original dynamic range you must multiply the data by the matching **gain** value from **[train/test]_adc_info.csv** and then add the **offset** value, also from **[train/test]_adc_info.csv**.","metadata":{}},{"cell_type":"code","source":"train_785834_FGS = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/785834/FGS1_signal.parquet')\nprint(train_785834_FGS.shape)\ntrain_785834_FGS.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"signal_array = train_785834_FGS.values\n\nreshaped_signal_data = signal_array.reshape(135000, 32, 32)\n\ndata = reshaped_signal_data[0]\nplt.matshow(data, cmap=\"coolwarm\")\nplt.colorbar(shrink=0.8, aspect=20, pad=0.1)\nplt.title('Restored Signal Image at 0th Time Step')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#(673, 5)\n#FGS1_adc_offset\tFGS1_adc_gain\tAIRS-CH0_adc_offset\tAIRS-CH0_adc_gain\t\n#star planet_id\t\n\nplanet_id = 785834\n\ngain = train_adc_info.loc[planet_id, 'FGS1_adc_gain']\noffset = train_adc_info.loc[planet_id, 'FGS1_adc_offset']\n\noriginal_signal_data = reshaped_signal_data * gain + offset\n\nprint(original_signal_data.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = reshaped_signal_data[0]\nplt.matshow(data, cmap=\"coolwarm\")\nplt.colorbar(shrink=0.8, aspect=20, pad=0.1)\nplt.title('Restored Original Signal Image at 0th Time Step')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"max_value = np.max(data)\nmax_position = np.unravel_index(np.argmax(data), data.shape)\nprint(max_value, max_position)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slice_index = 0  # Select the first time slice\n\n# Select the slice\nZ = reshaped_signal_data[slice_index, :, :]\n\n# Creating a grid\nX = np.arange(Z.shape[1])  # length\nY = np.arange(Z.shape[0])  # value\nX, Y = np.meshgrid(X, Y)  # Creating a grid\n\n# Draw a 3D thermal map\nfig = plt.figure(figsize=(12, 8), dpi=80)\nax = fig.add_subplot(111, projection='3d')\nsurf = ax.plot_surface(X, Y, Z, cmap='coolwarm')\n\n# Add a colorbar\ncbar = fig.colorbar(surf, ax=ax, shrink=0.8, aspect=20, pad=0.1)\n\nax.set_xlabel('Length')\nax.set_ylabel('Value')\nax.set_zlabel('Signal')\nax.set_title(f'Slice at Time Index {slice_index}')\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slice_index = 5000  # Select the 5000th time slice\n\n# Select the slice\nZ = reshaped_signal_data[slice_index, :, :]\n\n# Creating a grid\nX = np.arange(Z.shape[1])  # length\nY = np.arange(Z.shape[0])  # value\nX, Y = np.meshgrid(X, Y)  # Creating a grid\n\n# Draw a 3D thermal map\nfig = plt.figure(figsize=(12, 8), dpi=80)\nax = fig.add_subplot(111, projection='3d')\nsurf = ax.plot_surface(X, Y, Z, cmap='coolwarm')\n\n# Add a colorbar\ncbar = fig.colorbar(surf, ax=ax, shrink=0.8, aspect=20, pad=0.1)\n\nax.set_xlabel('Length')\nax.set_ylabel('Value')\nax.set_zlabel('Signal')\nax.set_title(f'Slice at Time Index {slice_index}')\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From the 3D figure, we can see that the FGS1 signal is tapered in each time slice. **Assuming that all planets are distributed in this way (to be verified later)**, we extract peaks as features of the training data.","metadata":{}},{"cell_type":"code","source":"time = reshaped_signal_data[0]\nvalue = reshaped_signal_data[1]\nlength = reshaped_signal_data[2]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"length_value_matrix_avg = np.mean(reshaped_signal_data, axis=0)\n\nplt.matshow(length_value_matrix_avg, cmap='coolwarm', origin='lower', fignum=1)\nplt.xlabel('Length')\nplt.ylabel('Value')\nplt.title('Length vs Value')\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"length_time_matrix_avg = np.mean(reshaped_signal_data, axis=1)\n\nplt.matshow(length_time_matrix_avg.T, cmap='coolwarm', origin='lower', aspect='auto')\nplt.xlabel('Time')\nplt.ylabel('Length')\nplt.title('Time vs Length')\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"value_time_matrix_avg = np.mean(reshaped_signal_data, axis=2)\nplt.matshow(value_time_matrix_avg.T, cmap='coolwarm', origin='lower', aspect='auto')\nplt.xlabel('Time')\nplt.ylabel('Value')\nplt.title('Time vs Value')\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"max_values = np.max(reshaped_signal_data, axis=(1, 2))\n\nprint(max_values, \"shape:\", max_values.shape)\nplt.figure(figsize=(12, 2))\nplt.plot(max_values)\nplt.xlabel(\"Time\")\nplt.ylabel('Maximum Value')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"haha, that's a lot like a bunch of white noise. To deal with these timing signals later, I first prepare a stationarity test.\n\n* **ADF test**: Used to test whether the time series is stationary. If the p-value is less than the significance level alpha, then the time series is considered stationary.\n* **KPSS test**: Used to test whether the time series is stationary. If the p-value is greater than the significance level alpha, then the time series is considered stationary.\n* **Stationarity judgment**: If the results of both the ADF test and the KPSS test indicate that the time series is stationary, then the comprehensive judgment is that the time series is stationary.","metadata":{}},{"cell_type":"code","source":"def check_stationarity(data, alpha=0.05):\n    \"\"\"\n    Check the stationarity of a time series.\n    \n    Parameters:\n    data: pandas Series or numpy array, the time series data.\n    alpha: significance level, default is 0.05.\n    \n    Returns:\n    dict: A dictionary containing the results of ADF and KPSS tests and stationarity judgement.\n    \"\"\"\n    \n    # ADF test\n    adf_result = adfuller(data)\n    adf_statistic = adf_result[0]\n    adf_p_value = adf_result[1]\n    adf_critical_values = adf_result[4]\n    \n    adf_result_str = (\n        f'ADF Statistic: {adf_statistic}\\n'\n        f'p-value: {adf_p_value}\\n'\n        f'Critical Values: {adf_critical_values}\\n'\n        f'ADF Test: {\"Stationary\" if adf_p_value < alpha else \"Non-stationary\"}'\n    )\n    \n    # KPSS test\n    with warnings.catch_warnings(record=True) as w:\n        warnings.simplefilter('always')\n        kpss_result = kpss(data, regression='c')\n        \n        kpss_statistic = kpss_result[0]\n        kpss_p_value = kpss_result[1]\n        kpss_critical_values = kpss_result[3]\n        \n        # Check for the specific warning\n        if any('InterpolationWarning' in str(warning.message) for warning in w):\n            kpss_result_str = (\n                f'KPSS Statistic: {kpss_statistic}\\n'\n                f'p-value: {kpss_p_value} (Warning: actual p-value is smaller than this)\\n'\n                f'Critical Values: {kpss_critical_values}\\n'\n                f'KPSS Test: {\"Stationary\" if kpss_p_value > alpha else \"Non-stationary\"}'\n            )\n        else:\n            kpss_result_str = (\n                f'KPSS Statistic: {kpss_statistic}\\n'\n                f'p-value: {kpss_p_value}\\n'\n                f'Critical Values: {kpss_critical_values}\\n'\n                f'KPSS Test: {\"Stationary\" if kpss_p_value > alpha else \"Non-stationary\"}'\n            )\n    \n    # Stationarity judgement\n    is_stationary = (\n        adf_p_value < alpha and kpss_p_value > alpha\n    )\n    \n    return {\n        'ADF Result': adf_result_str,\n        'KPSS Result': kpss_result_str,\n        'Overall Stationarity': 'Stationary' if is_stationary else 'Non-stationary'\n    }","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check stationarity\nresults = check_stationarity(max_values)\nfor key, value in results.items():\n    print(f'{key}:\\n{value}\\n')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Although the ADF test indicates that the time series is stationary, the KPSS test indicates that the time series is non-stationary. Since the KPSS test is particularly **sensitive to low-frequency trends**, it is more likely to reject the stationarity hypothesis in the presence of a trend or nonstationarity component. Let's take a look at the signal peaks at odd and even times separately.","metadata":{}},{"cell_type":"code","source":"max_values_odd = max_values[1::2]\nmax_values_even = max_values[0::2]\n\nplt.figure(figsize=(12, 4))\n\nplt.plot(range(1, len(max_values), 2), max_values_odd, label='Odd Time Indices', color='blue')\n\nplt.plot(range(0, len(max_values), 2), max_values_even, label='Even Time Indices', color='green')\n\nplt.title('Signal Peaks at Odd and Even Time Indices')\nplt.xlabel('Time')\nplt.ylabel('Maximum Value')\n\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check stationarity\nresults = check_stationarity(max_values_odd)\nfor key, value in results.items():\n    print(f'{key}:\\n{value}\\n')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check stationarity\nresults = check_stationarity(max_values_even)\nfor key, value in results.items():\n    print(f'{key}:\\n{value}\\n')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It's still like white noise, I don't know how to extract features. I will try **Filtering**. Because we have selected the peak value, I want to remove the high frequency noise and keep the low frequency component, that is, use the **low frequency filter**.","metadata":{}},{"cell_type":"code","source":"def butter_lowpass(cutoff, fs, order=5):\n    nyquist = 0.5 * fs\n    normal_cutoff = cutoff / nyquist\n    b, a = butter(order, normal_cutoff, btype='low', analog=False)\n    return b, a\n\ndef lowpass_filter(data, cutoff, fs, order=5):\n    b, a = butter_lowpass(cutoff, fs, order=order)\n    y = filtfilt(b, a, data)\n    return y","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:27:54.337366Z","iopub.execute_input":"2024-08-06T17:27:54.338153Z","iopub.status.idle":"2024-08-06T17:27:54.344611Z","shell.execute_reply.started":"2024-08-06T17:27:54.338116Z","shell.execute_reply":"2024-08-06T17:27:54.343229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fs = 1000.0  # Frequency of sampling\ncutoff = 1.0  # Frequency of cutoff\nfiltered_data_odd = lowpass_filter(max_values_odd, cutoff, fs)\nfiltered_data_even = lowpass_filter(max_values_even, cutoff, fs)\nfiltered_data_all = lowpass_filter(max_values, cutoff, fs)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 2))\nplt.plot(filtered_data_odd, label='Odd Time Indices', color='blue')\nplt.xlabel(\"Time\")\nplt.ylabel('Maximum Value')\nplt.legend()\nplt.show()\n\nplt.figure(figsize=(12, 2))\nplt.plot(filtered_data_even, label='Odd Time Indices', color='green')\nplt.xlabel(\"Time\")\nplt.ylabel('Maximum Value')\nplt.legend()\nplt.show()\n\nplt.figure(figsize=(12, 2))\nplt.plot(filtered_data_all, label='All Time Indices', color='red')\nplt.xlabel(\"Time\")\nplt.ylabel('Maximum Value')\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check stationarity\nresults = check_stationarity(filtered_data_all)\nfor key, value in results.items():\n    print(f'{key}:\\n{value}\\n')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check stationarity\nresults = check_stationarity(filtered_data_all[40000:])\nfor key, value in results.items():\n    print(f'{key}:\\n{value}\\n')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"haha, the fluctuation before 40000 steps is bound to make the data difficult to smooth, even after we delete this part still cannot pass KPSS. However, it is worth noting that the **follow-up work should look at more data to see if there is a similar shake before 40000 steps**. If it is widespread, it may be due to systematic errors caused by some climatic and equipment problems. Things seem to be getting better. Let me try the **Fourier Transform** and **Wavelet Transform**.\n\n**Fourier Transform**\n* Application: Fourier transform is suitable for analyzing stationary signals with fixed frequency components. It converts time domain signals into frequency domain signals to provide spectral information of the whole time series.\n* Limitations: For non-stationary signals, the Fourier transform is less effective because it fails to provide the local characteristics of the signal in time. Even if the frequency components in the signal change over time, the Fourier transform cannot detect these changes.\n\n**Wavelet Transform (Wavelet Transform)**\n* Application: Wavelet transform is suitable for analyzing signals with local non-stationary characteristics. It can provide multi-resolution analysis of the signal in the time-frequency domain and is able to capture the signal variation in time and frequency simultaneously.\n* Advantages: Wavelet transform is particularly suitable for dealing with non-stationary signals because it provides information on time-frequency localization and can detect transient features and changes in the signal.\n\nPassing the ADF test but not the KPSS test indicates that your time series may contain some complex, time-varying characteristics or trends. The Fourier transform may not be sufficient to reveal these characteristics over time, since it assumes that the signal is stationary, and it may be more appropriate to use the wavelet transform to analyze your time series. **Different levels can be tried in the future and the best level can be decided by analyzing the validity of the results.**","metadata":{}},{"cell_type":"code","source":"def perform_wavelet_transform(data, wavelet='db1'):\n    \"\"\"\n    Perform wavelet transform on the given data and plot the results.\n\n    Parameters:\n    - data: The input time series data (1D array or list).\n    - wavelet: The type of wavelet to use (default is 'db1' for Daubechies wavelet).\n    \"\"\"\n    # Perform wavelet transform\n    coeffs = pywt.wavedec(data, wavelet)\n    \n    # Plot the wavelet transform results\n    plt.figure(figsize=(12, 18))\n    for i, coeff in enumerate(coeffs):\n        plt.subplot(len(coeffs), 1, i + 1)\n        plt.plot(coeff)\n        plt.title(f'Wavelet Coefficient Level {i}')\n    \n    plt.tight_layout()\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perform_wavelet_transform(max_values)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perform_wavelet_transform(filtered_data_all)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here I want to try the wavelet coeddicient level 6 for original data first.","metadata":{}},{"cell_type":"code","source":"def get_wavelet_coefficients_at_level(data, wavelet='db1', level=8):\n    \"\"\"\n    Perform wavelet transform on the given data and return the coefficients for the specified level.\n\n    Parameters:\n    - data: The input time series data (1D array or list).\n    - wavelet: The type of wavelet to use (default is 'db1' for Daubechies wavelet).\n    - level: The level of wavelet decomposition (default is 6).\n\n    Returns:\n    - level_coeffs: The coefficients at the specified level.\n    \"\"\"\n    # Perform wavelet transform\n    coeffs = pywt.wavedec(data, wavelet, level)\n    \n    # The coefficients are returned in order: [cA_n, cD_n, cD_(n-1), ..., cD_1]\n    # where cA_n is the approximation coefficients at level n\n    # and cD_i are the detail coefficients at level i.\n    \n    # The approximation coefficients at the highest level\n    # The detail coefficients for each level are at indices 1 through level\n    level_coeffs = coeffs[-1]  # Coefficients at the specified level\n\n    return level_coeffs","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"w_data = get_wavelet_coefficients_at_level(max_values, wavelet='db1', level=8)\nwf_data = get_wavelet_coefficients_at_level(filtered_data_all, wavelet='db1', level=8)\n\nplt.figure(figsize=(12, 2))\n\nplt.plot(w_data, label='Wave', color='grey')\nplt.plot(filtered_data_all, label='Filter', color='red')\nplt.plot(wf_data, label='Wave+Filter', color='orange')\n\nplt.xlabel(\"Time\")\nplt.ylabel('Maximum Value')\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The wavelet transform is not so obvious at the moment, so I'll just use the filter data.","metadata":{}},{"cell_type":"code","source":"# # Calculate the Fourier transform\n# fs = 1000.0  # Sampling frequency (assuming 100 Hz)\n# T = 100.0 / fs  # Period of sampling\n# n = len(max_values_odd)  # Length of signal\n# yf = fft(max_values_odd)\n# xf = fftfreq(n, T)[:n//2]\n\n# plt.figure(figsize=(12, 6))\n# plt.plot(xf, 2.0/n * np.abs(yf[:n//2]))\n# plt.title('Fourier Transform of max_values_odd')\n# plt.xlabel('Frequency (Hz)')\n# plt.ylabel('Amplitude')\n# plt.grid()\n# plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## [train/test]/[planet_id]/AIRS-CH0_signal.parquet\nSignal data from the AIRS-CH0 instrument. Each file contains 11,250 rows of images captured at constant time steps noted in axis_info.parquet file for details of the time steps. Each 32 x 356 image has been flattened into 11392 columns. You can un-flatten the data with numpy.reshape(11250, 32, 356). The instruments generate data as uint16. To restore the full dynamic range you must multiply the data by the matching gain value from","metadata":{}},{"cell_type":"code","source":"train_785834_AIRS = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/train/785834/AIRS-CH0_signal.parquet')\nprint(train_785834_AIRS.shape)\ntrain_785834_AIRS.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"signal_array = train_785834_AIRS.values\n\nreshaped_signal_data = signal_array.reshape(11250, 32, 356)\n\ndata = reshaped_signal_data[0]\nplt.matshow(data, cmap=\"coolwarm\")\nplt.colorbar()\nplt.title('Restored Signal Image at 0th Time Step')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#(673, 5)\n#FGS1_adc_offset\tFGS1_adc_gain\tAIRS-CH0_adc_offset\tAIRS-CH0_adc_gain\t\n#star planet_id\t\n\nplanet_id = 785834\n\ngain = train_adc_info.loc[planet_id, 'AIRS-CH0_adc_gain']\noffset = train_adc_info.loc[planet_id, 'AIRS-CH0_adc_offset']\n\noriginal_signal_data = reshaped_signal_data * gain + offset\n\nprint(original_signal_data.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = reshaped_signal_data[0]\nplt.matshow(data, cmap=\"coolwarm\")\nplt.colorbar()\nplt.title('Restored Original Signal Image at 0th Time Step')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slice_index = 0  # Select the first time slice\n\n# Select the slice\nZ = reshaped_signal_data[slice_index, :, :]\n\n# Creating a grid\nX = np.arange(Z.shape[1])  # length \nY = np.arange(Z.shape[0])  # value \nX, Y = np.meshgrid(X, Y)  \n\n# Draw a 3D thermal map\nfig = plt.figure(figsize=(12, 8), dpi=80)\nax = fig.add_subplot(111, projection='3d')\nsurf = ax.plot_surface(X, Y, Z, cmap='coolwarm')\n\n# Add a colorbar\ncbar = fig.colorbar(surf, ax=ax, shrink=0.8, aspect=20, pad=0.1)\n\nax.set_xlabel('Length')\nax.set_ylabel('Value')\nax.set_zlabel('Signal')\nax.set_title(f'Slice at Time Index {slice_index}')\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"time = reshaped_signal_data[0]\nvalue = reshaped_signal_data[1]\nlength = reshaped_signal_data[2]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"length_value_matrix_avg = np.mean(reshaped_signal_data, axis=0)\n\nplt.matshow(length_value_matrix_avg, cmap='coolwarm', origin='lower', fignum=1)\nplt.xlabel('Length')\nplt.ylabel('Value')\nplt.title('Length vs Value')\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"length_time_matrix_avg = np.mean(reshaped_signal_data, axis=1)\n\nplt.matshow(length_time_matrix_avg.T, cmap='coolwarm', origin='lower', aspect='auto')\nplt.xlabel('Time')\nplt.ylabel('Length')\nplt.title('Time vs Length')\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"value_time_matrix_avg = np.mean(reshaped_signal_data, axis=2)\nplt.matshow(value_time_matrix_avg.T, cmap='coolwarm', origin='lower', aspect='auto')\nplt.xlabel('Time')\nplt.ylabel('Value')\nplt.title('Time vs Value')\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## [train/test]/[planet_id]/[AIRS-CH0/FGS1]_calibration/dark.parquet\nDark frames are exposures taken with the shutter closed, capturing the thermal noise and bias level of the sensor. These are used to subtract the dark current from science images.\n\n*--Details--*\n* **Thermal noise correction**: The dark frame records the thermal noise of the sensor in the absence of light. These noises are generated by the electronic components of the sensor when exposed to ambient temperature. Thermal noise in the scientific image needs to be subtracted from the observed image to improve the accuracy of the data.\n* **Bias correction**: The dark frame also captures the sensor's bias level (that is, the signal baseline under zero light conditions). In actual observations, this bias signal will have an impact on the image, so it needs to be subtracted from the scientific image to eliminate the effect of baseline drift on the data.\n* **Dark frame correction**: Using read dark frame data to correct scientific images. This typically involves subtracting dark frame data from scientific images to remove thermal noise and biased signals.","metadata":{}},{"cell_type":"code","source":"train_785834_dark = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/train/785834/AIRS-CH0_calibration/dark.parquet')\nprint(train_785834_dark.shape)\ntrain_785834_dark.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = train_785834_dark.values  \nplt.matshow(data, cmap=\"coolwarm\")\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## [train/test]/[planet_id]/[AIRS-CH0/FGS1]_calibration/dead.parquet\nIdentifies dead or hot pixels on the sensor. Dead pixels do not respond to light, while hot pixels consistently produce high signal levels regardless of incoming light.\n\n*--Details--*\n* **Dead Pixels**: are pixels that do not respond to light. In the image, these pixels usually appear as constant values, which may be zero or a fixed value. They cannot capture any actual optical signal when collecting data, thus affecting the quality of the image.\n* **Hot Pixels**: produce unusually high signal levels, even in the absence of light. These pixels are usually caused by electronic noise or malfunction of the sensor. Thermal pixels appear as bright spots in images and may interfere with data analysis and scientific measurements.\n* **Image correction**: In the actual image processing process, you can use this data to correct scientific images. For example, for dead pixels, the value can be replaced by the average of neighboring pixels; For hot pixels, an appropriate algorithm can be used to reduce their impact.","metadata":{}},{"cell_type":"code","source":"train_785834_dead = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/train/785834/AIRS-CH0_calibration/dead.parquet')\nprint(train_785834_dead.shape)\ntrain_785834_dead.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = train_785834_dead.values  \nplt.matshow(data, cmap=\"binary\")\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## [train/test]/[planet_id]/[AIRS-CH0/FGS1]_calibration/flat.parquet\nFlat field frames are created by imaging a uniformly illuminated surface. They are used to correct for variations in pixel-to-pixel sensitivity and optical system irregularities.\n\n--Details--\n* Uneven image brightness caused by sensors or optical components.\n* **The flat-field frame**: records the relative sensitivity of each pixel and the irregularity of the optical system.\n* **Image Correction**: Each pixel value of the scientific image was divided by the normalized flat field frame value at the corresponding position.","metadata":{}},{"cell_type":"code","source":"train_785834_flat = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/train/785834/AIRS-CH0_calibration/flat.parquet')\nprint(train_785834_flat.shape)\ntrain_785834_flat.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = train_785834_flat.values  \nplt.matshow(data, cmap=\"coolwarm\")\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##[train/test]/[planet_id]/[AIRS-CH0/FGS1]_calibration/linear_corr.parquet\nInformation about the linearity correction of the sensor. The response of the pixels in the detector becomes less linear as they fill with electrons, approaching the point of saturation, where the pixel can no longer collect additional electrons and its response to light becomes flat. For an accurate estimate of the signal, the instrument's response as a function of the received charge is calibrated, and the correction is calculated using a polynomial of degree n. This polynomial allows for the conversion of the number of electrons collected/measured by the pixel into the number of electrons that the detector would have generated with a linear response.\n\n*--Details--*\n* $P(x)=a^0+a^1x+a^2x^2+…+a^nx^n$ where $x$ is the number of electrons measured by the pixel, $P(x)$ is the corrected number of electrons.\n* **Image Correction**: For each pixel in the scientific image, the nonlinear response is converted to a linear response using the calibrated polynomial parameters corrected.","metadata":{}},{"cell_type":"code","source":"train_785834_linear_corr = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/train/785834/AIRS-CH0_calibration/linear_corr.parquet')\nprint(train_785834_linear_corr.shape)\ntrain_785834_linear_corr.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = train_785834_linear_corr.values  \nplt.matshow(data, cmap=\"coolwarm\")\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## [train/test]/[planet_id]/[AIRS-CH0/FGS1]_calibration/read.parquet\nRead noise frames capture the electronic noise introduced during the readout process of the sensor. This noise is present even when no light falls on the detector.\n\n*--Details--*\n* **Objective**: By reading noisy frames, the noise characteristics of sensors can be understood and corrected in actual scientific images.\n* **Methods**: In scientific image processing, the noise signal in the read noise frame can be subtracted from the actual image, so as to improve the signal-to-noise ratio (SNR) of the image signal.","metadata":{}},{"cell_type":"code","source":"train_785834_read = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/train/785834/AIRS-CH0_calibration/read.parquet')\nprint(train_785834_read.shape)\ntrain_785834_read.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = train_785834_read.values  \nplt.matshow(data, cmap=\"coolwarm\")\nplt.colorbar()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data processing: extract the signal sequence from the two-dimensional image\n\nOnly do for **FGS1_signals** as a try. We must have missed a lot of valid information here, and the two-dimensional information is simply extracted into a list of features. Don't worry, this is just the beginning and experiment:\n\n1.  ***Data directory structure and file reading***:\nThe data is stored in the /kaggle/input/ariel-data-challenge-2024/train folder.\nThe name of each subfolder is the ID of the planet, and the subfolder contains a file named FGS1_signal.parquet, which stores the signal data for the planet.\n1.  ***Deformation of data***:\nRead the FGS1_signal.parquet file and load the data as signal_array.\nDeform signal_array to rearrange it into a 3D array with the shape (135000, 32, 32) reshaped_signal_data. Here, 135000 is the time step and 32x32 is the two-dimensional spatial dimension of the signal.\n1. ***Maximum value extraction***:\nThe maximum value of each time step is extracted from reshaped_signal_data to obtain a one-dimensional array max_values with the shape (135000,). This operation is to calculate the maximum value on the (1, 2) axis, that is, to calculate the maximum value on the spatial dimension.\n1. ***Filtering of noise***:\nThe max_values are filtered for noise using a low-pass filter. The design parameters of low-pass filter include cutoff frequency, sampling frequency fs and filter order.\nDefine the lowpass_filter function lowpass_filter(data, cutoff, fs, order=5) and enter max_values into the function to obtain the filtered data filtered_data_all.\n1. ***Store filtered data***:\nCreate a new DataFrame output_df indexed by the ID of the planet, with each row corresponding to the filtered_data_all of the planet.\nIterate through the folders of all planets, perform the above steps, and store the filtered_data_all for each planet into output_df.\n1. ***Save the result***:\nSave the final DataFrame output_df to the /kaggle/working/ folder for later use and analysis.","metadata":{}},{"cell_type":"code","source":"# # try  \n####################### train_labels = train_labels.head(5)       ","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:29:27.875837Z","iopub.execute_input":"2024-08-06T17:29:27.876276Z","iopub.status.idle":"2024-08-06T17:29:27.88179Z","shell.execute_reply.started":"2024-08-06T17:29:27.876239Z","shell.execute_reply":"2024-08-06T17:29:27.88046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Define file paths and parameters\ninput_dir = '/kaggle/input/ariel-data-challenge-2024/train'\n\n# Load the train_labels to get planet_ids\nplanet_ids = train_labels.index  ## set index\nnum_columns = 135000\nprint(\"Build:\")\noutput_df = pd.DataFrame(index=planet_ids, columns=range(num_columns))\n\n# Parameters for low-pass filter\ncutoff = 1.0\nfs = 1000.0","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:29:37.105628Z","iopub.execute_input":"2024-08-06T17:29:37.106054Z","iopub.status.idle":"2024-08-06T17:29:42.389588Z","shell.execute_reply.started":"2024-08-06T17:29:37.106017Z","shell.execute_reply":"2024-08-06T17:29:42.388097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"####################### output_df = output_df.head(5)     \n####################### print(output_df.shape)        ","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:29:42.392084Z","iopub.execute_input":"2024-08-06T17:29:42.392636Z","iopub.status.idle":"2024-08-06T17:29:42.400298Z","shell.execute_reply.started":"2024-08-06T17:29:42.392591Z","shell.execute_reply":"2024-08-06T17:29:42.399056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Iterate over each planet_id\nprint(\"Start:\")\nfor planet_id in tqdm(planet_ids):\n    planet_dir = os.path.join(input_dir, str(planet_id))\n    signal_file = os.path.join(planet_dir, 'FGS1_signal.parquet')\n\n    if os.path.isfile(signal_file):\n        # Step 1: Load and reshape the data\n        signal_data = pd.read_parquet(signal_file)\n        signal_array = signal_data.values\n        reshaped_signal_data = signal_array.reshape(135000, 32, 32)\n\n        # Step 2: Extract the maximum values\n        max_values = np.max(reshaped_signal_data, axis=(1, 2))\n\n        # Step 3: Filter the noise\n        filtered_data_all = lowpass_filter(max_values, cutoff, fs)\n        #print(len(filtered_data_all))\n\n        # Step 4: Store the filtered data in the DataFrame\n        output_df.loc[planet_id] = filtered_data_all\n\n# # Reset index of the output DataFrame\n# output_df.reset_index(inplace=True)\n# output_df.rename(columns={'index': 'planet_id'}, inplace=True)\n\n# Display the resulting DataFrame\nprint(output_df.shape)\noutput_df.head(5)","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:29:50.716606Z","iopub.execute_input":"2024-08-06T17:29:50.717157Z","iopub.status.idle":"2024-08-06T17:29:58.434825Z","shell.execute_reply.started":"2024-08-06T17:29:50.717102Z","shell.execute_reply":"2024-08-06T17:29:58.43362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# %%time\n# # Save the DataFrame to /kaggle/working/\noutput_file_path = 'FGS1_time_max_filtered.csv'\noutput_df.to_csv(output_file_path, index=False)\nprint(f'Saved the output DataFrame to {output_file_path}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 4))\nfor i in range(min(5, len(output_df))):\n    plt.plot(output_df.columns, output_df.iloc[i], label=f'Row {output_df.index[i]}')\n\nplt.legend(title='Index')\nplt.title('Line Plots of the First Five Rows')\nplt.xlabel('Columns')\nplt.ylabel('Values')\nplt.grid(True)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-06T16:02:14.039663Z","iopub.execute_input":"2024-08-06T16:02:14.040901Z","iopub.status.idle":"2024-08-06T16:02:16.44135Z","shell.execute_reply.started":"2024-08-06T16:02:14.040835Z","shell.execute_reply":"2024-08-06T16:02:16.4396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"output_df = output_df.fillna(method='ffill')\ntrain_labels = train_labels.fillna(method='ffill')\noutput_df = output_df.astype(float)\ntrain_labels = train_labels.astype(float)","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:30:02.572134Z","iopub.execute_input":"2024-08-06T17:30:02.572605Z","iopub.status.idle":"2024-08-06T17:30:10.342834Z","shell.execute_reply.started":"2024-08-06T17:30:02.572567Z","shell.execute_reply":"2024-08-06T17:30:10.341801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data processing: extract time series features","metadata":{}},{"cell_type":"code","source":"def extract_features(df):\n    features = pd.DataFrame(index=df.index)\n    \n    features['mean'] = df.mean(axis=1)\n    features['std'] = df.std(axis=1)\n    features['max'] = df.max(axis=1)\n    features['min'] = df.min(axis=1)\n    features['median'] = df.median(axis=1)\n    features['skewness'] = df.skew(axis=1)\n    features['kurtosis'] = df.kurtosis(axis=1)\n    features['iqr'] = df.apply(lambda x: np.percentile(x, 75) - np.percentile(x, 25), axis=1)\n    features['cv'] = df.apply(lambda x: np.std(x) / np.mean(x) if np.mean(x) != 0 else 0, axis=1)\n    features['rms'] = df.apply(lambda x: np.sqrt(np.mean(np.square(x))), axis=1)\n    features['energy'] = df.apply(lambda x: np.sum(np.square(x)), axis=1)\n    \n    features['zero_crossing_rate'] = df.apply(lambda x: ((x[:-1] * x[1:]) < 0).sum(), axis=1)\n    \n    features['first_quartile'] = df.apply(lambda x: np.percentile(x, 25), axis=1)\n    features['third_quartile'] = df.apply(lambda x: np.percentile(x, 75), axis=1)\n    \n    features['autocorrelation'] = df.apply(lambda x: x.autocorr(lag=1), axis=1)\n\n    return features","metadata":{"execution":{"iopub.status.busy":"2024-08-06T17:30:10.344818Z","iopub.execute_input":"2024-08-06T17:30:10.345238Z","iopub.status.idle":"2024-08-06T17:30:10.362194Z","shell.execute_reply.started":"2024-08-06T17:30:10.345207Z","shell.execute_reply":"2024-08-06T17:30:10.360742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"extract_features(output_df)","metadata":{},"execution_count":null,"outputs":[]}]}