Adds new examples on using the Glitch search with one and two glitches

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()
+```
+

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()
+

examples/semi_coherent_glitch_search.py

+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()

examples/semi_coherent_two_glitch_search.py

+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()

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):
         """