{"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":"!pip install git+https://github.com/PyFstat/PyFstat@python37","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-21T14:23:10.4133Z","iopub.execute_input":"2022-12-21T14:23:10.413738Z","iopub.status.idle":"2022-12-21T14:23:25.438726Z","shell.execute_reply.started":"2022-12-21T14:23:10.413704Z","shell.execute_reply":"2022-12-21T14:23:25.43755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\nimport h5py","metadata":{"execution":{"iopub.status.busy":"2022-12-21T14:23:25.44117Z","iopub.execute_input":"2022-12-21T14:23:25.441558Z","iopub.status.idle":"2022-12-21T14:23:25.448805Z","shell.execute_reply.started":"2022-12-21T14:23:25.441524Z","shell.execute_reply":"2022-12-21T14:23:25.44712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random\nrandom.choice([-1,1])*random.uniform(0.1,0.15)","metadata":{"execution":{"iopub.status.busy":"2022-12-21T18:29:48.242485Z","iopub.execute_input":"2022-12-21T18:29:48.24318Z","iopub.status.idle":"2022-12-21T18:29:48.251702Z","shell.execute_reply.started":"2022-12-21T18:29:48.243129Z","shell.execute_reply":"2022-12-21T18:29:48.250078Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import shutil","metadata":{"execution":{"iopub.status.busy":"2022-12-21T14:23:25.450333Z","iopub.execute_input":"2022-12-21T14:23:25.450698Z","iopub.status.idle":"2022-12-21T14:23:25.463242Z","shell.execute_reply.started":"2022-12-21T14:23:25.450668Z","shell.execute_reply":"2022-12-21T14:23:25.462403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pyfstat\nfrom pyfstat.utils import get_sft_as_arrays","metadata":{"execution":{"iopub.status.busy":"2022-12-21T14:23:25.464813Z","iopub.execute_input":"2022-12-21T14:23:25.465121Z","iopub.status.idle":"2022-12-21T14:23:25.475851Z","shell.execute_reply.started":"2022-12-21T14:23:25.465094Z","shell.execute_reply":"2022-12-21T14:23:25.47482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nimport pyfstat\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom PIL import Image as Img\nimport random\nimport io\nimport os\nimport json\nfrom pyfstat.utils import get_sft_as_arrays\n\n\ndef json_filename_from_params(signal_v,tstart,f0,h0,f1,alpha,delta,cosi,psi,phi,sqrtSX):\n    savefname = [str(int(tstart))]\n    for value in [f0,h0,f1,alpha,delta,cosi,psi,phi,sqrtSX * 1e24]:\n        savefname.append(str(round(value, 2)))\n    savefname = \"_\".join(savefname) + f\"_{signal_v}.json\"\n    return savefname        \n\n\ndef full_json_from_data(g):\n    startpx = int(random.random() * (g[\"H1\"].shape[0] - 361))\n    \n    h1raw = g[\"H1\"][startpx:startpx+360]\n    l1raw = g[\"L1\"][startpx:startpx+360]\n    fullacross = min(h1raw.shape[1], l1raw.shape[1])\n    \n    img = np.empty((4, 360, fullacross), dtype=np.float32)\n\n    startpx = int(random.random() * (g[\"H1\"].shape[0] - 361))\n    a = g[\"H1\"][startpx:startpx+360][:, :fullacross] * 1e22\n    img[0] = a.real\n    img[1] = a.imag\n    \n    a = g[\"L1\"][startpx:startpx+360][:, :fullacross] * 1e22\n    img[2] = a.real\n    img[3] = a.imag\n    \n    return img\n\ndef make_h5py_file_slice(freqs,fourier_data,timestamps,i,path=\"/kaggle/working/noise/\",name=\"Gausian_Noise\"):\n#frequency,fourier_data,timestamps,path=\"/kaggle/working/noise/\",name=\"Gausian_Noise\",i        \n        \n        print(freqs.shape)\n        slice_min=freqs.shape[0]//2-180\n        slice_max=freqs.shape[0]//2+180\n        file_path=path+name\n        f = h5py.File(f'{file_path}{i}.hdf5','w')\n        g0 = f.create_group(f\"{name}{i}\")\n        g1 = g0.create_group(\"H1\")\n        g2 = g0.create_group(\"L1\")\n        print(freqs.shape[0])\n        frequency= freqs[slice_min:slice_max]\n        print(frequency.shape[0])\n        dset = g0.create_dataset(\"frequency_Hz\", data=frequency)\n        temp=fourier_data['H1']\n        temp= temp[slice_min:slice_max,:]\n        print(temp.shape)\n        temp=temp[0:360,:]\n        g1.create_dataset(\"SFTs\",data=temp)\n        g1.create_dataset(\"timestamps_GPS\",data=timestamps['H1'])\n        temp=fourier_data['L1']\n        temp= temp[slice_min:slice_max,:]\n        temp=temp[0:360,:]\n        g2.create_dataset(\"SFTs\",data=temp)\n        g2.create_dataset(\"timestamps_GPS\",data=timestamps['L1'])\n        f.close()\ndef make_h5py_file_bslice(freqs,fourier_data_l1,fourier_data_h1,timestamps_h1,timestamps_l1,i,path=\"/kaggle/working/noise/\",name=\"Gausian_Noise\"):\n#frequency,fourier_data,timestamps,path=\"/kaggle/working/noise/\",name=\"Gausian_Noise\",i        \n        \n        print(freqs.shape)\n        slice_min=freqs.shape[0]//2-180\n        slice_max=freqs.shape[0]//2+180\n        file_path=path+name\n        f = h5py.File(f'{file_path}{i}.hdf5','w')\n        g0 = f.create_group(f\"{name}{i}\")\n        g1 = g0.create_group(\"H1\")\n        g2 = g0.create_group(\"L1\")\n        print(freqs.shape[0])\n        frequency= freqs[slice_min:slice_max]\n        print(frequency.shape[0])\n        dset = g0.create_dataset(\"frequency_Hz\", data=frequency)\n        temp=fourier_data_h1\n        temp= temp[slice_min:slice_max,:]\n        print(temp.shape)\n        temp=temp[0:360,:]\n        g1.create_dataset(\"SFTs\",data=temp)\n        g1.create_dataset(\"timestamps_GPS\",data=timestamps_h1)\n        temp=fourier_data_l1\n        temp= temp[slice_min:slice_max,:]\n        temp=temp[0:360,:]\n        g2.create_dataset(\"SFTs\",data=temp)\n        g2.create_dataset(\"timestamps_GPS\",data=timestamps_l1)\n        f.close()\n\n\n# os.mkdir(\"./signals\")\n# os.mkdir(\"./signals_test\")\nnoise='/kaggle/working/noise/'\nif not os.path.exists(noise):\n    os.mkdir(noise)\nsignal='/kaggle/working/signal/'\nif not os.path.exists(signal):\n    os.mkdir(signal)\n\n\n\n\n# for foldername in [\"signal\", \"noise\"]:\n\n    print(\"GENERATING\")\n    \n    SIGNALS_TO_GENERATE = 150\n    \nfor i_images in range(SIGNALS_TO_GENERATE):\n    print(\"====================== IMAGE \", i_images)\n\n    sft_path = []\n\n#     try:\n    pi = 3.141592\n\n    tstart = 1238166018\n    totaldiff = 200000\n    tdiff = int(random.random() * totaldiff) - (totaldiff / 2)\n    tstart = tstart - tdiff\n\n    f0 = random.random() * (500 - 50) + 50\n    h0 = random.random() * (40 - 10) + 10\n\n    f1 = random.random() * (1.0e-8 - 1.0e-12) + 1.0e-12\n    alpha = random.random() * (2*pi - 0) + 0\n    delta = random.random() * ((pi/2) - (-pi/2)) + (-pi/2)\n    cosi = random.random() * (1 - (-1)) + (-1)\n    psi = random.random() * ((pi/4) - (-pi/4)) + (-pi/4)\n    phi = random.random() * ((2*pi) - (0)) + (0)\n\n    sqrtSX = 5e-24\n\n    band = 0.36\n    tdur = 12*7.5 * 86400\n\n#         signal_kwargs = {\n#             \"outdir\": \"gen_outdir\",\n#             \"tstart\": tstart,\n#             \"duration\": tdur,\n#             \"detectors\": \"H1,L1\",\n#             \"Tsft\": 1800,\n#             \"Band\": band,\n\n#             \"sqrtSX\": sqrtSX,\n#             \"SFTWindowType\": \"tukey\", \n#             \"SFTWindowBeta\": 0.01,\n\n#             \"F0\": f0, # 50 to 500,\n#             \"Alpha\": alpha, # 0 to 2pi\n#             \"Delta\": delta, # -pi/2 to pi/2\n#             \"cosi\": cosi, # -1.0 to 1.0\n#             \"psi\": psi, # any (mapped -pi/4 to pi/4)\n#             \"phi\": phi, # any (mapped 0 to 2pi)\n#         }\n\n#         signal_writer = pyfstat.Writer( **signal_kwargs )\n#         signal_writer.make_data()        \n#         sft_path.append(signal_writer.sftfilepath)\n\n#         freqs, times, sft_data = get_sft_as_arrays(\";\".join(sft_path))\n\n#         gen_signal = random.random() < 0.5\n#         #gen_signal = random.random() < -0.5\n\n#         writer_kwargs = {\n#             \"label\": \"single_detector_spectral_line\",\n#             \"outdir\": \"PyFstat_example_data\",\n#             \"tstart\": 1238166018,  # Starting time of the observation [GPS time]\n#             \"duration\": 5 * 86400,  # Duration [seconds]\n#             \"detectors\": \"H1\",  # Detector to simulate, in this case LIGO Hanford\n#             \"F0\": f0,  # Central frequency of the band to be generated [Hz]\n#             \"phi\": 1.3,  # Initial phase of the spectral line\n#             \"Band\": 1.0,  # Frequency band-width around F0 [Hz]                \"h0\": 1e-24,              # Amplitude of the spectral line\n#             \"sqrtSX\": 5e-24,  # Single-sided Amplitude Spectral Density of the noise\n#             \"Tsft\": 1800,  # Fourier transform time duration\n#             \"SFTWindowType\": \"tukey\",\n#             \"SFTWindowBeta\": 0.01,\n#         }\n\n#         writer = pyfstat.LineWriter(**writer_kwargs)\n#         writer.make_data()\n#         sft_path.append(writer.sftfilepath)\n# #         frequency, timestamps, fourier_data = get_sft_as_arrays(\";\".join(sft_path))\n#         frequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)\n\n#         if not gen_signal:\n#             print(i_images)\n#             print(sft_data['H1'].shape,np.min( freqs),np.max( freqs),np.max( freqs)-np.min( freqs),np.max( freqs)-np.min( freqs))\n#             fourier_data_l1=sft_data ['L1']\n#             fourier_data_h1=fourier_data['H1']\n#             timestamps_h1=timestamps['H1']\n#             timestamps_l1=times['L1']\n#             make_h5py_file_bslice(frequency,fourier_data_l1,fourier_data_h1,timestamps_h1,timestamps_l1,i_images,path=\"/kaggle/working/noise/\",name=\"Gausian_Noise\")\n\n#             #make_h5py_file_slice(freqs,sft_data,times,i_images,path=\"/kaggle/working/noise/\",name=\"SNRGausian_Noise\")\n# #             make_h5py_file_slice(frequency, timestamps, fourier_data,i_images,path=\"/kaggle/working/noise/\",name=\"SNRArtefact_Noise\")\n# #             gg = full_json_from_data(sft_data)    \n# #             fname = json_filename_from_params(\"0\",tstart,f0,h0,f1,alpha,delta,cosi,psi,phi,sqrtSX)\n# #             np.save(f\"./{foldername}/{fname}\", gg) \n\n\n#         if gen_signal:\n#             signal_kwargs2 = {\n#                 \"outdir\": \"gen_outdir\",\n#                 \"tstart\": tstart,\n#                 \"duration\": 12*7.5 * 86400,\n#                 \"detectors\": \"H1,L1\",\n#                 \"Tsft\": 1800,\n#                 \"Band\": band, \n\n#                 \"sqrtSX\": 0,\n#                 \"SFTWindowType\": \"tukey\", \n#                 \"SFTWindowBeta\": 0.01,\n\n#                 \"F0\": f0, # 50 to 500\n\n#                 \"h0\": sqrtSX / h0, # div by 10 to 1000\n#                 \"F1\": f1, # 1.0e-12 to 1.0e-8 (or so; could be to 0)\n#                 \"Alpha\": alpha, # 0 to 2pi\n#                 \"Delta\": delta, # -pi/2 to pi/2\n#                 \"cosi\": cosi, # -1.0 to 1.0\n#                 \"psi\": psi, # any (mapped -pi/4 to pi/4)\n#                 \"phi\": phi, # any (mapped 0 to 2pi)\n\n#                 \"noiseSFTs\": \";\".join(sft_path)\n#             }\n#             signal_writer = pyfstat.Writer( **signal_kwargs2 )\n#             signal_writer.make_data()\n#             signal_sft_path = signal_writer.sftfilepath\n\n\n#             freqs, times, sft_data = get_sft_as_arrays(signal_sft_path)\n#             print(sft_data['H1'].shape,np.min( freqs),np.max( freqs),np.max( freqs)-np.min( freqs),np.max( freqs)-np.min( freqs))\n#             make_h5py_file_slice(freqs,sft_data,times,i_images,path=\"/kaggle/working/signal/\",name=\"SignalSNRGausian_Noise\")\n\n#             #             gg = full_json_from_data(sft_data)\n\n# #             fname = json_filename_from_params(\"1\",tstart,f0,h0,f1,alpha,delta,cosi,psi,phi,sqrtSX)\n# #             np.save(f\"./{foldername}/{fname}\", gg)  \n\n\n####################\n    detector=random.choice(['H1','L1'])\n    def other_one(detector):\n        if detector=='H1':\n            return 'L1'\n        else:\n            return 'H1'\n    detector2=other_one(detector)\n    no_signal = random.random() < 0.5\n    #Artefacts in one detector\n    writer_kwargs2 = {\n        \"label\": \"single_detector_spectral_line\",\n        \"outdir\": \"gen_outdir2\",\n        \"tstart\": tstart,\n        \"duration\": tdur,  # Duration [seconds]\n        \"detectors\": detector,  # Detector to simulate, in this case LIGO Hanford\n        \"F0\":  f0-random.choice([-1,1])*random.uniform(0.1,0.15),  # Central frequency of the band to be generated [Hz]\n        \"phi\": phi,  # Initial phase of the spectral line\n        \"Band\": band,  # Frequency band-width around F0 [Hz]                \"h0\": 1e-24,              # Amplitude of the spectral line\n        \"sqrtSX\": sqrtSX,  # Single-sided Amplitude Spectral Density of the noise\n        \"Tsft\": 1800,  # Fourier transform time duration\n        \"SFTWindowType\": \"tukey\",\n        \"SFTWindowBeta\": 0.01,\n    }\n    writer2 = pyfstat.LineWriter(**writer_kwargs2)\n    writer2.make_data()\n    sft_path = []\n    frequency, timestamps, fourier_data = get_sft_as_arrays(writer2.sftfilepath)\n    print(frequency.shape)\n    fourier_data_1=fourier_data[detector]\n    timestamps_1=timestamps[detector]\n    signal_kwargs = {\n        \"outdir\": \"gen_outdir\",\n        \"tstart\": tstart,\n        \"duration\": tdur,\n        \"detectors\": detector2,\n        \"Tsft\": 1800,\n        \"Band\": band,\n\n        \"sqrtSX\": sqrtSX,\n        \"SFTWindowType\": \"tukey\", \n        \"SFTWindowBeta\": 0.01,\n\n        \"F0\": f0, # 50 to 500\n    }\n\n    signal_writer = pyfstat.Writer(**signal_kwargs)\n    sft_path.append(signal_writer.sftfilepath)\n    # Create SFTs\n    signal_writer.make_data()\n    frequency, timestamps, fourier_data = get_sft_as_arrays(signal_writer.sftfilepath)\n    fourier_data_2=fourier_data[detector2]\n    timestamps_2=timestamps[detector2]\n    make_h5py_file_bslice(frequency,fourier_data_1,fourier_data_2,timestamps_1,timestamps_2,i_images,path=\"/kaggle/working/noise/\",name=\"ArtefactandGaussianNoise\")\n            #Noise in tehe second detector\n            # Setup Writer\n\n    signal_kwargs2 = {\n            \"outdir\": \"gen_outdir\",\n            \"tstart\": tstart,\n            \"duration\": 12*7.5 * 86400,\n            \"detectors\": \"H1,L1\",\n            \"Tsft\": 1800,\n            \"Band\": band, \n            \"sqrtSX\": 0,\n            \"SFTWindowType\": \"tukey\", \n            \"SFTWindowBeta\": 0.01,\n            \"F0\": f0, # 50 to 500\n            \"h0\": sqrtSX / h0, # div by 10 to 1000\n            \"F1\": f1, # 1.0e-12 to 1.0e-8 (or so; could be to 0)\n            \"Alpha\": alpha, # 0 to 2pi\n            \"Delta\": delta, # -pi/2 to pi/2\n            \"cosi\": cosi, # -1.0 to 1.0\n            \"psi\": psi, # any (mapped -pi/4 to pi/4)\n            \"phi\": phi, # any (mapped 0 to 2pi)\n\n            \"noiseSFTs\": \";\".join(sft_path)\n        }\n\n    signal_writer = pyfstat.Writer( **signal_kwargs2 )\n    signal_writer.make_data()    \n    signal_sft_path = signal_writer.sftfilepath\n\n#         writer3 = pyfstat.Writer(**writer_kwargs3, **signal_parameters)\n#         writer3.make_data()\n    frequency, timestamps, fourier_data = get_sft_as_arrays(signal_sft_path)\n\n    # Create SFTs\n    fourier_data_3=fourier_data[detector2]\n    timestamps_3=timestamps[detector2]  \n    make_h5py_file_bslice(frequency,fourier_data_1,fourier_data_3,timestamps_1,timestamps_3,i_images,path=\"/kaggle/working/signal/\",name=\"NoisySignal\")\n    try:\n        shutil.rmtree(\"/kaggle/working/gen_outdir\")\n    except:\n        print('directory not found')\n#                 os.chdir(\"/kaggle/working/gen_outdir/\")\n#                 for item in clean:\n#                     if item.endswith(\".sft\"):\n#                         try:\n#                             os.remove(item)\n#                         except:\n#                             print(item)\n#                     if item.endswith(\".cff\"):\n#                         try:\n#                             os.remove(item)\n#                         except:\n#                             print(item)\n#####################    \n\n#         for sft_f in signal_sft_path.split(\";\"):\n#             try:\n#                 os.remove(sft_f)\n#             except:\n#                 print(\"NO FILE 1: \", sft_f)\n#         for sft_f in sft_path:\n#             for sft_f_f in sft_f.split(\";\"):\n#                 try:\n#                     os.remove(sft_f_f)\n#                 except:\n#                     print(\"NO FILE 2: \", sft_f_f)\n\n\nprint(\"DONE\")\n","metadata":{"execution":{"iopub.status.busy":"2022-12-21T14:24:27.756202Z","iopub.execute_input":"2022-12-21T14:24:27.756648Z","iopub.status.idle":"2022-12-21T14:25:03.062442Z","shell.execute_reply.started":"2022-12-21T14:24:27.756614Z","shell.execute_reply":"2022-12-21T14:25:03.061017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lista= [\"/kaggle/working/signal/\",\"/kaggle/working/noise/\",\"/kaggle/working/gen_outdir\"]\nprint(lista)\ntry:\n    for elem in lista:\n        clean=os.listdir(elem)\n        os.chdir(elem)\n        for item in clean:\n            if item.endswith(\".sft\"):\n                try:\n                    os.remove(item)\n                except:\n                    print(item)\n            if item.endswith(\".cff\"):\n                try:\n                    os.remove(item)\n                except:\n                    print(item)\nexcept:\n    print('directorul PyFstat_example_data nu exista')","metadata":{"execution":{"iopub.status.busy":"2022-12-21T14:23:29.015759Z","iopub.status.idle":"2022-12-21T14:23:29.016219Z","shell.execute_reply.started":"2022-12-21T14:23:29.016019Z","shell.execute_reply":"2022-12-21T14:23:29.016039Z"},"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":"","metadata":{},"execution_count":null,"outputs":[]}]}