diff --git a/docs/img/semi_coherent_glitch_search_corner.png b/docs/img/semi_coherent_glitch_search_corner.png new file mode 100644 index 0000000000000000000000000000000000000000..870925ceb1983be9010fdd6c6245ec6a5ae557b2 Binary files /dev/null and b/docs/img/semi_coherent_glitch_search_corner.png differ diff --git a/docs/img/semi_coherent_two_glitch_search_corner.png b/docs/img/semi_coherent_two_glitch_search_corner.png new file mode 100644 index 0000000000000000000000000000000000000000..2f03f7a926d0712d50ab7b7a8aaffb5f9e5938fb Binary files /dev/null and b/docs/img/semi_coherent_two_glitch_search_corner.png differ diff --git a/docs/make_fake_data.md b/docs/make_fake_data.md index 0869431aeffa2182457384c6bee3e751edc7f75c..53536009cfef0b753b7c9db7672759fbff8ab86e 100644 --- a/docs/make_fake_data.md +++ b/docs/make_fake_data.md @@ -142,3 +142,29 @@ So calling ``` Notice that the predicted value will be the same for both sets of data. + +## Making data with multiple glitches + +Finally, one can also use the `Writer` to generate data with multiple glitches. +To do this, simply pass in `dtglitch`, `delta_phi`, `delta_F0`, `delta_F1`, and +`delta_F2` as arrays (with a length equal to the number of glitches). Note +that all these must be of equal length. Moreover, the glitches are applied +sequentially and additively as implemented +`pyfstat.BaseSearchClass.calculate_thetas`. Here is an example with two +glitches, one a quarter of the way through and the other a fifth from the end. + +``` +dtglitch = [duration/4.0, 4*duration/5.0] +delta_phi = [0, 0] +delta_F0 = [0.4e-5, 0.3e-6] +delta_F1 = [0, 0] +delta_F2 = [0, 0] + +two_glitch_data = Writer( + label='two_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_phi=delta_phi, delta_F0=delta_F0, + delta_F1=delta_F1, delta_F2=delta_F2) +two_glitch_data.make_data() +``` + diff --git a/examples/make_fake_data.py b/examples/make_fake_data.py index b6315e1140817fc92e581fcbef39e63572dd2e51..f402907378c2674d9c641076216f6826dc6dcbec 100644 --- a/examples/make_fake_data.py +++ b/examples/make_fake_data.py @@ -40,3 +40,19 @@ glitch_data.make_data() # The predicted twoF, given by lalapps_predictFstat can be accsessed by print data.predict_fstat() + +# Making data with two glitches + +dtglitch = [duration/4.0, 4*duration/5.0] +delta_phi = [0, 0] +delta_F0 = [0.4e-5, 0.3e-6] +delta_F1 = [0, 0] +delta_F2 = [0, 0] + +two_glitch_data = Writer( + label='two_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_phi=delta_phi, delta_F0=delta_F0, + delta_F1=delta_F1, delta_F2=delta_F2) +two_glitch_data.make_data() + diff --git a/examples/semi_coherent_glitch_search.py b/examples/semi_coherent_glitch_search.py new file mode 100644 index 0000000000000000000000000000000000000000..b4dc3b727de34433042c3e85bbc8cf0b9520c538 --- /dev/null +++ b/examples/semi_coherent_glitch_search.py @@ -0,0 +1,37 @@ +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': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)}, + 'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)}, + 'F2': F2, + 'Alpha': Alpha, + 'Delta': Delta, + 'delta_F0': {'type': 'halfnorm', 'loc': 0, + 'scale': 1e-5*F0}, + 'delta_F1': 0, + 'tglitch': {'type': 'unif', + 'lower': tstart+0.1*duration, + 'upper': tstart+0.9*duration}, + } + +nwalkers = 500 +nsteps = [1000, 1000, 1000] + +mcmc = pyfstat.MCMCGlitchSearch( + 'semi_coherent_glitch_search', 'data', sftlabel='glitch', sftdir='data', + theta_prior=theta_prior, tref=tref, tstart=tstart, tend=tend, + nsteps=nsteps, nwalkers=nwalkers, scatter_val=1e-10, nglitch=1) + +mcmc.run() +mcmc.plot_corner(add_prior=True) +mcmc.print_summary() diff --git a/examples/semi_coherent_two_glitch_search.py b/examples/semi_coherent_two_glitch_search.py new file mode 100644 index 0000000000000000000000000000000000000000..8968d815c7c2051ec33b0a41c772ae44e25184bd --- /dev/null +++ b/examples/semi_coherent_two_glitch_search.py @@ -0,0 +1,37 @@ +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': 'norm', 'loc': F0, 'scale': abs(1e-6*F0)}, + 'F1': {'type': 'norm', 'loc': F1, 'scale': abs(1e-6*F1)}, + 'F2': F2, + 'Alpha': Alpha, + 'Delta': Delta, + 'delta_F0': {'type': 'halfnorm', 'loc': 0, + 'scale': 1e-7*F0}, + 'delta_F1': 0, + 'tglitch': {'type': 'unif', + 'lower': tstart+0.01*duration, + 'upper': tstart+0.99*duration}, + } + +nwalkers = 100 +nsteps = [500, 500, 500] + +mcmc = pyfstat.MCMCGlitchSearch( + 'semi_coherent_two_glitch_search', 'data', sftlabel='two_glitch', + sftdir='data', theta_prior=theta_prior, tref=tref, tstart=tstart, + tend=tend, nsteps=nsteps, nwalkers=nwalkers, scatter_val=1e-10, nglitch=2) + +mcmc.run() +mcmc.plot_corner(add_prior=True, tglitch_ratio=True, figsize=(10, 10)) +mcmc.print_summary() diff --git a/pyfstat.py b/pyfstat.py index 15417cfd1c55540b124f6c628d7072d4e606582d..592c7029249955c051be4c131c7256a69fd75ca5 100755 --- a/pyfstat.py +++ b/pyfstat.py @@ -600,21 +600,24 @@ class MCMCSearch(BaseSearchClass): self.lnlikes = lnlikes self.save_data(sampler, samples, lnprobs, lnlikes) - def plot_corner(self, corner_figsize=(7, 7), deltat=False, + def plot_corner(self, figsize=(7, 7), tglitch_ratio=False, add_prior=False, nstds=None, label_offset=0.4, **kwargs): fig, axes = plt.subplots(self.ndim, self.ndim, - figsize=corner_figsize) + figsize=figsize) samples_plt = copy.copy(self.samples) theta_symbols_plt = copy.copy(self.theta_symbols) theta_symbols_plt = [s.replace('_{glitch}', r'_\textrm{glitch}') for s in theta_symbols_plt] - if deltat: - samples_plt[:, self.theta_keys.index('tglitch')] -= self.tref - theta_symbols_plt[self.theta_keys.index('tglitch')] = ( - r'$t_{\textrm{glitch}} - t_{\textrm{ref}}$') + if tglitch_ratio: + for j, k in enumerate(self.theta_keys): + if k == 'tglitch': + s = samples_plt[:, j] + samples_plt[:, j] = (s - self.tstart)/( + self.tend - self.tstart) + theta_symbols_plt[j] = r'$R_{\textrm{glitch}}$' if type(nstds) is int and 'range' not in kwargs: _range = [] @@ -976,7 +979,7 @@ class MCMCGlitchSearch(MCMCSearch): def __init__(self, label, outdir, sftlabel, sftdir, theta_prior, tref, tstart, tend, nglitch=1, nsteps=[100, 100, 100], nwalkers=100, ntemps=1, log10temperature_min=-5, theta_initial=None, - scatter_val=1e-4, dtglitchmin=20*86400, detector=None, + scatter_val=1e-4, dtglitchmin=1*86400, detector=None, minCoverFreq=None, maxCoverFreq=None, earth_ephem=None, sun_ephem=None): """