{"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":"# Here is a simple script to generate gravitational waves.\n\nSeted up two gravitational wave detectors 90 degrees apart.\n\nYou can simply set the preprocessing method.\n\nThe default preprocessing method is to take the absolute value, normalize it, and compress it.\n\nAlso, include the method of generating noise.\n\nOutput to /kaggle/working/signal and /kaggle/working/noise\n","metadata":{}},{"cell_type":"markdown","source":"# Install module and import","metadata":{}},{"cell_type":"code","source":"!pip install git+https://github.com/PyFstat/PyFstat@python37","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2022-11-05T04:49:40.846368Z","iopub.execute_input":"2022-11-05T04:49:40.847109Z","iopub.status.idle":"2022-11-05T04:50:34.550798Z","shell.execute_reply.started":"2022-11-05T04:49:40.847009Z","shell.execute_reply":"2022-11-05T04:50:34.549453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport torch as tc\nimport pyfstat,sys,os,io,cv2,warnings\nfrom pyfstat.utils import get_sft_as_arrays\nfrom typing import TYPE_CHECKING, Iterable, Optional\nimport logging,shutil,tqdm\nimport matplotlib.pyplot as plt\nwarnings.filterwarnings(\"ignore\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-05T04:50:34.553101Z","iopub.execute_input":"2022-11-05T04:50:34.554126Z","iopub.status.idle":"2022-11-05T04:50:38.157542Z","shell.execute_reply.started":"2022-11-05T04:50:34.554087Z","shell.execute_reply":"2022-11-05T04:50:38.156483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Clear unnecessary log","metadata":{}},{"cell_type":"code","source":"\"\"\"fc:\"set_up_logger\" is come from https://www.kaggle.com/code/crischir/pyfstat-tutorial-adapted-to-kaggle\"\"\"\ndef set_up_logger(\n    outdir: Optional[str] = None,\n    label: Optional[str] = \"pyfstat\",\n    log_level: str = \"INFO\",  # FIXME: Requires Python 3.8 Literal[\"CRITICAL\", \"ERROR\", \"WARNING\", \"INFO\", \"DEBUG\"] = \"INFO\",\n    streams: Optional[Iterable[\"io.TextIOWrapper\"]] = (sys.stdout,),\n    append: bool = True,\n) -> logging.Logger:\n#     \"\"\"Add file and stream handlers to the `pyfstat` logger.\n#     Handler names generated from ``streams`` and ``outdir, label``\n#     must be unique and no duplicated handler will be attached by\n#     this function.\n#     Parameters\n#     ----------\n#     outdir:\n#         Path to outdir directory. If ``None``, no file handler will be added.\n#     label:\n#         Label for the file output handler, i.e.\n#         the log file will be called `label.log`.\n#         Required, in conjunction with ``outdir``, to add a file handler.\n#         Ignored otherwise.\n#     log_level:\n#         Level of logging. This level is imposed on the logger itself and\n#         *every single handler* attached to it.\n#     streams:\n#         Stream to which logging messages will be passed using a\n#         StreamHandler object. By default, log to ``sys.stdout``.\n#         Other common streams include e.g. ``sys.stderr``.\n#     append:\n#         If ``False``, removes all handlers from the `pyfstat` logger\n#         before adding new ones. This removal is not propagated to\n#         handlers on the `root` logger.\n#     Returns\n#     -------\n#     obj:\n#         Configured instance of the ``logging.Logger`` class.\n#     \"\"\"\n    logger = logging.getLogger(\"pyfstat\")\n    logger.setLevel(log_level)\n\n    if not append:\n        for handler in logger.handlers:\n            logger.removeHandler(handler)\n    else:\n        for handler in logger.handlers:\n            handler.setLevel(log_level)\n\n    stream_names = [\n        handler.stream.name\n        for handler in logger.handlers\n        if type(handler) == logging.StreamHandler\n    ]\n    file_names = [\n        handler.baseFilename\n        for handler in logger.handlers\n        if type(handler) == logging.FileHandler\n    ]\n\n    common_formatter = logging.Formatter(\n        \"%(asctime)s.%(msecs)03d %(name)s %(levelname)-8s: %(message)s\",\n        datefmt=\"%y-%m-%d %H:%M:%S\",  # intended to match LALSuite's format\n    )\n\n    for stream in streams or []:\n        if stream.name in stream_names:\n            continue\n        stream_handler = logging.StreamHandler(stream)\n        stream_handler.setFormatter(common_formatter)\n        stream_handler.setLevel(log_level)\n        logger.addHandler(stream_handler)\n\n    if label and outdir:\n        os.makedirs(outdir, exist_ok=True)\n        log_file = os.path.join(outdir, f\"{label}.log\")\n\n        if log_file not in file_names:\n\n            file_handler = logging.FileHandler(log_file)\n            file_handler.setFormatter(common_formatter)\n            file_handler.setLevel(log_level)\n            logger.addHandler(file_handler)\n\n    return logger\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-05T04:50:38.15886Z","iopub.execute_input":"2022-11-05T04:50:38.159999Z","iopub.status.idle":"2022-11-05T04:50:38.174487Z","shell.execute_reply.started":"2022-11-05T04:50:38.159932Z","shell.execute_reply":"2022-11-05T04:50:38.173168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def Normalization(x:tc.Tensor)->tc.Tensor:\n    \"\"\"input.shape=(batch,f1,f2,...)\"\"\"\n    #[batch,f1,f2]->dim[1,2]\n    dim=list(range(1,x.ndim))\n    mean=x.mean(dim=dim,keepdim=True)\n    var=x.std(dim=dim,keepdim=True)\n    return (x-mean)/var","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-05T04:50:38.176711Z","iopub.execute_input":"2022-11-05T04:50:38.177519Z","iopub.status.idle":"2022-11-05T04:50:38.189804Z","shell.execute_reply.started":"2022-11-05T04:50:38.177479Z","shell.execute_reply":"2022-11-05T04:50:38.189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Generate gravitational waves","metadata":{}},{"cell_type":"code","source":"def generate(signal_rate=0.01)->np.ndarray:\n    \"\"\"signal_rate must >=0\n        output=[H1_signal,L1_signal]\n        \"\"\"\n    writer_kwargs = {\n        \"label\": \"single_detector_gaussian_noise\",\n        \"outdir\": \"PyFstat_example_data\",\n        \"tstart\": 1238166018,\n        \"duration\": 365 * 86400,\n        \"detectors\": \"H1\",\n        \"sqrtSX\": 1e-22,\n        \"Tsft\": 1800,\n        \"SFTWindowType\": \"tukey\",\n        \"SFTWindowBeta\": 0.01,\n    }\n    signal_parameters = {\n        \"F0\": np.random.randint(50,500),#meta_data_arange=(50,400)\n        \"F1\": -1e-9*(1+np.random.randn()*0.01),#[-1e-9,0]\n        \"Alpha\": 2*np.pi*np.random.rand(),#(0,2pi)\n        \"Delta\": np.pi*(0.5-np.random.rand()),#[-pi/2,pi/2]\n        \"h0\": 1e-22*signal_rate,\n        \"cosi\": 1-2*np.random.rand()*0,#(-1,1)Affects signal strength\n        \"psi\": 0.5*np.pi*(0.5-np.random.rand()),#(-pi/4,pi/4)\n        \"phi\": 2*np.pi*np.random.rand(),#(0,2pi)\n        \"tref\": writer_kwargs[\"tstart\"],\n    }\n\n    writer = pyfstat.Writer(**writer_kwargs, **signal_parameters)\n    writer.make_data()# Create SFTs\n\n    frequency, timestamps, fourier_data_H1 = get_sft_as_arrays(writer.sftfilepath)#H1_detector\n\n    signal_parameters[\"psi\"]=signal_parameters[\"psi\"]+np.pi/2#90 degree\n    writer = pyfstat.Writer(**writer_kwargs, **signal_parameters)\n    writer.make_data()\n    \n    frequency, timestamps, fourier_data_L1 = get_sft_as_arrays(writer.sftfilepath)#L1_detector\n\n    \n    ####preprocessing method\n    fourier_data_H1[\"H1\"]/=fourier_data_H1[\"H1\"].mean()\n    fourier_data_L1[\"H1\"]/=fourier_data_L1[\"H1\"].mean()\n    \n    output_H=fourier_data_H1[\"H1\"].real**2+fourier_data_H1[\"H1\"].imag**2\n    #L1 is another detector with a 90 degree difference\n    output_L=fourier_data_L1[\"H1\"].real**2+fourier_data_L1[\"H1\"].imag**2\n    \n    ####compress gravitational waves to speed up training\n    output=np.stack((output_H,output_L),axis=0)\n\n    output=output.transpose(1,2,0)\n    #output.shape=(200~360(perhaps),17520, 2)\n    output=cv2.resize(output,(3500,360)).transpose(2,0,1)#Compression_ratio=17520/3500~=5\n    \n    return frequency,Normalization(tc.from_numpy(output))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-11-05T04:50:38.191227Z","iopub.execute_input":"2022-11-05T04:50:38.191842Z","iopub.status.idle":"2022-11-05T04:50:38.20603Z","shell.execute_reply.started":"2022-11-05T04:50:38.19181Z","shell.execute_reply":"2022-11-05T04:50:38.204738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The output path","metadata":{}},{"cell_type":"code","source":"\"\"\"path settings\"\"\"\n\nlogger = set_up_logger(label=\"0_generating_noise\", log_level=\"WARNING\")\nfile_path=\"/kaggle/working/\"\ntry:\n    os.mkdir(path=file_path+\"signal\")\n    os.mkdir(path=file_path+\"noise\")\nexcept:\n    shutil.rmtree(file_path+\"signal\")#delete folder\n    shutil.rmtree(file_path+\"noise\")#delete folder\n    os.mkdir(path=file_path+\"signal\")\n    os.mkdir(path=file_path+\"noise\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-05T04:50:38.207456Z","iopub.execute_input":"2022-11-05T04:50:38.207888Z","iopub.status.idle":"2022-11-05T04:50:38.221425Z","shell.execute_reply.started":"2022-11-05T04:50:38.207855Z","shell.execute_reply":"2022-11-05T04:50:38.220383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Show gravitational waves","metadata":{}},{"cell_type":"code","source":"\"\"\"show Generate noisy gravitational waves signal_rate=0.1\"\"\"\nfrequency,output=generate(signal_rate=0.1)\nprint(output.shape)\noutput=output.numpy().transpose(1,2,0)\noutput=cv2.resize(output,(500,360)).transpose(2,0,1)\n\nfig, ax = plt.subplots(1,1,figsize=(8, 8))\nax.set(xlabel=\"SFT index\")\nc = ax.pcolorfast(np.arange(output.shape[1]),np.arange(output.shape[2]), output[0])\nfig.colorbar(c, ax=ax, orientation=\"horizontal\", label=\"Value\")\nfig","metadata":{"execution":{"iopub.status.busy":"2022-11-05T04:50:38.222779Z","iopub.execute_input":"2022-11-05T04:50:38.224271Z","iopub.status.idle":"2022-11-05T04:50:44.432532Z","shell.execute_reply.started":"2022-11-05T04:50:38.224223Z","shell.execute_reply":"2022-11-05T04:50:44.431207Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Output gravitational waves","metadata":{}},{"cell_type":"code","source":"\"\"\"In this script, noisy gravitational waves (signal_rate=0.01) are generated in /kaggle/working/signal\n     Generate noise without gravitational waves in /kaggle/working/noise\"\"\"\n\ntotal=100#noise_generate_time+signal_generate_time\nsignal_rate=100\n    \ntime=tqdm.tqdm(total=total)\nfor i in range(total//2):\n    frequency,output=generate(signal_rate=signal_rate)#Generate noisy gravitational waves \n    #output:tc.Tensor to np.ndarray ->output.numpy()\n    tc.save(output,file_path+f\"signal/_{i}.pth\")\n    frequency,output=generate(0)#Generate noisy\n    tc.save(output,file_path+f\"noise/_{i}.pth\")\n    time.update(2)","metadata":{"execution":{"iopub.status.busy":"2022-11-05T04:50:44.434361Z","iopub.execute_input":"2022-11-05T04:50:44.434733Z","iopub.status.idle":"2022-11-05T04:50:55.795799Z","shell.execute_reply.started":"2022-11-05T04:50:44.434701Z","shell.execute_reply":"2022-11-05T04:50:55.793818Z"},"trusted":true},"execution_count":null,"outputs":[]}]}