{"cells":[{"metadata":{"_uuid":"5829966057ac4dfb9536083174567bb31dfd7e05"},"cell_type":"markdown","source":"# Strategies for Flux Time Series Preprocessing\n### From simple to advanced methods"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport warnings\nfrom itertools import chain\nsns.set_style('whitegrid')\nwarnings.simplefilter('ignore', FutureWarning)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"One of the first major challenges in this competition is to properly preprocess the light curve time series data. Unlike many other \"well-behaved\" time series, the light curves here are not only irregular in terms of observation intervals, but also unsynchronised across different light bands. This makes applying many existing tools and techniques difficult because they often require a regular time series, preferably without missing values. Let us first look at what we have:"},{"metadata":{"trusted":true,"_uuid":"bab1c32e45589927cf5f2b5b8c6d22503013f0d5"},"cell_type":"code","source":"train_series = pd.read_csv('../input/training_set.csv')\ntrain_metadata = pd.read_csv('../input/training_set_metadata.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ea2cebe13bb8c9f55108371545c8f7fff8e7fa5d"},"cell_type":"code","source":"train_series.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"703d1e3fa4f3a25d41f962147b596b343d0bcc02"},"cell_type":"markdown","source":"There are actually 3 time series (flux, flux_err, detected) per band per object. However, as a first step, here we will only focus on the \"flux\" time series. Even so, we still have 6 passbands for each object."},{"metadata":{"_uuid":"ff0d32326a26ef32c414da2f130f99cbc199077b"},"cell_type":"markdown","source":"We first look at the distribution of time series lengths."},{"metadata":{"trusted":true,"_uuid":"5542d5016c2b7300ab25410280063dfc4955669a"},"cell_type":"code","source":"ts_lens = train_series.groupby(['object_id', 'passband']).size()\nf, ax = plt.subplots(figsize=(12, 6))\nsns.distplot(ts_lens, ax=ax)\nax.set_title('distribution of time series lengths')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a1b147f71e15b123ec560356283356f1af41f1c6"},"cell_type":"markdown","source":"As we see here, the lengths of light curve time series have a large range, so we may face difficulties trying to convert them all into the same length."},{"metadata":{"_uuid":"b23049ea27eb242726bafac07bd2e3ff9e452593"},"cell_type":"markdown","source":"Let us also look at the times at which observations happen, and count the number of observations made at each point in time:"},{"metadata":{"trusted":true,"_uuid":"d755a3eb284d3b309c0e0c00465fb737ee766516"},"cell_type":"code","source":"f, ax = plt.subplots(figsize=(12, 6))\nsns.distplot(train_series['mjd'], ax=ax, bins=200)\nax.set_title('number of observations made at each time point')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"57c377af9d6a59a78baa08141d1240b3582e612c"},"cell_type":"markdown","source":"What we see here shows that the observations are not taken evenly through the entire sampling period. There are periods where more samples than usual are taken, but also periods with few observations.\n\nAs for each individual object:"},{"metadata":{"trusted":true,"_uuid":"95e73d35dfc78193f25948712178ac6a36561720"},"cell_type":"code","source":"f, ax = plt.subplots(figsize=(12, 6))\nsns.distplot(train_series[train_series['object_id'] == 615]['mjd'], ax=ax, bins=200)\nax.set_title('number of observations made at each time point')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"47ae311530ca0162078d1a978da8ce5ed54ee1c4"},"cell_type":"markdown","source":"What about for a single light band?"},{"metadata":{"trusted":true,"_uuid":"cb41eabdaa29c879536d33e694480165d9dccabe"},"cell_type":"code","source":"f, ax = plt.subplots(figsize=(12, 6))\nsns.distplot(train_series[(train_series['object_id'] == 615) \n                          & (train_series['passband'] == 2)]['mjd'], ax=ax, bins=200)\nax.set_title('number of observations made at each time point')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6a57677fcf3eedb153247d01966f5350c5feb4e5"},"cell_type":"markdown","source":"We see that even for individual objects, observations are not taken evenly. There can be large gaps between periods of rich observations."},{"metadata":{"_uuid":"28a0eb78f56c5c49bef6dd2637f0991cf25de95e"},"cell_type":"markdown","source":"Let us also look at the number of observations for each object:"},{"metadata":{"trusted":true,"_uuid":"9dad88fcfbcae3bf4cb157d0a43e6884b79332ab"},"cell_type":"code","source":"obj_obs_count = train_series['object_id'].value_counts().reset_index()\nobj_obs_count.columns = ['object_id', 'count']\n\nobj_obs_count_w_ddf = pd.merge(\n    obj_obs_count, train_metadata[['object_id', 'ddf']], on='object_id')\n\nselected = obj_obs_count_w_ddf.groupby('ddf')['count'].value_counts()\nselected.index.names = ['ddf', 'count_val']\nselected = selected.reset_index().pivot('count_val', 'ddf',\n                                        'count').rename(columns={\n                                            0: 'nonddf',\n                                            1: 'ddf'\n                                        }).fillna(0)\nselected['total'] = selected['nonddf'] + selected['ddf']\n\nf, ax = plt.subplots(figsize=(12, 6))\nax.vlines(x=selected.index, ymin=0, ymax=selected['total'], \n          colors='red', label='ddf')\nax.vlines(x=selected.index, ymin=0, ymax=selected['nonddf'], \n          colors='blue', label='non-ddf')\nax.set_title('number of observations per object')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ed40c2ba5e0109f2b1a5f37a759a9063b6e29506"},"cell_type":"markdown","source":"As can be seen above, there is a large gap between the number of observations for an average DDF object and non-DDF object."},{"metadata":{"_uuid":"1307ca1ef66c2e3021f2f75a27727cee6128652f"},"cell_type":"markdown","source":"Judging by the observations we made so far, the processing of the light curve time series is indeed a challenging task. We have the following difficulties to overcome:\n1. The observations, even for one object in one band, are not evenly spaced, so adjacent observations in the time series are not equally \"close\" to each other.\n2. The length of the time series vary quite a lot, even among the colour bands of the same object. At no point do we have multiple colour band observations simultaneously, so we cannot naively treat the multiple light bands' time series as a multivariate series.\n3. There are some large gaps between periods where data is avaiable for an object, so we may miss characteristics of the time series at a certain period if we only look at global features.\n4. The time series for different objects are not synchronised, making comparisons between objects difficult. We will have to disocover features that are invariant with time shifting or slight stretching."},{"metadata":{"_uuid":"7171764829d0c1f2b1383e64ead2266e06890d89"},"cell_type":"markdown","source":"#### Ignore time values"},{"metadata":{"_uuid":"e63a385c3c33ad8a6f8e68b1dda1c2484f063bb2"},"cell_type":"markdown","source":"There are multiple strategies to approach this problem, with varying degrees of complexity.\n* We may completely ignore the time values and simply treat the time series as sequences. Many sequence features such as min, max, range, std, monotonicity etc. are invariant with the loss of time interval information, so we can still build many useful features. However, it comes with two major disadvantages:\n    1. We will not able to align observations that occur close together, so we will have to analyse the time series for each passband individually without exploiting their relationships.\n    2. We will lose access to many useful feature extraction techniques, or at least the form of these techniques that are typically implemented in common tools if we throw out time interval information. Features like autocorrelation do not make sense if we do not know the actual time gaps between observations.\n    \nPerforming time-less analysis is simple. We do not even need much preprocessing, which makes it a good starting point. In the following example, we extract four very basic features from the time-less time series:"},{"metadata":{"trusted":true,"_uuid":"083197d2abf306cc32006e47d444a8f5be0d7842"},"cell_type":"code","source":"simple_features = train_series.groupby(\n    ['object_id', 'passband'])['flux'].agg(\n    ['mean', 'max', 'min', 'std']).unstack('passband')\nsimple_features.head().T","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d0f40e7f07b5f50678a8adddfb99302e3fbc09f0"},"cell_type":"markdown","source":"Even if we cannot extract interaction features from multiple bands directly, we may still be able to discover their relationships from the features we have in downstream analysis and learning tasks.\n\n**What not to do here:** grab cesium or tsfresh, spin up a few cores and extract all the fancy features you can get. Be aware that many of the features in these libraries do not work well with unevenly-spaced time series. I wasted many CPU hours here myself :("},{"metadata":{"_uuid":"42445da7488bde4a1fabf3a7b6e6a2f0b8cde747"},"cell_type":"markdown","source":"However, we cannot simply ignore the rich information we just threw out. The time-less approach won't even allow us to discover periodicity! Plus, we know from some astronomy background that it is the composition of light at different frequencies (and by extension, bands) that allow us to identify objects in the first place. We really do not want to ignore the relationship between passbands. So we can come up with a different approach:"},{"metadata":{"_uuid":"f36e6f5cf66ecf37ace6bd5ec3687ff1afd68c83"},"cell_type":"markdown","source":"#### Convert to regular TS"},{"metadata":{"_uuid":"3a5d6c153543596b44951ef75ee1f0edbe4ebe1e"},"cell_type":"markdown","source":"* We may take all **mjd** values in the dataset and declare them 'valid observation timestamps', then construct a time series for each object and each passband that include all these observation timestamps, filling in NA if obseration is not available. Now that all objects and all passband have time series of the same length, they are comparable, and many time series analysis techniques can be applied as long as they can deal with missing values. One immediate drawback of this method is that none of the objects will have multiple passband data at a given timestamp, and there will be a huge number of NAs. The time series will also be really long due to the large number of possible timestamps. Also, we will face more issues when we want to build a model and apply it to new data, as the new timestamps will not match the timestamps in our current dataset. We can make a compromise by binning **mjd** and averaging the observations within a bin, and use all bins within the range of **mjd** as possible time steps in the time series. By doing so, we will be able to obtain a collection of time series that are equally spaced, not overwhemingly NA, and have the same length."},{"metadata":{"trusted":true,"_uuid":"ce10139e65ecb7a60686be29845bf2300feb83ba"},"cell_type":"code","source":"print(f'mjd unique values: {train_series[\"mjd\"].nunique()}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c2fd4f1362b75564beb538444b0d563e09cb4376"},"cell_type":"code","source":"print(f'int (1day) mjd unique values: {train_series[\"mjd\"].astype(int).nunique()}')\nprint(f'''int mjd bins: {train_series[\"mjd\"].astype(int).max()\n      - train_series[\"mjd\"].astype(int).min() + 1}''')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ec7a2e1409943b5d7ace338f214ba6698d94a3fd"},"cell_type":"code","source":"print(f'5day mjd unique values: {(train_series[\"mjd\"]/5).astype(int).nunique()}')\nprint(f'''5day mjd bins: {int((train_series[\"mjd\"].astype(int).max()\n      - train_series[\"mjd\"].astype(int).min()) / 5 + 1)}''')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5452661d6488b4b669a4066c19b575d8233eb4ba"},"cell_type":"markdown","source":"As we see above, while using each unique **mjd** value as time steps is really bad, binning by day or even by 5 days greatly reduces time series length. Here is an example of constructing the resulting time series using binned observations:"},{"metadata":{"trusted":true,"_uuid":"a8159c6777d0e4baa3d64c5f98ffd8bcc7b94cb5"},"cell_type":"code","source":"# binning\nts_mod = train_series[['object_id', 'mjd', 'passband', 'flux']].copy()\nts_mod['mjd_d5'] = (ts_mod['mjd'] / 5).astype(int)\nts_mod = ts_mod.groupby(['object_id', 'mjd_d5', 'passband'])['flux'].mean().reset_index()\nts_mod.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c4e6f3fc6d52f64839bdefd472b506d4b2a2c97e"},"cell_type":"code","source":"# pivotting\nts_piv = pd.pivot_table(ts_mod, \n                        index='object_id', \n                        columns=['mjd_d5', 'passband'], \n                        values='flux',\n                        dropna=False)\nts_piv.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9208e742a5ae59a02e7c10165b1413b744ba3a0b"},"cell_type":"code","source":"# resetting column index to fill mjd_d5 gaps \nt_min, t_max = ts_piv.columns.levels[0].min(), ts_piv.columns.levels[0].max()\nt_range = range(t_min, t_max + 1)\nmux = pd.MultiIndex.from_product([list(t_range), list(range(6))], \n                                 names=['mjd_d5', 'passband'])\nts_piv = ts_piv.reindex(columns=mux).stack('passband')\nts_piv.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ef9a4d95b9d173a87a0f1d387ee11d5f5242c4c2"},"cell_type":"markdown","source":"We have now converted the unaligned time series into aligned, equally spaced time series using binning-and-averaging. We can see below that despite our efforts, nearly 90% of the values in our new time series are NaN, so there will be major challenges trying to make use of the converted data with tools that cannot cope well with NAs. Normal procedures like mean imputation or interpolation won't work well here due to the prevalance of missing values. However, some algorithms like LSTM may learn to ignore padded values and data gaps and still manage to extract valuable information from the time series. It is even possible to build multiple time series collections with different bins to look for features at different time scales."},{"metadata":{"trusted":true,"_uuid":"4420f30c9098c92196988e99ce7820e8ed7422fc"},"cell_type":"code","source":"np.mean(np.ravel(pd.isna(ts_piv).values))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c1925a410c655d9e9d956a84a330e5456ab4ef00"},"cell_type":"markdown","source":"One of the drawbacks of the method above is that binning appears to be quite arbitrary. What we are effectively doing is sampling at different time points and averaging flux values around these points. So why not do it properly?"},{"metadata":{"_uuid":"73cfeecc39ffff7d31c582ba99f8a0b9877e9995"},"cell_type":"markdown","source":"#### Sample with a distance kernel"},{"metadata":{"_uuid":"06ce7f495cec36066fe17f02af5dfb4d8425cf5f"},"cell_type":"markdown","source":"* We may sample the time series data at evenly-spaced time steps, using a time kernel to determine how much weight we put on each observation when taking the moving average."},{"metadata":{"trusted":true,"_uuid":"cbdc51e621fc882e83d659d75ef7966b69a89ba4"},"cell_type":"code","source":"def time_kernel(diff, tau):\n    return np.exp(-diff ** 2 / (2 * tau ** 2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5e138e18129b61392f8543644d01267aebaa1f29"},"cell_type":"code","source":"t_min, t_max = train_series['mjd'].min(), train_series['mjd'].max()\nt_min, t_max","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3c4dadd6df784520475913b150c2ca85443d3359"},"cell_type":"code","source":"sample_points = np.array(np.arange(t_min, t_max, 20))\nsample_points","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4ef9bfefba08e3d169b60adc36b780b54e99ef67"},"cell_type":"code","source":"weights = time_kernel(np.expand_dims(sample_points, 0) \n                      - np.expand_dims(train_series['mjd'].values, 1), 5)\nts_mod = train_series[['object_id', 'mjd', 'passband', 'flux']].copy()\nfor i in range(len(sample_points)):\n    ts_mod[f'sw_{i}'] = weights[:, i]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"59ede54c0cbf62ebf6b3ea360c465fab5f2d28ed"},"cell_type":"code","source":"def group_transform(chunk):\n    sample_weights = chunk[[f'sw_{i}' for i in range(len(sample_points))]]\n    sample_weights /= np.sum(sample_weights, axis=0)\n    weighted_flux = np.expand_dims(chunk['flux'].values, 1) * sample_weights.fillna(0)\n    return np.sum(weighted_flux, axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"62c85ebb3260823c9a68b31e0e876734d10c6587"},"cell_type":"code","source":"# only using a small sample as this step is slower than the other steps\nts_samp = ts_mod[ts_mod['object_id'].isin([615, 713])].groupby(\n    ['object_id', 'passband']).apply(group_transform)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3a5ac08dd672c8f5bb284b552b27416afc3382ad"},"cell_type":"code","source":"ts_samp","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a1994a7c8631a4dede992012e7e2074439761256"},"cell_type":"markdown","source":"We can see above that we have managed to sample from the time series, using the time difference from the sample point to determine moving average weights.\n\nThere are still many issues with this method. For example, it still does not solve the problem of long gaps between available data, and the use of 0 as a filler value for NA is questionable (however we cannot retain NA, as that will make calculating moving average difficult). Overall, this is still a better approach than simple binning, and still relatively cheap to compute."},{"metadata":{"_uuid":"8885d8648606d0b5f761e5aecde38f8dd830573e"},"cell_type":"markdown","source":"#### Use time / time difference as a feature"},{"metadata":{"_uuid":"00c0025a1f59c77b1841a3b49489205e44c94f6d"},"cell_type":"markdown","source":"This idea is quite simple. Why not just concatenate the flux values with the timestamp at which they are observed, or the time elapsed since last observation, and pass the whole sequence to a learning algorithm? Indeed, some algorithms have been shown to be able to cope with that. This is not so different from \"given (x, y), fit the curve and extract features from the curve\". An LSTM or even a regular MLP might be able to deal with that. There is an [example](https://arxiv.org/pdf/1711.10609.pdf) of this being applied to a very similar problem.\n\nAlso, some time series feature extractors accept (t, x) pairs instead of (x) sequences.\n\n(WIP)"},{"metadata":{"_uuid":"963bda5a526e6cb7fd95d9ae112212b5df8fd411"},"cell_type":"markdown","source":"#### Convert time series to phase series"},{"metadata":{"_uuid":"e578121e488648d094c49d5b6dd3665858744223"},"cell_type":"markdown","source":"If we know that some of the time series will be periodic, we can attempt to find the optimal period and convert the time series to phase domain, so that we can avoid the issues of observation gap, alleviate the problem of sparse observations and make better use of the periodic property. We may also assume that for most objects that exhibit period behaviour, the period is the same for all colour bands and the phases should all be in sync, because they are likely driven by the same physical event.\nTo discover the period of the time series, we first normalise all series then extract features about the most prominent frequencies / periods."},{"metadata":{"trusted":true,"_uuid":"a79c416e18b9cf27707b1189c5ef9bc3c4a72793"},"cell_type":"code","source":"groups = train_series.groupby(['object_id', 'passband'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"519d8d02b498769b72c742138edcdf02fcce0a19"},"cell_type":"code","source":"def normalise(ts):\n    return (ts - ts.mean()) / ts.std()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"59c625e9457e442d8d1b384944cff573fecef4a6"},"cell_type":"code","source":"times = groups.apply(\n    lambda block: block['mjd'].values).reset_index().rename(columns={0: 'seq'})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"66a0d8a0dcb14d7010140e8eb6e11c624ca96b70"},"cell_type":"code","source":"flux = groups.apply(\n    lambda block: normalise(block['flux']).values\n).reset_index().rename(columns={0: 'seq'})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cd32fe6d05f7f5c8ef4b87c1b99ce2c067e96126"},"cell_type":"code","source":"times_list = times.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()\nflux_list = flux.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a79b1049a6f7d108fa3f7734122e45521c626657"},"cell_type":"code","source":"import cesium.featurize as featurize\nfrom scipy import signal\nimport warnings","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"43ded3de908cef55df1fa8627e637d75dafb3e14"},"cell_type":"markdown","source":"To save some kernel time, we will only extract features from a small subset of objects.\n\n([**FIXED - See below**]For some reason, the frequencies features of the cesium package are not working properly for me, so I have to revert to the method used by Michal Haltuf in his notebook [Feature extraction using period analysis](https://www.kaggle.com/rejpalcz/feature-extraction-using-period-analysis). It takes quite some time even for just 20 objects as you can see below:\n"},{"metadata":{"trusted":true,"_uuid":"9b60234c7b6ae2a37136bdd27565eb3b9d689575"},"cell_type":"code","source":"# def extract_freq(t, m, e):\n#     fs = np.linspace(2*np.pi/0.1, 2*np.pi/500, 10000)\n#     pgram = signal.lombscargle(t, m, fs, normalize=True)\n#     return fs[np.argmax(pgram)]\n\n# N = 20\n# warnings.simplefilter('ignore', RuntimeWarning)\n# feats = featurize.featurize_time_series(times=times_list[:N],\n#                                         values=flux_list[:N],\n#                                         features_to_use=['freq1'],\n#                                         custom_functions={'freq1': extract_freq},\n#                                         scheduler=None)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1af560d13b3fafe77308c5507af9e1c2448ba7fe"},"cell_type":"markdown","source":"**EDIT:** I found why the features from the cesium package is acting weird. The frequency features extracted by cesium are oscilation frequencies (cycles per unit time) rather than angular frequencies (rad per unit time), so they differ by $2\\pi$. We can now use cesium's 'freqN_freq' feature extractor to extract features slightly faster:"},{"metadata":{"trusted":true,"_uuid":"bf6d8999ffcd12704995fa0a24d8727840376b2b"},"cell_type":"code","source":"warnings.simplefilter('ignore', RuntimeWarning)\nN = 100\ncfeats = featurize.featurize_time_series(times=times_list[:N],\n                                        values=flux_list[:N],\n                                        features_to_use=['freq1_freq',\n                                                        'freq1_signif',\n                                                        'freq1_amplitude1'],\n                                        scheduler=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"463bc1851e3d81253d1c0afbac7f36caa34fbaba"},"cell_type":"code","source":"cfeats.stack('channel').iloc[:24]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7a536ac3cbb5f2ecd66a029cfeacd2bac69b3c7d"},"cell_type":"markdown","source":"We can see that for some objects, there appears to be a common frequency / period for all its observed passbands, whereas for others the calculated frequencies are all different and seemingly random. We may assume that the objects with completely mismatching frequencies or very low frequencies (very long periods) actually do not exhibit periodic behaviour. Let us look at a few examples where the periods do match up between most bands:"},{"metadata":{"trusted":true,"_uuid":"b224e1f1ee83b215233c5d4870c66fb8f63852c5"},"cell_type":"code","source":"def plot_phase(n, fr):\n    selected_times = times_list[n]\n    selected_flux = flux_list[n]\n    colors = ['red', 'orange', 'yellow', 'green', 'blue', 'purple']\n    f, ax = plt.subplots(figsize=(12, 6))\n    for band in range(6):\n        ax.scatter(x=(selected_times[band] * fr) % 1, \n                   y=selected_flux[band], \n                   c=colors[band])\n    ax.set_xlabel('phase')\n    ax.set_ylabel('relative flux')\n    ax.set_title(\n        f'object {train_metadata[\"object_id\"][n]}, class {train_metadata[\"target\"][n]}')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8faff52441f8d6e577fedcc07f22c77874b14ee9","scrolled":true},"cell_type":"code","source":"plot_phase(0, 3.081631)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"34533ba26c7f6c25d5745d3c74413ee221ebff22"},"cell_type":"code","source":"plot_phase(3, 0.001921)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fa31d2508163bdf4ec562b278e114befed8bc0d7"},"cell_type":"code","source":"plot_phase(6, 1.005547)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e04e3d6e5365b1a30d3aa9f6c3f7728df45ef066"},"cell_type":"markdown","source":"For periodic objects, observations made in the same phase should have similar relative flux values. As we see, Object 615 is clearly periodically changing its magnitude,  whereas object 745 and 1598's 'periodicity' are more like staying stationary most of the time with sudden bursts.\n\nWe kind of found a pattern here. If all colour bands agree in frequency, an object is likely to have periodic behaviours. When the colour bands have wildly different frequencies, the object is likely aperiodic. When most bands have very small frequencies (longer period than the observation window) and a few bands have larger detected frequencies, it is likely a burst event.\n\nIt appears that calculated periods alone cannot tell us the full picture, as there are plenty of objects in the dataset without clear periodic patterns. However, determining whether an object has a periodic pattern can be very important in determining the object class, and for objects that are actually periodic, the time series in phase space can tell us a lot more than in time space. Therefore, for a large portion of the dataset, conversion to phase space is actually a very powerful preprocessing step to use.\n\nWe will plot below a few more examples of where phase conversion works well:"},{"metadata":{"trusted":true,"_uuid":"24f62f7b33ec23b74856a01f898c431d590ddb67"},"cell_type":"code","source":"plot_phase(46, 1.781711)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c0928a08a845fbde0044a2d37941d7ce91ff1ffa"},"cell_type":"code","source":"plot_phase(62, 4.771011)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5701d06e34e319fd4df16eea33688262345173df"},"cell_type":"code","source":"plot_phase(69, 1.597659)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d049d51c7ca33a6344b177a0238badf6518ee4b9"},"cell_type":"code","source":"plot_phase(91, 2.668811)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1eeb02000a9d895f46ee5b185f8c6c2d99d867dc"},"cell_type":"markdown","source":"It appears that we are onto something. The classes 16 and 92 appear to be quite periodic where class 90 is what appears to be a group of burst events. So by just calculating a few frequency features and comparing these features between bands, we can already roughly tell them apart. We might still have trouble distinguishing between class 16 and 92, but with frequency divided out through phase-space conversion, we essentially have higher density data for these two classes to work with, and building a model on top of that should be easier."},{"metadata":{"_uuid":"95a3b31f2945ff672b1093effb899d568ff2cd52"},"cell_type":"markdown","source":"One more thing... We have to figure out if the same procedure can work well on non-DDF objects."},{"metadata":{"trusted":true,"_uuid":"9b69e1a9eddce7ae7d7b13d995b6dea61336a3d1"},"cell_type":"code","source":"nonddf_pos = train_metadata[train_metadata['ddf'] == 0].index\nnonddf_times_list = [v for i, v in enumerate(times_list) if i in set(nonddf_pos)]\nnonddf_flux_list = [v for i, v in enumerate(flux_list) if i in set(nonddf_pos)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5233e394101f0b816c1dd1cfbb63d9d6414992c6"},"cell_type":"code","source":"warnings.simplefilter('ignore', RuntimeWarning)\nN = 50\ncfeats = featurize.featurize_time_series(times=nonddf_times_list[:N],\n                                        values=nonddf_flux_list[:N],\n                                        features_to_use=['freq1_freq',\n                                                        'freq1_signif',\n                                                        'freq1_amplitude1'],\n                                        scheduler=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1b4a26d3d9dabcf0362c87a3f264f132032a583b"},"cell_type":"code","source":"cfeats.stack('channel').iloc[:24]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5f7f81ecf14cc2435ce366087c30ba9d40e3271e"},"cell_type":"code","source":"plot_phase(nonddf_pos[15], 4.440506)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dcaadbf6821a08d02007a6605df71cb6c70707d3"},"cell_type":"code","source":"plot_phase(nonddf_pos[20], 0.831637)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"865706fe5ce8d2d99314dc062e9e51144d2356e8"},"cell_type":"code","source":"plot_phase(nonddf_pos[39], 0.817696)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5a9ac3c6f50e4a854ea79270e68ecd47386e7e74"},"cell_type":"code","source":"plot_phase(nonddf_pos[46], 2.514430)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6462ca4672dd99ce81125d7bb8558eb67a2b2f7c"},"cell_type":"markdown","source":"Generally, we found that we are still able to use frequency analysis to find the most apparent periodic classes like 16 and 92, but we have to be much more lenient on how different bands agree with each other (sometimes only 3-4 bands will strongly agree in terms of frequency), and we might have to pay much more attention to the significance level of these frequencies. This could be evidence to support the case that DDF and non-DDF objects should be dealt with using different models, or at least treated differently in model."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}