{"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":"\"\"\"\nMCMC search: Fully coherent F-statistic\n=======================================\n\nDirected MCMC search for an isolated CW signal using the\nfully coherent F-statistic.\n\"\"\"\n\nimport os\n\nimport numpy as np\n\nimport pyfstat\nfrom pyfstat.utils import get_predict_fstat_parameters_from_dict\n\n#idx = '9f5acbd94'\n#idx = '00f226552'\nidx = 'fb1e0df09'\n\nlabel = \"PyFstat_example_fully_coherent_MCMC_search_%s\" % idx\noutdir = os.path.join(\"PyFstat_example_data2\", label)\nlogger = pyfstat.set_up_logger(label=label, outdir=outdir)\n\nsftfilepattern = '/home/jon/kaggle/contgrav/sfts_%s/*H1*.sft' % idx\nwith np.load('/home/jon/kaggle/contgrav/sfts_%s/H1.meta.npy.npz' % idx) as data:\n    H1_duration = int(data['duration1'])\n    f1b = float(data['fb'])\n    f1e = float(data['fe'])\n    t1b = int(data['tb'])\n    t1e = int(data['te'])\nwith np.load('/home/jon/kaggle/contgrav/sfts_%s/L1.meta.npy.npz' % idx) as data:\n    L1_duration = int(data['duration1'])\n    f2b = float(data['fb'])\n    f2e = float(data['fe'])\n    t2b = int(data['tb'])\n    t2e = int(data['te'])\n\ntbmin = int(min(t1b, t2b))\ntbmax = int(min(t1e, t2e))\n\nmid_time = 0.5 * (tbmin + tbmax)\nmid_time = int(mid_time)\nF0mid1 = 0.5 * (f1b + f1e)\nF0mid2 = 0.5 * (f2b + f2e)\nF0mid = 0.5 * (F0mid1 + F0mid2)\nF1guess = 1e-11\ntstart = tbmin + 10 * 1800\ntend = tbmax - 10 * 1800\n\ndata_parameters = dict(duration=max(H1_duration, L1_duration))\n#fudge = 0.0002\n#lower = float(min(f1b, f2b)) * (1 + fudge)\n#higher = float(min(f1e, f2e)) * (1 - fudge)\n\nDeltaF0 = 0.2/360.0\nF0mid = int(F0mid / DeltaF0) * DeltaF0\n\n# signal priors\nsignal_parameters = dict(F0=F0mid,\n                         F1=F1guess,\n                         F2=0,\n                         Alpha=np.radians(83.6292),\n                         Delta=np.radians(22.0144),\n                         )\n\n#DeltaF0 = 1e-7  # what's this?\n#DeltaF0 = 0.2/360.0\n#DeltaF0 = f1b + 2 * 0.2/360.0\nlower = f1b\nhigher = f1e\n\nlower = 406.49666666666667\nhigher = lower + 30 * 0.2/360\n\n#lower = int(lower/DeltaF0) * DeltaF0\n#higher = int(higher/DeltaF0) * DeltaF0\n\n#DeltaF0 =\n#lower = signal_parameters[\"F0\"] - DeltaF0 / 2.0\n#higher = signal_parameters[\"F0\"] + DeltaF0 / 2.0\nprint(\"F0: %s lower: %s higher: %s\" % (signal_parameters[\"F0\"], lower, higher))\n\nF1lower = -1e-10\nF1higher = -1e-12\n\n#DeltaF1 = 1e-13\n#VF0 = (np.pi * data_parameters[\"duration\"] * DeltaF0) ** 2 / 3.0\n#VF1 = (np.pi * data_parameters[\"duration\"] ** 2 * DeltaF1) ** 2 * 4 / 45.0\n#logger.info(\"\\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\\n\".format(VF0 * VF1, VF0, VF1))\n\n\n\n# signal priors\n# [\"F0\", \"F1\", \"F2\", \"Alpha\", \"Delta\"]\ntheta_prior = {\n    \"F0\": {\n        \"type\": \"unif\",\n        \"lower\": lower,\n        \"upper\": higher,\n    },\n    \"F1\": {\n        \"type\": \"unif\",\n        \"lower\": F1lower,\n        \"upper\": F1higher,\n    },\n    \"Alpha\": {\n        \"type\": \"unif\",\n        \"lower\": 0,\n        \"upper\": 2.0*np.pi,\n    },\n    \"Delta\": {\n        \"type\": \"unif\",\n        \"lower\": -1,\n        \"upper\": +1,\n    },\n    #\"cosi\": {\n    #    \"type\": \"unif\",\n    #    \"lower\": -1,\n    #    \"upper\": +1,\n    #},\n    #\"psi\": {\n    #    \"type\": \"unif\",\n    #    \"lower\": -.25*np.pi,\n    #    \"upper\": +0.25*np.pi,\n    #},\n    #\"phi\": {\n    #    \"type\": \"unif\",\n    #    \"lower\": 0,\n    #    \"upper\": 2.0*np.pi,\n    #},\n}\nfor key in [\"F2\"]:\n    theta_prior[key] = signal_parameters[key]\n\nntemps = 2\nlog10beta_min = -0.5\nnwalkers = 100\nnsteps = [300, 300]\n\nmcmc = pyfstat.MCMCSearch(\n    label=label,\n    outdir=outdir,\n    sftfilepattern=sftfilepattern,\n    theta_prior=theta_prior,\n    tref=mid_time,\n    #tref=tstart,\n    minStartTime=tstart,\n    maxStartTime=tend,\n    #minCoverFreq=lower,\n    #maxCoverFreq=higher,\n    nsteps=nsteps,\n    nwalkers=nwalkers,\n    ntemps=ntemps,\n    log10beta_min=log10beta_min,\n)\nmcmc.transform_dictionary = dict(\n    F0=dict(subtractor=F0mid, symbol=\"$f-f^\\\\mathrm{s}$\"),\n    F1=dict(\n        subtractor=signal_parameters[\"F1\"], symbol=\"$\\\\dot{f}-\\\\dot{f}^\\\\mathrm{s}$\"\n    ),\n)\nmcmc.run(\n    walker_plot_args={\"plot_det_stat\": True, \"injection_parameters\": signal_parameters}\n)\nmcmc.print_summary()\nmcmc.plot_corner(add_prior=True, truths=signal_parameters)\nmcmc.plot_prior_posterior(injection_parameters=signal_parameters)\n\nmcmc.generate_loudest()\n\n# plot cumulative 2F, first building a dict as required for PredictFStat\nd, maxtwoF = mcmc.get_max_twoF()\nfor key, val in mcmc.theta_prior.items():\n    if key not in d:\n        d[key] = val\nd[\"h0\"] = 1e-23\nd[\"cosi\"] = 0.0\nd[\"psi\"] = 0.0\nPFS_input = get_predict_fstat_parameters_from_dict(d)\nmcmc.plot_cumulative_max(PFS_input=PFS_input)\n\nprint(\"done\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]}]}