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
  • 72-improve-docs-for_optimal_setup
  • add-higher-spindown-components
  • develop-GA
  • master
  • os-path-join
  • v1.0.1
  • v1.1.0
  • v1.1.2
  • v1.2
  • v1.3
10 results

Target

Select target project
No results found
Select Git revision
  • 72-improve-docs-for_optimal_setup
  • add-higher-spindown-components
  • develop-GA
  • master
  • os-path-join
  • v1.0.1
  • v1.1.0
  • v1.1.2
  • v1.2
  • v1.3
10 results
Show changes

Commits on Source 294

194 additional commits have been omitted to prevent performance issues.
104 files
+ 8461
− 6846
Compare changes
  • Side-by-side
  • Inline

Files

.gitignore

0 → 100644
+2 −0
Original line number Original line Diff line number Diff line 
*pyc

*swp

.gitlab-ci.yml

+39 −6
Original line number Original line Diff line number Diff line 
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:

  script:

    - . /home/greg/lalsuite-install/etc/lalapps-user-env.sh

    - pip install -e $CI_PROJECT_DIR

    - /home/greg/anaconda2/bin/python tests.py TestWriter

      # make sure to test *installed* version of pyFstat

    - /home/greg/anaconda2/bin/python tests.py TestBaseSearchClass

    - (cd .. && pytest $CI_PROJECT_DIR/tests.py --log-file=$CI_PROJECT_DIR/tests.log)

    - /home/greg/anaconda2/bin/python tests.py TestComputeFstat



    - /home/greg/anaconda2/bin/python tests.py TestSemiCoherentGlitchSearch
 
  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

Paper/AllSkyMC/AllSkyMC_repeat.sh

deleted100755 → 0
+0 −11
Original line number Original line Diff line number Diff line 
#!/bin/bash



. /home/gregory.ashton/lalsuite-install/etc/lalapps-user-env.sh

export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"

export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib



for ((n=0;n<10;n++))

do

/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive

done

cp /local/user/gregory.ashton/MCResults_"$1".txt $(pwd)/CollectedOutput

Paper/AllSkyMC/generate_data.py

deleted100644 → 0
+0 −100
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import os

import sys

import time



ID = sys.argv[1]

outdir = sys.argv[2]



label = 'run_{}'.format(ID)

data_label = '{}_data'.format(label)

results_file_name = '{}/MCResults_{}.txt'.format(outdir, ID)



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = -1e-10

F2 = 0

tref = .5*(tstart+tend)





VF0 = VF1 = 100

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)



DeltaAlpha = 0.02

DeltaDelta = 0.02



depths = np.linspace(100, 400, 9)

depths = [118.75, 156.25]



nsteps = 50

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 11, False),

             ((nsteps, 0), 6, False),

             ((nsteps, 0), 3, False),

             ((nsteps, nsteps), 1, False)]





for depth in depths:

    h0 = sqrtSX / float(depth)

    F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0

    F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1

    Alpha_center = np.random.uniform(DeltaAlpha, 2*np.pi-DeltaAlpha)

    Delta_center = np.arccos(2*np.random.uniform(0, 1)-1)-np.pi/2

    Alpha = Alpha_center + np.random.uniform(-0.5, 0.5)*DeltaAlpha

    Delta = Delta_center + np.random.uniform(-0.5, 0.5)*DeltaDelta

    psi = np.random.uniform(-np.pi/4, np.pi/4)

    phi = np.random.uniform(0, 2*np.pi)

    cosi = np.random.uniform(-1, 1)



    data = pyfstat.Writer(

        label=data_label, outdir=outdir, tref=tref,

        tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,

        Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,

        detector='H1,L1')

    data.make_data()



    startTime = time.time()

    theta_prior = {'F0': {'type': 'unif',

                          'lower': F0_center-DeltaF0,

                          'upper': F0_center+DeltaF0},

                   'F1': {'type': 'unif',

                          'lower': F1_center-DeltaF1,

                          'upper': F1_center+DeltaF1},

                   'F2': F2,

                   'Alpha': {'type': 'unif',

                             'lower': Alpha_center-DeltaAlpha,

                             'upper': Alpha_center+DeltaAlpha},

                   'Delta': {'type': 'unif',

                             'lower': Delta_center-DeltaDelta,

                             'upper': Delta_center+DeltaDelta},

                   }



    ntemps = 2

    log10temperature_min = -1

    nwalkers = 100



    mcmc = pyfstat.MCMCFollowUpSearch(

        label=label, outdir=outdir,

        sftfilepath='{}/*{}*sft'.format(outdir, data_label),

        theta_prior=theta_prior,

        tref=tref, minStartTime=tstart, maxStartTime=tend,

        nwalkers=nwalkers, ntemps=ntemps,

        log10temperature_min=log10temperature_min)

    mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,

             gen_tex_table=False)

    d, maxtwoF = mcmc.get_max_twoF()

    dF0 = F0 - d['F0']

    dF1 = F1 - d['F1']

    runTime = time.time() - startTime

    with open(results_file_name, 'a') as f:

        f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {}\n'

                .format(depth, h0, dF0, dF1, maxtwoF, runTime))

    os.system('rm {}/*{}*'.format(outdir, label))

Paper/AllSkyMC/generate_failures.py

deleted100644 → 0
+0 −92
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import os

import time



outdir = 'data'



label = 'run_failures'

data_label = '{}_data'.format(label)

results_file_name = '{}/MCResults_failures.txt'.format(outdir)



# Properties of the GW data

sqrtSX = 2e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = 1e-10

F2 = 0

tref = .5*(tstart+tend)





VF0 = VF1 = 100

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)



DeltaAlpha = 0.02

DeltaDelta = 0.02



depths = [140]



nsteps = 50

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 11, False),

             ((nsteps, 0), 6, False),

             ((nsteps, 0), 3, False),

             ((nsteps, nsteps), 1, False)]





for depth in depths:

    h0 = sqrtSX / float(depth)

    F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0

    F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1

    Alpha_center = np.random.uniform(0, 2*np.pi)

    Delta_center = np.arccos(2*np.random.uniform(0, 1)-1)-np.pi/2

    Alpha = Alpha_center + np.random.uniform(-0.5, 0.5)*DeltaAlpha

    Delta = Delta_center + np.random.uniform(-0.5, 0.5)*DeltaDelta

    psi = np.random.uniform(-np.pi/4, np.pi/4)

    phi = np.random.uniform(0, 2*np.pi)

    cosi = np.random.uniform(-1, 1)



    data = pyfstat.Writer(

        label=data_label, outdir=outdir, tref=tref,

        tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,

        Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,

        detector='H1,L1')

    data.make_data()

    predicted_twoF = data.predict_fstat()



    startTime = time.time()

    theta_prior = {'F0': {'type': 'unif',

                          'lower': F0_center-DeltaF0,

                          'upper': F0_center+DeltaF0},

                   'F1': {'type': 'unif',

                          'lower': F1_center-DeltaF1,

                          'upper': F1_center+DeltaF1},

                   'F2': F2,

                   'Alpha': {'type': 'unif',

                             'lower': Alpha_center-DeltaAlpha,

                             'upper': Alpha_center+DeltaAlpha},

                   'Delta': {'type': 'unif',

                             'lower': Delta_center-DeltaDelta,

                             'upper': Delta_center+DeltaDelta},

                   }



    ntemps = 2

    log10temperature_min = -1

    nwalkers = 100



    mcmc = pyfstat.MCMCFollowUpSearch(

        label=label, outdir=outdir,

        sftfilepath='{}/*{}*sft'.format(outdir, data_label),

        theta_prior=theta_prior,

        tref=tref, minStartTime=tstart, maxStartTime=tend,

        nwalkers=nwalkers, ntemps=ntemps,

        log10temperature_min=log10temperature_min)

    mcmc.run(run_setup=run_setup, create_plots=True, log_table=False,

             gen_tex_table=False)

    d, maxtwoF = mcmc.get_max_twoF()

    print 'MaxtwoF = {}'.format(maxtwoF)

Paper/AllSkyMC/generate_table.py

deleted100644 → 0
+0 −88
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np



outdir = 'data'



label = 'allsky_setup'

data_label = '{}_data'.format(label)



# Properties of the GW data

sqrtSX = 2e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = 1e-10

F2 = 0

tref = .5*(tstart+tend)





VF0 = VF1 = 100

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)



DeltaAlpha = 0.02

DeltaDelta = 0.02



depth = 100



nsteps = 50

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 11, False),

             ((nsteps, 0), 6, False),

             ((nsteps, 0), 3, False),

             ((nsteps, nsteps), 1, False)]



h0 = sqrtSX / float(depth)

r = np.random.uniform(0, 1)

theta = np.random.uniform(0, 2*np.pi)

F0 = F0_center + 3*np.sqrt(r)*np.cos(theta)/(np.pi**2 * Tspan**2)

F1 = F1_center + 45*np.sqrt(r)*np.sin(theta)/(4*np.pi**2 * Tspan**4)



Alpha_center = 0

Delta_center = 0

Alpha = Alpha_center + np.random.uniform(-0.5, 0.5)*DeltaAlpha

Delta = Delta_center + np.random.uniform(-0.5, 0.5)*DeltaDelta



psi = np.random.uniform(-np.pi/4, np.pi/4)

phi = np.random.uniform(0, 2*np.pi)

cosi = np.random.uniform(-1, 1)



data = pyfstat.Writer(

    label=data_label, outdir=outdir, tref=tref,

    tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,

    Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,

    detector='H1,L1')

data.make_data()

predicted_twoF = data.predict_fstat()



theta_prior = {'F0': {'type': 'unif',

                      'lower': F0_center-DeltaF0,

                      'upper': F0_center+DeltaF0},

               'F1': {'type': 'unif',

                      'lower': F1_center-DeltaF1,

                      'upper': F1_center+DeltaF1},

               'F2': F2,

               'Alpha': {'type': 'unif',

                         'lower': Alpha_center-DeltaAlpha,

                         'upper': Alpha_center+DeltaAlpha},

               'Delta': {'type': 'unif',

                         'lower': Delta_center-DeltaDelta,

                         'upper': Delta_center+DeltaDelta},

               }



ntemps = 2

log10temperature_min = -1

nwalkers = 100



mcmc = pyfstat.MCMCFollowUpSearch(

    label=label, outdir=outdir,

    sftfilepath='{}/*{}*sft'.format(outdir, data_label),

    theta_prior=theta_prior,

    tref=tref, minStartTime=tstart, maxStartTime=tend,

    nwalkers=nwalkers, ntemps=ntemps, nsteps=[nsteps, nsteps],

    log10temperature_min=log10temperature_min)

mcmc.run(Nsegs0=20, R=10)

#mcmc.run(run_setup)

Paper/AllSkyMC/plot_data.py

deleted100644 → 0
+0 −87
Original line number Original line Diff line number Diff line 
import matplotlib.pyplot as plt

import pandas as pd

import numpy as np

import os

from tqdm import tqdm

from oct2py import octave

import glob



filenames = glob.glob("CollectedOutput/*.txt")



plt.style.use('paper')



Tspan = 100 * 86400





def Recovery(Tspan, Depth, twoFstar=60, detectors='H1,L1'):

    numDetectors = len(detectors.split(','))

    cmd = ("DetectionProbabilityStackSlide('Nseg', 1, 'Tdata', {},"

           "'misHist', createDeltaHist(0), 'avg2Fth', {}, 'detectors', '{}',"

           "'Depth', {})"

           ).format(numDetectors*Tspan, twoFstar, detectors, Depth)

    return octave.eval(cmd, verbose=False)





def binomialConfidenceInterval(N, K, confidence=0.95):

    cmd = '[fLow, fUpper] = binomialConfidenceInterval({}, {}, {})'.format(

        N, K, confidence)

    [l, u] =  octave.eval(cmd, verbose=False, return_both=True)[0].split('\n')

    return float(l.split('=')[1]), float(u.split('=')[1])



df_list = []

for fn in filenames:

    df = pd.read_csv(

        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF', 'runTime'])

    df['CLUSTER_ID'] = fn.split('_')[1]

    df_list.append(df)

df = pd.concat(df_list)



twoFstar = 70

depths = np.unique(df.depth.values)

recovery_fraction = []

recovery_fraction_CI = []

for d in depths:

    twoFs = df[df.depth == d].twoF.values

    N = len(twoFs)

    K = np.sum(twoFs > twoFstar)

    print d, N, K

    recovery_fraction.append(K/float(N))

    [fLower, fUpper] = binomialConfidenceInterval(N, K)

    recovery_fraction_CI.append([fLower, fUpper])



yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)

fig, ax = plt.subplots()

ax.errorbar(depths, recovery_fraction, yerr=yerr, fmt='sr', marker='s', ms=2,

            capsize=1, capthick=0.5, elinewidth=0.5,

            label='Monte-Carlo result', zorder=10)



fname = 'analytic_data_{}.txt'.format(twoFstar)

if os.path.isfile(fname):

    depths_smooth, recovery_analytic = np.loadtxt(fname)

else:

    depths_smooth = np.linspace(10, 550, 100)

    recovery_analytic = []

    for d in tqdm(depths_smooth):

        recovery_analytic.append(Recovery(Tspan, d, twoFstar))

    np.savetxt(fname, np.array([depths_smooth, recovery_analytic]))

depths_smooth = np.concatenate(([0], depths_smooth))

recovery_analytic = np.concatenate(([1], recovery_analytic))

ax.plot(depths_smooth, recovery_analytic, '-k', label='Theoretical maximum')





ax.set_ylim(0, 1.05)

ax.set_xlabel(r'Sensitivity depth', size=10)

ax.set_ylabel(r'Recovered fraction', size=10)

ax.legend(loc=1, frameon=False)



fig.tight_layout()

fig.savefig('allsky_recovery.png')





total_number_steps = 6 * 50.

