{"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":"import os\nimport sys\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\npd.options.display.max_columns = 100\n\n#from skimage.filters import difference_of_gaussians\nfrom sklearn.model_selection import StratifiedKFold, GroupKFold\nfrom sklearn.metrics import f1_score\nimport random\nimport time\n\n# import torch\n# import torch.nn as nn\n# import torch.nn.functional as F\n# from torch.utils.data import DataLoader, Dataset, WeightedRandomSampler\n# from torch.cuda.amp import autocast, GradScaler\n\nimport librosa\nimport librosa.display\n\nfrom scipy.special import logit, expit\n\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import Normalize\n%matplotlib inline\n\nimport cv2\nfrom scipy.interpolate import interp1d\nimport pywt\n\ndef sigmoid(var):\n    return 1/(1+np.exp(-var))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-09-06T16:17:02.098897Z","iopub.execute_input":"2021-09-06T16:17:02.099484Z","iopub.status.idle":"2021-09-06T16:17:03.223328Z","shell.execute_reply.started":"2021-09-06T16:17:02.099445Z","shell.execute_reply":"2021-09-06T16:17:03.221687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install pycbc\nimport pycbc\nfrom pycbc.waveform import td_approximants, fd_approximants, get_td_waveform\nfrom pycbc.detector import Detector","metadata":{"execution":{"iopub.status.busy":"2021-09-06T16:17:03.225482Z","iopub.execute_input":"2021-09-06T16:17:03.225883Z","iopub.status.idle":"2021-09-06T16:17:05.401988Z","shell.execute_reply.started":"2021-09-06T16:17:03.225840Z","shell.execute_reply":"2021-09-06T16:17:05.400536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# List of td approximants that are available\nprint(td_approximants())","metadata":{"execution":{"iopub.status.busy":"2021-09-06T16:17:05.404302Z","iopub.execute_input":"2021-09-06T16:17:05.404726Z","iopub.status.idle":"2021-09-06T16:17:05.410842Z","shell.execute_reply.started":"2021-09-06T16:17:05.404682Z","shell.execute_reply":"2021-09-06T16:17:05.409532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# List of fd approximants that are currently available\nprint(fd_approximants())","metadata":{"execution":{"iopub.status.busy":"2021-09-06T16:17:05.412553Z","iopub.execute_input":"2021-09-06T16:17:05.412992Z","iopub.status.idle":"2021-09-06T16:17:05.428580Z","shell.execute_reply.started":"2021-09-06T16:17:05.412953Z","shell.execute_reply":"2021-09-06T16:17:05.427234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gwlist = ['SEOBNRv2', 'SEOBNRv2_opt', 'SEOBNRv4', 'SEOBNRv4_opt', 'SEOBNRv2T', 'SEOBNRv4T', ]\ngwlist = ['SEOBNRv2', 'SEOBNRv4' ]","metadata":{"execution":{"iopub.status.busy":"2021-09-06T16:17:05.430751Z","iopub.execute_input":"2021-09-06T16:17:05.431193Z","iopub.status.idle":"2021-09-06T16:17:05.439219Z","shell.execute_reply.started":"2021-09-06T16:17:05.431151Z","shell.execute_reply":"2021-09-06T16:17:05.438154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.random.seed(1)\n\n#Define the Detectors\ndet_h1 = Detector('H1')\ndet_l1 = Detector('L1')\ndet_v1 = Detector('V1')\n\nfor n in range(100):\n    #Define the GW params\n    gwapprox = np.random.choice( gwlist )\n    print(gwapprox)\n    hp, hc = get_td_waveform(approximant=gwapprox,\n                             mass1=16 + np.random.randint(0,10),\n                             mass2=16 + np.random.randint(0,10),\n                             delta_t=1.0/4096,\n                             spin1z=0.5 + np.random.rand()*0.5,\n                             spin2z=0.25 + np.random.rand()*0.5,\n                             inclination= 2 * np.pi * np.random.rand(),\n                             coa_phase= 2 * np.pi * np.random.rand(),\n                             phase_order = np.random.randint(2,8),\n                             f_lower=np.random.randint(24,64),\n                             distance=int(np.random.randint(1,1000)),\n                            )\n\n    \n    # Choose a GPS end time, sky location, and polarization phase for the merger\n    # NOTE: Right ascension and polarization phase runs from 0 to 2pi\n    #       Declination runs from pi/2. to -pi/2 with the poles at pi/2. and -pi/2.\n    end_time = 1192529720 + np.random.randint(1192529720//100000)\n    declination = np.pi * np.random.rand() - np.pi/2\n    right_ascension = 2 * np.pi * np.random.rand()\n    polarization = 2 * np.pi * np.random.rand()\n    hp.start_time += end_time\n    hc.start_time += end_time\n\n    signal_h1 = det_h1.project_wave(hp, hc,  right_ascension, declination, polarization)\n    signal_l1 = det_l1.project_wave(hp, hc,  right_ascension, declination, polarization)\n    signal_v1 = det_v1.project_wave(hp, hc,  right_ascension, declination, polarization)    \n    minlen = np.min( [len(signal_h1), len(signal_l1), len(signal_v1)] )\n    data = np.stack( (signal_h1[:minlen], signal_l1[:minlen], signal_v1[:minlen]),  ) * 1e19\n    print(data.shape)\n    \n    if data.shape[1]>4096:\n        data = data[:,data.shape[1]-4096:]\n        for N in range(80):\n            data[:,N] *= 1./(N+1)\n    \n    if len(hp)<4096:\n        for N in range(80):\n            data[:,N] *= 1./(N+1)\n        data = np.pad(data, ((0,0),(4096-data.shape[1],0)) )\n    \n    \n    plt.plot(data[0])\n    plt.plot(data[1])\n    plt.plot(data[2])\n    plt.show()\n    \n    cwt, freqs = pywt.cwt(data, scales=np.arange(1, 95, 0.62), wavelet='cmor1.5-0.95', sampling_period=1/2048, method='fft')\n    cwt = cwt.transpose(0,2,1)\n    print(cwt.shape)\n    cwt = np.log1p( np.abs(cwt) )\n    print( cwt.min(), cwt.max())\n    cwt -= cwt.min()\n    cwt /= cwt.max()\n    plt.imshow(cv2.resize(cwt,(256, 256)) )\n    plt.title(gwapprox)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-06T16:39:54.574437Z","iopub.execute_input":"2021-09-06T16:39:54.574870Z","iopub.status.idle":"2021-09-06T16:41:10.612815Z","shell.execute_reply.started":"2021-09-06T16:39:54.574836Z","shell.execute_reply":"2021-09-06T16:41:10.608357Z"},"trusted":true},"execution_count":null,"outputs":[]}]}