{"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":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-10-14T09:12:33.011414Z","iopub.execute_input":"2022-10-14T09:12:33.011978Z","iopub.status.idle":"2022-10-14T09:12:39.260794Z","shell.execute_reply.started":"2022-10-14T09:12:33.011865Z","shell.execute_reply":"2022-10-14T09:12:39.259753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"! pip install giotto-tda","metadata":{"execution":{"iopub.status.busy":"2022-10-14T10:20:53.075704Z","iopub.execute_input":"2022-10-14T10:20:53.076145Z","iopub.status.idle":"2022-10-14T10:21:08.395538Z","shell.execute_reply.started":"2022-10-14T10:20:53.076112Z","shell.execute_reply":"2022-10-14T10:21:08.393721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Generating gravitational-wave signals #\nthis code was not adapted for the competition\n","metadata":{}},{"cell_type":"markdown","source":"## Giotto-Tda ## \nGiotto-Tda is a high performance topological machine learning toolbox in Python built on top of scikit-learn and is distributed under the GNU AGPLv3 license. It is part of the Giotto family of open-source projects.\n\ngiotto-tda: A Topological Data Analysis Toolkit for Machine Learning and Data Exploration, Tauzin et al, J. Mach. Learn. Res. 22.39 (2021): 1-6.\n","metadata":{}},{"cell_type":"markdown","source":"## Giotto-TDA example ##\nThis a light customisation of Giotto-TDA prebuild example.","metadata":{}},{"cell_type":"code","source":"from IPython.display import YouTubeVideo\n\nYouTubeVideo(\"Y3eR49ogsF0\", width=600, height=400)","metadata":{"execution":{"iopub.status.busy":"2022-10-14T09:23:28.491598Z","iopub.execute_input":"2022-10-14T09:23:28.492579Z","iopub.status.idle":"2022-10-14T09:23:28.668286Z","shell.execute_reply.started":"2022-10-14T09:23:28.492541Z","shell.execute_reply":"2022-10-14T09:23:28.66723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The Sound of Two Black Holes Colliding ##","metadata":{}},{"cell_type":"markdown","source":"Gravitational waves sent out from a pair of colliding black holes have been converted to sound waves, as heard in this animation. On September 14, 2015, LIGO observed gravitational waves from the merger of two black holes, each about 30 times the mass of our sun. The incredibly powerful event, which released 50 times more energy than all the stars in the observable universe, lasted only fractions of a second.\nIn the first two runs of the animation, the sound-wave frequencies exactly match the frequencies of the gravitational waves. The second two runs of the animation play the sounds again at higher frequencies that better fit the human hearing range. The animation ends by playing the original frequencies again twice.\nAs the black holes spiral closer and closer in together, the frequency of the gravitational waves increases. Scientists call these sounds \"chirps,\" because some events that generate gravitation waves would sound like a bird's chirp.\nImage credit: LIGO [Source Wikipedia]","metadata":{}},{"cell_type":"code","source":"YouTubeVideo(\"QyDcTbR-kEA\", width=600, height=400)","metadata":{"execution":{"iopub.status.busy":"2022-10-14T09:24:46.75022Z","iopub.execute_input":"2022-10-14T09:24:46.751689Z","iopub.status.idle":"2022-10-14T09:24:46.85357Z","shell.execute_reply.started":"2022-10-14T09:24:46.751634Z","shell.execute_reply":"2022-10-14T09:24:46.852353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Generate the data ###\n\nChrisopher Bresten and Jae-Hun Jung have proposed a synthetic training set as follows:\n\nGenerate gravitational wave signals that correspond to __non-spinning binary black hole__ mergers\nGenerate a noisy time series and embed a gravitational wave signal with probability 0.5 at a random time.\nThe result is a set of time series of the form\n\n$$ s = g + \\epsilon \\frac{1}{R}\\xi $$\nwhere g is a gravitational wave signal from the reference set,\n$ \\epsilon $ is Gaussian noise,\n$$ \\epsilon=10^{-19} $$\nscales the noise amplitude to the signal,\nand\n$ R \\in (0.075, 0.65)$\nis a parameter that controls the signal-to-noise-ratio (SNR).","metadata":{}},{"cell_type":"markdown","source":"### Constant signal-to-noise ratio###","metadata":{}},{"cell_type":"markdown","source":"make gravitational waves function","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\ndef make_gravitational_waves(\n    path_to_data: Path,\n    n_signals: int = 30,\n    downsample_factor: int = 2,\n    r_min: float = 0.075,\n    r_max: float = 0.65,\n    n_snr_values: int = 10,\n        ):\n    def padrand(V, n, kr):\n        cut = np.random.randint(n)\n        rand1 = np.random.randn(cut)\n        rand2 = np.random.randn(n - cut)\n        out = np.concatenate((rand1 * kr, V, rand2 * kr))\n        return out\n\n    Rcoef = np.linspace(r_min, r_max, n_snr_values)\n    Npad = 500  # number of padding points on either side of the vector\n    #gw = np.load(path_to_data / \"gravitational_wave_signals.npy\")\n    #../input/gravitational-wave-signals/gravitational_wave_signals.npy\n    gw = np.load(\"../input/gravitational-wave-signals/gravitational_wave_signals.npy\")\n    Norig = len(gw[\"data\"][0])\n    Ndat = len(gw[\"signal_present\"])\n    N = int(Norig / downsample_factor)\n\n    ncoeff = []\n    Rcoeflist = []\n\n    for j in range(n_signals):\n        ncoeff.append(10 ** (-19) * (1 / Rcoef[j % n_snr_values]))\n        Rcoeflist.append(Rcoef[j % n_snr_values])\n\n    noisy_signals = []\n    gw_signals = []\n    k = 0\n    labels = np.zeros(n_signals)\n\n    for j in range(n_signals):\n        signal = gw[\"data\"][j % Ndat][range(0, Norig, downsample_factor)]\n        sigp = int((np.random.randn() < 0))\n        noise = ncoeff[j] * np.random.randn(N)\n        labels[j] = sigp\n        if sigp == 1:\n            rawsig = padrand(signal + noise, Npad, ncoeff[j])\n            if k == 0:\n                k = 1\n        else:\n            rawsig = padrand(noise, Npad, ncoeff[j])\n        noisy_signals.append(rawsig.copy())\n        gw_signals.append(signal)\n\n    return noisy_signals, gw_signals, labels","metadata":{"execution":{"iopub.status.busy":"2022-10-14T10:49:16.57483Z","iopub.execute_input":"2022-10-14T10:49:16.575342Z","iopub.status.idle":"2022-10-14T10:49:16.59586Z","shell.execute_reply.started":"2022-10-14T10:49:16.575302Z","shell.execute_reply":"2022-10-14T10:49:16.59472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#from data.generate_datasets import make_gravitational_waves\nfrom pathlib import Path\n\nR = 0.65\nn_signals = 100\nDATA = Path(\"./data\")\n\nnoisy_signals, gw_signals, labels = make_gravitational_waves(\n    path_to_data=DATA, n_signals=n_signals, r_min=R, r_max=R, n_snr_values=1\n)\n\nprint(f\"Number of noisy signals: {len(noisy_signals)}\")\nprint(f\"Number of timesteps per series: {len(noisy_signals[0])}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-14T10:49:20.528554Z","iopub.execute_input":"2022-10-14T10:49:20.529533Z","iopub.status.idle":"2022-10-14T10:49:20.700697Z","shell.execute_reply.started":"2022-10-14T10:49:20.529499Z","shell.execute_reply":"2022-10-14T10:49:20.699549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Two different types of time series visualisation: \n* pure noise and: \n* noise plus an embedded gravitational wave signal:","metadata":{}},{"cell_type":"code","source":"import numpy as np\nfrom plotly.subplots import make_subplots\nimport plotly.graph_objects as go\n\n# get the index corresponding to the first pure noise time series\nbackground_idx = np.argmin(labels)\n# get the index corresponding to the first noise + gravitational wave time series\nsignal_idx = np.argmax(labels)\n\nts_noise = noisy_signals[background_idx]\nts_background = noisy_signals[signal_idx]\nts_signal = gw_signals[signal_idx]\n\nfig = make_subplots(rows=1, cols=2)\n\nfig.add_trace(\n    go.Scatter(x=list(range(len(ts_noise))), y=ts_noise, mode=\"lines\", name=\"noise\"),\n    row=1,\n    col=1,\n)\n\nfig.add_trace(\n    go.Scatter(\n        x=list(range(len(ts_background))),\n        y=ts_background,\n        mode=\"lines\",\n        name=\"background\",\n    ),\n    row=1,\n    col=2,\n)\n\nfig.add_trace(\n    go.Scatter(x=list(range(len(ts_signal))), y=ts_signal, mode=\"lines\", name=\"signal\"),\n    row=1,\n    col=2,\n)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-14T10:50:06.70008Z","iopub.execute_input":"2022-10-14T10:50:06.700528Z","iopub.status.idle":"2022-10-14T10:50:07.247431Z","shell.execute_reply.started":"2022-10-14T10:50:06.700496Z","shell.execute_reply":"2022-10-14T10:50:07.245892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"https://giotto-ai.github.io/gtda-docs/0.3.1/notebooks/gravitational_waves_detection.html","metadata":{}},{"cell_type":"markdown","source":"1. It is hard to distinguish the signal by eye, \n2. The signal features some regularity or periodicity.\n\nBoth observations lead us to examining the Takens embedding of the signal \n, in order to pick up the recurrent structure. \nIndeed, if \nis sampled from a dynamical system with a non-trivial recurrent structure, then, for appropriate parameters, the image by the embedding will have non-trivial topology.","metadata":{}},{"cell_type":"markdown","source":"More formally,, we extract a sequence of vectors in $\\mathbb{R}^{d}$ of the form\n\n$$\\begin{split}TD_{d,\\tau} s : \\mathbb{R} \\to \\mathbb{R}^{d}\\,, \\qquad t \\to \\begin{bmatrix}\n           s(t) \\\\\n           s(t + \\tau) \\\\\n           s(t + 2\\tau) \\\\\n           \\vdots \\\\\n           s(t + (d-1)\\tau)\n         \\end{bmatrix},\\end{split} $$\n \n \nwhere d\n is the embedding dimension and $\\tau$\n is the time delay. The quantity $(d-1)\\tau$\n is known as the “window size” and the difference between $t_{i+1}$\n and \n is called the stride.\n\nLet’s examine what the time delay embedding of a pure gravitational wave signal looks like:","metadata":{}},{"cell_type":"code","source":"from gtda.time_series import SingleTakensEmbedding\nembedding_dimension = 30\nembedding_time_delay = 30\nstride = 5\n\nembedder = SingleTakensEmbedding(\n    parameters_type=\"search\", n_jobs=6, time_delay=embedding_time_delay, dimension=embedding_dimension, stride=stride\n)\n\ny_gw_embedded = embedder.fit_transform(gw_signals[0])","metadata":{"execution":{"iopub.status.busy":"2022-10-14T11:22:04.68328Z","iopub.execute_input":"2022-10-14T11:22:04.683775Z","iopub.status.idle":"2022-10-14T11:22:08.134804Z","shell.execute_reply.started":"2022-10-14T11:22:04.683741Z","shell.execute_reply":"2022-10-14T11:22:08.133314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.decomposition import PCA\nfrom gtda.plotting import plot_point_cloud\n\npca = PCA(n_components=3)\ny_gw_embedded_pca = pca.fit_transform(y_gw_embedded)\n\nplot_point_cloud(y_gw_embedded_pca)","metadata":{"execution":{"iopub.status.busy":"2022-10-14T11:22:34.908044Z","iopub.execute_input":"2022-10-14T11:22:34.908808Z","iopub.status.idle":"2022-10-14T11:22:35.202782Z","shell.execute_reply.started":"2022-10-14T11:22:34.908747Z","shell.execute_reply":"2022-10-14T11:22:35.201818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From the plot we can see that the decaying periodic signal generated by a black hole merger emerges as a spiral in the time delay embedding space! For contrast, let’s compare this to one of the pure noise time series in our sample:","metadata":{}},{"cell_type":"code","source":"embedding_dimension = 30\nembedding_time_delay = 30\nstride = 5\n\nembedder = SingleTakensEmbedding(\n    parameters_type=\"search\", n_jobs=6, time_delay=embedding_time_delay, dimension=embedding_dimension, stride=stride\n)\n\ny_noise_embedded = embedder.fit_transform(noisy_signals[background_idx])\n\npca = PCA(n_components=3)\ny_noise_embedded_pca = pca.fit_transform(y_noise_embedded)\n\nplot_point_cloud(y_noise_embedded_pca)","metadata":{"execution":{"iopub.status.busy":"2022-10-14T11:23:47.052168Z","iopub.execute_input":"2022-10-14T11:23:47.052672Z","iopub.status.idle":"2022-10-14T11:23:47.868315Z","shell.execute_reply.started":"2022-10-14T11:23:47.052602Z","shell.execute_reply":"2022-10-14T11:23:47.866893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Evidently, pure noise resembles a high-dimensional ball in the time delay embedding space. Let’s see if we can use persistent homology to tease apart which time series contain a gravitational wave signal versus those that don’t. To do so we will adapt the strategy from the original article:\n\nGenerate 200-dimensional time delay embeddings of each time series\nUse PCA to reduce the time delay embeddings to 3-dimensions\nUse the Vietoris-Rips construction to calculate persistence diagrams of \n and \n generators\nExtract feature vectors using persistence entropy\nTrain a binary classifier on the topological features","metadata":{}},{"cell_type":"markdown","source":"## Define the topological feature generation pipeline\n\nWe can do steps 1 and 2 by using the following giotto-tda tools:\n\nThe TakensEmbedding transformer – instead of SingleTakensEmbedding – which will transform each time series in noisy_signals separately and return a collection of point clouds;\nCollectionTransformer, which is a convenience “meta-estimator” for applying the same PCA to each point cloud resulting from step 1.\nUsing the Pipeline class from giotto-tda, we can chain all operations up to and including step 4 as follows:","metadata":{}},{"cell_type":"code","source":"from gtda.diagrams import PersistenceEntropy, Scaler\nfrom gtda.homology import VietorisRipsPersistence\nfrom gtda.metaestimators import CollectionTransformer\nfrom gtda.pipeline import Pipeline\nfrom gtda.time_series import TakensEmbedding\n\nembedding_dimension = 200\nembedding_time_delay = 10\nstride = 10\n\nembedder = TakensEmbedding(time_delay=embedding_time_delay,\n                           dimension=embedding_dimension,\n                           stride=stride)\n\nbatch_pca = CollectionTransformer(PCA(n_components=3), n_jobs=-1)\n\npersistence = VietorisRipsPersistence(homology_dimensions=[0, 1], n_jobs=-1)\n\nscaling = Scaler()\n\nentropy = PersistenceEntropy(normalize=True, nan_fill_value=-10)\n\n\nsteps = [(\"embedder\", embedder),\n         (\"pca\", batch_pca),\n         (\"persistence\", persistence),\n         (\"scaling\", scaling),\n         (\"entropy\", entropy)]\ntopological_transfomer = Pipeline(steps)","metadata":{"execution":{"iopub.status.busy":"2022-10-14T11:25:39.295029Z","iopub.execute_input":"2022-10-14T11:25:39.296688Z","iopub.status.idle":"2022-10-14T11:25:39.344727Z","shell.execute_reply.started":"2022-10-14T11:25:39.296614Z","shell.execute_reply":"2022-10-14T11:25:39.343478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features = topological_transfomer.fit_transform(noisy_signals)","metadata":{"execution":{"iopub.status.busy":"2022-10-14T11:25:58.696123Z","iopub.execute_input":"2022-10-14T11:25:58.696581Z","iopub.status.idle":"2022-10-14T11:26:14.214851Z","shell.execute_reply.started":"2022-10-14T11:25:58.696545Z","shell.execute_reply":"2022-10-14T11:26:14.213458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train and evaluate a model \n\nFor the final step, let’s train a simple classifier on our topological features. As usual we create training and validation sets","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX_train, X_valid, y_train, y_valid = train_test_split(\n    features, labels, test_size=0.1, random_state=42\n)","metadata":{"execution":{"iopub.status.busy":"2022-10-14T11:27:14.171583Z","iopub.execute_input":"2022-10-14T11:27:14.1721Z","iopub.status.idle":"2022-10-14T11:27:14.180882Z","shell.execute_reply.started":"2022-10-14T11:27:14.172056Z","shell.execute_reply":"2022-10-14T11:27:14.17932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.metrics import accuracy_score, roc_auc_score\n\n\ndef print_scores(fitted_model):\n    res = {\n        \"Accuracy on train:\": accuracy_score(fitted_model.predict(X_train), y_train),\n        \"ROC AUC on train:\": roc_auc_score(\n            y_train, fitted_model.predict_proba(X_train)[:, 1]\n        ),\n        \"Accuracy on valid:\": accuracy_score(fitted_model.predict(X_valid), y_valid),\n        \"ROC AUC on valid:\": roc_auc_score(\n            y_valid, fitted_model.predict_proba(X_valid)[:, 1]\n        ),\n    }\n\n    for k, v in res.items():\n        print(k, round(v, 3))","metadata":{"execution":{"iopub.status.busy":"2022-10-14T11:27:37.916069Z","iopub.execute_input":"2022-10-14T11:27:37.916518Z","iopub.status.idle":"2022-10-14T11:27:37.925121Z","shell.execute_reply.started":"2022-10-14T11:27:37.916483Z","shell.execute_reply":"2022-10-14T11:27:37.923862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"from sklearn.linear_model import LogisticRegression\n\nmodel = LogisticRegression()\nmodel.fit(X_train, y_train)\nprint_scores(model)","metadata":{"execution":{"iopub.status.busy":"2022-10-14T11:27:59.727309Z","iopub.execute_input":"2022-10-14T11:27:59.727846Z","iopub.status.idle":"2022-10-14T11:27:59.749057Z","shell.execute_reply.started":"2022-10-14T11:27:59.727795Z","shell.execute_reply":"2022-10-14T11:27:59.747382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"https://www.kaggle.com/code/mistag/reverse-engineering-create-clean-gw-signals/notebook","metadata":{}}]}