Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • GregAshton/PyFstat
  • Yuanhao.Zhang/PyFstat
  • dkeitel/PyFstat
3 results
Select Git revision
Show changes
Commits on Source (48)
Show whitespace changes
Inline Side-by-side
Showing
with 780 additions and 518 deletions
.gitlab-ci.yml
View file @ d8571a7d
test_app:
stages:
- Test
- Static Analysis
variables:
VENV_DIR: $CI_PROJECT_DIR/../venv-pyFstat
INSTALLER_DIR: $CI_PROJECT_DIR/install-cw-software
pytest:
stage: Test
tags: [ pyFstat ]
before_script:
- python3 -m venv $VENV_DIR
- source ${VENV_DIR}/bin/activate
- pip install --upgrade pip
- pip install -r requirements.txt
- pip install lalsuite
- pip install pytest
- export LAL_DATA_PATH=$HOME/ephemeris
- export LALPULSAR_DATADIR=$LAL_DATA_PATH
script:
- . /home/user1/lalsuite-install/etc/lalapps-user-env.sh
- /home/user1/anaconda2/bin/python tests.py Writer
- /home/user1/anaconda2/bin/python tests.py par
- /home/user1/anaconda2/bin/python tests.py BaseSearchClass
- /home/user1/anaconda2/bin/python tests.py ComputeFstat
- /home/user1/anaconda2/bin/python tests.py SemiCoherentSearch
- /home/user1/anaconda2/bin/python tests.py SemiCoherentGlitchSearch
- /home/user1/anaconda2/bin/python tests.py MCMCSearch
- /home/user1/anaconda2/bin/python tests.py GridSearch
- pip install -e $CI_PROJECT_DIR
# make sure to test *installed* version of pyFstat
- (cd .. && pytest $CI_PROJECT_DIR/tests.py --log-file=$CI_PROJECT_DIR/tests.log)
artifacts:
paths:
- ./*.log
name: testlogs
when: always
expire_in: 24h
static:
stage: Static Analysis
tags: [ pyFstat ]
script:
- source ${VENV_DIR}/bin/activate
- black --check --diff .
# - flake8 . ## not ready
README.md
View file @ d8571a7d
......@@ -3,110 +3,13 @@
This is a python package providing an interface to perform F-statistic based
continuous gravitational wave (CW) searches.
For documentation, please use the [wiki](https://gitlab.aei.uni-hannover.de/GregAshton/PyFstat/wikis/home).
This repository is now for archival purposes only,
active development has been moved to:
In the
[examples](https://gitlab.aei.uni-hannover.de/GregAshton/PyFstat/tree/master/examples),
we have a number of scripts demonstrating different use cases.
**===> https://github.com/PyFstat/PyFstat <===**
## Installation
### `python` installation
The scripts are written in `python 2.7+` and therefore require a working
`python` installation. While many systems come with a system wide python
installation, it can often be easier to manage a user-specific python
installation. This way one does not require root access to install or remove
modules. One method to do this, is to use the `conda` system, either through
the stripped down [miniconda](http://conda.pydata.org/miniconda.html)
installation, or the full-featured
[anaconda](https://www.continuum.io/downloads) (these are essentially the
same, but the `anaconda` version installs a variety of useful packages such as
`numpy` and `scipy` by default).
The fastest/easiest method is to follow your OS instructions
[here](https://conda.io/docs/install/quick.html) which will install Miniconda.
For the rest of this tutorial, we will make use of `pip` to install modules (
not all packages can be installed with `conda` and for those using alternatives
to `conda`, `pip` is more universal).
This can be installed with
```
$ conda install pip
```
### Clone the repository
In a terminal, clone the directory:
```
$ git clone https://gitlab.aei.uni-hannover.de/GregAshton/PyFstat.git
```
### Dependencies
`pyfstat` makes uses the following external python modules:
* [numpy](http://www.numpy.org/)
* [matplotlib](http://matplotlib.org/) >= 1.4
* [scipy](https://www.scipy.org/)
* [ptemcee](https://github.com/willvousden/ptemcee)
* [corner](https://pypi.python.org/pypi/corner/)
* [dill](https://pypi.python.org/pypi/dill)
* [peakutils](https://pypi.python.org/pypi/PeakUtils)
*Optional*
* [tqdm](https://pypi.python.org/pypi/tqdm)(optional), if installed, this
provides a useful progress bar and estimate of the remaining run-time.
* [bashplotlib](https://github.com/glamp/bashplotlib), if installed, presents
a histogram of the loaded SFT data
* [pathos](https://pypi.python.org/pypi/pathos), if installed, this provides
support for multiprocessing some functions.
For an introduction to installing modules see
[here](https://docs.python.org/3.5/installing/index.html). If you are using
`pip`, to install all of these modules, run
```
$ pip install -r /PATH/TO/THIS/DIRECTORY/requirements.txt
```
In addition to these modules, you also need a working **swig-enabled**
[`lalapps`](http://software.ligo.org/docs/lalsuite/lalsuite/) with
at least `lalpulsar`. A minimal confuration line to use when installing
`lalapps` is
```
$ ./configure --prefix=${HOME}/lalsuite-install --disable-all-lal --enable-lalpulsar --enable-lalapps --enable-swig
```
### `pyfstat` installation
The script can be installed system wide, assuming you are in the source directory, via
```
$ python setup.py install
```
or simply add this directory to your python path. To check that the installation
was successful, run
```
$ python -c 'import pyfstat'
```
if no error message is output, then you have installed `pyfstat`. Note that
the module will be installed to whichever python executable you call it from.
### Ephemeris installation
The scripts require a path to ephemeris files in order to use the
`lalpulsar.ComputeFstat` module. This can either be specified when initialising
each search (as one of the arguments), or simply by placing a file
`~/.pyfstat.conf` into your home directory which looks like
```
earth_ephem = '/home/<USER>/lalsuite-install/share/lalpulsar/earth00-19-DE421.dat.gz'
sun_ephem = '/home/<USER>/lalsuite-install/share/lalpulsar/sun00-19-DE421.dat.gz'
```
here, we use the default ephemeris files provided with `lalsuite`.
Please also refer to the github repository for installation and usage instructions
and to report issues.
### Contributors
......@@ -116,25 +19,29 @@ here, we use the default ephemeris files provided with `lalsuite`.
* Karl Wette
* Sylvia Zhu
This project is open to development, please feel free to contact us
for advice or just jump in and submit a pull request.
## Citing this work
If you use `PyFstat` in a publication we would appreciate if you cite the
original paper introducing the code, the [ADS page can be found
here](http://adsabs.harvard.edu/abs/2018arXiv180205450A) and the version
release:
If you use `PyFstat` in a publication we would appreciate if you cite both the
original paper introducing the code
(the [ADS page can be found here](http://adsabs.harvard.edu/abs/2018arXiv180205450A))
and the DOI for the software itself.
The latest version released on Zenodo from this repository here was v1.3,
see the following Bibtex entry.
Or see the [new github repository](https://github.com/PyFstat/PyFstat)
for the latest information.
```
@misc{pyfstat,
author = {{Ashton}, G. and {Keitel}, D.},
title = {{PyFstat-v1.2}},
month = may,
year = 2018,
doi = {10.5281/zenodo.1243931},
url = {https://doi.org/10.5281/zenodo.1243931},
note= {\url{https://doi.org/10.5281/zenodo.1243931}}
author = {Ashton, Gregory and
Keitel, David and
Prix, Reinhard},
title = {PyFstat-v1.3},
month = jan,
year = 2020,
publisher = {Zenodo},
version = {1.3},
doi = {10.5281/zenodo.3620861},
url = {https://doi.org/10.5281/zenodo.3620861}
note = {\url{https://doi.org/10.5281/zenodo.3620861}}
}
```
......
examples/MCMC_examples/fully_coherent_search_using_MCMC.py
View file @ d8571a7d
......@@ -13,38 +13,45 @@ F1 = -1e-10
F2 = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = .5*(tstart+tend)
tref = 0.5 * (tstart + tend)
depth = 10
h0 = sqrtSX / depth
label = 'fully_coherent_search_using_MCMC'
outdir = 'data'
label = "fully_coherent_search_using_MCMC"
outdir = "data"
data = pyfstat.Writer(
label=label, outdir=outdir, tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX)
label=label,
outdir=outdir,
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
)
data.make_data()
# The predicted twoF, given by lalapps_predictFstat can be accessed by
twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
print("Predicted twoF value: {}\n".format(twoF))
DeltaF0 = 1e-7
DeltaF1 = 1e-13
VF0 = (np.pi * duration * DeltaF0) ** 2 / 3.0
VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
VF1 = (np.pi * duration ** 2 * DeltaF1) ** 2 * 4 / 45.0
print("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))
theta_prior = {'F0': {'type': 'unif',
'lower': F0-DeltaF0/2.,
'upper': F0+DeltaF0/2.},
'F1': {'type': 'unif',
'lower': F1-DeltaF1/2.,
'upper': F1+DeltaF1/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
theta_prior = {
"F0": {"type": "unif", "lower": F0 - DeltaF0 / 2.0, "upper": F0 + DeltaF0 / 2.0},
"F1": {"type": "unif", "lower": F1 - DeltaF1 / 2.0, "upper": F1 + DeltaF1 / 2.0},
"F2": F2,
"Alpha": Alpha,
"Delta": Delta,
}
ntemps = 2
......@@ -53,13 +60,22 @@ nwalkers = 100
nsteps = [300, 300]
mcmc = pyfstat.MCMCSearch(
label=label, outdir=outdir,
sftfilepattern='{}/*{}*sft'.format(outdir, label), theta_prior=theta_prior,
tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,
nwalkers=nwalkers, ntemps=ntemps, log10beta_min=log10beta_min)
label=label,
outdir=outdir,
sftfilepattern="{}/*{}*sft".format(outdir, label),
theta_prior=theta_prior,
tref=tref,
minStartTime=tstart,
maxStartTime=tend,
nsteps=nsteps,
nwalkers=nwalkers,
ntemps=ntemps,
log10beta_min=log10beta_min,
)
mcmc.transform_dictionary = dict(
F0=dict(subtractor=F0, symbol='$f-f^\mathrm{s}$'),
F1=dict(subtractor=F1, symbol='$\dot{f}-\dot{f}^\mathrm{s}$'))
F0=dict(subtractor=F0, symbol="$f-f^\mathrm{s}$"),
F1=dict(subtractor=F1, symbol="$\dot{f}-\dot{f}^\mathrm{s}$"),
)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
examples/MCMC_examples/semi_coherent_search_using_MCMC.py
View file @ d8571a7d
......@@ -13,38 +13,45 @@ F1 = -1e-10
F2 = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = .5*(tstart+tend)
tref = 0.5 * (tstart + tend)
depth = 10
h0 = sqrtSX / depth
label = 'semicoherent_search_using_MCMC'
outdir = 'data'
label = "semicoherent_search_using_MCMC"
outdir = "data"
data = pyfstat.Writer(
label=label, outdir=outdir, tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX)
label=label,
outdir=outdir,
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
)
data.make_data()
# The predicted twoF, given by lalapps_predictFstat can be accessed by
twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
print("Predicted twoF value: {}\n".format(twoF))
DeltaF0 = 1e-7
DeltaF1 = 1e-13
VF0 = (np.pi * duration * DeltaF0) ** 2 / 3.0
VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
VF1 = (np.pi * duration ** 2 * DeltaF1) ** 2 * 4 / 45.0
print("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))
theta_prior = {'F0': {'type': 'unif',
'lower': F0-DeltaF0/2.,
'upper': F0+DeltaF0/2.},
'F1': {'type': 'unif',
'lower': F1-DeltaF1/2.,
'upper': F1+DeltaF1/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
theta_prior = {
"F0": {"type": "unif", "lower": F0 - DeltaF0 / 2.0, "upper": F0 + DeltaF0 / 2.0},
"F1": {"type": "unif", "lower": F1 - DeltaF1 / 2.0, "upper": F1 + DeltaF1 / 2.0},
"F2": F2,
"Alpha": Alpha,
"Delta": Delta,
}
ntemps = 1
......@@ -53,14 +60,23 @@ nwalkers = 100
nsteps = [300, 300]
mcmc = pyfstat.MCMCSemiCoherentSearch(
label=label, outdir=outdir, nsegs=10,
sftfilepattern='{}/*{}*sft'.format(outdir, label),
theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
label=label,
outdir=outdir,
nsegs=10,
sftfilepattern="{}/*{}*sft".format(outdir, label),
theta_prior=theta_prior,
tref=tref,
minStartTime=tstart,
maxStartTime=tend,
nsteps=nsteps,
nwalkers=nwalkers,
ntemps=ntemps,
log10beta_min=log10beta_min,
)
mcmc.transform_dictionary = dict(
F0=dict(subtractor=F0, symbol='$f-f^\mathrm{s}$'),
F1=dict(subtractor=F1, symbol='$\dot{f}-\dot{f}^\mathrm{s}$'))
F0=dict(subtractor=F0, symbol="$f-f^\mathrm{s}$"),
F1=dict(subtractor=F1, symbol="$\dot{f}-\dot{f}^\mathrm{s}$"),
)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
examples/followup_examples/semi_coherent_directed_follow_up.py
View file @ d8571a7d
......@@ -13,36 +13,44 @@ sqrtSX = 1e-23
tstart = 1000000000
duration = 100 * 86400
tend = tstart + duration
tref = .5*(tstart+tend)
tref = 0.5 * (tstart + tend)
depth = 40
label = 'semicoherent_directed_follow_up'
outdir = 'data'
label = "semicoherent_directed_follow_up"
outdir = "data"
h0 = sqrtSX / depth
data = pyfstat.Writer(
label=label, outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX)
label=label,
outdir=outdir,
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
)
data.make_data()
# The predicted twoF, given by lalapps_predictFstat can be accessed by
twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
print("Predicted twoF value: {}\n".format(twoF))
# Search
VF0 = VF1 = 1e5
DeltaF0 = np.sqrt(VF0) * np.sqrt(3) / (np.pi * duration)
DeltaF1 = np.sqrt(VF1) * np.sqrt(180) / (np.pi * duration ** 2)
theta_prior = {'F0': {'type': 'unif',
'lower': F0-DeltaF0/2.,
'upper': F0+DeltaF0/2},
'F1': {'type': 'unif',
'lower': F1-DeltaF1/2.,
'upper': F1+DeltaF1/2},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
theta_prior = {
"F0": {"type": "unif", "lower": F0 - DeltaF0 / 2.0, "upper": F0 + DeltaF0 / 2},
"F1": {"type": "unif", "lower": F1 - DeltaF1 / 2.0, "upper": F1 + DeltaF1 / 2},
"F2": F2,
"Alpha": Alpha,
"Delta": Delta,
}
ntemps = 3
......@@ -51,23 +59,35 @@ nwalkers = 100
nsteps = [100, 100]
mcmc = pyfstat.MCMCFollowUpSearch(
label=label, outdir=outdir,
sftfilepattern='{}/*{}*sft'.format(outdir, label),
theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
nwalkers=nwalkers, nsteps=nsteps, ntemps=ntemps,
log10beta_min=log10beta_min)
label=label,
outdir=outdir,
sftfilepattern="{}/*{}*sft".format(outdir, label),
theta_prior=theta_prior,
tref=tref,
minStartTime=tstart,
maxStartTime=tend,
nwalkers=nwalkers,
nsteps=nsteps,
ntemps=ntemps,
log10beta_min=log10beta_min,
)
NstarMax = 1000
Nsegs0 = 100
fig, axes = plt.subplots(nrows=2, figsize=(3.4, 3.5))
fig, axes = mcmc.run(
NstarMax=NstarMax, Nsegs0=Nsegs0, labelpad=0.01,
plot_det_stat=False, return_fig=True, fig=fig,
axes=axes)
NstarMax=NstarMax,
Nsegs0=Nsegs0,
labelpad=0.01,
plot_det_stat=False,
return_fig=True,
fig=fig,
axes=axes,
)
for ax in axes:
ax.grid()
ax.set_xticks(np.arange(0, 600, 100))
ax.set_xticklabels([str(s) for s in np.arange(0, 700, 100)])
axes[-1].set_xlabel(r'$\textrm{Number of steps}$', labelpad=0.1)
axes[-1].set_xlabel(r"$\textrm{Number of steps}$", labelpad=0.1)
fig.tight_layout()
fig.savefig('{}/{}_walkers.png'.format(mcmc.outdir, mcmc.label), dpi=400)
fig.savefig("{}/{}_walkers.png".format(mcmc.outdir, mcmc.label), dpi=400)
examples/glitch_examples/make_simulated_data.py
View file @ d8571a7d
from pyfstat import Writer, GlitchWriter
import numpy as np
outdir = 'data'
outdir = "data"
# First, we generate data with a reasonably strong smooth signal
# Define parameters of the Crab pulsar as an example
......@@ -22,8 +22,19 @@ tend = tstart+duration
tref = tstart + 0.5 * duration
data = Writer(
label='0_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX)
label="0_glitch",
outdir=outdir,
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
)
data.make_data()
# Next, taking the same signal parameters, we include a glitch half way through
......@@ -32,9 +43,22 @@ delta_F0 = 5e-6
delta_F1 = 0
glitch_data = GlitchWriter(
label='1_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)
label="1_glitch",
outdir=outdir,
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
dtglitch=dtglitch,
delta_F0=delta_F0,
delta_F1=delta_F1,
)
glitch_data.make_data()
# Making data with two glitches
......@@ -46,8 +70,22 @@ delta_F1_2 = [0, 0]
delta_F2_2 = [0, 0]
two_glitch_data = GlitchWriter(
label='2_glitch', outdir=outdir, tref=tref, tstart=tstart, F0=F0, F1=F1,
F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,
dtglitch=dtglitch_2, delta_phi=delta_phi_2, delta_F0=delta_F0_2,
delta_F1=delta_F1_2, delta_F2=delta_F2_2)
label="2_glitch",
outdir=outdir,
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
dtglitch=dtglitch_2,
delta_phi=delta_phi_2,
delta_F0=delta_F0_2,
delta_F1=delta_F1_2,
delta_F2=delta_F2_2,
)
two_glitch_data.make_data()
examples/glitch_examples/semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch.py
View file @ d8571a7d
......@@ -3,33 +3,41 @@ import matplotlib.pyplot as plt
import pyfstat
import gridcorner
import time
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, delta_F0, dtglitch, outdir
from make_simulated_data import (
tstart,
duration,
tref,
F0,
F1,
F2,
Alpha,
Delta,
delta_F0,
dtglitch,
outdir,
)
plt.style.use('./paper.mplstyle')
plt.style.use("./paper.mplstyle")
label = 'semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch'
label = "semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch"
Nstar = 1000
F0_width = np.sqrt(Nstar) * np.sqrt(12) / (np.pi * duration)
F1_width = np.sqrt(Nstar) * np.sqrt(180) / (np.pi * duration ** 2)
theta_prior = {
'F0': {'type': 'unif',
'lower': F0-F0_width/2.,
'upper': F0+F0_width/2.},
'F1': {'type': 'unif',
'lower': F1-F1_width/2.,
'upper': F1+F1_width/2.},
'F2': F2,
'delta_F0': {'type': 'unif',
'lower': 0,
'upper': 1e-5},
'delta_F1': 0,
'tglitch': {'type': 'unif',
'lower': tstart+0.1*duration,
'upper': tstart+0.9*duration},
'Alpha': Alpha,
'Delta': Delta,
"F0": {"type": "unif", "lower": F0 - F0_width / 2.0, "upper": F0 + F0_width / 2.0},
"F1": {"type": "unif", "lower": F1 - F1_width / 2.0, "upper": F1 + F1_width / 2.0},
"F2": F2,
"delta_F0": {"type": "unif", "lower": 0, "upper": 1e-5},
"delta_F1": 0,
"tglitch": {
"type": "unif",
"lower": tstart + 0.1 * duration,
"upper": tstart + 0.9 * duration,
},
"Alpha": Alpha,
"Delta": Delta,
}
ntemps = 3
......@@ -38,33 +46,49 @@ nwalkers = 100
nsteps = [250, 250]
mcmc = pyfstat.MCMCGlitchSearch(
label=label, sftfilepattern='data/*1_glitch*sft', theta_prior=theta_prior,
tref=tref, minStartTime=tstart, maxStartTime=tstart+duration,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min, nglitch=1)
mcmc.transform_dictionary['F0'] = dict(
subtractor=F0, multiplier=1e6, symbol='$f-f_\mathrm{s}$')
mcmc.unit_dictionary['F0'] = '$\mu$Hz'
mcmc.transform_dictionary['F1'] = dict(
subtractor=F1, multiplier=1e12, symbol='$\dot{f}-\dot{f}_\mathrm{s}$')
mcmc.unit_dictionary['F1'] = '$p$Hz/s'
mcmc.transform_dictionary['delta_F0'] = dict(
multiplier=1e6, subtractor=delta_F0,
symbol='$\delta f-\delta f_\mathrm{s}$')
mcmc.unit_dictionary['delta_F0'] = '$\mu$Hz/s'
mcmc.transform_dictionary['tglitch']['subtractor'] = tstart + dtglitch
mcmc.transform_dictionary['tglitch']['label'] ='$t^\mathrm{g}-t^\mathrm{g}_\mathrm{s}$\n[d]'
label=label,
sftfilepattern="data/*1_glitch*sft",
theta_prior=theta_prior,
tref=tref,
minStartTime=tstart,
maxStartTime=tstart + duration,
nsteps=nsteps,
nwalkers=nwalkers,
ntemps=ntemps,
log10beta_min=log10beta_min,
nglitch=1,
)
mcmc.transform_dictionary["F0"] = dict(
subtractor=F0, multiplier=1e6, symbol="$f-f_\mathrm{s}$"
)
mcmc.unit_dictionary["F0"] = "$\mu$Hz"
mcmc.transform_dictionary["F1"] = dict(
subtractor=F1, multiplier=1e12, symbol="$\dot{f}-\dot{f}_\mathrm{s}$"
)
mcmc.unit_dictionary["F1"] = "$p$Hz/s"
mcmc.transform_dictionary["delta_F0"] = dict(
multiplier=1e6, subtractor=delta_F0, symbol="$\delta f-\delta f_\mathrm{s}$"
)
mcmc.unit_dictionary["delta_F0"] = "$\mu$Hz/s"
mcmc.transform_dictionary["tglitch"]["subtractor"] = tstart + dtglitch
mcmc.transform_dictionary["tglitch"][
"label"
] = "$t^\mathrm{g}-t^\mathrm{g}_\mathrm{s}$\n[d]"
t1 = time.time()
mcmc.run()
dT = time.time() - t1
fig_and_axes = gridcorner._get_fig_and_axes(4, 2, 0.05)
mcmc.plot_corner(label_offset=0.25, truths=[0, 0, 0, 0],
fig_and_axes=fig_and_axes, quantiles=(0.16, 0.84),
mcmc.plot_corner(
label_offset=0.25,
truths=[0, 0, 0, 0],
fig_and_axes=fig_and_axes,
quantiles=(0.16, 0.84),
hist_kwargs=dict(lw=1.5, zorder=-1),
truth_color='C3')
truth_color="C3",
)
mcmc.print_summary()
print('Prior widths =', F0_width, F1_width)
print("Actual run time = {}".format(dT))
print(("Prior widths =", F0_width, F1_width))
print(("Actual run time = {}".format(dT)))
examples/glitch_examples/semicoherent_glitch_robust_directed_grid_search_on_1_glitch.py
View file @ d8571a7d
import pyfstat
import numpy as np
import matplotlib.pyplot as plt
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, delta_F0, outdir, dtglitch
from make_simulated_data import (
tstart,
duration,
tref,
F0,
F1,
F2,
Alpha,
Delta,
delta_F0,
outdir,
dtglitch,
)
import time
try:
......@@ -9,18 +21,19 @@ try:
except ImportError:
raise ImportError(
"Python module 'gridcorner' not found, please install from "
"https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner")
"https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner"
)
label = 'semicoherent_glitch_robust_directed_grid_search_on_1_glitch'
label = "semicoherent_glitch_robust_directed_grid_search_on_1_glitch"
plt.style.use('./paper.mplstyle')
plt.style.use("./paper.mplstyle")
Nstar = 1000
F0_width = np.sqrt(Nstar) * np.sqrt(12) / (np.pi * duration)
F1_width = np.sqrt(Nstar) * np.sqrt(180) / (np.pi * duration ** 2)
N = 20
F0s = [F0-F0_width/2., F0+F0_width/2., F0_width/N]
F1s = [F1-F1_width/2., F1+F1_width/2., F1_width/N]
F0s = [F0 - F0_width / 2.0, F0 + F0_width / 2.0, F0_width / N]
F1s = [F1 - F1_width / 2.0, F1 + F1_width / 2.0, F1_width / N]
F2s = [F2]
Alphas = [Alpha]
Deltas = [Delta]
......@@ -33,10 +46,21 @@ delta_F1s = [0]
t1 = time.time()
search = pyfstat.GridGlitchSearch(
label, outdir, 'data/*1_glitch*sft', F0s=F0s, F1s=F1s, F2s=F2s,
Alphas=Alphas, Deltas=Deltas, tref=tref, minStartTime=tstart,
maxStartTime=tstart+duration, tglitchs=tglitchs, delta_F0s=delta_F0s,
delta_F1s=delta_F1s)
label,
outdir,
"data/*1_glitch*sft",
F0s=F0s,
F1s=F1s,
F2s=F2s,
Alphas=Alphas,
Deltas=Deltas,
tref=tref,
minStartTime=tstart,
maxStartTime=tstart + duration,
tglitchs=tglitchs,
delta_F0s=delta_F0s,
delta_F1s=delta_F1s,
)
search.run()
dT = time.time() - t1
......@@ -44,24 +68,32 @@ F0_vals = np.unique(search.data[:, 0]) - F0
F1_vals = np.unique(search.data[:, 1]) - F1
delta_F0s_vals = np.unique(search.data[:, 5]) - delta_F0
tglitch_vals = np.unique(search.data[:, 7])
tglitch_vals_days = (tglitch_vals-tstart) / 86400. - dtglitch/86400.
tglitch_vals_days = (tglitch_vals - tstart) / 86400.0 - dtglitch / 86400.0
twoF = search.data[:, -1].reshape((len(F0_vals), len(F1_vals),
len(delta_F0s_vals), len(tglitch_vals)))
twoF = search.data[:, -1].reshape(
(len(F0_vals), len(F1_vals), len(delta_F0s_vals), len(tglitch_vals))
)
xyz = [F0_vals * 1e6, F1_vals * 1e12, delta_F0s_vals * 1e6, tglitch_vals_days]
labels = ['$f - f_\mathrm{s}$\n[$\mu$Hz]',
'$\dot{f} - \dot{f}_\mathrm{s}$\n[$p$Hz/s]',
'$\delta f-\delta f_\mathrm{s}$\n[$\mu$Hz]',
'$t^\mathrm{g} - t^\mathrm{g}_\mathrm{s}$\n[d]',
'$\widehat{2\mathcal{F}}$']
labels = [
"$f - f_\mathrm{s}$\n[$\mu$Hz]",
"$\dot{f} - \dot{f}_\mathrm{s}$\n[$p$Hz/s]",
"$\delta f-\delta f_\mathrm{s}$\n[$\mu$Hz]",
"$t^\mathrm{g} - t^\mathrm{g}_\mathrm{s}$\n[d]",
"$\widehat{2\mathcal{F}}$",
]
fig, axes = gridcorner(
twoF, xyz, projection='log_mean', labels=labels,
showDvals=False, lines=[0, 0, 0, 0], label_offset=0.25,
max_n_ticks=4)
fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label),
bbox_inches='tight')
twoF,
xyz,
projection="log_mean",
labels=labels,
showDvals=False,
lines=[0, 0, 0, 0],
label_offset=0.25,
max_n_ticks=4,
)
fig.savefig("{}/{}_projection_matrix.png".format(outdir, label), bbox_inches="tight")
print('Prior widths =', F0_width, F1_width)
print("Actual run time = {}".format(dT))
print("Actual number of grid points = {}".format(search.data.shape[0]))
print(("Prior widths =", F0_width, F1_width))
print(("Actual run time = {}".format(dT)))
print(("Actual number of grid points = {}".format(search.data.shape[0])))
examples/glitch_examples/standard_directed_MCMC_search_on_1_glitch.py
View file @ d8571a7d
......@@ -3,24 +3,20 @@ import matplotlib.pyplot as plt
import pyfstat
from make_simulated_data import tstart, duration, tref, F0, F1, F2, Alpha, Delta, outdir
plt.style.use('paper')
plt.style.use("paper")
label = 'standard_directed_MCMC_search_on_1_glitch'
label = "standard_directed_MCMC_search_on_1_glitch"
Nstar = 10000
F0_width = np.sqrt(Nstar) * np.sqrt(12) / (np.pi * duration)
F1_width = np.sqrt(Nstar) * np.sqrt(180) / (np.pi * duration ** 2)
theta_prior = {
'F0': {'type': 'unif',
'lower': F0-F0_width/2.,
'upper': F0+F0_width/2.},
'F1': {'type': 'unif',
'lower': F1-F1_width/2.,
'upper': F1+F1_width/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
"F0": {"type": "unif", "lower": F0 - F0_width / 2.0, "upper": F0 + F0_width / 2.0},
"F1": {"type": "unif", "lower": F1 - F1_width / 2.0, "upper": F1 + F1_width / 2.0},
"F2": F2,
"Alpha": Alpha,
"Delta": Delta,
}
ntemps = 2
......@@ -29,15 +25,22 @@ nwalkers = 100
nsteps = [500, 2000]
mcmc = pyfstat.MCMCSearch(
label=label, sftfilepattern='data/*1_glitch*sft', theta_prior=theta_prior,
tref=tref, minStartTime=tstart, maxStartTime=tstart+duration,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
mcmc.transform_dictionary['F0'] = dict(
subtractor=F0, symbol='$f-f^\mathrm{s}$')
mcmc.transform_dictionary['F1'] = dict(
subtractor=F1, symbol='$\dot{f}-\dot{f}^\mathrm{s}$')
label=label,
sftfilepattern="data/*1_glitch*sft",
theta_prior=theta_prior,
tref=tref,
minStartTime=tstart,
maxStartTime=tstart + duration,
nsteps=nsteps,
nwalkers=nwalkers,
ntemps=ntemps,
log10beta_min=log10beta_min,
)
mcmc.transform_dictionary["F0"] = dict(subtractor=F0, symbol="$f-f^\mathrm{s}$")
mcmc.transform_dictionary["F1"] = dict(
subtractor=F1, symbol="$\dot{f}-\dot{f}^\mathrm{s}$"
)
mcmc.run()
mcmc.plot_corner()
......
examples/grid_examples/grid_F0F1F2.py
View file @ d8571a7d
......@@ -7,7 +7,8 @@ try:
except ImportError:
raise ImportError(
"Python module 'gridcorner' not found, please install from "
"https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner")
"https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner"
)
F0 = 30.0
F1 = 1e-10
......@@ -20,18 +21,28 @@ sqrtSX = 1e-23
tstart = 1000000000
duration = 10 * 86400
tend = tstart + duration
tref = .5*(tstart+tend)
tref = 0.5 * (tstart + tend)
depth = 20
label = 'grid_F0F1F2'
outdir = 'data'
label = "grid_F0F1F2"
outdir = "data"
h0 = sqrtSX / depth
data = pyfstat.Writer(
label=label, outdir=outdir, tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX)
label=label,
outdir=outdir,
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
)
data.make_data()
m = 0.01
......@@ -42,14 +53,24 @@ N = 100
DeltaF0 = N * dF0
DeltaF1 = N * dF1
DeltaF2 = N * dF2
F0s = [F0-DeltaF0/2., F0+DeltaF0/2., dF0]
F1s = [F1-DeltaF1/2., F1+DeltaF1/2., dF1]
F2s = [F2-DeltaF2/2., F2+DeltaF2/2., dF2]
F0s = [F0 - DeltaF0 / 2.0, F0 + DeltaF0 / 2.0, dF0]
F1s = [F1 - DeltaF1 / 2.0, F1 + DeltaF1 / 2.0, dF1]
F2s = [F2 - DeltaF2 / 2.0, F2 + DeltaF2 / 2.0, dF2]
Alphas = [Alpha]
Deltas = [Delta]
search = pyfstat.GridSearch(
'grid_F0F1F2', 'data', data.sftfilepath, F0s, F1s,
F2s, Alphas, Deltas, tref, tstart, tend)
"grid_F0F1F2",
"data",
data.sftfilepath,
F0s,
F1s,
F2s,
Alphas,
Deltas,
tref,
tstart,
tend,
)
search.run()
F0_vals = np.unique(search.data[:, 2]) - F0
......@@ -57,8 +78,13 @@ F1_vals = np.unique(search.data[:, 3]) - F1
F2_vals = np.unique(search.data[:, 4]) - F2
twoF = search.data[:, -1].reshape((len(F0_vals), len(F1_vals), len(F2_vals)))
xyz = [F0_vals, F1_vals, F2_vals]
labels = ['$f - f_0$', '$\dot{f} - \dot{f}_0$', '$\ddot{f} - \ddot{f}_0$',
'$\widetilde{2\mathcal{F}}$']
labels = [
"$f - f_0$",
"$\dot{f} - \dot{f}_0$",
"$\ddot{f} - \ddot{f}_0$",
"$\widetilde{2\mathcal{F}}$",
]
fig, axes = gridcorner(
twoF, xyz, projection='log_mean', labels=labels, whspace=0.1, factor=1.8)
fig.savefig('{}/{}_projection_matrix.png'.format(outdir, label))
twoF, xyz, projection="log_mean", labels=labels, whspace=0.1, factor=1.8
)
fig.savefig("{}/{}_projection_matrix.png".format(outdir, label))
examples/grid_examples/grided_frequency_search.py
View file @ d8571a7d
......@@ -13,50 +13,70 @@ sqrtSX = 1e-23
tstart = 1000000000
duration = 100 * 86400
tend = tstart + duration
tref = .5*(tstart+tend)
tref = 0.5 * (tstart + tend)
depth = 70
data_label = 'grided_frequency_depth_{:1.0f}'.format(depth)
data_label = "grided_frequency_depth_{:1.0f}".format(depth)
h0 = sqrtSX / depth
data = pyfstat.Writer(
label=data_label, outdir='data', tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX)
label=data_label,
outdir="data",
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
)
data.make_data()
m = 0.001
dF0 = np.sqrt(12 * m) / (np.pi * duration)
DeltaF0 = 800 * dF0
F0s = [F0-DeltaF0/2., F0+DeltaF0/2., dF0]
F0s = [F0 - DeltaF0 / 2.0, F0 + DeltaF0 / 2.0, dF0]
F1s = [F1]
F2s = [F2]
Alphas = [Alpha]
Deltas = [Delta]
search = pyfstat.GridSearch(
'grided_frequency_search', 'data', 'data/*'+data_label+'*sft', F0s, F1s,
F2s, Alphas, Deltas, tref, tstart, tend)
"grided_frequency_search",
"data",
"data/*" + data_label + "*sft",
F0s,
F1s,
F2s,
Alphas,
Deltas,
tref,
tstart,
tend,
)
search.run()
fig, ax = plt.subplots()
xidx = search.keys.index('F0')
xidx = search.keys.index("F0")
frequencies = np.unique(search.data[:, xidx])
twoF = search.data[:, -1]
# mismatch = np.sign(x-F0)*(duration * np.pi * (x - F0))**2 / 12.0
ax.plot(frequencies, twoF, 'k', lw=1)
ax.plot(frequencies, twoF, "k", lw=1)
DeltaF = frequencies - F0
sinc = np.sin(np.pi * DeltaF * duration) / (np.pi * DeltaF * duration)
A = np.abs((np.max(twoF) - 4) * sinc ** 2 + 4)
ax.plot(frequencies, A, '-r', lw=1)
ax.set_ylabel('$\widetilde{2\mathcal{F}}$')
ax.set_xlabel('Frequency')
ax.plot(frequencies, A, "-r", lw=1)
ax.set_ylabel("$\widetilde{2\mathcal{F}}$")
ax.set_xlabel("Frequency")
ax.set_xlim(F0s[0], F0s[1])
dF0 = np.sqrt(12 * 1) / (np.pi * duration)
xticks = [F0 - 10 * dF0, F0, F0 + 10 * dF0]
ax.set_xticks(xticks)
xticklabels = ['$f_0 {-} 10\Delta f$', '$f_0$', '$f_0 {+} 10\Delta f$']
xticklabels = ["$f_0 {-} 10\Delta f$", "$f_0$", "$f_0 {+} 10\Delta f$"]
ax.set_xticklabels(xticklabels)
plt.tight_layout()
fig.savefig('{}/{}_1D.png'.format(search.outdir, search.label), dpi=300)
fig.savefig("{}/{}_1D.png".format(search.outdir, search.label), dpi=300)
examples/other_examples/sliding_window.py
View file @ d8571a7d
......@@ -13,23 +13,43 @@ F1 = -1e-10
F2 = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = .5*(tstart+tend)
tref = 0.5 * (tstart + tend)
depth = 60
h0 = sqrtSX / depth
data_label = 'sliding_window'
data_label = "sliding_window"
data = pyfstat.Writer(
label=data_label, outdir='data', tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX)
label=data_label,
outdir="data",
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
)
data.make_data()
DeltaF0 = 1e-5
search = pyfstat.FrequencySlidingWindow(
label='sliding_window', outdir='data', sftfilepattern='data/*sliding_window*sft',
F0s=[F0-DeltaF0, F0+DeltaF0, DeltaF0/100.], F1=F1, F2=0,
Alpha=Alpha, Delta=Delta, tref=tref, minStartTime=tstart,
maxStartTime=tend, window_size=25*86400, window_delta=1*86400)
label="sliding_window",
outdir="data",
sftfilepattern="data/*sliding_window*sft",
F0s=[F0 - DeltaF0, F0 + DeltaF0, DeltaF0 / 100.0],
F1=F1,
F2=0,
Alpha=Alpha,
Delta=Delta,
tref=tref,
minStartTime=tstart,
maxStartTime=tend,
window_size=25 * 86400,
window_delta=1 * 86400,
)
search.run()
search.plot_sliding_window()
examples/other_examples/twoF_cumulative.py
View file @ d8571a7d
......@@ -13,37 +13,45 @@ F1 = -1e-10
F2 = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = .5*(tstart+tend)
tref = 0.5 * (tstart + tend)
depth = 100
h0 = sqrtSX / depth
data_label = 'twoF_cumulative'
data_label = "twoF_cumulative"
data = pyfstat.Writer(
label=data_label, outdir='data', tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX, detectors='H1,L1')
label=data_label,
outdir="data",
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
detectors="H1,L1",
)
data.make_data()
# The predicted twoF, given by lalapps_predictFstat can be accessed by
twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
print("Predicted twoF value: {}\n".format(twoF))
DeltaF0 = 1e-7
DeltaF1 = 1e-13
VF0 = (np.pi * duration * DeltaF0) ** 2 / 3.0
VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
VF1 = (np.pi * duration ** 2 * DeltaF1) ** 2 * 4 / 45.0
print("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))
theta_prior = {'F0': {'type': 'unif',
'lower': F0-DeltaF0/2.,
'upper': F0+DeltaF0/2.},
'F1': {'type': 'unif',
'lower': F1-DeltaF1/2.,
'upper': F1+DeltaF1/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
theta_prior = {
"F0": {"type": "unif", "lower": F0 - DeltaF0 / 2.0, "upper": F0 + DeltaF0 / 2.0},
"F1": {"type": "unif", "lower": F1 - DeltaF1 / 2.0, "upper": F1 + DeltaF1 / 2.0},
"F2": F2,
"Alpha": Alpha,
"Delta": Delta,
}
ntemps = 1
......@@ -52,10 +60,18 @@ nwalkers = 100
nsteps = [50, 50]
mcmc = pyfstat.MCMCSearch(
label='twoF_cumulative', outdir='data',
sftfilepattern='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,
minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,
ntemps=ntemps, log10beta_min=log10beta_min)
label="twoF_cumulative",
outdir="data",
sftfilepattern="data/*" + data_label + "*sft",
theta_prior=theta_prior,
tref=tref,
minStartTime=tstart,
maxStartTime=tend,
nsteps=nsteps,
nwalkers=nwalkers,
ntemps=ntemps,
log10beta_min=log10beta_min,
)
mcmc.run()
mcmc.plot_corner(add_prior=True)
mcmc.print_summary()
......
examples/other_examples/using_initialisation.py
View file @ d8571a7d
......@@ -13,38 +13,45 @@ F1 = -1e-10
F2 = 0
Alpha = np.radians(83.6292)
Delta = np.radians(22.0144)
tref = .5*(tstart+tend)
tref = 0.5 * (tstart + tend)
depth = 10
h0 = sqrtSX / depth
label = 'using_initialisation'
outdir = 'data'
label = "using_initialisation"
outdir = "data"
data = pyfstat.Writer(
label=label, outdir=outdir, tref=tref,
tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
Delta=Delta, h0=h0, sqrtSX=sqrtSX)
label=label,
outdir=outdir,
tref=tref,
tstart=tstart,
F0=F0,
F1=F1,
F2=F2,
duration=duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
)
data.make_data()
# The predicted twoF, given by lalapps_predictFstat can be accessed by
twoF = data.predict_fstat()
print 'Predicted twoF value: {}\n'.format(twoF)
print("Predicted twoF value: {}\n".format(twoF))
DeltaF0 = 1e-7
DeltaF1 = 1e-13
VF0 = (np.pi * duration * DeltaF0) ** 2 / 3.0
VF1 = (np.pi * duration**2 * DeltaF1)**2 * 4/45.
print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
VF1 = (np.pi * duration ** 2 * DeltaF1) ** 2 * 4 / 45.0
print("\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n".format(VF0 * VF1, VF0, VF1))
theta_prior = {'F0': {'type': 'unif',
'lower': F0-DeltaF0/2.,
'upper': F0+DeltaF0/2.},
'F1': {'type': 'unif',
'lower': F1-DeltaF1/2.,
'upper': F1+DeltaF1/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta
theta_prior = {
"F0": {"type": "unif", "lower": F0 - DeltaF0 / 2.0, "upper": F0 + DeltaF0 / 2.0},
"F1": {"type": "unif", "lower": F1 - DeltaF1 / 2.0, "upper": F1 + DeltaF1 / 2.0},
"F2": F2,
"Alpha": Alpha,
"Delta": Delta,
}
ntemps = 1
......@@ -53,11 +60,18 @@ nwalkers = 100
nsteps = [100, 100]
mcmc = pyfstat.MCMCSearch(
label=label, outdir=outdir,
sftfilepattern='{}/*{}*sft'.format(outdir, label),
theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
log10beta_min=log10beta_min)
label=label,
outdir=outdir,
sftfilepattern="{}/*{}*sft".format(outdir, label),
theta_prior=theta_prior,
tref=tref,
minStartTime=tstart,
maxStartTime=tend,
nsteps=nsteps,
nwalkers=nwalkers,
ntemps=ntemps,
log10beta_min=log10beta_min,
)
mcmc.setup_initialisation(100, scatter_val=1e-10)
mcmc.run()
mcmc.plot_corner(add_prior=True)
......
examples/transient_examples/long_transient_search_MCMC.py
View file @ d8571a7d
......@@ -16,19 +16,18 @@ tref = minStartTime
DeltaF0 = 6e-7
DeltaF1 = 1e-13
theta_prior = {'F0': {'type': 'unif',
'lower': F0-DeltaF0/2.,
'upper': F0+DeltaF0/2.},
'F1': {'type': 'unif',
'lower': F1-DeltaF1/2.,
'upper': F1+DeltaF1/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
'transient_tstart': minStartTime,
'transient_duration': {'type': 'halfnorm',
'loc': 0.001*Tspan,
'scale': 0.5*Tspan}
theta_prior = {
"F0": {"type": "unif", "lower": F0 - DeltaF0 / 2.0, "upper": F0 + DeltaF0 / 2.0},
"F1": {"type": "unif", "lower": F1 - DeltaF1 / 2.0, "upper": F1 + DeltaF1 / 2.0},
"F2": F2,
"Alpha": Alpha,
"Delta": Delta,
"transient_tstart": minStartTime,
"transient_duration": {
"type": "halfnorm",
"loc": 0.001 * Tspan,
"scale": 0.5 * Tspan,
},
}
ntemps = 2
......@@ -37,12 +36,19 @@ nwalkers = 100
nsteps = [100, 100]
mcmc = pyfstat.MCMCTransientSearch(
label='transient_search', outdir='data_l',
sftfilepattern='data_l/*simulated_transient_signal*sft',
theta_prior=theta_prior, tref=tref, minStartTime=minStartTime,
maxStartTime=maxStartTime, nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
label="transient_search",
outdir="data_l",
sftfilepattern="data_l/*simulated_transient_signal*sft",
theta_prior=theta_prior,
tref=tref,
minStartTime=minStartTime,
maxStartTime=maxStartTime,
nsteps=nsteps,
nwalkers=nwalkers,
ntemps=ntemps,
log10beta_min=log10beta_min,
transientWindowType='rect')
transientWindowType="rect",
)
mcmc.run()
mcmc.plot_corner(label_offset=0.7)
mcmc.print_summary()
examples/transient_examples/long_transient_search_make_simulated_data.py
View file @ d8571a7d
......@@ -19,8 +19,20 @@ h0 = 1e-23
sqrtSX = 1e-22
transient = pyfstat.Writer(
label='simulated_transient_signal', outdir='data_l', tref=tref,
tstart=transient_tstart, F0=F0, F1=F1, F2=F2, duration=transient_duration,
Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX, minStartTime=minStartTime,
maxStartTime=maxStartTime, transientWindowType='rect')
label="simulated_transient_signal",
outdir="data_l",
tref=tref,
tstart=transient_tstart,
F0=F0,
F1=F1,
F2=F2,
duration=transient_duration,
Alpha=Alpha,
Delta=Delta,
h0=h0,
sqrtSX=sqrtSX,
minStartTime=minStartTime,
maxStartTime=maxStartTime,
transientWindowType="rect",
)
transient.make_data()
examples/transient_examples/short_transient_search_MCMC.py
View file @ d8571a7d
......@@ -18,21 +18,22 @@ Tsft = 1800
DeltaF0 = 1e-2
DeltaF1 = 1e-9
theta_prior = {'F0': {'type': 'unif',
'lower': F0-DeltaF0/2.,
'upper': F0+DeltaF0/2.},
'F1': {'type': 'unif',
'lower': F1-DeltaF1/2.,
'upper': F1+DeltaF1/2.},
'F2': F2,
'Alpha': Alpha,
'Delta': Delta,
'transient_tstart': {'type': 'unif',
'lower': minStartTime,
'upper': maxStartTime-2*Tsft},
'transient_duration': {'type': 'unif',
'lower': 2*Tsft,
'upper': Tspan-2*Tsft}
theta_prior = {
"F0": {"type": "unif", "lower": F0 - DeltaF0 / 2.0, "upper": F0 + DeltaF0 / 2.0},
"F1": {"type": "unif", "lower": F1 - DeltaF1 / 2.0, "upper": F1 + DeltaF1 / 2.0},
"F2": F2,
"Alpha": Alpha,
"Delta": Delta,
"transient_tstart": {
"type": "unif",
"lower": minStartTime,
"upper": maxStartTime - 2 * Tsft,
},
"transient_duration": {
"type": "unif",
"lower": 2 * Tsft,
"upper": Tspan - 2 * Tsft,
},
}
ntemps = 2
......@@ -41,12 +42,19 @@ nwalkers = 100
nsteps = [100, 100]
mcmc = pyfstat.MCMCTransientSearch(
label='transient_search', outdir='data_s',
sftfilepattern='data_s/*simulated_transient_signal*sft',
theta_prior=theta_prior, tref=tref, minStartTime=minStartTime,
maxStartTime=maxStartTime, nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
label="transient_search",
outdir="data_s",
sftfilepattern="data_s/*simulated_transient_signal*sft",
theta_prior=theta_prior,
tref=tref,
minStartTime=minStartTime,
maxStartTime=maxStartTime,
nsteps=nsteps,
nwalkers=nwalkers,
ntemps=ntemps,
log10beta_min=log10beta_min,
transientWindowType='rect')
transientWindowType="rect",
)
mcmc.run()
mcmc.plot_corner(label_offset=0.7)
mcmc.print_summary()
examples/transient_examples/short_transient_search_gridded.py
View file @ d8571a7d
......@@ -5,7 +5,7 @@ import os
import numpy as np
import matplotlib.pyplot as plt
datadir = 'data_s'
datadir = "data_s"
F0 = 30.0
F1 = -1e-10
......@@ -23,37 +23,53 @@ Tsft = 1800
m = 0.001
dF0 = np.sqrt(12 * m) / (np.pi * Tspan)
DeltaF0 = 100 * dF0
F0s = [F0-DeltaF0/2., F0+DeltaF0/2., dF0]
F0s = [F0 - DeltaF0 / 2.0, F0 + DeltaF0 / 2.0, dF0]
F1s = [F1]
F2s = [F2]
Alphas = [Alpha]
Deltas = [Delta]
print('Standard CW search:')
print("Standard CW search:")
search1 = pyfstat.GridSearch(
label='CW', outdir=datadir,
sftfilepattern=os.path.join(datadir,'*simulated_transient_signal*sft'),
F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref,
minStartTime=minStartTime, maxStartTime=maxStartTime,
BSGL=False)
label="CW",
outdir=datadir,
sftfilepattern=os.path.join(datadir, "*simulated_transient_signal*sft"),
F0s=F0s,
F1s=F1s,
F2s=F2s,
Alphas=Alphas,
Deltas=Deltas,
tref=tref,
minStartTime=minStartTime,
maxStartTime=maxStartTime,
BSGL=False,
)
search1.run()
search1.print_max_twoF()
search1.plot_1D(xkey='F0',
xlabel='freq [Hz]', ylabel='$2\mathcal{F}$')
search1.plot_1D(xkey="F0", xlabel="freq [Hz]", ylabel="$2\mathcal{F}$")
print('with t0,tau bands:')
print("with t0,tau bands:")
search2 = pyfstat.TransientGridSearch(
label='tCW', outdir=datadir,
sftfilepattern=os.path.join(datadir,'*simulated_transient_signal*sft'),
F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref,
minStartTime=minStartTime, maxStartTime=maxStartTime,
transientWindowType='rect', t0Band=Tspan-2*Tsft, tauBand=Tspan,
label="tCW",
outdir=datadir,
sftfilepattern=os.path.join(datadir, "*simulated_transient_signal*sft"),
F0s=F0s,
F1s=F1s,
F2s=F2s,
Alphas=Alphas,
Deltas=Deltas,
tref=tref,
minStartTime=minStartTime,
maxStartTime=maxStartTime,
transientWindowType="rect",
t0Band=Tspan - 2 * Tsft,
tauBand=Tspan,
BSGL=False,
outputTransientFstatMap=True,
tCWFstatMapVersion='lal')
tCWFstatMapVersion="lal",
)
search2.run()
search2.print_max_twoF()
search2.plot_1D(xkey='F0',
xlabel='freq [Hz]', ylabel='$2\mathcal{F}$')
search2.plot_1D(xkey="F0", xlabel="freq [Hz]", ylabel="$2\mathcal{F}$")
examples/transient_examples/short_transient_search_make_simulated_data.py
View file @ d8571a7d
......@@ -17,15 +17,27 @@ tref = minStartTime
h0 = 1e-23
sqrtSX = 1e-22
detectors = 'H1,L1'
detectors = "H1,L1"
Tsft = 1800
transient = pyfstat.Writer(
label='simulated_transient_signal', outdir='data_s',
tref=tref, tstart=transient_tstart, duration=transient_duration,
F0=F0, F1=F1, F2=F2, Alpha=Alpha, Delta=Delta, h0=h0,
detectors=detectors,sqrtSX=sqrtSX,
minStartTime=minStartTime, maxStartTime=maxStartTime,
transientWindowType='rect', Tsft=Tsft)
label="simulated_transient_signal",
outdir="data_s",
tref=tref,
tstart=transient_tstart,
duration=transient_duration,
F0=F0,
F1=F1,
F2=F2,
Alpha=Alpha,
Delta=Delta,
h0=h0,
detectors=detectors,
sqrtSX=sqrtSX,
minStartTime=minStartTime,
maxStartTime=maxStartTime,
transientWindowType="rect",
Tsft=Tsft,
)
transient.make_data()
pyfstat/__init__.py
View file @ d8571a7d
from __future__ import division as _division
from .core import (
BaseSearchClass,
ComputeFstat,
SemiCoherentSearch,
SemiCoherentGlitchSearch,
)
from .make_sfts import (
Writer,
GlitchWriter,
FrequencyModulatedArtifactWriter,
FrequencyAmplitudeModulatedArtifactWriter,
)
from .mcmc_based_searches import (
MCMCSearch,
MCMCGlitchSearch,
MCMCSemiCoherentSearch,
MCMCFollowUpSearch,
MCMCTransientSearch,
)
from .grid_based_searches import (
GridSearch,
GridUniformPriorSearch,
GridGlitchSearch,
FrequencySlidingWindow,
DMoff_NO_SPIN,
SliceGridSearch,
TransientGridSearch,
)
from .core import BaseSearchClass, ComputeFstat, SemiCoherentSearch, SemiCoherentGlitchSearch
from .make_sfts import Writer, GlitchWriter, FrequencyModulatedArtifactWriter, FrequencyAmplitudeModulatedArtifactWriter
from .mcmc_based_searches import MCMCSearch, MCMCGlitchSearch, MCMCSemiCoherentSearch, MCMCFollowUpSearch, MCMCTransientSearch
from .grid_based_searches import GridSearch, GridUniformPriorSearch, GridGlitchSearch, FrequencySlidingWindow, DMoff_NO_SPIN, SliceGridSearch, TransientGridSearch
from .helper_functions import get_version_information
__version__ = get_version_information()

AEI Privacy Policy Imprint