{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"In this notebook we'll try to analyze, what we have in the competition by the simplest possible methods. Nevertheless, it could help to understand the data and split it into some subclasses for futher deeper analysis.\n\nLet's start with importing nessesary libraries.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport seaborn as sns\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport h5py\n\nfrom scipy import stats","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:49.487829Z","iopub.execute_input":"2022-11-10T11:21:49.488436Z","iopub.status.idle":"2022-11-10T11:21:50.508537Z","shell.execute_reply.started":"2022-11-10T11:21:49.488324Z","shell.execute_reply":"2022-11-10T11:21:50.507456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1) Load and concat the labels\n\nHere we load the available labels and put train and test data together. It will be more comfortable to work on both datasets simultaniously. We will differ them by the label: in test dataset we have `target`=0.5, and in the train it is 0 or 1.\n\n`dataset_file_mask` will be used for easy extraction of file path based on its `id`","metadata":{}},{"cell_type":"code","source":"df1 = pd.read_csv('../input/g2net-detecting-continuous-gravitational-waves/sample_submission.csv')\ndf2 = pd.read_csv('../input/g2net-detecting-continuous-gravitational-waves/train_labels.csv')\n\ndf = pd.concat((df1, df2[df2.target!=-1]), ignore_index=True)\n\ndataset_file_mask = ('../input/g2net-detecting-continuous-gravitational-waves/test/%s.hdf5',\n                     '../input/g2net-detecting-continuous-gravitational-waves/train/%s.hdf5')","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:50.513738Z","iopub.execute_input":"2022-11-10T11:21:50.51423Z","iopub.status.idle":"2022-11-10T11:21:50.562997Z","shell.execute_reply.started":"2022-11-10T11:21:50.514185Z","shell.execute_reply":"2022-11-10T11:21:50.56153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2) Load the hdf5 file\nThe function below loads the data and also performs some moving average (MA) to reduce the noise","metadata":{}},{"cell_type":"code","source":"def get_MA_SFT(n, MA=12):\n    # Load `id` and `target` from dataframe\n    id = df.iloc[n]['id']\n    target = df.iloc[n]['target']\n\n    # Get the hdf5 file name\n    fname = dataset_file_mask[int(target!=0.5)] % id\n\n    # Read the data from the file\n    data = h5py.File(fname, 'r')[id]\n\n    freqs = np.array(data['frequency_Hz'])\n    ma_sft = {}\n    for detector in ['L1', 'H1']:\n        # Get the SFT values\n        sft = np.array(data[detector]['SFTs'])*1e22\n\n        # Transform them into power spectrum values\n        power = np.power(np.abs(sft), 2)\n\n        # Now calculate the moving average in every `MA` time points\n        # For this we get the `cumsum` along the time axis and substract for every `MA` time points\n        # The last block may not contain `MA` points, and it is ignored\n        power_ma_cumsum = np.cumsum(power, axis=1)\n        power_ma_cumsum_0 = np.concatenate((np.zeros((power.shape[0],1)), power_ma_cumsum), axis=1)[:,::MA]\n        power_ma = np.diff(power_ma_cumsum_0, axis=1)[:,:-1]/MA\n        ma_sft[detector] = power_ma\n        \n    return ma_sft, freqs","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:50.565065Z","iopub.execute_input":"2022-11-10T11:21:50.565921Z","iopub.status.idle":"2022-11-10T11:21:50.578946Z","shell.execute_reply.started":"2022-11-10T11:21:50.565859Z","shell.execute_reply":"2022-11-10T11:21:50.577689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3) Testing some of the records\n\nThe idea of the testing is to plot the average values of SFT with respect both to the frequency or to the time axis. \nThen we do the statistical test, whether or not the achieved values satisfy the normal distribution.\n\nOur null hypothesis will be 'the record is a white gaussian noise'. We calculate the corresponding p-values. If they are suffisiently low, we reject the null hypothesis and say, that the signal is not gaussian. So we must pay attention to it.","metadata":{}},{"cell_type":"code","source":"def test_record_view(id):\n    # SFT with MA\n    ma_sft, freqs = get_MA_SFT(id)\n\n    fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(18,10))\n    ax_names = ('Interval No', 'Frequency, Hz')\n    for col, detector in enumerate(['L1', 'H1']):\n        # Prepare the axes\n        axs[0, col].set_title(f'Detector={detector}')\n        axs[0, col].set_xlabel(ax_names[0])\n        axs[0, col].set_ylabel(ax_names[1])\n        axs[0, col].pcolormesh(range(ma_sft[detector].shape[1]), freqs, ma_sft[detector])\n\n        for row in range(2):\n            # Calculate mean value for each frequency or for each time interval\n            x = np.mean(ma_sft[detector], axis=row)\n            k2, p = stats.normaltest(x)\n\n            axs[row+1, col].set_title(f'Detector={detector}; p-value={p:.3}')\n            axs[row+1, col].set_xlabel(ax_names[row])\n            if row==0:\n                axs[row+1, col].plot(x)\n            else:\n                axs[row+1, col].plot(freqs, x)\n                \n    fig.tight_layout(pad=1.5)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:50.582387Z","iopub.execute_input":"2022-11-10T11:21:50.583792Z","iopub.status.idle":"2022-11-10T11:21:50.598152Z","shell.execute_reply.started":"2022-11-10T11:21:50.583723Z","shell.execute_reply":"2022-11-10T11:21:50.596563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.1) `id`=609","metadata":{}},{"cell_type":"code","source":"test_record_view(id=609)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:50.600451Z","iopub.execute_input":"2022-11-10T11:21:50.601429Z","iopub.status.idle":"2022-11-10T11:21:52.576085Z","shell.execute_reply.started":"2022-11-10T11:21:50.601362Z","shell.execute_reply":"2022-11-10T11:21:52.574883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* It can be seen from the upper two graphs, that this record does contain the signal. At the same time the noise is not stationary.\n* From the middle two graphs we can also see, that the noise is not stationary.\n  P-value is low (2e-28) only for the left graph. For the right it is not so low (0.155), that we can reject the null hypothesis.\n* From the bottom graphs we can say, it looks like white noise, and the corresponding p-values (0.287 and 0.731) are high.\n\nConclusion: **the record has non-gaussian noise background and may or may not contain a signal. We must pay attention to it.**\n\nOf course, we can see the signal on the upper graphs by our eyes, but the used statistical criteria is rather weak for this))","metadata":{}},{"cell_type":"markdown","source":"## 3.2) `id`=7850","metadata":{}},{"cell_type":"code","source":"test_record_view(id=7850)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:52.57731Z","iopub.execute_input":"2022-11-10T11:21:52.577738Z","iopub.status.idle":"2022-11-10T11:21:54.266377Z","shell.execute_reply.started":"2022-11-10T11:21:52.577677Z","shell.execute_reply":"2022-11-10T11:21:54.265146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* It can be seen from the upper two graphs, that this record does contain the signal (line at the bottom).\n* From the middle two graphs we can see, that the noise is likely to be stationary at L1 and non-stationary at H1 (p-value is about 0.006).\n* From the bottom graphs we can say, it is completly non-gaussian. The record may contain a signal.\n\nConclusion: **the record may contain a signal. We must pay attention to it.**","metadata":{}},{"cell_type":"markdown","source":"## 3.3) `id`=3283","metadata":{}},{"cell_type":"code","source":"test_record_view(id=3283)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:54.268298Z","iopub.execute_input":"2022-11-10T11:21:54.269088Z","iopub.status.idle":"2022-11-10T11:21:55.997886Z","shell.execute_reply.started":"2022-11-10T11:21:54.269041Z","shell.execute_reply":"2022-11-10T11:21:55.99655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* We can't see the signal on the left upper graph. The right upper graph is something strange, because the signal is to strong. Also it doesn't have frequency gradient, so it may be just (electrical?) interference. If we had such strong signal on H1, we must had definitely registered it at L1, too. But this is not the case.\n* From the middle left graph we can see, that the noise is not stationary at L1. The p-value is also low fow H1, but we need to exclude the interference first to make a conclusion here.\n* From the bottom graphs we can clearly see the interference at H1.\n\nConclusion: **the record may or may not contain a signal. But it definitely contains interference that should not be confused with it.**","metadata":{}},{"cell_type":"markdown","source":"## 3.4) `id`=8100","metadata":{}},{"cell_type":"code","source":"test_record_view(id=8100)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:55.999337Z","iopub.execute_input":"2022-11-10T11:21:55.999754Z","iopub.status.idle":"2022-11-10T11:21:57.732718Z","shell.execute_reply.started":"2022-11-10T11:21:55.999695Z","shell.execute_reply":"2022-11-10T11:21:57.731223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is the record from the train set. \nAll the graphs say the same thing: **the record looks like white noise, and me have to investigate it further**.","metadata":{}},{"cell_type":"markdown","source":"# 4) Process all records\n\nNow we do the same, but with all records and without plots. Then we'll compare the results for train and test sets.","metadata":{}},{"cell_type":"code","source":"def test_record(id):\n    ma_sft, freqs = get_MA_SFT(id)\n\n    results = {'freq' : freqs[0]}\n    \n    for detector in ['L1', 'H1']:\n        results[detector] = []\n        for row in range(2):\n            # Calculate mean value for each frequency or for each time interval\n            x = np.mean(ma_sft[detector], axis=row)\n            k2, p = stats.normaltest(x)\n            results[detector].append(p)\n            \n    return [results['freq']] + results['L1'] + results['H1']","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:57.735182Z","iopub.execute_input":"2022-11-10T11:21:57.735867Z","iopub.status.idle":"2022-11-10T11:21:57.748259Z","shell.execute_reply.started":"2022-11-10T11:21:57.735811Z","shell.execute_reply":"2022-11-10T11:21:57.746571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"columns = ['freq', 'L1_time', 'L1_freq', 'H1_time', 'H1_freq']\nif True:\n    !wget http://93.180.51.46/exploration/df.csv -O ./df.csv\n    df = pd.read_csv('df.csv', index_col=0)\nelse:\n    # Add zero-filled columns to our dataframe\n    df[columns] = 0\n\n    # Now calculate p-values and put it inside the dataframe\n    for id in range(len(df)):\n        if id % 100 == 0: \n            print(f'Processing #{id}')\n        df.loc[id, columns] = test_record(id)\n    \ndf.head()","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:21:57.753843Z","iopub.execute_input":"2022-11-10T11:21:57.754816Z","iopub.status.idle":"2022-11-10T11:22:00.419132Z","shell.execute_reply.started":"2022-11-10T11:21:57.754759Z","shell.execute_reply":"2022-11-10T11:22:00.417296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5) Finally, let's display some statistics\n\nWe introduse some `treshold` and a set of criteria, when we reject the particulra null hypothesis. All is put into one function","metadata":{}},{"cell_type":"code","source":"def print_criteria_statistics(ds, treshold=0.01):\n    # A set of criteria for every column\n    criteria = [df[col]<treshold for col in columns[1:]]\n    \n    # A dict of criteria for datasets\n    criteria_ds = {'train' : df.target != 0.5, 'test': df.target == 0.5}\n\n    # Total length\n    total = len(df[criteria_ds[ds]])\n\n    print(f'Dataset = {ds}; length = {total}')\n    print('-'*80)\n\n    for n, col in enumerate(columns[1:]):\n        l = len(df[criteria_ds[ds] & criteria[n]])\n        print(f'{col}  is not gauss for\\t{l}/{total},\\twhich is {100*l/total:.3}%')\n\n    print('-'*80)\n\n    l = len(df[criteria_ds[ds] & criteria[0] & criteria[2]])\n    print(f'Time L&H is not gauss for\\t{l}/{total},\\twhich is {100*l/total:.3}%')\n\n    l = len(df[criteria_ds[ds] & criteria[1] & criteria[3]])\n    print(f'Freq L&H is not gauss for\\t{l}/{total},\\twhich is {100*l/total:.3}%')\n\n    print('-'*80)\n\n    l = len(df[criteria_ds[ds] & (criteria[0] | criteria[2])])\n    print(f'Time L|H is not gauss for\\t{l}/{total},\\twhich is {100*l/total:.3}%')\n\n    l = len(df[criteria_ds[ds] & (criteria[1] | criteria[3])])\n    print(f'Freq L|H is not gauss for\\t{l}/{total},\\twhich is {100*l/total:.3}%')","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:22:00.420994Z","iopub.execute_input":"2022-11-10T11:22:00.42145Z","iopub.status.idle":"2022-11-10T11:22:00.435071Z","shell.execute_reply.started":"2022-11-10T11:22:00.421403Z","shell.execute_reply":"2022-11-10T11:22:00.433377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, display the results","metadata":{}},{"cell_type":"code","source":"print_criteria_statistics(ds = 'train')\n\nprint('\\n')\n\nprint_criteria_statistics(ds = 'test')","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:22:00.437256Z","iopub.execute_input":"2022-11-10T11:22:00.437914Z","iopub.status.idle":"2022-11-10T11:22:00.466368Z","shell.execute_reply.started":"2022-11-10T11:22:00.437855Z","shell.execute_reply":"2022-11-10T11:22:00.464799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It's interesting! So...\n* **Stationarity.**\n\nIn the train set we never have non stationary noise on both detectors simultaniously, and only in 3.5% we have it on one of the detectors.\nIn the test set we have it on both detectors in 11.2% and on at least one of them in 21.3%.\n\n* **Frequency \"artifacts\"**\n\nThese artifacts could be potential signals or a glitches, or an interferences.\nIn the train set we have it in 7% on the both detectors and in 15% on at least one of them. In the test the corresponding values are 1.9% and 8.8%. \n\nNow let's set even less value for the `treshold` and look, what happens.","metadata":{}},{"cell_type":"code","source":"print_criteria_statistics(ds = 'train', treshold=1e-5)\n\nprint('\\n')\n\nprint_criteria_statistics(ds = 'test', treshold=1e-5)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:22:00.468853Z","iopub.execute_input":"2022-11-10T11:22:00.469394Z","iopub.status.idle":"2022-11-10T11:22:00.490938Z","shell.execute_reply.started":"2022-11-10T11:22:00.469339Z","shell.execute_reply":"2022-11-10T11:22:00.48991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see, that now all train signals look stationary (speaking corectly, we don't reject our null hypothesis), and 18.5% of test signals not. But the most interesting is that frequency artifacts are not independent on both detectors. If it was so, than, for test set, we'll have 0.0176*0.0406 = 0.07% of not gauss 'Freq L&H', but we have 1.17%. So, almost all of the records with not gauss 'Freq L&H' (we have selected 93 of them) are likely to contain a signal.\n\nLet's get their exact ids and look at them to see, if it's true.","metadata":{}},{"cell_type":"markdown","source":"# 6) Investigating not gaussian 'Freq L&H'","metadata":{}},{"cell_type":"markdown","source":"Let's look at, say, first 5 candidates","metadata":{}},{"cell_type":"code","source":"treshold = 1e-5\n\ncandidates = df[(df.target==0.5) & (df.L1_freq<treshold) & (df.H1_freq<treshold)].index.tolist()\n\nfor n in candidates[:5]:\n    print(f'Candidate {n}\\t Target={df.target[n]}')\n    test_record_view(id=n)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T11:22:00.492192Z","iopub.execute_input":"2022-11-10T11:22:00.49264Z","iopub.status.idle":"2022-11-10T11:22:09.659161Z","shell.execute_reply.started":"2022-11-10T11:22:00.492599Z","shell.execute_reply":"2022-11-10T11:22:09.657675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The first one (#50) is very interesting! On L1 detector we see the increasing frequency, and on H1 detector it is different and approximately constant. So, we may deal with the signal and the interference at the same time))\nWhen looking carefully we can also see the small peak near 250.1Hz on H1 detector, which may correspond to the signal on L1.\n\nOther 4 records have approximately the same frequencies of peaks on L1 and H1.","metadata":{}},{"cell_type":"markdown","source":"You may notice, that if the signal frequency varies a lot, it doesn't affect the average power spectrum much. This signal will be definitely lost with the proposed approach. So, to conclude, there's much work, and we must be very carefull with the signals we have. ","metadata":{}},{"cell_type":"markdown","source":"# 7) Getting rid of artifacts\n\nWe are facing at least 3 types of artifacts.\n* **Not stationary background noise**\nThis may be the easiest artifact to suppress. To do so, we can just normalize the record by subtracting the mean value from the `ma_sft` and divide it by its standart deviation. \nThe disadvantage is that if we have a glitch, its amplitude affects the mean value and standart deviation. So, the signal (if it exists in the record) is affected by the procedure. We have to take care about it.\n\n* **(Elecromagnetic?) interferencies** on the single frequency, which are active on the whole record.\nThe features of these artifacts are, first, the constant frequency, second, the presence of the interference on the only detector and, third, relatively large amplitude. The idea is to put all the values on this frequency to zero.\n\n* **Glitches**, which are active only a small amount of time.\nAs was said before, they affect both time and frequency amplitude distributions. Like electromagnetic interferencies, they present on the only detector and have relatively high amplitude. We have to suppress not the single frequency, but some area of `ma_sft` array. Then, when calculating statistical properties of the remainder, we must take into account the absence of these points. Otherwise we may get wrong result.","metadata":{}},{"cell_type":"code","source":"class G2Preprocessor:\n    detectors = ['L1', 'H1']\n    other_detector ={'H1':'L1', 'L1':'H1'}\n    \n    def __init__(self, df, dataset_file_mask):\n        self.df = df\n\n    def load(self, id):\n        \"\"\"\n            Loads the record from dataframe and hdf5 file\n            \n            id: the raw index of the record (i.e. 0, 1, 2 ...)\n            \n            Returns: power spectrograms on each of the detectors\n        \"\"\"\n        \n        # Load `id` and `target` from dataframe\n        self.id_ = id\n        self.id = df.iloc[id]['id']\n        self.target = df.iloc[id]['target']\n\n        # Get the hdf5 file name and read the data from it\n        fname = dataset_file_mask[int(self.target!=0.5)] % self.id\n        data = h5py.File(fname, 'r')[self.id]\n\n        # Get the power values and frequencies\n        self.freqs = np.array(data['frequency_Hz'])\n        self.sft = {}\n        for detector in self.detectors:\n            # We multiply here by 1e22 not to deal with small numbers\n            self.sft[detector] = np.array(data[detector]['SFTs'])*1e22\n            \n        return self.sft\n        \n    def get_power(self, sft=None):\n        \"\"\"\n            Calculates the power of SFT, i.e. the square of its absolute value\n            \n            sft: the SFTs of the record. If it is None (default) uses the sft stored in the class variable\n            \n            Returns: power spectrograms on each of the detectors\n        \"\"\"\n        if sft is None:\n            sft = self.sft\n        power = {}\n        for detector in self.sft.keys():\n            power[detector] = np.power(np.abs(sft[detector]), 2)\n            \n        return power\n        \n    def exclude_time_outliers(self, power, time_treshold=5):\n        \"\"\"\n            Excludes some of the samples on the time axis, which seem to be outliers\n            \n            power: the power spectrograms\n            time_treshold: the treshold used to find outliers. We call the sample outlier, if its value is bigger, than q75 + treshold*(q75-q25),\n                           where q75 and q25 are corresponding 75% and 25% quantiles\n            \n            Returns: power spectrograms on each of the detectors\n        \"\"\"\n            \n        for detector in power.keys():\n            # Here we take median, as it is more robust to the interference, which can exist on certain frequencies\n            power_time_median = np.median(power[detector], axis=0)  \n            # If this median falls outside some interval, we exclude the correspondin point\n            q25, q75 = np.quantile(power_time_median, [0.25, 0.75])\n            power_time_mask = np.where(power_time_median < q75+time_treshold*(q75-q25))[0]\n            power[detector] = power[detector][:, power_time_mask]\n        \n        return power\n    \n    def _get_intervals(self, stack, needle):\n        \"\"\"\n            For every item in the `needle` we return the closest points to it from the `stack` (to the left and to the right)\n            So, these two points represent thesmallest interval containing the needle.\n        \"\"\"\n        stack = np.array(stack).reshape((1,-1))\n        needle = np.array(needle).reshape((-1,1))\n        result = np.array((np.max(np.where(stack < needle, stack, -1), axis=1),\n                np.min(np.where(stack > needle, stack, needle), axis=1))).T\n        return result\n\n    def _points_to_intervals(self, points):\n        \"\"\"\n            Given a set of points (which are integers), this function rearanges them and returns the\n            set of intervals, both ends are included\n        \"\"\"\n        # 1) Sort the points\n        points_sorted = np.sort(points)\n\n        # 2) The distance d[n] is 1 if the point n has a neighbour to the right\n        #    The distance d[n-1] is 1 if the point n has a neighbour to the left\n        points_dist = np.diff(points_sorted) != 1\n        \n        # 3) We include the first and the last point, as they are definitely ends\n        points_right = np.append(points_dist, True)\n        points_left = np.insert(points_dist, 0, True)\n        #is_end = np.logical_or(points_left, points_right)\n        #return  points_sorted[is_end]\n        \n        return np.vstack((points_sorted[points_left], points_sorted[points_right]))\n\n    def _intervals_containing_needle(self, intervals, needle):\n        \"\"\"\n            Returns only those intervals, which contain at least one point from the needle array\n        \"\"\"\n        int_mtr = np.logical_and(intervals[0].reshape((-1,1))<=needle, intervals[1].reshape((-1,1))>=needle)\n        int_idx = np.sum(int_mtr, axis=1)\n        return intervals[:, np.where(int_idx)[0]]\n    \n    def get_time_ma(self, power, window=12):\n        \"\"\"\n            Calculates the simplest moving average on the time axis. The samples are averaged inside not intersecting rectangular windows\n            \n            power: the power spectrograms\n            window: the window width\n            \n            Returns: averaged power spectrograms on each of the detectors\n        \"\"\"\n            \n        for detector in power.keys():\n            # For this we get the `cumsum` along the time axis and substract for every `window` time points\n            # The last block may not contain `window` points, and it is ignored\n            power_ma_cumsum = np.cumsum(power[detector], axis=1)\n            power_ma_cumsum_0 = np.concatenate((np.zeros((power[detector].shape[0],1)), power_ma_cumsum), axis=1)[:,::window]\n            power[detector] = np.diff(power_ma_cumsum_0, axis=1)[:,:-1]/window\n        \n        return power\n\n    def norm_time(self, power):\n        \"\"\"\n            Normalize the data along the time axis by subtracting its mean value and dividing by standart deviation\n            \n            power: the power spectrograms\n            \n            Returns: normilized power spectrograms on each of the detectors\n        \"\"\"\n            \n        for detector in power.keys():\n            power_mean = np.mean(power[detector], axis=0)\n            power_std = np.std(power[detector], axis=0)\n\n            power[detector] = (power[detector]-power_mean)/power_std\n        return power\n\n    def exclude_freq_outliers(self, power, freq_treshold1=5., freq_treshold2=1., max_interval=5):\n        \"\"\"\n            Excludes some of the samples on the frequency axis, which seem to be outliers, i.e. interferences\n            \n            power: the power spectrograms.\n            freq_treshold1, freq_treshold2: \n                   the tresholds used to find outliers. We call the sample outlier, if its value is bigger, than q75 + treshold*(q75-q25),\n                   where q75 and q25 are corresponding 75% and 25% quantiles\n                   The second condition is that we must see the outlier on the only detector. If it present on the both, it is a signal!\n                   So we use 2 treshold: to be an outlier the signal must be over freq_treshold1 on one of the detector and below treshold2 on the other.\n            max_interval: the maximum possible excluded interval length (in points). If this length is exeeded, the set of failed detectors is returned.\n            \n            Returns: power spectrograms on each of the detectors on success\n                     the set of failed detectors on failure\n        \"\"\"\n\n        freq_treshold = [freq_treshold1, freq_treshold2, 0]\n        power_freq_neg_mask = [{}, {}, {}]\n        \n        # Now we only get the frequencies, which are candidates to exclude\n        # Notice: these are negative masks, i.e. they contain the frequencies to exclude, not to remain\n        power_freq_mean = {}\n        power_freq_mean_median = {}\n        for detector in power.keys():\n            power_freq_mean[detector] = np.mean(power[detector], axis=1)\n            power_freq_mean_median[detector] = np.median(power_freq_mean[detector])\n            q25, q75 = np.quantile(power_freq_mean[detector], [0.25, 0.75])\n            for n in range(len(freq_treshold)):\n                power_freq_neg_mask[n][detector] = np.where(power_freq_mean[detector] > q75+freq_treshold[n]*(q75-q25))[0]\n                print(power_freq_neg_mask[n][detector])\n                print(q75+freq_treshold[n]*(q75-q25), '=>', power_freq_neg_mask[n][detector])\n                #plt.figure()\n                #plt.plot(power_freq_mean)\n            \n            #mm = np.mean(power_freq_mean)\n            #print(f'Detector={detector}, mean={mm}')\n\n        # Now we have 3 sets of points for each of the detectors, each of them corresponds to a treshold\n        # For each point higher, than treshold1, we select the continious interval of points higher than treshold2,\n        # which includes this point. This interval is a candidate for removal\n        exclusion_masks = {}\n        for detector in power.keys():\n            intervals = self._points_to_intervals(power_freq_neg_mask[2][detector])\n            #print(intervals)\n            intervals_candidates = self._intervals_containing_needle(intervals, power_freq_neg_mask[0][detector])\n            \n            # The outliers of other detector (if the other detector data exists)\n            if self.other_detector[detector] in power_freq_neg_mask[1]:\n                other_detector_outliers = set(power_freq_neg_mask[1][self.other_detector[detector]])\n            else:\n                other_detector_outliers = set()\n\n            # Now check each interval candidate\n            exclusion_mask = set()\n            detectors_failed = []\n            for left, right in intervals_candidates.T:\n                # Check, if the interval is small enough\n                if right-left+1 > max_interval:\n                    print(f'!!! The interval of [{left}; {right}] is too big to delete from {self.id_}/{self.id} record, detector={detector}')\n                    detectors_failed.append(detector)\n                    break\n\n                # Notice `right+1`, which means we do include the right end\n                interval = set(range(left,right+1))\n\n                # If the selected interval has no intersections with outliers of the other detector,\n                # and the signal amplitudes on two detectors are comparable,\n                # we add it to our exclusion mask\n                intersection = list(interval & other_detector_outliers)\n                if len(intersection) == 0:\n                    exclusion_mask = exclusion_mask | interval\n                else:\n                    # Get the amplitudes of current and other detector\n                    a_cur   = power_freq_mean[detector] - power_freq_mean_median[detector]\n                    a_other = power_freq_mean[self.other_detector[detector]] - power_freq_mean_median[self.other_detector[detector]]\n                    # Their ratio\n                    a_rel   = np.max(a_cur[intersection])/np.max(a_other[intersection])\n                    #print(f'a_rel={a_rel} at {detector}:{interval}')\n                    if a_rel > 2:\n                        exclusion_mask = exclusion_mask | interval\n\n            power[detector] = np.delete(power[detector], list(exclusion_mask), axis=0)\n            print(f'detector={detector}')\n            print('intervals=', intervals_candidates.T)\n            print('other_detector_outliers', other_detector_outliers)\n            print(f'finally exclude from {detector}:', exclusion_mask)\n            exclusion_masks[detector] = exclusion_mask\n        \n        detectors_failed = set(detectors_failed)\n        for detector in detectors_failed:\n            del power[detector]\n        return power, detectors_failed, exclusion_masks\n    \n    def check_gauss(self, signal, treshold=1e-10, axis=None):\n        \"\"\"\n            Calculates the amount of non-gaussian samples in the provided signal.\n        \"\"\"\n        if axis is None:\n            return self.check_gauss(signal, treshold=treshold, axis=0) + self.check_gauss(signal, treshold=treshold, axis=1)\n        k2, p = stats.normaltest(np.sqrt(signal), axis=axis)\n        return np.sum(p < treshold) / len(p)\n\n    def delete_bad_detectors(self, power, treshold=0.5):\n        \"\"\"\n            Deletes the data from detectors, which is considered bad\n        \"\"\"\n        result = {}\n        for detector in power.keys():\n            if(self.check_gauss(power[detector]) < treshold):\n                result[detector] = power[detector]\n        return result\n\n    def subtract_median(self, power, axis=1):\n        \"\"\"\n            Subtracts the median values from the power data to make the signal more homogenious\n            \n            power: the power spectrograms\n            axis: the axis for which to prform the subtraction (0 for time and 1 for freq, default freq)\n            \n            Returns: power spectrograms on each of the detectors\n        \"\"\"\n        result = {}\n        for detector in power.keys():\n            med = np.median(power[detector], axis=axis)\n            #medf = np.flip(med)\n            #med3 = np.concatenate([medf, med, medf])\n            #N = 60\n            #r_med3 = np.convolve(med3, np.ones(N)/N, mode='same').reshape((-1,1))[len(med):2*len(med)]\n            gp = gaussian_filter1d(med, sigma=60).reshape((-1,1))\n            result[detector] = power[detector] - gp + np.mean(gp)\n        return result\n\ng2p = G2Preprocessor(df = df, dataset_file_mask = dataset_file_mask)\n\"\"\"    \nfrom scipy.ndimage import gaussian_filter1d\n\nfrom numpy.fft import fft, ifft\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\n\nclass Resonance(nn.Module):\n\n    def __init__(self, len=360):\n        super(Resonance, self).__init__()\n        \n        self.a =  nn.Parameter(torch.randn(1))\n        self.b =  nn.Parameter(torch.randn(1))\n        self.c =  nn.Parameter(torch.randn(1))\n        \n        self.x = torch.arange(0, len)\n\n    def forward(self):\n        return self.a*self.x**2 + self.b*self.x + self.c\n    #/(self.a*torch.pow(self.x-self.x0, 2)+1)\n        \"\"\"\nfor n in [128,6489,6809,7474,7602,335,6220,5095,2866,2645]:#range(0):#[\n          #1389,8262,8146,8327,8033,\n          #3235,8380,5366,8363,7024,\n          #5863,8553,361,5730,8314,\n          #1147,1782,4068,8293,1070,\n          #2866,5164,8414,8072,4330,\n          #2394,5014,7850,8407,8421,\n          #8509,7240,\n        #]:\n    print('-'*80)\n    print(f'id={n}')\n\n    sft = g2p.load(n)\n    ps = g2p.get_power(sft)\n    ps = g2p.exclude_time_outliers(ps)\n    ps = g2p.get_time_ma(ps)\n    #ps = g2p.subtract_median(ps)\n    #ps = g2p.delete_bad_detectors(ps)\n    #ps = g2p.norm_time(ps)\n\n    print(f'Considering detectors: {ps.keys()}')\n\n    ps, _, _= g2p.exclude_freq_outliers(ps, freq_treshold1=5., freq_treshold2=1., max_interval=50)\n\n    axes = ('time', 'freq')\n    for detector in ps.keys():\n        for axis in [0,1]:\n            res = g2p.check_gauss(ps[detector], axis=axis)\n            print(f'{detector} has {res:.3} non-gaussian samples in axis {axes[axis]}')\n\n        fig, axs = plt.subplots(ncols=2, figsize=(18,5))\n        fig.suptitle(f'id={n}/{detector}')\n        x = np.mean(ps[detector], axis=1)\n\n        \"\"\"res = Resonance()\n        loss = nn.MSELoss()\n        optimizer = optim.Adam(res.parameters(), lr=1e-3)\n\n        for step in range(10000):\n            y=res()\n            loss_value = loss(y, 1/torch.tensor(x).type(torch.FloatTensor))\n            optimizer.zero_grad()\n            loss_value.backward()\n            optimizer.step()\n            if step % 100 ==0 :\n                print(f'\\t{step} => {loss_value};')\n             \n            #print(res.a, res.b, res.x0)\"\"\"\n            \n        \n        axs[0].plot(x);\n        #axs[0].plot(1/y.detach().numpy() );\n        axs[1].pcolormesh(abs(ps[detector]));\n        fig.savefig(f'{n}_{detector}.png')","metadata":{"execution":{"iopub.status.busy":"2022-11-10T14:45:03.882263Z","iopub.execute_input":"2022-11-10T14:45:03.882762Z","iopub.status.idle":"2022-11-10T14:45:27.92057Z","shell.execute_reply.started":"2022-11-10T14:45:03.882723Z","shell.execute_reply":"2022-11-10T14:45:27.919085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"check_df = pd.DataFrame(columns=['id','L1_exclude', 'H1_exclude', 'failed_detectors', 'L1_time', 'L1_freq', 'H1_time', 'H1_freq'])\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"g2p = G2Preprocessor(df = df, dataset_file_mask = dataset_file_mask)\n\n\nfor n in range(6971, len(df)):\n    print('-'*80)\n    print(f'id={n}')\n    \n    sft = g2p.load(n)\n    ps = g2p.get_power(sft)\n    ps = g2p.exclude_time_outliers(ps)\n    ps = g2p.get_time_ma(ps)\n    #ps = g2p.subtract_median(ps)\n    ps = g2p.delete_bad_detectors(ps)\n\n    print(f'Considering detectors: {ps.keys()}')\n\n    ps, detectors_failed, exclusion_masks = g2p.exclude_freq_outliers(ps, freq_treshold1=5., freq_treshold2=1., max_interval=50)\n    #ps = g2p.norm_time()\n\n    axes = ('time', 'freq')\n    res = [{'L1':'-','H1':'-'}] * 2\n    for detector in ps.keys():\n        for axis in [0,1]:\n            res[axis][detector] = g2p.check_gauss(ps[detector], axis=axis)\n            print(f'{detector} has {res[axis][detector]:.3} non-gaussian samples in axis {axes[axis]}')\n    \n    for detector in g2p.detectors:\n        if detector not in exclusion_masks.keys():\n            exclusion_masks[detector]='-'\n            \n    check_df.loc[len(check_df.index)] = [n, exclusion_masks['L1'], exclusion_masks['H1'], detectors_failed,\n                                         res[0]['L1'], res[1]['L1'], res[0]['H1'], res[1]['H1']]\n\ncheck_df.to_csv('check.csv')","metadata":{"_kg_hide-output":false,"execution":{"iopub.status.busy":"2022-11-10T13:08:27.080525Z","iopub.execute_input":"2022-11-10T13:08:27.081104Z","iopub.status.idle":"2022-11-10T13:08:33.569962Z","shell.execute_reply.started":"2022-11-10T13:08:27.081064Z","shell.execute_reply":"2022-11-10T13:08:33.567718Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"check_df","metadata":{"execution":{"iopub.status.busy":"2022-11-10T13:08:40.128869Z","iopub.execute_input":"2022-11-10T13:08:40.130692Z","iopub.status.idle":"2022-11-10T13:08:40.164893Z","shell.execute_reply.started":"2022-11-10T13:08:40.130628Z","shell.execute_reply":"2022-11-10T13:08:40.163365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"check_df.head(20)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T12:22:19.841393Z","iopub.status.idle":"2022-11-10T12:22:19.841959Z","shell.execute_reply.started":"2022-11-10T12:22:19.841728Z","shell.execute_reply":"2022-11-10T12:22:19.841752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"check_df.to_csv('check_6970.csv')","metadata":{"execution":{"iopub.status.busy":"2022-11-10T12:50:15.795926Z","iopub.execute_input":"2022-11-10T12:50:15.796515Z","iopub.status.idle":"2022-11-10T12:50:15.861432Z","shell.execute_reply.started":"2022-11-10T12:50:15.796475Z","shell.execute_reply":"2022-11-10T12:50:15.859786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"g2p = G2Preprocessor(df = df, dataset_file_mask = dataset_file_mask)\nfor n in range(0, len(df)):\n    print('-'*80)\n    print(f'id={n}')\n    \n    g2p.load(n)\n    g2p.get_power()\n    g2p.exclude_time_outliers()\n    g2p.get_time_ma()\n    ps = g2p.exclude_freq_outliers(freq_treshold1=5., freq_treshold2=1.)\n","metadata":{"execution":{"iopub.status.busy":"2022-11-10T12:22:19.843557Z","iopub.status.idle":"2022-11-10T12:22:19.844102Z","shell.execute_reply.started":"2022-11-10T12:22:19.843851Z","shell.execute_reply":"2022-11-10T12:22:19.843877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"g2p = G2Preprocessor(df = df, dataset_file_mask = dataset_file_mask)\n#for n in [1147,1658,1834,2394,3235,3593,3977]:  # These are VERY BAD records\nfor n in [1147,1658,1834,2394,3235,3593,3977,             25,50,361,464,610,831,886,915,1249,1595,2059,\n          2608,2693,3020,3083,3132,3211,3265,3908,4029,4061,4068,4194,4243]: # Just BAD, but to much)\n#for n in range(50):  # Good?? Or not??\n    print('-'*80)\n    print(f'id={n}')\n    \n    g2p.load(n)\n    g2p.get_power()\n    g2p.exclude_time_outliers()\n    g2p.get_time_ma()\n    g2p.exclude_freq_outliers(freq_treshold1=5., freq_treshold2=1.)\n    ps = g2p.norm_time()\n    plt.figure(figsize=(24,18))\n    plt.pcolormesh(ps['H1'])\n    test_record_view(id=n)","metadata":{"execution":{"iopub.status.busy":"2022-11-06T07:10:31.507913Z","iopub.status.idle":"2022-11-06T07:10:31.508222Z","shell.execute_reply.started":"2022-11-06T07:10:31.508068Z","shell.execute_reply":"2022-11-06T07:10:31.508086Z"}}},{"cell_type":"code","source":"len(df)","metadata":{"execution":{"iopub.status.busy":"2022-11-10T13:18:08.124775Z","iopub.execute_input":"2022-11-10T13:18:08.125432Z","iopub.status.idle":"2022-11-10T13:18:08.137463Z","shell.execute_reply.started":"2022-11-10T13:18:08.125385Z","shell.execute_reply":"2022-11-10T13:18:08.135423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"def test_record_view2(id):\n    # SFT with MA\n    ma_sft, freqs = get_MA_SFT(id, MA=1)\n    \n    \"\"\"    for col, detector in enumerate(['L1', 'H1']):\n            sft=np.random.randn(*ma_sft[detector].shape)\n            ma_sft[detector] = np.power(sft, 2)\n    \"\"\"\n    fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(18,24))\n    ax_names = ('Interval No', 'Frequency, Hz')\n    for col, detector in enumerate(['L1', 'H1']):\n        # Calculate mean & std value for each time interval\n        mean_time = np.mean(ma_sft[detector], axis=0)\n        std_time  = np.std(ma_sft[detector], axis=0)\n        rem_time = (ma_sft[detector] - mean_time)/std_time\n        \n        # Prepare the axes\n        axs[0, col].set_title(f'Detector={detector}')\n        axs[0, col].set_xlabel(ax_names[0])\n        axs[0, col].set_ylabel(ax_names[1])\n        c = axs[0, col].pcolormesh(range(rem_time.shape[1]), freqs, rem_time)\n        plt.colorbar(c, ax=axs[0, col])\n\n        # Calculate mean value for each frequency and variation for each time interval\n        x = np.mean(rem_time, axis=1)\n        k2, p = stats.normaltest(x)\n\n        axs[1, col].set_title(f'Detector={detector}; p-value={p:.3}')\n        axs[1, col].set_xlabel(ax_names[1])\n        axs[1, col].plot(freqs, x)\n                \n    fig.tight_layout(pad=1.5)\n    plt.show()\n    \n    \ntest_record_view2(6988)","metadata":{"execution":{"iopub.status.busy":"2022-11-06T07:10:31.511107Z","iopub.status.idle":"2022-11-06T07:10:31.511381Z","shell.execute_reply.started":"2022-11-06T07:10:31.511239Z","shell.execute_reply":"2022-11-06T07:10:31.511253Z"}}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"            \n            \n            \n        \"\"\"# We have to exclude not only the peaks themselves, but also some area near the peak, which is affected by interference\n        # We use the values from `power_freq_neg_mask[1][other_detector]` to check if the considered peak is an outlier for it (with freq_treshold2)\n        # In this case we have two outliers at the same frequency, so they are likely to be a signal, not an interference\n        detectors_failed = []\n        for detector in self.detectors:\n            if len(power_freq_neg_mask[0][detector]) > 0:\n                # 1st derivative sign: if 1 is folowed by -1, it is maximum, and \n                #                      if -1 is folowed by 1, it is minimum\n                power_freq_mean_ds = np.sign(np.diff(power_freq_mean))\n\n                # The indices of all minima of power_freq_median. \n                # Note '+1', which is used because of the shift made by taking the derivative\n                power_freq_mean_min_idx = np.where(np.diff(power_freq_mean_ds)==2)[0] + 1\n\n                # We get the intervals 'minimum--outlier--minimum', which we will try to exclude as a whole\n                intervals = self._get_intervals(power_freq_mean_min_idx, power_freq_neg_mask[0][detector]).astype(int)\n                #print(detector, intervals)\n                # The outliers of other detector\n                other_detector_outliers = set(power_freq_neg_mask[1][self.other_detector[detector]])\n\n                exclusion_mask = set()\n                for left, right in intervals:\n                    # Check, if the interval is small\n                    if right-left-1 > max_interval:\n                        print(f'!!! The interval of [{left}; {right}] is to big to delete from {self.id_}/{self.id} record, detector={detector}')\n                        detectors_failed.append(detector)\n                        break\n                    \n                    # Notice `left+1`, which means we don't include the left minimum itself\n                    interval = set(range(left+1,right))\n                    \n                    # If the selected interval has no intersections with outliers of the other detector,\n                    # We add it to our exclusion mask\n                    if len(interval & other_detector_outliers) == 0:\n                        exclusion_mask = exclusion_mask | interval\n\n                power[detector] = np.delete(power[detector], list(exclusion_mask), axis=0)\n                print(f'detector={detector}')\n                print('intervals=', intervals)\n                print('other_detector_outliers', other_detector_outliers)\n                print('finally exclude=', exclusion_mask)\n        \n        if len(detectors_failed)>0:\n            detectors_failed = set(detectors_failed)\n            return detectors_failed\n        else:\n            self.power = power\n            return self.power\"\"\"","metadata":{}}]}