fig, ax = plt.subplots()

ax.hist(df.runTime/total_number_steps, bins=50)

ax.set_xlabel('run-time per step [s]')

fig.tight_layout()

fig.savefig('runTimeHist.png')

Paper/AllSkyMC/submitfile

deleted100644 → 0
+0 −12
Original line number Original line Diff line number Diff line 
Executable=AllSkyMC_repeat.sh

Arguments=$(Cluster)_$(Process)

Universe=vanilla

Input=/dev/null

accounting_group = ligo.dev.o2.cw.explore.test

Output=CollectedOutput/out.$(Cluster).$(Process)

Error=CollectedOutput/err.$(Cluster).$(Process)

Log=CollectedOutput/log.$(Cluster).$(Process)

request_cpus = 1

request_memory = 16 GB



Queue 1

Paper/AllSkyMCNoiseOnly/AllSkyMCNoiseOnly_repeat.sh

deleted100755 → 0
+0 −11
Original line number Original line Diff line number Diff line 
#!/bin/bash



. /home/gregory.ashton/lalsuite-install/etc/lalapps-user-env.sh

export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"

export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib



for ((n=0;n<90;n++))

do

/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive

done

cp /local/user/gregory.ashton/NoiseOnlyMCResults_"$1".txt $(pwd)/CollectedOutput

Paper/AllSkyMCNoiseOnly/generate_data.py

deleted100644 → 0
+0 −96
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import os

import sys

import time



ID = sys.argv[1]

outdir = sys.argv[2]



label = 'run_{}'.format(ID)

data_label = '{}_data'.format(label)

results_file_name = '{}/NoiseOnlyMCResults_{}.txt'.format(outdir, ID)



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = -1e-10

F2 = 0

tref = .5*(tstart+tend)





VF0 = VF1 = 100

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)



DeltaAlpha = 0.02

DeltaDelta = 0.02



nsteps = 50

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 11, False),

             ((nsteps, 0), 6, False),

             ((nsteps, 0), 3, False),

             ((nsteps, nsteps), 1, False)]





h0 = 0

F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0

F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1

Alpha_center = np.random.uniform(DeltaAlpha, 2*np.pi-DeltaAlpha)

Delta_center = np.arccos(2*np.random.uniform(0, 1)-1)-np.pi/2

Alpha = Alpha_center + np.random.uniform(-0.5, 0.5)*DeltaAlpha

Delta = Delta_center + np.random.uniform(-0.5, 0.5)*DeltaDelta

psi = np.random.uniform(-np.pi/4, np.pi/4)

phi = np.random.uniform(0, 2*np.pi)

cosi = np.random.uniform(-1, 1)



data = pyfstat.Writer(

    label=data_label, outdir=outdir, tref=tref,

    tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,

    Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,

    detector='H1,L1')

data.make_data()



startTime = time.time()

theta_prior = {'F0': {'type': 'unif',

                      'lower': F0_center-DeltaF0,

                      'upper': F0_center+DeltaF0},

               'F1': {'type': 'unif',

                      'lower': F1_center-DeltaF1,

                      'upper': F1_center+DeltaF1},

               'F2': F2,

               'Alpha': {'type': 'unif',

                         'lower': Alpha_center-DeltaAlpha,

                         'upper': Alpha_center+DeltaAlpha},

               'Delta': {'type': 'unif',

                         'lower': Delta_center-DeltaDelta,

                         'upper': Delta_center+DeltaDelta},

               }



ntemps = 2

log10temperature_min = -1

nwalkers = 100



mcmc = pyfstat.MCMCFollowUpSearch(

    label=label, outdir=outdir,

    sftfilepath='{}/*{}*sft'.format(outdir, data_label),

    theta_prior=theta_prior,

    tref=tref, minStartTime=tstart, maxStartTime=tend,

    nwalkers=nwalkers, ntemps=ntemps,

    log10temperature_min=log10temperature_min)

mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,

         gen_tex_table=False)

d, maxtwoF = mcmc.get_max_twoF()

dF0 = F0 - d['F0']

dF1 = F1 - d['F1']

runTime = time.time() - startTime

with open(results_file_name, 'a') as f:

    f.write('{:1.8e} {:1.8e} {:1.8e} {}\n'

            .format(dF0, dF1, maxtwoF, runTime))

os.system('rm {}/*{}*'.format(outdir, label))

Paper/AllSkyMCNoiseOnly/generate_table.py

deleted100644 → 0
+0 −85
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np



outdir = 'data'



label = 'allsky_setup'

data_label = '{}_data'.format(label)



# Properties of the GW data

sqrtSX = 2e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = 1e-10

F2 = 0

tref = .5*(tstart+tend)





VF0 = VF1 = 100

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)



DeltaAlpha = 0.05

DeltaDelta = 0.05



depth = 100



nsteps = 50

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 11, False),

             ((nsteps, 0), 6, False),

             ((nsteps, 0), 3, False),

             ((nsteps, nsteps), 1, False)]



h0 = sqrtSX / float(depth)

r = np.random.uniform(0, 1)

theta = np.random.uniform(0, 2*np.pi)

F0 = F0_center + 3*np.sqrt(r)*np.cos(theta)/(np.pi**2 * Tspan**2)

F1 = F1_center + 45*np.sqrt(r)*np.sin(theta)/(4*np.pi**2 * Tspan**4)



Alpha = 0

Delta = 0



psi = np.random.uniform(-np.pi/4, np.pi/4)

phi = np.random.uniform(0, 2*np.pi)

cosi = np.random.uniform(-1, 1)



data = pyfstat.Writer(

    label=data_label, outdir=outdir, tref=tref,

    tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,

    Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,

    detector='H1,L1')

data.make_data()

predicted_twoF = data.predict_fstat()



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': {'type': 'unif',

                         'lower': Alpha-DeltaAlpha/2.,

                         'upper': Alpha+DeltaAlpha/2.},

               'Delta': {'type': 'unif',

                         'lower': Delta-DeltaDelta/2.,

                         'upper': Delta+DeltaDelta/2.},

               }



ntemps = 1

log10temperature_min = -1

nwalkers = 100



mcmc = pyfstat.MCMCFollowUpSearch(

    label=label, outdir=outdir,

    sftfilepath='{}/*{}*sft'.format(outdir, data_label),

    theta_prior=theta_prior,

    tref=tref, minStartTime=tstart, maxStartTime=tend,

    nwalkers=nwalkers, ntemps=ntemps,

    log10temperature_min=log10temperature_min)

mcmc.run(run_setup)

Paper/AllSkyMCNoiseOnly/plot_data.py

deleted100644 → 0
+0 −66
Original line number Original line Diff line number Diff line 
import pyfstat

import matplotlib.pyplot as plt

import pandas as pd

import numpy as np

import os

from tqdm import tqdm

from oct2py import octave

import glob

from scipy.stats import rv_continuous, chi2



filenames = glob.glob("CollectedOutput/*.txt")



plt.style.use('paper')



Tspan = 100 * 86400





class maxtwoFinNoise_gen(rv_continuous):

    def _pdf(self, twoF, Ntrials):

        F = twoF/2.0

        alpha = (1 + F)*np.exp(-F)

        a = Ntrials/2.0*F*np.exp(-F)

        b = (1 - alpha)**(Ntrials-1)

        return a*b





df_list = []

for fn in filenames:

    df = pd.read_csv(

        fn, sep=' ', names=['dF0', 'dF1', 'twoF', 'runTime'])

    df['CLUSTER_ID'] = fn.split('_')[1]

    df_list.append(df)

df = pd.concat(df_list)

print 'Number of samples = ', len(df)



fig, ax = plt.subplots()



ax.hist(df.twoF, bins=50, histtype='step', color='k', normed=True, linewidth=1,

        label='Monte-Carlo histogram')



maxtwoFinNoise = maxtwoFinNoise_gen(a=0)

Ntrials_effective, loc, scale = maxtwoFinNoise.fit(df.twoF.values, floc=0, fscale=1)

print 'Ntrials effective = {:1.2e}'.format(Ntrials_effective)

twoFsmooth = np.linspace(0, df.twoF.max(), 1000)

best_fit_pdf = maxtwoFinNoise.pdf(twoFsmooth, Ntrials_effective)

ax.plot(twoFsmooth, best_fit_pdf, '-r',

        label=r'$p(2\mathcal{{F}}_{{\rm max}})$ for {} $N_{{\rm trials}}$'

        .format(pyfstat.texify_float(Ntrials_effective, d=2)))



pval = 1e-6

twoFsmooth_HD = np.linspace(

    twoFsmooth[np.argmax(best_fit_pdf)], df.twoF.max(), 100000)

best_fit_pdf_HD = maxtwoFinNoise.pdf(twoFsmooth_HD, Ntrials_effective)

spacing = twoFsmooth_HD[1]-twoFsmooth_HD[0]

print twoFsmooth_HD[np.argmin(np.abs(best_fit_pdf_HD - pval))], spacing



ax.set_xlabel('$\widetilde{2\mathcal{F}}$')

ax.set_xlim(0, 60)

ax.legend(frameon=False, fontsize=6, loc=2)

fig.tight_layout()

fig.savefig('allsky_noise_twoF_histogram.png')



from latex_macro_generator import write_to_macro

write_to_macro('AllSkyMCNoiseOnlyMaximum', '{:1.1f}'.format(np.max(df.twoF)),

               '../macros.tex')

write_to_macro('AllSkyMCNoiseN', len(df), '../macros.tex')

Paper/AllSkyMCNoiseOnly/submitfile

deleted100644 → 0
+0 −12
Original line number Original line Diff line number Diff line 
Executable= AllSkyMCNoiseOnly_repeat.sh

Arguments= $(Cluster)_$(Process)

Universe=vanilla

Input=/dev/null

accounting_group = ligo.dev.o2.cw.explore.test

Output=CollectedOutput/out.$(Cluster).$(Process)

Error=CollectedOutput/err.$(Cluster).$(Process)

Log=CollectedOutput/log.$(Cluster).$(Process)

request_cpus = 1

request_memory = 16 GB



Queue 100

Paper/DirectedMC/DirectedMC_repeat.sh

deleted100755 → 0
+0 −11
Original line number Original line Diff line number Diff line 
#!/bin/bash



. /home/gregory.ashton/lalsuite-install/etc/lalapps-user-env.sh

export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"

export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib



for ((n=0;n<10;n++))

do

/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive

done

cp /local/user/gregory.ashton/MCResults_"$1".txt $(pwd)/CollectedOutput

Paper/DirectedMC/generate_data.py

deleted100644 → 0
+0 −91
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import os

import sys

import time



ID = sys.argv[1]

outdir = sys.argv[2]



label = 'run_{}'.format(ID)

data_label = '{}_data'.format(label)

results_file_name = '{}/MCResults_{}.txt'.format(outdir, ID)



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = -1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)

tref = .5*(tstart+tend)





VF0 = VF1 = 100

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)



depths = np.linspace(100, 400, 9)

depths = [118.75, 156.25]



nsteps = 25

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 7, False),

             ((nsteps, 0), 2, False),

             ((nsteps, nsteps), 1, False)]



for depth in depths:

    h0 = sqrtSX / float(depth)

    F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0

    F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1



    psi = np.random.uniform(-np.pi/4, np.pi/4)

    phi = np.random.uniform(0, 2*np.pi)

    cosi = np.random.uniform(-1, 1)



    data = pyfstat.Writer(

        label=data_label, outdir=outdir, tref=tref,

        tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,

        Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,

        detector='H1,L1')

    data.make_data()



    startTime = time.time()

    theta_prior = {'F0': {'type': 'unif',

                          'lower': F0_center-DeltaF0,

                          'upper': F0_center+DeltaF0},

                   'F1': {'type': 'unif',

                          'lower': F1_center-DeltaF1,

                          'upper': F1_center+DeltaF1},

                   'F2': F2,

                   'Alpha': Alpha,

                   'Delta': Delta

                   }



    ntemps = 2

    log10temperature_min = -1

    nwalkers = 100



    mcmc = pyfstat.MCMCFollowUpSearch(

        label=label, outdir=outdir,

        sftfilepath='{}/*{}*sft'.format(outdir, data_label),

        theta_prior=theta_prior,

        tref=tref, minStartTime=tstart, maxStartTime=tend,

        nwalkers=nwalkers, ntemps=ntemps,

        log10temperature_min=log10temperature_min)

    mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,

             gen_tex_table=False)

    mcmc.print_summary()

    d, maxtwoF = mcmc.get_max_twoF()

    dF0 = F0 - d['F0']

    dF1 = F1 - d['F1']

    runTime = time.time() - startTime

    with open(results_file_name, 'a') as f:

        f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.8e} {}\n'

                .format(depth, h0, dF0, dF1, maxtwoF, runTime))

    os.system('rm {}/*{}*'.format(outdir, label))

Paper/DirectedMC/generate_table.py

deleted100644 → 0
+0 −77
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np



outdir = 'data'



label = 'directed_setup'

data_label = '{}_data'.format(label)



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = 1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)

tref = .5*(tstart+tend)





VF0 = VF1 = 100

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)



depth = 100



nsteps = 25

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 7, False),

             ((nsteps, 0), 2, False),

             ((nsteps, nsteps), 1, False)]



h0 = sqrtSX / float(depth)

r = np.random.uniform(0, 1)

theta = np.random.uniform(0, 2*np.pi)

F0 = F0_center + 3*np.sqrt(r)*np.cos(theta)/(np.pi**2 * Tspan**2)

F1 = F1_center + 45*np.sqrt(r)*np.sin(theta)/(4*np.pi**2 * Tspan**4)



psi = np.random.uniform(-np.pi/4, np.pi/4)

phi = np.random.uniform(0, 2*np.pi)

cosi = np.random.uniform(-1, 1)



data = pyfstat.Writer(

    label=data_label, outdir=outdir, tref=tref,

    tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,

    Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,

    detector='H1,L1')

data.make_data()

predicted_twoF = data.predict_fstat()



theta_prior = {'F0': {'type': 'unif',

                      'lower': F0-DeltaF0,

                      'upper': F0+DeltaF0},

               'F1': {'type': 'unif',

                      'lower': F1-DeltaF1,

                      'upper': F1+DeltaF1},

               'F2': F2,

               'Alpha': Alpha,

               'Delta': Delta

               }



ntemps = 1

log10temperature_min = -1

nwalkers = 100



mcmc = pyfstat.MCMCFollowUpSearch(

    label=label, outdir=outdir,

    sftfilepath='{}/*{}*sft'.format(outdir, data_label),

    theta_prior=theta_prior,

    tref=tref, minStartTime=tstart, maxStartTime=tend,

    nwalkers=nwalkers, ntemps=ntemps, nsteps=[nsteps, nsteps],

    log10temperature_min=log10temperature_min)

#mcmc.run(Nsegs0=20, R=10)

mcmc.run(run_setup)

Paper/DirectedMC/plot_data.py

deleted100644 → 0
+0 −87
Original line number Original line Diff line number Diff line 
import matplotlib.pyplot as plt

import pandas as pd

import numpy as np

import os

from tqdm import tqdm

from oct2py import octave

import glob



filenames = glob.glob("CollectedOutput/*.txt")



plt.style.use('paper')



Tspan = 100 * 86400





def Recovery(Tspan, Depth, twoFstar=60, detectors='H1,L1'):

    numDetectors = len(detectors.split(','))

    cmd = ("DetectionProbabilityStackSlide('Nseg', 1, 'Tdata', {},"

           "'misHist', createDeltaHist(0), 'avg2Fth', {}, 'detectors', '{}',"

           "'Depth', {})"

           ).format(numDetectors*Tspan, twoFstar, detectors, Depth)

    return octave.eval(cmd, verbose=False)





def binomialConfidenceInterval(N, K, confidence=0.95):

    cmd = '[fLow, fUpper] = binomialConfidenceInterval({}, {}, {})'.format(

        N, K, confidence)

    [l, u] = octave.eval(cmd, verbose=False, return_both=True)[0].split('\n')

    return float(l.split('=')[1]), float(u.split('=')[1])



df_list = []

for fn in filenames:

    df = pd.read_csv(

        fn, sep=' ', names=['depth', 'h0', 'dF0', 'dF1', 'twoF', 'runTime'])

    df['CLUSTER_ID'] = fn.split('_')[1]

    df_list.append(df)

df = pd.concat(df_list)



twoFstar = 60

depths = np.unique(df.depth.values)

recovery_fraction = []

recovery_fraction_CI = []

for d in depths:

    twoFs = df[df.depth == d].twoF.values

    N = len(twoFs)

    K = np.sum(twoFs > twoFstar)

    print d, N, K

    recovery_fraction.append(K/float(N))

    [fLower, fUpper] = binomialConfidenceInterval(N, K)

    recovery_fraction_CI.append([fLower, fUpper])



yerr = np.abs(recovery_fraction - np.array(recovery_fraction_CI).T)

fig, ax = plt.subplots()

ax.errorbar(depths, recovery_fraction, yerr=yerr, fmt='sr', marker='s', ms=2,

            capsize=1, capthick=0.5, elinewidth=0.5,

            label='Monte-Carlo result', zorder=10)



fname = 'analytic_data_{}.txt'.format(twoFstar)

if os.path.isfile(fname):

    depths_smooth, recovery_analytic = np.loadtxt(fname)

else:

    depths_smooth = np.linspace(10, 550, 100)

    recovery_analytic = []

    for d in tqdm(depths_smooth):

        recovery_analytic.append(Recovery(Tspan, d, twoFstar))

    np.savetxt(fname, np.array([depths_smooth, recovery_analytic]))

depths_smooth = np.concatenate(([0], depths_smooth))

recovery_analytic = np.concatenate(([1], recovery_analytic))

ax.plot(depths_smooth, recovery_analytic, '-k', label='Theoretical maximum')





ax.set_ylim(0, 1.05)

ax.set_xlabel(r'Sensitivity depth', size=10)

ax.set_ylabel(r'Recovered fraction', size=10)

ax.legend(loc=1, frameon=False)



fig.tight_layout()

fig.savefig('directed_recovery.png')





total_number_steps = 5*25.

fig, ax = plt.subplots()

ax.hist(df.runTime/total_number_steps, bins=50)

ax.set_xlabel('run-time per step [s]')

fig.tight_layout()

fig.savefig('runTimeHist.png')

Paper/DirectedMC/plot_recovery.py

deleted100644 → 0
+0 −30
Original line number Original line Diff line number Diff line 
import matplotlib.pyplot as plt

import numpy as np

import scipy.stats





def Recovery(Tspan, Depth, twoFstar=60):

    rho2 = 4*Tspan/25./Depth**2

    twoF_Hs = scipy.stats.distributions.ncx2(df=4, nc=rho2)

    return 1 - twoF_Hs.cdf(twoFstar)



N = 500

Tspan = np.linspace(0.1, 365*86400, N)

Depth = np.linspace(10, 300, N)



X, Y = np.meshgrid(Tspan, Depth)

X = X / 86400



Z = [[Recovery(t, d) for t in Tspan] for d in Depth]



fig, ax = plt.subplots()

pax = ax.pcolormesh(X, Y, Z, cmap=plt.cm.viridis)

CS = ax.contour(X, Y, Z, [0.95])

plt.clabel(CS, inline=1, fontsize=12, fmt='%s', manual=[(200, 180)])

plt.colorbar(pax, label='Recovery fraction')

ax.set_xlabel(r'$T_{\rm span}$ [days]', size=16)

ax.set_ylabel(r'Depth=$\frac{\sqrt{S_{\rm n}}}{h_0}$', size=14)

ax.set_xlim(min(Tspan)/86400., max(Tspan)/86400.)

ax.set_ylim(min(Depth), max(Depth))



fig.savefig('recovery.png')

Paper/DirectedMC/submitfile

deleted100644 → 0
+0 −12
Original line number Original line Diff line number Diff line 
Executable=DirectedMC_repeat.sh

Arguments=$(Cluster)_$(Process)

Universe=vanilla

Input=/dev/null

accounting_group = ligo.dev.o2.cw.explore.test

Output=CollectedOutput/out.$(Cluster).$(Process)

Error=CollectedOutput/err.$(Cluster).$(Process)

Log=CollectedOutput/log.$(Cluster).$(Process)

request_cpus = 1

request_memory = 16 GB



Queue 1

Paper/DirectedMCNoiseOnly/DirectedMCNoiseOnly_repeat.sh

deleted100755 → 0
+0 −11
Original line number Original line Diff line number Diff line 
#!/bin/bash



. /home/gregory.ashton/lalsuite-install/etc/lalapps-user-env.sh

export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"

export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib



for ((n=0;n<100;n++))

do

/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-template-counting --no-interactive

done

cp /local/user/gregory.ashton/NoiseOnlyMCResults_"$1".txt $(pwd)/CollectedOutput

Paper/DirectedMCNoiseOnly/generate_data.py

deleted100644 → 0
+0 −87
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import os

import sys

import time





ID = sys.argv[1]

outdir = sys.argv[2]



label = 'run_{}'.format(ID)

data_label = '{}_data'.format(label)

results_file_name = '{}/NoiseOnlyMCResults_{}.txt'.format(outdir, ID)



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = -1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)

tref = .5*(tstart+tend)



VF0 = VF1 = 100

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*Tspan)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*Tspan**2)



nsteps = 25

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 7, False),

             ((nsteps, 0), 2, False),

             ((nsteps, nsteps), 1, False)]



h0 = 0

F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0

F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1



psi = np.random.uniform(-np.pi/4, np.pi/4)

phi = np.random.uniform(0, 2*np.pi)

cosi = np.random.uniform(-1, 1)



data = pyfstat.Writer(

    label=data_label, outdir=outdir, tref=tref,

    tstart=tstart, F0=F0, F1=F1, F2=F2, duration=Tspan, Alpha=Alpha,

    Delta=Delta, h0=h0, sqrtSX=sqrtSX, psi=psi, phi=phi, cosi=cosi,

    detector='H1,L1')

data.make_data()



startTime = time.time()

theta_prior = {'F0': {'type': 'unif',

                      'lower': F0_center-DeltaF0,

                      'upper': F0_center+DeltaF0},

               'F1': {'type': 'unif',

                      'lower': F1_center-DeltaF1,

                      'upper': F1_center+DeltaF1},

               'F2': F2,

               'Alpha': Alpha,

               'Delta': Delta

               }



ntemps = 2

log10temperature_min = -1

nwalkers = 100



mcmc = pyfstat.MCMCFollowUpSearch(

    label=label, outdir=outdir,

    sftfilepath='{}/*{}*sft'.format(outdir, data_label),

    theta_prior=theta_prior,

    tref=tref, minStartTime=tstart, maxStartTime=tend,

    nwalkers=nwalkers, ntemps=ntemps,

    log10temperature_min=log10temperature_min)

mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,

         gen_tex_table=False)

mcmc.print_summary()

d, maxtwoF = mcmc.get_max_twoF()

dF0 = F0 - d['F0']

dF1 = F1 - d['F1']

runTime = time.time() - startTime

with open(results_file_name, 'a') as f:

    f.write('{:1.8e} {:1.8e} {:1.8e} {}\n'

            .format(dF0, dF1, maxtwoF, runTime))

os.system('rm {}/*{}*'.format(outdir, label))

Paper/DirectedMCNoiseOnly/plot_data.py

deleted100644 → 0
+0 −65
Original line number Original line Diff line number Diff line 
import pyfstat

import matplotlib.pyplot as plt

import pandas as pd

import numpy as np

import os

from tqdm import tqdm

from oct2py import octave

import glob

from scipy.stats import rv_continuous, chi2



filenames = glob.glob("CollectedOutput/*.txt")



plt.style.use('paper')



Tspan = 100 * 86400





class maxtwoFinNoise_gen(rv_continuous):

    def _pdf(self, twoF, Ntrials):

        F = twoF/2.0

        alpha = (1 + F)*np.exp(-F)

        a = Ntrials/2.0*F*np.exp(-F)

        b = (1 - alpha)**(Ntrials-1)

        return a*b





df_list = []

for fn in filenames:

    df = pd.read_csv(

        fn, sep=' ', names=['dF0', 'dF1', 'twoF', 'runTime'])

    df['CLUSTER_ID'] = fn.split('_')[1]

    df_list.append(df)

df = pd.concat(df_list)

print 'Number of samples = ', len(df)



fig, ax = plt.subplots()

ax.hist(df.twoF, bins=50, histtype='step', color='k', normed=True, linewidth=1,

        label='Monte-Carlo histogram')



maxtwoFinNoise = maxtwoFinNoise_gen(a=0)

Ntrials_effective, loc, scale = maxtwoFinNoise.fit(df.twoF.values, floc=0, fscale=1)

print 'Ntrials effective = {:1.2e}'.format(Ntrials_effective)

twoFsmooth = np.linspace(0, df.twoF.max(), 1000)

best_fit_pdf = maxtwoFinNoise.pdf(twoFsmooth, Ntrials_effective)

ax.plot(twoFsmooth, best_fit_pdf, '-r',

        label=r'$p(2\mathcal{{F}}_{{\rm max}})$ for {} $N_{{\rm trials}}$'

        .format(pyfstat.texify_float(Ntrials_effective, d=2)))



pval = 1e-6

twoFsmooth_HD = np.linspace(

    twoFsmooth[np.argmax(best_fit_pdf)], df.twoF.max(), 100000)

best_fit_pdf_HD = maxtwoFinNoise.pdf(twoFsmooth_HD, Ntrials_effective)

spacing = twoFsmooth_HD[1]-twoFsmooth_HD[0]

print twoFsmooth_HD[np.argmin(np.abs(best_fit_pdf_HD - pval))], spacing



ax.set_xlabel('$\widetilde{2\mathcal{F}}$')

ax.set_xlim(0, 60)

ax.legend(frameon=False, fontsize=6)

fig.tight_layout()

fig.savefig('directed_noise_twoF_histogram.png')



from latex_macro_generator import write_to_macro

write_to_macro('DirectedMCNoiseOnlyMaximum', '{:1.1f}'.format(np.max(df.twoF)),

               '../macros.tex')

write_to_macro('DirectedMCNoiseN', len(df), '../macros.tex')

Paper/DirectedMCNoiseOnly/submitfile

deleted100644 → 0
+0 −12
Original line number Original line Diff line number Diff line 
Executable= DirectedMCNoiseOnly_repeat.sh

Arguments= $(Cluster)_$(Process)

Universe=vanilla

Input=/dev/null

accounting_group = ligo.dev.o2.cw.explore.test

Output=CollectedOutput/out.$(Cluster).$(Process)

Error=CollectedOutput/err.$(Cluster).$(Process)

Log=CollectedOutput/log.$(Cluster).$(Process)

request_cpus = 1

request_memory = 16 GB



Queue 1

Paper/Examples/Makefile

deleted100644 → 0
+0 −10
Original line number Original line Diff line number Diff line 
frequency_grid_files = grided_frequency_search_1D.png

fully_coherent_files = fully_coherent_search_using_MCMC_walkers.png

follow_up_files = follow_up_run_setup.tex follow_up_walkers.png

transient_files = transient_search_initial_stage_twoFcumulative.png transient_search_corner.png

glitch_files = single_glitch_F0F1_grid_2D.png single_glitch_glitchSearch_corner.png



all_files = $(frequency_grid_files) $(fully_coherent_files) $(follow_up_files) $(transient_files) $(glitch_files)



copyfiles:

	cd data; cp $(all_files) ../../

Paper/Examples/follow_up.py

deleted100644 → 0
+0 −76
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import matplotlib.pyplot as plt

import matplotlib



F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = 1.0

Delta = 0.5



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

duration = 100*86400

tend = tstart+duration

tref = .5*(tstart+tend)



depth = 50

data_label = 'follow_up'



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)

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)



# Search

VF0 = VF1 = 500

DeltaF0 = VF0 * np.sqrt(3)/(np.pi*duration)

DeltaF1 = VF1 * np.sqrt(45/4.)/(np.pi*duration**2)

DeltaAlpha = 1e-1

DeltaDelta = 1e-1

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': {'type': 'unif', 'lower': Alpha-DeltaAlpha,

                         'upper': Alpha+DeltaAlpha},

               'Delta': {'type': 'unif', 'lower': Delta-DeltaDelta,

                         'upper': Delta+DeltaDelta},

               }



ntemps = 3

log10temperature_min = -0.5

nwalkers = 100

scatter_val = 1e-10

nsteps = [200, 200]



mcmc = pyfstat.MCMCFollowUpSearch(

    label='follow_up', outdir='data',

    sftfilepath='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,

    minStartTime=tstart, maxStartTime=tend, nwalkers=nwalkers, nsteps=nsteps,

    ntemps=ntemps, log10temperature_min=log10temperature_min,

    scatter_val=scatter_val)



fig, axes = plt.subplots(nrows=2, ncols=2)

fig, axes = mcmc.run(

    R=10, Nsegs0=100, subtractions=[F0, F1, Alpha, Delta], labelpad=0.01,

    fig=fig, axes=axes, plot_det_stat=False, return_fig=True)

axes[3].set_xlabel(r'$\textrm{Number of steps}$', labelpad=0.1)

for ax in axes:

    ax.set_xlim(0, axes[0].get_xlim()[-1])

    ax.xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(5))

fig.tight_layout()

fig.savefig('{}/{}_walkers.png'.format(mcmc.outdir, mcmc.label), dpi=400)



mcmc.plot_corner(add_prior=True)

mcmc.print_summary()

Paper/Examples/fully_coherent_search_using_MCMC.py

deleted100644 → 0
+0 −76
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

duration = 100*86400

tend = tstart + duration



# Properties of the signal

F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = 5e-3

Delta = 6e-2

tref = .5*(tstart+tend)



depth = 10

h0 = sqrtSX / depth

data_label = 'fully_coherent_search_using_MCMC'



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)

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)



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)



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

               }



ntemps = 1

log10temperature_min = -1

nwalkers = 1000

nsteps = [50, 50]



mcmc = pyfstat.MCMCSearch(

    label='fully_coherent_search_using_MCMC', outdir='data',

    sftfilepath='data/*'+data_label+'-*sft', theta_prior=theta_prior, tref=tref,

    minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,

    ntemps=ntemps, log10temperature_min=log10temperature_min)

mcmc.run(context='paper', subtractions=[30, -1e-10])

mcmc.plot_corner(add_prior=True)

mcmc.print_summary()



from latex_macro_generator import write_to_macro

write_to_macro('BasicExampleF0', '{:1.0f}'.format(F0), '../macros.tex')

write_to_macro('BasicExampleF1', F1, '../macros.tex')

write_to_macro('BasicExampleh0', h0, '../macros.tex')

write_to_macro('BasicExampleSqrtSn', sqrtSX, '../macros.tex')

write_to_macro('BasicExampleDepth', depth, '../macros.tex')

write_to_macro('BasicExampleDeltaF0', DeltaF0, '../macros.tex')

write_to_macro('BasicExampleDeltaF1', DeltaF1, '../macros.tex')

write_to_macro('BasicExampleVF0', VF0, '../macros.tex')

write_to_macro('BasicExampleVF1', VF1, '../macros.tex')

write_to_macro('BasicExampleV', VF0*VF1, '../macros.tex')

write_to_macro('BasicExamplenburn', nsteps[0], '../macros.tex')

write_to_macro('BasicExamplenprod', nsteps[1], '../macros.tex')

Paper/Examples/fully_coherent_search_using_MCMC_convergence.py

deleted100644 → 0
+0 −64
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

duration = 100*86400

tend = tstart + duration



# Properties of the signal

F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = 5e-3

Delta = 6e-2

tref = .5*(tstart+tend)



depth = 10

h0 = sqrtSX / depth

data_label = 'fully_coherent_search_using_MCMC_convergence'



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)

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)



DeltaF0 = 5e-7

DeltaF1 = 1e-12

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)



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

               }



ntemps = 1

log10temperature_min = -1

nwalkers = 100

nsteps = [900, 100]



mcmc = pyfstat.MCMCSearch(

    label='fully_coherent_search_using_MCMC_convergence', outdir='data',

    sftfilepath='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,

    minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,

    ntemps=ntemps, log10temperature_min=log10temperature_min)

mcmc.setup_convergence_testing(

    convergence_threshold_number=5, convergence_plot_upper_lim=10,

    convergence_burnin_fraction=0.1)

mcmc.run(context='paper', subtractions=[30, -1e-10])

mcmc.plot_corner(add_prior=True)

mcmc.print_summary()

Paper/Examples/grided_frequency_search.py

deleted100644 → 0
+0 −65
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import matplotlib.pyplot as plt



plt.style.use('paper')



F0 = 30.0

F1 = 0

F2 = 0

Alpha = 1.0

Delta = 1.5



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

duration = 100*86400

tend = tstart+duration

tref = .5*(tstart+tend)



depth = 70

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)

data.make_data()



m = 1

dF0 = np.sqrt(12*m)/(np.pi*duration)

DeltaF0 = 30*dF0

F0s = [F0-DeltaF0/2., F0+DeltaF0/2., 1e-2*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)

search.run()



fig, ax = plt.subplots()

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=0.8)

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.2, dashes=(0.2, 0.2))

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 = [r'$f_0 {-} 10\Delta f(\tilde{\mu}=1)$', '$f_0$',

               r'$f_0 {+} 10\Delta f(\tilde{\mu}=1)$']

ax.set_xticklabels(xticklabels)

plt.tight_layout()

fig.savefig('{}/{}_1D.png'.format(search.outdir, search.label), dpi=300)

Paper/Examples/matplotlibrc

deleted100644 → 0
+0 −36
Original line number Original line Diff line number Diff line 
figure.figsize: 3.4, 2.5

figure.facecolor: w



axes.labelsize: 6

axes.titlesize: 6

xtick.labelsize: 4

ytick.labelsize: 4

legend.fontsize: 6

font.size: 6



grid.linewidth: 0.8

lines.linewidth: 1.0

patch.linewidth: 0.24

lines.markersize: 3.6

lines.markeredgewidth: 0



xtick.major.size: 2.0

ytick.major.size: 2.0

xtick.minor.size: 1.5

xtick.minor.visible: True

ytick.minor.size: 1.5

ytick.minor.visible: True



xtick.major.pad: 2.5

ytick.major.pad: 2.5



pdf.compression: 9

savefig.format: png

savefig.dpi: 300



font.family: serif

font.serif: Computer Modern

text.usetex: True

axes.formatter.use_mathtext: True

axes.formatter.useoffset: False

axes.formatter.limits : -3, 4

Paper/Examples/single_glitch.py

deleted100644 → 0
+0 −115
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import matplotlib.pyplot as plt

import matplotlib



sqrtSX = 1e-22

tstart = 1000000000

duration = 100*86400

tend = tstart+duration



# Define parameters of the Crab pulsar as an example

tref = .5*(tstart + tend)

F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)



# Signal strength

depth = 10

h0 = sqrtSX / depth



VF0 = VF1 = 200

dF0 = np.sqrt(3)/(np.pi*duration)

dF1 = np.sqrt(45/4.)/(np.pi*duration**2)

DeltaF0 = VF0 * dF0

DeltaF1 = VF1 * dF1



# Next, taking the same signal parameters, we include a glitch half way through

R = 0.5

dtglitch = R*duration

delta_F0 = 0.25*DeltaF0

delta_F1 = -0.1*DeltaF1



glitch_data = pyfstat.Writer(

    label='single_glitch', outdir='data', 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()

FCtwoFPM = glitch_data.predict_fstat()



F0s = [F0-DeltaF0/2., F0+DeltaF0/2., 1*dF0]

F1s = [F1-DeltaF1/2., F1+DeltaF1/2., 1*dF1]

F2s = [F2]

Alphas = [Alpha]

Deltas = [Delta]

search = pyfstat.GridSearch(

    'single_glitch_F0F1_grid', 'data', 'data/*single_glitch*sft', F0s, F1s,

    F2s, Alphas, Deltas, tref, tstart, tend)

search.run()

ax = search.plot_2D('F0', 'F1', save=False)

ax.xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(6))

plt.tight_layout()

plt.savefig('data/single_glitch_F0F1_grid_2D.png', dpi=500)

FCtwoF = search.get_max_twoF()['twoF']





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.,

                      'upper': F0+DeltaF0/2},

               'F1': {'type': 'unif', 'lower': F1-DeltaF1/2.,

                      'upper': F1+DeltaF1/2},

               'F2': F2,

               'Alpha': Alpha,

               'Delta': Delta,

               'tglitch': {'type': 'unif', 'lower': tstart+0.1*duration,

                           'upper': tend-0.1*duration},

               'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': 1e-7*F0},

               'delta_F1': {'type': 'norm', 'loc': 0, 'scale': abs(1e-3*F1)},

               }

ntemps = 3

log10temperature_min = -0.1

nwalkers = 100

nsteps = [1000, 1000]

glitch_mcmc = pyfstat.MCMCGlitchSearch(

    'single_glitch_glitchSearch', 'data',

    sftfilepath='data/*_single_glitch*.sft', theta_prior=theta_prior,

    tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,

    nwalkers=nwalkers, ntemps=ntemps,

    log10temperature_min=log10temperature_min)

glitch_mcmc.write_prior_table()

glitch_mcmc.run()

glitch_mcmc.plot_corner(figsize=(6, 6))

glitch_mcmc.print_summary()

GlitchtwoF = glitch_mcmc.get_max_twoF()[1]



from latex_macro_generator import write_to_macro

write_to_macro('SingleGlitchTspan', '{:1.1f}'.format(duration/86400),

               '../macros.tex')

write_to_macro('SingleGlitchR', '{:1.1f}'.format(R), '../macros.tex')

write_to_macro('SingleGlitchDepth',

               '{}'.format(pyfstat.texify_float(sqrtSX/h0)), '../macros.tex')

write_to_macro('SingleGlitchF0', '{}'.format(F0), '../macros.tex')

write_to_macro('SingleGlitchF1', '{}'.format(F1), '../macros.tex')

write_to_macro('SingleGlitchdeltaF0',

               '{}'.format(pyfstat.texify_float(delta_F0)), '../macros.tex')

write_to_macro('SingleGlitchdeltaF1',

               '{}'.format(pyfstat.texify_float(delta_F1)), '../macros.tex')

write_to_macro('SingleGlitchFCtwoFPM', '{:1.1f}'.format(FCtwoFPM),

               '../macros.tex')

write_to_macro('SingleGlitchFCtwoF', '{:1.1f}'.format(FCtwoF),

               '../macros.tex')

write_to_macro('SingleGlitch_FCMismatch',

               '{:1.1f}'.format((FCtwoFPM-FCtwoF)/FCtwoFPM), '../macros.tex')

write_to_macro('SingleGlitchGlitchtwoF', '{:1.1f}'.format(GlitchtwoF),

               '../macros.tex')

Paper/Examples/transient_search_using_MCMC.py

deleted100644 → 0
+0 −95
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import matplotlib.pyplot as plt



plt.style.use('thesis')



F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = 5e-3

Delta = 6e-2



tstart = 1000000000

duration = 100*86400

data_tstart = tstart - duration

data_tend = data_tstart + 3*duration

tref = .5*(data_tstart+data_tend)



h0 = 4e-24

sqrtSX = 1e-22



transient = pyfstat.Writer(

    label='transient', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,

    F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,

    minStartTime=data_tstart, maxStartTime=data_tend)

transient.make_data()

print transient.predict_fstat()



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)



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

               }



ntemps = 3

log10temperature_min = -1

nwalkers = 100

nsteps = [100, 100]



mcmc = pyfstat.MCMCSearch(

    label='transient_search_initial_stage', outdir='data',

    sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref,

    minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps,

    nwalkers=nwalkers, ntemps=ntemps,

    log10temperature_min=log10temperature_min)

mcmc.run()

fig, ax = plt.subplots()

mcmc.write_par()

mcmc.generate_loudest()

mcmc.plot_cumulative_max(ax=ax)

ax.set_xlabel('Days from $t_\mathrm{start}$')

ax.legend_.remove()

fig.savefig('data/transient_search_initial_stage_twoFcumulative')

mcmc.print_summary()



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': data_tstart,

                                    'upper': data_tend-0.2*duration},

               'transient_duration': {'type': 'halfnorm',

                                      'loc': 0.01*duration,

                                      'scale': 0.5*duration}

               }



nwalkers = 500

nsteps = [200, 200]



mcmc = pyfstat.MCMCTransientSearch(

    label='transient_search', outdir='data',

    sftfilepath='data/*transient*sft', theta_prior=theta_prior, tref=tref,

    minStartTime=data_tstart, maxStartTime=data_tend, nsteps=nsteps,

    nwalkers=nwalkers, ntemps=ntemps,

    log10temperature_min=log10temperature_min)

mcmc.run()

mcmc.plot_corner()

mcmc.print_summary()

Paper/GlitchMCNoiseOnly/SingleGlitchMCNoiseOnly.sh

deleted100755 → 0
+0 −11
Original line number Original line Diff line number Diff line 
#!/bin/bash



. /home/gregory.ashton/lalsuite-install/etc/lalapps-user-env.sh

export PATH="/home/gregory.ashton/anaconda2/bin:$PATH"

export MPLCONFIGDIR=/home/gregory.ashton/.config/matplotlib



for ((n=0;n<5;n++))

do

/home/gregory.ashton/anaconda2/bin/python generate_data.py "$1" /local/user/gregory.ashton --no-interactive

done

cp /local/user/gregory.ashton/NoiseOnlyMCResults_"$1".txt $(pwd)/CollectedOutput

Paper/GlitchMCNoiseOnly/generate_data.py

deleted100644 → 0
+0 −100
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import os

import sys

import time



nglitch = 2



ID = sys.argv[1]

outdir = sys.argv[2]



label = 'run_{}'.format(ID)

data_label = '{}_data'.format(label)

results_file_name = '{}/NoiseOnlyMCResults_{}.txt'.format(outdir, ID)



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

Tspan = 100*86400

tend = tstart + Tspan



# Fixed properties of the signal

F0_center = 30

F1_center = -1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)

tref = .5*(tstart+tend)



VF0 = VF1 = 200

dF0 = np.sqrt(3)/(np.pi*Tspan)

dF1 = np.sqrt(45/4.)/(np.pi*Tspan**2)

DeltaF0 = VF0 * dF0

DeltaF1 = VF1 * dF1



nsteps = 25

run_setup = [((nsteps, 0), 20, False),

             ((nsteps, 0), 7, False),

             ((nsteps, 0), 2, False),

             ((nsteps, nsteps), 1, False)]



h0 = 0

F0 = F0_center + np.random.uniform(-0.5, 0.5)*DeltaF0

F1 = F1_center + np.random.uniform(-0.5, 0.5)*DeltaF1



psi = np.random.uniform(-np.pi/4, np.pi/4)

phi = np.random.uniform(0, 2*np.pi)

cosi = np.random.uniform(-1, 1)



# Next, taking the same signal parameters, we include a glitch half way through

dtglitch = Tspan/2.0

delta_F0 = 0.25*DeltaF0

delta_F1 = -0.1*DeltaF1



glitch_data = pyfstat.Writer(

    label=data_label, outdir=outdir, tref=tref, tstart=tstart, F0=F0,

    F1=F1, F2=F2, duration=Tspan, Alpha=Alpha, Delta=Delta, h0=h0,

    sqrtSX=sqrtSX, dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)

glitch_data.make_data()





startTime = time.time()

theta_prior = {'F0': {'type': 'unif', 'lower': F0-DeltaF0/2., # PROBLEM

                      'upper': F0+DeltaF0/2},

               'F1': {'type': 'unif', 'lower': F1-DeltaF1/2.,

                      'upper': F1+DeltaF1/2},

               'F2': F2,

               'Alpha': Alpha,

               'Delta': Delta,

               'tglitch': {'type': 'unif', 'lower': tstart+0.1*Tspan,

                           'upper': tend-0.1*Tspan},

               'delta_F0': {'type': 'halfnorm', 'loc': 0, 'scale': DeltaF0},

               'delta_F1': {'type': 'norm', 'loc': 0, 'scale': DeltaF1},

               }

ntemps = 2

log10temperature_min = -0.1

nwalkers = 100

nsteps = [500, 500]

glitch_mcmc = pyfstat.MCMCGlitchSearch(

    label=label, outdir=outdir,

    sftfilepath='{}/*{}*sft'.format(outdir, data_label),

    theta_prior=theta_prior,

    tref=tref, minStartTime=tstart, maxStartTime=tend, nsteps=nsteps,

    nwalkers=nwalkers, ntemps=ntemps, nglitch=nglitch,

    log10temperature_min=log10temperature_min)

glitch_mcmc.run(run_setup=run_setup, create_plots=False, log_table=False,

                gen_tex_table=False)

glitch_mcmc.print_summary()

d, maxtwoF = glitch_mcmc.get_max_twoF()

dF0 = F0 - d['F0']

dF1 = F1 - d['F1']

#tglitch = d['tglitch']

#R = (tglitch - tstart) / Tspan

#delta_F0 = d['delta_F0']

#delta_F1 = d['delta_F1']

runTime = time.time() - startTime

with open(results_file_name, 'a') as f:

    f.write('{} {:1.8e} {:1.8e} {:1.8e} {:1.1f}\n'

            .format(nglitch, dF0, dF1, maxtwoF, runTime))

os.system('rm {}/*{}*'.format(outdir, label))

Paper/GlitchMCNoiseOnly/plot_data.py

deleted100644 → 0
+0 −91
Original line number Original line Diff line number Diff line 
import pyfstat

import matplotlib.pyplot as plt

import pandas as pd

import numpy as np

import os

from tqdm import tqdm

from oct2py import octave

import glob

from scipy.stats import rv_continuous, chi2

from scipy.special import gammaincc

from latex_macro_generator import write_to_macro





def CF_twoFmax_integrand(theta, twoFmax, Nt):

    Fmax = twoFmax/2.0

    return np.exp(1j*theta*twoFmax)*Nt/2.0*Fmax*np.exp(-Fmax)*(1-(1+Fmax)*np.exp(-Fmax))**(Nt-1.)





def pdf_twoFhat(twoFhat, Ng, Nts, twoFtildemax=100, dtwoF=0.05):

    if np.ndim(Nts) == 0:

        Nts = np.zeros(Ng+1) + Nts



    twoFtildemax_int = np.arange(0, twoFtildemax, dtwoF)

    theta_int = np.arange(-1./dtwoF, 1./dtwoF, 1./twoFtildemax)



    CF_twoFtildemax_theta = np.array(

        [[np.trapz(CF_twoFmax_integrand(t, twoFtildemax_int, Nt), twoFtildemax_int)

          for t in theta_int]

         for Nt in Nts])



    CF_twoFhat_theta = np.prod(CF_twoFtildemax_theta, axis=0)

    print CF_twoFhat_theta.shape, theta_int.shape

    pdf = (1/(2*np.pi)) * np.array(

        [np.trapz(np.exp(-1j*theta_int*twoFhat_val)*CF_twoFhat_theta,

                  theta_int) for twoFhat_val in twoFhat])

    print np.trapz(pdf.real, x=twoFhat)

    return pdf



filenames = glob.glob("CollectedOutput/*.txt")



plt.style.use('paper')



Tspan = 100 * 86400



df_list = []

for fn in filenames:

    df = pd.read_csv(

        fn, sep=' ', names=['nglitches', 'dF0', 'dF1', 'twoF', 'runTime'])

    df['CLUSTER_ID'] = fn.split('_')[1]

    df_list.append(df)

df = pd.concat(df_list)



colors = ['C0', 'C1']

fig, ax = plt.subplots()

handles = []

labels = []

for ng, c in zip(np.unique(df.nglitches.values), colors):

    print 'ng={}'.format(ng)

    df_temp = df[df.nglitches == ng]

    #df_temp = df_temp[[str(x).isalpha() for x in df_temp.CLUSTER_ID.values]]

    print df_temp.tail()

    Nsamples = len(df_temp)

    MaxtwoF = df_temp.twoF.max()

    print 'Number of samples = ', Nsamples

    print 'Max twoF', MaxtwoF

    print np.any(np.isnan(df_temp.twoF.values))

    ax.hist(df_temp.twoF, bins=40, histtype='stepfilled',

            normed=True, align='mid', alpha=0.5,

            linewidth=1, label='$N_\mathrm{{glitches}}={}$'.format(ng),

            color=c)



    write_to_macro('DirectedMC{}GlitchNoiseOnlyMaximum'.format(ng),

                   '{:1.1f}'.format(MaxtwoF), '../macros.tex')

    write_to_macro('DirectedMC{}GlitchNoiseN'.format(ng),

                   '{:1.0f}'.format(Nsamples), '../macros.tex')



    twoFmax = np.linspace(0, 100, 200)

    ax.plot(twoFmax, pdf_twoFhat(twoFmax, ng, Nsamples,

                                 twoFtildemax=2*MaxtwoF, dtwoF=0.1),

            color=c, label='$N_\mathrm{{glitches}}={}$ predicted'.format(ng))



ax.set_xlabel('$\widehat{2\mathcal{F}}$')

ax.set_xlim(0, 90)

handles, labels = ax.get_legend_handles_labels()

idxs = np.argsort(labels)

ax.legend(np.array(handles)[idxs], np.array(labels)[idxs], frameon=False,

          fontsize=6)

fig.tight_layout()

fig.savefig('glitch_noise_twoF_histogram.png')

Paper/GlitchMCNoiseOnly/submitfile

deleted100644 → 0
+0 −12
Original line number Original line Diff line number Diff line 
Executable= SingleGlitchMCNoiseOnly.sh

Arguments= $(Cluster)_$(Process)

Universe=vanilla

Input=/dev/null

accounting_group=ligo.dev.o1.cw.directedisolatedother.semicoherent 

Output=CollectedOutput/out.$(Cluster).$(Process)

Error=CollectedOutput/err.$(Cluster).$(Process)

Log=CollectedOutput/log.$(Cluster).$(Process)

request_cpus = 1

request_memory = 16 GB



Queue 1806

Paper/Makefile

deleted100644 → 0
+0 −21
Original line number Original line Diff line number Diff line 
SHELL = /bin/bash



main = paper_cw_mcmc.pdf



.PHONY : main clean figclean arxiv jones



main : $(main)



git_tag.tex :

	./git-tag.sh $@



$(main): *.tex git_tag.tex

	pdflatex ${@:.pdf=} && bibtex ${@:.pdf=} && pdflatex ${@:.pdf=} && pdflatex ${@:.pdf=}



clean :

	rm -f git_tag.tex $(main:.pdf=){.aux,.bbl,.blg,.log,.out,.pdf,Notes.bib} $(texfigs) $(texonlyfigs)



arxiv :

	rm -f $(main)

	rm -f git_tag.tex

	tar -cvf arxiv_submission.tar *tex *bib bibliography.bib *pdf *.bbl

Paper/Todo.txt

deleted100644 → 0
+0 −8
Original line number Original line Diff line number Diff line 
- Add macros for the directed MC results

- Add macros for the all-sky MC results

- Adds macros for the examples values

- Investigate cause of all-sky losses again

- Run the transient search for a weaker signal

- Add results from Messenger et al for the transient - is there a fit?

- Weaken basyes factor definitions since they are not used, just put it in the context

Paper/allsky_setup_run_setup.tex

deleted100644 → 0
+0 −8
Original line number Original line Diff line number Diff line 
\begin{tabular}{c|cccccc}

Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline

0 & 20 & 5.0 & 50 & $4.2{\times}10^{2}$ & 7.5 & 56.0 \\

1 & 11 & 9.1 & 50 & $4.4{\times}10^{3}$ & 24.0 & $1.8{\times}10^{2}$ \\

2 & 6 & 16.7 & 50 & $4.2{\times}10^{4}$ & 69.0 & $6.1{\times}10^{2}$ \\

3 & 3 & 33.3 & 50 & $2.8{\times}10^{5}$ & $1.2{\times}10^{2}$ & $2.4{\times}10^{3}$ \\

4 & 1 & 100.0 & 50,50 & $2.0{\times}10^{6}$ & $2.0{\times}10^{2}$ & $1.0{\times}10^{4}$ \\

\end{tabular}

Paper/bibliography.bib

deleted100644 → 0
+0 −555
Original line number Original line Diff line number Diff line 
@ARTICLE{jks1998,

   author = {{Jaranowski}, P. and {Kr{\'o}lak}, A. and {Schutz}, B.~F.},

    title = "{Data analysis of gravitational-wave signals from spinning neutron stars: The signal and its detection}",

  journal = {\prd},

   eprint = {gr-qc/9804014},

 keywords = {Gravitational radiation detectors, mass spectrometers, and other instrumentation and techniques, Gravitational wave detectors and experiments, Mathematical procedures and computer techniques, Pulsars},

     year = 1998,

    month = sep,

   volume = 58,

   number = 6,

      eid = {063001},

    pages = {063001},

      doi = {10.1103/PhysRevD.58.063001},

   adsurl = {http://adsabs.harvard.edu/abs/1998PhRvD..58f3001J},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@ARTICLE{brady1998,

   author = {{Brady}, P.~R. and {Creighton}, T. and {Cutler}, C. and {Schutz}, B.~F.

	},

    title = "{Searching for periodic sources with LIGO}",

  journal = {\prd},

   eprint = {gr-qc/9702050},

 keywords = {Gravitational radiation detectors, mass spectrometers, and other instrumentation and techniques, Gravitational wave detectors and experiments, Mathematical procedures and computer techniques, Pulsars},

     year = 1998,

    month = feb,

   volume = 57,

    pages = {2101-2116},

      doi = {10.1103/PhysRevD.57.2101},

   adsurl = {http://adsabs.harvard.edu/abs/1998PhRvD..57.2101B},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@ARTICLE{brady2000,

   author = {{Brady}, P.~R. and {Creighton}, T.},

    title = "{Searching for periodic sources with LIGO. II. Hierarchical searches}",

  journal = {\prd},

   eprint = {gr-qc/9812014},

 keywords = {Gravitational wave detectors and experiments, Gravitational radiation detectors, mass spectrometers, and other instrumentation and techniques, Mathematical procedures and computer techniques, Pulsars},

     year = 2000,

    month = apr,

   volume = 61,

   number = 8,

      eid = {082001},

    pages = {082001},

      doi = {10.1103/PhysRevD.61.082001},

   adsurl = {http://adsabs.harvard.edu/abs/2000PhRvD..61h2001B},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@ARTICLE{shaltev2013,

   author = {{Shaltev}, M. and {Prix}, R.},

    title = "{Fully coherent follow-up of continuous gravitational-wave candidates}",

  journal = {\prd},

archivePrefix = "arXiv",

   eprint = {1303.2471},

 primaryClass = "gr-qc",

 keywords = {Gravitational waves: theory},

     year = 2013,

    month = apr,

   volume = 87,

   number = 8,

      eid = {084057},

    pages = {084057},

      doi = {10.1103/PhysRevD.87.084057},

   adsurl = {http://adsabs.harvard.edu/abs/2013PhRvD..87h4057S},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}





@ARTICLE{prix2012,

   author = {{Prix}, R. and {Shaltev}, M.},

    title = "{Search for continuous gravitational waves: Optimal StackSlide method at fixed computing cost}",

  journal = {\prd},

archivePrefix = "arXiv",

   eprint = {1201.4321},

 primaryClass = "gr-qc",

 keywords = {Gravitational waves: theory, Gravitational-wave astrophysics, Gravitational wave detectors and experiments, Data analysis: algorithms and implementation, data management},

     year = 2012,

    month = apr,

   volume = 85,

   number = 8,

      eid = {084010},

    pages = {084010},

      doi = {10.1103/PhysRevD.85.084010},

   adsurl = {http://adsabs.harvard.edu/abs/2012PhRvD..85h4010P},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@ARTICLE{krishnan2004,

   author = {{Krishnan}, B. and {Sintes}, A.~M. and {Papa}, M.~A. and {Schutz}, B.~F. and 

	{Frasca}, S. and {Palomba}, C.},

    title = "{Hough transform search for continuous gravitational waves}",

  journal = {\prd},

   eprint = {gr-qc/0407001},

 keywords = {Gravitational wave detectors and experiments, Wave generation and sources, Data analysis: algorithms and implementation, data management, Gravitational radiation detectors, mass spectrometers, and other instrumentation and techniques},

     year = 2004,

    month = oct,

   volume = 70,

   number = 8,

      eid = {082001},

    pages = {082001},

      doi = {10.1103/PhysRevD.70.082001},

   adsurl = {http://adsabs.harvard.edu/abs/2004PhRvD..70h2001K},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}





@ARTICLE{astone2014,

   author = {{Astone}, P. and {Colla}, A. and {D'Antonio}, S. and {Frasca}, S. and 

	{Palomba}, C.},

    title = "{Method for all-sky searches of continuous gravitational wave signals using the frequency-Hough transform}",

  journal = {\prd},

archivePrefix = "arXiv",

   eprint = {1407.8333},

 primaryClass = "astro-ph.IM",

 keywords = {Gravitational waves: theory, Gravitational wave detectors and experiments, Neutron stars},

     year = 2014,

    month = aug,

   volume = 90,

   number = 4,

      eid = {042002},

    pages = {042002},

      doi = {10.1103/PhysRevD.90.042002},

   adsurl = {http://adsabs.harvard.edu/abs/2014PhRvD..90d2002A},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@ARTICLE{cutler2005,

   author = {{Cutler}, C. and {Gholami}, I. and {Krishnan}, B.},

    title = "{Improved stack-slide searches for gravitational-wave pulsars}",

  journal = {\prd},

   eprint = {gr-qc/0505082},

 keywords = {Gravitational wave detectors and experiments, Mathematical procedures and computer techniques, Pulsars},

     year = 2005,

    month = aug,

   volume = 72,

   number = 4,

      eid = {042004},

    pages = {042004},

      doi = {10.1103/PhysRevD.72.042004},

   adsurl = {http://adsabs.harvard.edu/abs/2005PhRvD..72d2004C},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@ARTICLE{prix2009,

   author = {{Prix}, R. and {Krishnan}, B.},

    title = "{Targeted search for continuous gravitational waves: Bayesian versus maximum-likelihood statistics}",

  journal = {Classical and Quantum Gravity},

archivePrefix = "arXiv",

   eprint = {0907.2569},

 primaryClass = "gr-qc",

     year = 2009,

    month = oct,

   volume = 26,

   number = 20,

      eid = {204013},

    pages = {204013},

      doi = {10.1088/0264-9381/26/20/204013},

   adsurl = {http://adsabs.harvard.edu/abs/2009CQGra..26t4013P},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@Article{	  foreman-mackay2013,

  author	= {{Foreman-Mackey}, D. and {Hogg}, D.~W. and {Lang}, D. and

		  {Goodman}, J. },

  title		= "{emcee: The MCMC Hammer}",

  journal	= {PASP},

  archiveprefix	= "arXiv",

  eprint	= {1202.3665},

  primaryclass	= "astro-ph.IM",

  keywords	= {Data Analysis and Techniques},

  year		= 2013,

  month		= mar,

  volume	= 125,

  pages		= {306-312},

  doi		= {10.1086/670067},

  adsurl	= {http://adsabs.harvard.edu/abs/2013PASP..125..306F},

  adsnote	= {Provided by the SAO/NASA Astrophysics Data System}

}



@Article{	  goodman2010,

  author	= {{Goodman}, J. and {Weare}, J.},

  title		= {Ensemble samplers with affine invariance},

  journal	= {Comm. App. Math. Comp. Sci.},

  year		= 2010,

  volume	= 5,

  pages     = {65-80},

}



@article{     swendsen1986,

  title = {Replica Monte Carlo Simulation of Spin-Glasses},

  author = {Swendsen, Robert H. and Wang, Jian-Sheng},

  journal = {Phys. Rev. Lett.},

  volume = {57},

  issue = {21},

  pages = {2607--2609},

  numpages = {0},

  year = {1986},

  month = {Nov},

  publisher = {American Physical Society},

  doi = {10.1103/PhysRevLett.57.2607},

  url = {http://link.aps.org/doi/10.1103/PhysRevLett.57.2607}

}



@Article{goggans2004,

  author = "Goggans, Paul M. and Chi, Ying",

  title = "Using Thermodynamic Integration to Calculate the Posterior Probability in Bayesian Model Selection Problems",

  journal = "AIP Conference Proceedings",

  year = "2004",

  volume = "707",

  number = "1", 

  pages = "59-66",

  url = "http://scitation.aip.org/content/aip/proceeding/aipcp/10.1063/1.1751356",

  doi = "http://dx.doi.org/10.1063/1.1751356" 

}



@ARTICLE{wette2013,

   author = {{Wette}, K. and {Prix}, R.},

    title = "{Flat parameter-space metric for all-sky searches for gravitational-wave pulsars}",

  journal = {\prd},

archivePrefix = "arXiv",

   eprint = {1310.5587},

 primaryClass = "gr-qc",

 keywords = {Gravitational wave detectors and experiments, Gravitational radiation detectors, mass spectrometers, and other instrumentation and techniques, Mathematical procedures and computer techniques, Neutron stars},

     year = 2013,

    month = dec,

   volume = 88,

   number = 12,

      eid = {123005},

    pages = {123005},

      doi = {10.1103/PhysRevD.88.123005},

   adsurl = {http://adsabs.harvard.edu/abs/2013PhRvD..88l3005W},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@ARTICLE{wette2015,

   author = {{Wette}, K.},

    title = "{Parameter-space metric for all-sky semicoherent searches for gravitational-wave pulsars}",

  journal = {\prd},

archivePrefix = "arXiv",

   eprint = {1508.02372},

 primaryClass = "gr-qc",

 keywords = {Gravitational wave detectors and experiments, Gravitational radiation detectors, mass spectrometers, and other instrumentation and techniques, Mathematical procedures and computer techniques, Neutron stars},

     year = 2015,

    month = oct,

   volume = 92,

   number = 8,

      eid = {082003},

    pages = {082003},

      doi = {10.1103/PhysRevD.92.082003},

   adsurl = {http://adsabs.harvard.edu/abs/2015PhRvD..92h2003W},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@ARTICLE{prix2007,

   author = {{Prix}, R.},

    title = "{Template-based searches for gravitational waves: efficient lattice covering of flat parameter spaces}",

  journal = {Classical and Quantum Gravity},

archivePrefix = "arXiv",

   eprint = {0707.0428},

 primaryClass = "gr-qc",

     year = 2007,

    month = oct,

   volume = 24,

    pages = {S481-S490},

      doi = {10.1088/0264-9381/24/19/S11},

   adsurl = {http://adsabs.harvard.edu/abs/2007CQGra..24S.481P},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@article{prix2007metric,

    archivePrefix = {arXiv},

    arxivId = {gr-qc/0606088},

    author = {Prix, Reinhard},

    doi = {10.1103/PhysRevD.75.023004},

    eprint = {0606088},

    file = {:home/greg/Dropbox/Papers/Prix{\_}2007.pdf:pdf},

    issn = {15507998},

    journal = {Physical Review D - Particles, Fields, Gravitation and Cosmology},

    number = {2},

    pages = {1--20},

    primaryClass = {gr-qc},

    title = {{Search for continuous gravitational waves: Metric of the multidetector F-statistic}},

    volume = {75},

    year = {2007}

}





@ARTICLE{cutlershutz2005,

   author = {{Cutler}, C. and {Schutz}, B.~F.},

    title = "{Generalized F-statistic: Multiple detectors and multiple gravitational wave pulsars}",

  journal = {\prd},

   eprint = {gr-qc/0504011},

 keywords = {Gravitational radiation detectors, mass spectrometers, and other instrumentation and techniques, Gravitational wave detectors and experiments, Mathematical procedures and computer techniques, Pulsars},

     year = 2005,

    month = sep,

   volume = 72,

   number = 6,

      eid = {063006},

    pages = {063006},

      doi = {10.1103/PhysRevD.72.063006},

   adsurl = {http://adsabs.harvard.edu/abs/2005PhRvD..72f3006C},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@Article{	  prix2005,

  author	= {{Prix}, R. and {Itoh}, Y.},

  title		= "{Global parameter-space correlations of coherent searches

		  for continuous gravitational waves}",

  journal	= {Classical and Quantum Gravity},

  eprint	= {gr-qc/0504006},

  year		= 2005,

  month		= sep,

  volume	= 22,

  pages		= {1003},

  doi		= {10.1088/0264-9381/22/18/S14},

  adsurl	= {http://adsabs.harvard.edu/abs/2005CQGra..22S1003P},

  adsnote	= {Provided by the SAO/NASA Astrophysics Data System}

}



@misc{lalsuite,

  title = {{LALSuite: FreeSoftware (GPL) Tools for Data-Analysis}},

  author = {{LIGO Scientific Collaboration}},

  year   = 2014,

  note    = {\url{https://www.lsc-group.phys.uwm.edu/daswg/projects/lalsuite.html}}

}



@ARTICLE{allyskyS42008,

   author = {{Abbott}, B. and {Abbott}, R. and {Adhikari}, R. and {Agresti}, J. and 

	{Ajith}, P. and {Allen}, B. and {Amin}, R. and {Anderson}, S.~B. and 

	{Anderson}, W.~G. and {Arain}, M. and et al.},

    title = "{All-sky search for periodic gravitational waves in LIGO S4 data}",

  journal = {\prd},

archivePrefix = "arXiv",

   eprint = {0708.3818},

 primaryClass = "gr-qc",

 keywords = {Gravitational wave detectors and experiments, Data analysis: algorithms and implementation, data management, Gravitational radiation detectors, mass spectrometers, and other instrumentation and techniques, Pulsars},

     year = 2008,

    month = jan,

   volume = 77,

   number = 2,

      eid = {022001},

    pages = {022001},

      doi = {10.1103/PhysRevD.77.022001},

   adsurl = {http://adsabs.harvard.edu/abs/2008PhRvD..77b2001A},

  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

}



@article{pletsch2010,

    archivePrefix = {arXiv},

    arxivId = {1005.0395},

    author = {Pletsch, Holger J.},

    doi = {10.1103/PhysRevD.82.042002},

    eprint = {1005.0395},

    file = {:home/greg/Dropbox/Papers/Pletsch{\_}2014.pdf:pdf},

    issn = {15507998},

    journal = {Physical Review D},

    number = {4},

    title = {{Parameter-space metric of semicoherent searches for continuous gravitational waves}},

    volume = {82},

    year = {2010}

}



@article{leaci2015,

    archivePrefix = {arXiv},

    arxivId = {1502.00914},

    author = {Leaci, Paola and Prix, Reinhard},

    doi = {10.1103/PhysRevD.91.102003},

    eprint = {1502.00914},

    file = {:home/greg/Dropbox/Papers/Leaci{\_}Prix{\_}2015.pdf:pdf},

    issn = {15502368},

    journal = {Physical Review D},

    number = {10},

    pages = {1--25},

    title = {{Directed searches for continuous gravitational waves from binary systems: Parameter-space metrics and optimal scorpius X-1 sensitivity}},

    volume = {91},

    year = {2015}

}



@phdthesis{veitch2007,

  title={Applications of Markov Chain Monte Carlo methods to continuous gravitational wave data analysis},

  author={Veitch, John D},

  year={2007},

  school={University of Glasgow}

}



@article{prix2011,

author = {Prix, R. and Giampanis, S. and Messenger, C.},

doi = {10.1103/PhysRevD.84.023007},

isbn = {0556-2821},

issn = {15507998},

journal = {Physical Review D},

number = {2},

pages = {1--20},

title = {{Search method for long-duration gravitational-wave transients from neutron stars}},

volume = {84},

year = {2011}

}



@article{ashton2015,

archivePrefix = {arXiv},

arxivId = {1410.8044},

author = {Ashton, G. and Jones, D. I. and Prix, R.},

doi = {10.1103/PhysRevD.91.062009},

eprint = {1410.8044},

issn = {15502368},

journal = {Physical Review D - Particles, Fields, Gravitation and Cosmology},

number = {6},

pages = {1--9},

title = {{Effect of timing noise on targeted and narrow-band coherent searches for continuous gravitational waves from pulsars}},

volume = {91},

year = {2015}

}



@article{wette2012,

    arxivId = {1111.5650},

    author = {Wette, Karl},

    doi = {10.1103/PhysRevD.85.042003},

    eprint = {1111.5650},

    file = {:home/greg/Dropbox/Papers/Wette{\_}2012.pdf:pdf},

    isbn = {0556-2821},

    issn = {15507998},

    journal = {Physical Review D},

    number = {4},

    title = {{Estimating the sensitivity of wide-parameter-space searches for gravitational-wave pulsars}},

    volume = {85},

    year = {2012}

}



@article{melatos2008,

    archivePrefix = {arXiv},

    arxivId = {0710.1021},

    author = {Melatos, A. and Peralta, C. and Wyithe, J. S. B.},

    doi = {10.1086/523349},

    eprint = {0710.1021},

    file = {:home/greg/Dropbox/Papers/Mealtos{\_}2007.pdf:pdf},

    issn = {0004-637X},

    journal = {The Astrophysical Journal},

    keywords = {dense matter — pulsars: general — stars: interiors},

    number = {Jensen 1998},

    pages = {1103},

    title = {{Avalanche dynamics of radio pulsar glitches}},

    url = {http://arxiv.org/abs/0710.1021},

    volume = {672},

    year = {2008}

}



@article{espinoza2011,

    archivePrefix = {arXiv},

    arxivId = {1102.1743},

    author = {Espinoza, Crist{\'{o}}bal and Lyne, Andrew and Stappers, Ben and Kramer, Michael},

    doi = {10.1063/1.3615093},

    eprint = {1102.1743},

    file = {:home/greg/Dropbox/Papers/Espinoza{\_}2011.pdf:pdf},

    isbn = {9780735409156},

    issn = {0094243X},

    journal = {AIP Conference Proceedings},

    keywords = {general,glitches,pulsars,stars: neutron},

    number = {March},

    pages = {117--120},

    title = {{Glitches in the rotation of pulsars}},

    volume = {1357},

    year = {2011}

}



@article{cowles1996markov,

  title={Markov chain Monte Carlo convergence diagnostics: a comparative review},

  author={Cowles, Mary Kathryn and Carlin, Bradley P},

  journal={Journal of the American Statistical Association},

  volume={91},

  number={434},

  pages={883--904},

  year={1996},

  publisher={Taylor \& Francis}

}



@article{gelman1992,

  title={Inference from iterative simulation using multiple sequences},

  author={Gelman, Andrew and Rubin, Donald B},

  journal={Statistical science},

  pages={457--472},

  year={1992},

  publisher={JSTOR}

}



@article{brooks1998,

  title={General methods for monitoring convergence of iterative simulations},

  author={Brooks, Stephen P and Gelman, Andrew},

  journal={Journal of computational and graphical statistics},

  volume={7},

  number={4},

  pages={434--455},

  year={1998},

  publisher={Taylor \& Francis}

}



@Article{	  hobbs2010,

  author	= {Hobbs, G. and Lyne, A.~G. and Kramer, M.},

  doi		= {10.1111/j.1365-2966.2009.15938.x},

  issn		= {00358711},

  journal	= {Monthly Notices of the Royal Astronomical Society},

  keywords	= {1 i n t,6000 years of pulsar,archive of pulsar

		  observations,contains over,general,pulsars,ro d u

		  c,rotational history,t i o n,the jodrell bank data,the

		  pulsar timing method},

  month		= feb,

  number	= {2},

  pages		= {1027-1048},

  title		= {{An analysis of the timing irregularities for 366

		  pulsars}},

  url		= {http://doi.wiley.com/10.1111/j.1365-2966.2009.15938.x},

  volume	= {402},

  year		= {2010}

}



@article{ashton2017,

  title     = {{Statistical characterization of pulsar glitches and their

                potential impact on searches for continuous gravitational

                 waves}},

  author = {Ashton, G. and Prix, R. and Jones, D. I.},

  note = {In prep},

  year = {2017}

}



@article{Papa2016,

  title = {{Hierarchical follow-up of subthreshold candidates of an all-sky Einstein@Home search for continuous gravitational waves on LIGO sixth science run data}},

  author = {{Papa}, M.~A. and {Eggenstein}, H-B. and {Walsh}, S. and {Di Palma}, I. and {Allen}, B. and {Astone}, P. and {Bock}, O. and {Creighton}, T.~D. and {Keitel}, D. and {Machenschalk}, B. and {Prix}, R. and {Siemens}, X. and {Singh}, A. and {Zhu}, S.~J. and {Schutz}, B.~F.},

  journal = {Phys. Rev. D},

  volume = {94},

  issue = {12},

  pages = {122006},

  numpages = {14},

  year = {2016},

  month = {Dec},

  publisher = {American Physical Society},

  doi = {10.1103/PhysRevD.94.122006},

  url = {http://link.aps.org/doi/10.1103/PhysRevD.94.122006}

}



@article{dupuis2005,

  title = {Bayesian estimation of pulsar parameters from gravitational wave data},

  author = {Dupuis, R\'ejean J. and Woan, Graham},

  journal = {Phys. Rev. D},

  volume = {72},

  issue = {10},

  pages = {102002},

  numpages = {9},

  year = {2005},

  month = {Nov},

  publisher = {American Physical Society},

  doi = {10.1103/PhysRevD.72.102002},

  url = {http://link.aps.org/doi/10.1103/PhysRevD.72.102002}

}

Paper/definitions.tex

deleted100644 → 0
+0 −34
Original line number Original line Diff line number Diff line 
\newcommand{\fdot}{\dot{f}}

\newcommand{\F}{{\mathcal{F}}}

\newcommand{\twoFhat}{\widehat{2\F}}

\newcommand{\twoFtilde}{\widetilde{2\F}}

\newcommand{\A}{\boldsymbol{\mathcal{A}}}

\newcommand{\blambda}{\boldsymbol{\mathbf{\lambda}}}

\newcommand{\blambdaSignal}{\boldsymbol{\mathbf{\lambda}}^{\rm s}}

\newcommand{\tglitch}{t_{\rm glitch}}

\newcommand{\tstart}{t_{\rm start}}

\newcommand{\tend}{t_{\rm end}}

\newcommand{\Nglitches}{N_{\rm g}}

\newcommand{\Tspan}{T}

\newcommand{\Tcoh}{T_{\rm coh}}

\newcommand{\tref}{t_{\rm ref}}

\newcommand{\Nseg}{N_{\rm seg}}

\newcommand{\Nsteps}{N_{\rm steps}}

\newcommand{\Ntemps}{N_{\rm temps}}

\newcommand{\Nstages}{N_{\rm stages}}

\newcommand{\Nspindown}{N_{\rm spindowns}}

\renewcommand{\H}{\mathcal{H}}

\newcommand{\Hs}{\H_{\rm s}}

\newcommand{\Hn}{\H_{\rm n}}

\newcommand{\ho}{h_0}

\newcommand{\homax}{\ho^{\rm max}}

\newcommand{\Bsn}{B_{\rm S/N}}

\newcommand{\AmplitudePrior}{\Pi_\mathcal{A}}

\newcommand{\mutilde}{\tilde{\mu}}

\newcommand{\Sn}{S_{\rm n}}

\newcommand{\V}{\mathcal{V}}

\newcommand{\Vsky}{\V_{\rm Sky}}

\newcommand{\Vpe}{\V_{\rm PE}}

\newcommand{\smax}{s_{\textrm{max}}}

\newcommand{\fmax}{f_{\textrm{max}}}

Paper/directed_setup_run_setup.tex

deleted100644 → 0
+0 −7
Original line number Original line Diff line number Diff line 
\begin{tabular}{c|cccc}

Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\Vpe$ \\ \hline

0 & 20 & 5.0 & 25 & 56.0 \\

1 & 7 & 14.3 & 25 & $4.5{\times}10^{2}$ \\

2 & 2 & 50.0 & 25 & $5.0{\times}10^{3}$ \\

3 & 1 & 100.0 & 25,25 & $1.0{\times}10^{4}$ \\

\end{tabular}

Paper/follow_up_run_setup.tex

deleted100644 → 0
+0 −11
Original line number Original line Diff line number Diff line 
\begin{tabular}{c|cccccc}

Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline

0 & 100 & 1.0 & 200 & 95.0 & 6.8 & 14.0 \\

1 & 56 & 1.8 & 200 & $9.6{\times}10^{2}$ & 22.0 & 45.0 \\

2 & 31 & 3.2 & 200 & $1.0{\times}10^{4}$ & 69.0 & $1.5{\times}10^{2}$ \\

3 & 17 & 5.9 & 200 & $1.1{\times}10^{5}$ & $2.3{\times}10^{2}$ & $4.8{\times}10^{2}$ \\

4 & 9 & 11.1 & 200 & $1.3{\times}10^{6}$ & $7.5{\times}10^{2}$ & $1.7{\times}10^{3}$ \\

5 & 5 & 20.0 & 200 & $1.1{\times}10^{7}$ & $2.0{\times}10^{3}$ & $5.5{\times}10^{3}$ \\

6 & 2 & 50.0 & 200 & $1.0{\times}10^{8}$ & $3.3{\times}10^{3}$ & $3.1{\times}10^{4}$ \\

7 & 1 & 100.0 & 200,200 & $2.7{\times}10^{8}$ & $4.4{\times}10^{3}$ & $6.2{\times}10^{4}$ \\

\end{tabular}

Paper/git-tag.sh

deleted100755 → 0
+0 −29
Original line number Original line Diff line number Diff line 
#!/bin/bash



FILE="$1"

if [ "x$FILE" = x ]; then

    echo "usage: $0 <file.tex>" >&2

    exit 1

fi

TMPFILE="$FILE.tmp"



gitLOG=`git log -1 --pretty="format:%% generated by $0

\\newcommand{\\commitDATE}{%ai}

\\newcommand{\\commitID}{commitID: %H}

\\newcommand{\\commitIDshort}{commitID: %h}

"`

echo "$gitLOG" > $TMPFILE



diffcmd="git diff -- .";

diff=`$diffcmd`;

if test -n "$diff"; then

    echo "\\newcommand{\\commitSTATUS}{UNCLEAN}" >> $TMPFILE

else

    echo "\\newcommand{\\commitSTATUS}{CLEAN}" >> $TMPFILE

fi



if cmp -s $TMPFILE $FILE; then

    rm -f $TMPFILE

else

    mv -f $TMPFILE $FILE

fi

Paper/macros.tex

deleted100644 → 0
+0 −31
Original line number Original line Diff line number Diff line 
\def\AllSkyMCNoiseN{1.0{\times}10^{4}}

\def\AllSkyMCNoiseOnlyMaximum{59.8}

\def\BasicExampleDeltaFone{1.0{\times}10^{-13}}

\def\BasicExampleDeltaFzero{1.0{\times}10^{-7}}

\def\BasicExampleDepth{10.0}

\def\BasicExampleFone{-1.0{\times}10^{-10}}

\def\BasicExampleFzero{30}

\def\BasicExampleSqrtSn{1.0{\times}10^{-23}}

\def\BasicExampleV{120.0}

\def\BasicExampleVFone{49.0}

\def\BasicExampleVFzero{2.5}

\def\BasicExamplehzero{1.0{\times}10^{-24}}

\def\BasicExamplenburn{50.0}

\def\BasicExamplenprod{50.0}

\def\DirectedMCNoiseN{1.0{\times}10^{4}}

\def\DirectedMCNoiseOnlyMaximum{52.4}

\def\DirectedMConeGlitchNoiseN{10000}

\def\DirectedMConeGlitchNoiseOnlyMaximum{82.6}

\def\DirectedMCtwoGlitchNoiseN{9625}

\def\DirectedMCtwoGlitchNoiseOnlyMaximum{87.8}

\def\SingleGlitchDepth{10.0}

\def\SingleGlitchFCMismatch{0.7}

\def\SingleGlitchFCtwoF{718.5}

\def\SingleGlitchFCtwoFPM{2645.7}

\def\SingleGlitchFone{-1e-10}

\def\SingleGlitchFzero{30.0}

\def\SingleGlitchGlitchtwoF{2622.1}

\def\SingleGlitchR{0.5}

\def\SingleGlitchTspan{100.0}

\def\SingleGlitchdeltaFone{$-2.9{\times}10^{-13}$}

\def\SingleGlitchdeltaFzero{$3.2{\times}10^{-6}$}

Paper/paper_cw_mcmc.tex

deleted100644 → 0
+0 −1417

File deleted.

Preview size limit exceeded, changes collapsed.

README.md

+40 −95
Original line number Original line Diff line number Diff line
@@ -3,101 +3,46 @@ 
This is a python package providing an interface to perform F-statistic based

This is a python package providing an interface to perform F-statistic based

continuous gravitational wave (CW) searches.

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:

## Installation





**===> https://github.com/PyFstat/PyFstat <===**

### `python` installation



The scripts are written in `python 2.7+` and therefore require a working

Please also refer to the github repository for installation and usage instructions

`python` installation. While many systems come with a system wide python

and to report issues.

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

### Contributors

modules. One method to do this, is to use the `conda` system, either through



the stripped down [miniconda](http://conda.pydata.org/miniconda.html)

* Greg Ashton

installation, or the full-featured

* David Keitel

[anaconda](https://www.continuum.io/downloads) (these are essentially the

* Reinhard Prix

same, but the `anaconda` version installs a variety of useful packages such as

* Karl Wette

`numpy` and `scipy` by default).

* Sylvia Zhu





The fastest/easiest method is to follow your OS instructions

## Citing this work

[here](https://conda.io/docs/install/quick.html) which will install Miniconda.





If you use `PyFstat` in a publication we would appreciate if you cite both the

For the rest of this tutorial, we will make use of `pip` to install modules (

original paper introducing the code

not all packages can be installed with `conda` and for those using alternatives

(the [ADS page can be found here](http://adsabs.harvard.edu/abs/2018arXiv180205450A))

to `conda`, `pip` is more universal).

and the DOI for the software itself.



The latest version released on Zenodo from this repository here was v1.3,

This can be installed with

see the following Bibtex entry.

```

Or see the [new github repository](https://github.com/PyFstat/PyFstat)

$ conda install pip

for the latest information.

```

```



@misc{pyfstat,

### Clone the repository

  author       = {Ashton, Gregory and



                  Keitel, David and

In a terminal, to clone the directory:

                  Prix, Reinhard},



  title        = {PyFstat-v1.3},

```

  month        = jan,

$ git clone git@gitlab.aei.uni-hannover.de:GregAshton/PyFstat.git

  year         = 2020,

```

  publisher    = {Zenodo},



  version      = {1.3},

### Dependencies

  doi          = {10.5281/zenodo.3620861},



  url          = {https://doi.org/10.5281/zenodo.3620861}

`pyfstat` makes uses the following external python modules:

  note         = {\url{https://doi.org/10.5281/zenodo.3620861}}



}

* [numpy](http://www.numpy.org/)

* [matplotlib](http://matplotlib.org/) >= 1.4

* [scipy](https://www.scipy.org/)

* [emcee](http://dan.iel.fm/emcee/current/)

* [corner](https://pypi.python.org/pypi/corner/)

* [dill](https://pypi.python.org/pypi/dill)



*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



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`.

examples/MCMC_examples/fully_coherent_search_using_MCMC.py

0 → 100644
+81 −0
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

duration = 100 * 86400

tend = tstart + duration



# Properties of the signal

F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)

tref = 0.5 * (tstart + tend)



depth = 10

h0 = sqrtSX / depth

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,

)

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))



DeltaF0 = 1e-7

DeltaF1 = 1e-13

VF0 = (np.pi * duration * DeltaF0) ** 2 / 3.0

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.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

log10beta_min = -0.5

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,

)

mcmc.transform_dictionary = dict(

    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

0 → 100644
+82 −0
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

duration = 100 * 86400

tend = tstart + duration



# Properties of the signal

F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)

tref = 0.5 * (tstart + tend)



depth = 10

h0 = sqrtSX / depth

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,

)

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))



DeltaF0 = 1e-7

DeltaF1 = 1e-13

VF0 = (np.pi * duration * DeltaF0) ** 2 / 3.0

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.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

log10beta_min = -1

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,

)

mcmc.transform_dictionary = dict(

    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/Makefile

deleted100644 → 0
+0 −8
Original line number Original line Diff line number Diff line 
basic_sft = data/H-4800_H1_1800SFT_basic-1000000000-8640000.sft

glitch_sft = data/H-4800_H1_1800SFT_glitch-1000000000-8640000.sft



data/fully_coherent_corner.png : $(basic_sft) fully_coherent_search.py

	python fully_coherent_search.py



$(basic_sft) $(glitch_sft): make_fake_data.py

	python make_fake_data.py

examples/computing_the_Bayes_factor.py

deleted100644 → 0
+0 −34
Original line number Original line Diff line number Diff line 
from pyfstat import MCMCSearch



F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = 5e-3

Delta = 6e-2

tref = 362750407.0



tstart = 1000000000

duration = 100*86400

tend = tstart + duration



theta_prior = {'F0': {'type': 'unif', 'lower': F0*(1-1e-6), 'upper': F0*(1+1e-6)},

               'F1': {'type': 'unif', 'lower': F1*(1+1e-2), 'upper': F1*(1-1e-2)},

               'F2': F2,

               'Alpha': Alpha,

               'Delta': Delta

               }



ntemps = 20

log10temperature_min = -2

nwalkers = 100

nsteps = [500, 500]



mcmc = MCMCSearch(label='computing_the_Bayes_factor', outdir='data', 

                  sftfilepath='data/*basic*sft', theta_prior=theta_prior,

                  tref=tref, tstart=tstart, tend=tend, nsteps=nsteps,

                  nwalkers=nwalkers, ntemps=ntemps,

                  log10temperature_min=log10temperature_min)

mcmc.run()

mcmc.plot_corner(add_prior=True)

mcmc.print_summary()

mcmc.compute_evidence()

examples/follow_up.py

deleted100644 → 0
+0 −34
Original line number Original line Diff line number Diff line 
import pyfstat



F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = 5e-3

Delta = 6e-2

tref = 362750407.0



tstart = 1000000000

duration = 100*86400

tend = tstart + duration



theta_prior = {'F0': {'type': 'unif', 'lower': F0*(1-1e-6), 'upper': F0*(1+1e-5)},

               'F1': {'type': 'unif', 'lower': F1*(1+1e-2), 'upper': F1*(1-1e-2)},

               'F2': F2,

               'Alpha': Alpha,

               'Delta': Delta

               }



ntemps = 1

log10temperature_min = -1

nwalkers = 100

run_setup = [(1000, 50), (1000, 25), (1000, 1, False), 

             ((500, 500), 1, True)]



mcmc = pyfstat.MCMCFollowUpSearch(

    label='follow_up', outdir='data',

    sftfilepath='data/*basic*sft', theta_prior=theta_prior, tref=tref,

    minStartTime=tstart, maxStartTime=tend, nwalkers=nwalkers,

    ntemps=ntemps, log10temperature_min=log10temperature_min)

mcmc.run(run_setup)

mcmc.plot_corner(add_prior=True)

mcmc.print_summary()

examples/followup_examples/semi_coherent_directed_follow_up.py

0 → 100644
+93 −0
Original line number Original line Diff line number Diff line 
import pyfstat

import numpy as np

import matplotlib.pyplot as plt



F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)



# Properties of the GW data

sqrtSX = 1e-23

tstart = 1000000000

duration = 100 * 86400

tend = tstart + duration

tref = 0.5 * (tstart + tend)



depth = 40

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,

)

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))



# 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.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

log10beta_min = -0.5

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,

)



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,

)

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)

fig.tight_layout()

fig.savefig("{}/{}_walkers.png".format(mcmc.outdir, mcmc.label), dpi=400)

examples/fully_coherent_search_using_MCMC_on_glitching_data.py

deleted100644 → 0
+0 −35
Original line number Original line Diff line number Diff line 
import numpy as np

from pyfstat import MCMCSearch



F0 = 30.0

F1 = -1e-10

F2 = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)



tref = 362750407.0



tstart = 1000000000

duration = 100*86400

tend = tstart + duration



theta_prior = {'F0': {'type': 'unif', 'lower': F0-1e-4, 'upper': F0+1e-4},

               'F1': {'type': 'unif', 'lower': F1*(1+1e-3), 'upper': F1*(1-1e-3)},

               'F2': F2,

               'Alpha': Alpha,

               'Delta': Delta

               }



ntemps = 2

log10temperature_min = -0.01

nwalkers = 100

nsteps = [500, 500]



mcmc = MCMCSearch('fully_coherent_search_using_MCMC_on_glitching_data', 'data',

                  sftfilepath='data/*_glitch*.sft',

                  theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,

                  nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,

                  log10temperature_min=log10temperature_min)

mcmc.run()

mcmc.plot_corner(add_prior=True)

mcmc.print_summary()

examples/make_fake_data.py→examples/glitch_examples/make_simulated_data.py

+91 −0
Original line number Original line Diff line number Diff line 
from pyfstat import Writer

from pyfstat import Writer, GlitchWriter

import numpy as np

import numpy as np





outdir = "data"

# First, we generate data with a reasonably strong smooth signal

# First, we generate data with a reasonably strong smooth signal





# Define parameters of the Crab pulsar as an example

# Define parameters of the Crab pulsar as an example
@@ -9,48 +10,82 @@ F1 = -1e-10 
F2 = 0

F2 = 0

Alpha = np.radians(83.6292)

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)

Delta = np.radians(22.0144)

tref = 362750407.0





# Signal strength

# Signal strength

h0 = 1e-23

h0 = 5e-24





# Properties of the GW data

# Properties of the GW data

sqrtSX = 1e-22

sqrtSX = 1e-22

tstart = 1000000000

tstart = 1000000000

duration = 100*86400

duration = 50 * 86400

tend = tstart + duration

tend = tstart + duration

tref = tstart + 0.5 * duration





data = Writer(

data = Writer(

    label='basic', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,

    label="0_glitch",

    F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX)

    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()

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)



# Next, taking the same signal parameters, we include a glitch half way through

# Next, taking the same signal parameters, we include a glitch half way through

dtglitch = duration / 2.0

dtglitch = duration / 2.0

delta_F0 = 4e-5

delta_F0 = 5e-6

delta_F1 = 0

delta_F1 = 0





glitch_data = Writer(

glitch_data = GlitchWriter(

    label='glitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,

    label="1_glitch",

    F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,

    outdir=outdir,

    dtglitch=dtglitch, delta_F0=delta_F0, delta_F1=delta_F1)

    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()

glitch_data.make_data()





# Making data with two glitches

# Making data with two glitches





dtglitch = [duration/4.0, 4*duration/5.0]

dtglitch_2 = [duration / 4.0, 4 * duration / 5.0]

delta_phi = [0, 0]

delta_phi_2 = [0, 0]

delta_F0 = [4e-6, 3e-7]

delta_F0_2 = [4e-6, 3e-7]

delta_F1 = [0, 0]

delta_F1_2 = [0, 0]

delta_F2 = [0, 0]

delta_F2_2 = [0, 0]





two_glitch_data = Writer(

two_glitch_data = GlitchWriter(

    label='twoglitch', outdir='data', tref=tref, tstart=tstart, F0=F0, F1=F1,

    label="2_glitch",

    F2=F2, duration=duration, Alpha=Alpha, Delta=Delta, h0=h0, sqrtSX=sqrtSX,

    outdir=outdir,

    dtglitch=dtglitch, delta_phi=delta_phi, delta_F0=delta_F0,

    tref=tref,

    delta_F1=delta_F1, delta_F2=delta_F2)

    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()
 
two_glitch_data.make_data()

examples/glitch_examples/semicoherent_glitch_robust_directed_MCMC_search_on_1_glitch.py

0 → 100644
+94 −0
Original line number Original line Diff line number Diff line 
import numpy as np

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,

)



plt.style.use("./paper.mplstyle")



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.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

log10beta_min = -0.5

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]"



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),

    hist_kwargs=dict(lw=1.5, zorder=-1),

    truth_color="C3",

)



mcmc.print_summary()



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

0 → 100644
+99 −0
Original line number Original line Diff line number Diff line 
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,

)

import time



try:

    from gridcorner import gridcorner

except ImportError:

    raise ImportError(

        "Python module 'gridcorner' not found, please install from "

        "https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner"

    )



label = "semicoherent_glitch_robust_directed_grid_search_on_1_glitch"



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.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]



max_delta_F0 = 1e-5

tglitchs = [tstart + 0.1 * duration, tstart + 0.9 * duration, 0.8 * float(duration) / N]

delta_F0s = [0, max_delta_F0, max_delta_F0 / N]

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,

)

search.run()

dT = time.time() - t1



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.0 - dtglitch / 86400.0



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}}$",

]

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")





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

0 → 100644
+47 −0
Original line number Original line Diff line number Diff line 
import numpy as np

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")



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.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

log10beta_min = -0.5

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}$"

)



mcmc.run()

mcmc.plot_corner()

mcmc.print_summary()

examples/glitch_robust_search.py

deleted100644 → 0
+0 −54
Original line number Original line Diff line number Diff line 
import numpy as np

import pyfstat



outdir = 'data'

label = 'glitch_robust_search'



# Properties of the GW data

tstart = 1000000000

Tspan = 60 * 86400



# Fixed properties of the signal

F0s = 30

F1s = -1e-8

F2s = 0

Alpha = np.radians(83.6292)

Delta = np.radians(22.0144)



tref = tstart + .5 * Tspan



sftfilepath = 'data/*glitching_signal*sft'



F0_width = np.sqrt(3)/(np.pi*Tspan)

F1_width = np.sqrt(45/4.)/(np.pi*Tspan**2)

DeltaF0 = 50 * F0_width

DeltaF1 = 50 * F1_width



theta_prior = {'F0': {'type': 'unif',

                      'lower': F0s-DeltaF0,

                      'upper': F0s+DeltaF0},

               'F1': {'type': 'unif',

                      'lower': F1s-DeltaF1,

                      'upper': F1s+DeltaF1},

               'F2': F2s,

               'delta_F0': {'type': 'unif',

                            'lower': 0,

                            'upper': 1e-5},

               'delta_F1': {'type': 'unif',

                            'lower': -1e-11,

                            'upper': 1e-11},

               'tglitch': {'type': 'unif',

                           'lower': tstart+0.1*Tspan,

                           'upper': tstart+0.9*Tspan},

               'Alpha': Alpha,

               'Delta': Delta,

               }



search = pyfstat.MCMCGlitchSearch(

    label=label, outdir=outdir, sftfilepath=sftfilepath,

    theta_prior=theta_prior, nglitch=1, tref=tref, nsteps=[500, 500],

    ntemps=3, log10temperature_min=-0.5, minStartTime=tstart,

    maxStartTime=tstart+Tspan)

search.run()

search.plot_corner(label_offset=0.8, add_prior=True)

search.print_summary()

pyfstat/__init__.py

+31 −5

File changed.

Preview size limit exceeded, changes collapsed.

pyfstat/core.py

+1037 −684

File changed.

Preview size limit exceeded, changes collapsed.

pyfstat/make_sfts.py

0 → 100644
+810 −0

File added.

Preview size limit exceeded, changes collapsed.

requirements.txt

+4 −3

File changed.

Preview size limit exceeded, changes collapsed.

setup.cfg

0 → 100644
+8 −0

File added.

Preview size limit exceeded, changes collapsed.

setup.py

+97 −8

File changed.

Preview size limit exceeded, changes collapsed.

tests.py

+478 −174

File changed.

Preview size limit exceeded, changes collapsed.

AEI Privacy Policy Imprint