diff --git a/examples/Makefile b/examples/Makefile
deleted file mode 100644
index 888f4dc85e2b88bef574605d288fd4b926a09506..0000000000000000000000000000000000000000
--- a/examples/Makefile
+++ /dev/null
@@ -1,8 +0,0 @@
-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
diff --git a/examples/computing_the_Bayes_factor.py b/examples/computing_the_Bayes_factor.py
deleted file mode 100644
index 7ddead59a4bdcb635cb533a92ef87486ac52cb78..0000000000000000000000000000000000000000
--- a/examples/computing_the_Bayes_factor.py
+++ /dev/null
@@ -1,34 +0,0 @@
-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
-log10beta_min = -2
-nwalkers = 100
-nsteps = [500, 500]
-
-mcmc = MCMCSearch(label='computing_the_Bayes_factor', outdir='data', 
-                  sftfilepattern='data/*basic*sft', theta_prior=theta_prior,
-                  tref=tref, tstart=tstart, tend=tend, nsteps=nsteps,
-                  nwalkers=nwalkers, ntemps=ntemps,
-                  log10beta_min=log10beta_min)
-mcmc.run()
-mcmc.plot_corner(add_prior=True)
-mcmc.print_summary()
-mcmc.compute_evidence()
diff --git a/examples/follow_up.py b/examples/follow_up.py
deleted file mode 100644
index b9efbaaa58ad69a6b26892423566ca874e0359e3..0000000000000000000000000000000000000000
--- a/examples/follow_up.py
+++ /dev/null
@@ -1,35 +0,0 @@
-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
-log10beta_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',
-    sftfilepattern='data/*basic*sft', theta_prior=theta_prior, tref=tref,
-    minStartTime=tstart, maxStartTime=tend, nwalkers=nwalkers,
-    ntemps=ntemps, log10beta_min=log10beta_min)
-mcmc.run(run_setup, gen_tex_table=True)
-#mcmc.run(Nsegs0=50)
-mcmc.plot_corner(add_prior=True)
-mcmc.print_summary()
diff --git a/examples/fully_coherent_search_using_MCMC.py b/examples/fully_coherent_search_using_MCMC.py
index 994c8cfa97aaff66a4cc741c0ccbdcc28369acac..82fd4dd37425573a127c03e7580e7a2a2d956d12 100644
--- a/examples/fully_coherent_search_using_MCMC.py
+++ b/examples/fully_coherent_search_using_MCMC.py
@@ -17,10 +17,11 @@ tref = .5*(tstart+tend)
 
 depth = 10
 h0 = sqrtSX / depth
-data_label = 'fully_coherent_search_using_MCMC'
+label = 'fully_coherent_search_using_MCMC'
+outdir = 'data'
 
 data = pyfstat.Writer(
-    label=data_label, outdir='data', tref=tref,
+    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()
@@ -49,13 +50,13 @@ theta_prior = {'F0': {'type': 'unif',
 ntemps = 1
 log10beta_min = -1
 nwalkers = 100
-nsteps = [1000, 1000]
+nsteps = [300, 300]
 
 mcmc = pyfstat.MCMCSearch(
-    label='fully_coherent_search_using_MCMC', outdir='data',
-    sftfilepattern='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,
-    minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,
-    ntemps=ntemps, log10beta_min=log10beta_min)
-mcmc.run(context='paper', subtractions=[30, -1e-10], c=2)
+    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.run(subtractions=[F0, F1])
 mcmc.plot_corner(add_prior=True)
 mcmc.print_summary()
diff --git a/examples/fully_coherent_search_using_MCMC_on_glitching_data.py b/examples/glitch_examples/fully_coherent_search_using_MCMC_on_glitching_data.py
similarity index 100%
rename from examples/fully_coherent_search_using_MCMC_on_glitching_data.py
rename to examples/glitch_examples/fully_coherent_search_using_MCMC_on_glitching_data.py
diff --git a/examples/glitch_robust_search.py b/examples/glitch_examples/glitch_robust_search.py
similarity index 100%
rename from examples/glitch_robust_search.py
rename to examples/glitch_examples/glitch_robust_search.py
diff --git a/examples/glitch_robust_search_make_simulated_data.py b/examples/glitch_examples/glitch_robust_search_make_simulated_data.py
similarity index 100%
rename from examples/glitch_robust_search_make_simulated_data.py
rename to examples/glitch_examples/glitch_robust_search_make_simulated_data.py
diff --git a/examples/make_fake_data.py b/examples/glitch_examples/make_fake_data.py
similarity index 100%
rename from examples/make_fake_data.py
rename to examples/glitch_examples/make_fake_data.py
diff --git a/examples/semi_coherent_glitch_search_using_MCMC.py b/examples/glitch_examples/semi_coherent_glitch_search_using_MCMC.py
similarity index 100%
rename from examples/semi_coherent_glitch_search_using_MCMC.py
rename to examples/glitch_examples/semi_coherent_glitch_search_using_MCMC.py
diff --git a/examples/semi_coherent_twoglitch_search_using_MCMC.py b/examples/glitch_examples/semi_coherent_twoglitch_search_using_MCMC.py
similarity index 100%
rename from examples/semi_coherent_twoglitch_search_using_MCMC.py
rename to examples/glitch_examples/semi_coherent_twoglitch_search_using_MCMC.py
diff --git a/examples/grided_frequency_search.py b/examples/grid_examples/grided_frequency_search.py
similarity index 100%
rename from examples/grided_frequency_search.py
rename to examples/grid_examples/grided_frequency_search.py
diff --git a/examples/sliding_window.py b/examples/other_examples/sliding_window.py
similarity index 100%
rename from examples/sliding_window.py
rename to examples/other_examples/sliding_window.py
diff --git a/examples/transient_search_using_MCMC.py b/examples/other_examples/transient_search_using_MCMC.py
similarity index 100%
rename from examples/transient_search_using_MCMC.py
rename to examples/other_examples/transient_search_using_MCMC.py
diff --git a/examples/transient_search_using_MCMC_make_simulated_data.py b/examples/other_examples/transient_search_using_MCMC_make_simulated_data.py
similarity index 100%
rename from examples/transient_search_using_MCMC_make_simulated_data.py
rename to examples/other_examples/transient_search_using_MCMC_make_simulated_data.py
diff --git a/examples/twoF_cumulative.py b/examples/other_examples/twoF_cumulative.py
similarity index 100%
rename from examples/twoF_cumulative.py
rename to examples/other_examples/twoF_cumulative.py
diff --git a/examples/semi_coherent_directed_follow_up.py b/examples/semi_coherent_directed_follow_up.py
new file mode 100644
index 0000000000000000000000000000000000000000..a60d9ce949b0fd84ca92042becd7d4a1f6854125
--- /dev/null
+++ b/examples/semi_coherent_directed_follow_up.py
@@ -0,0 +1,73 @@
+import pyfstat
+import numpy as np
+import matplotlib.pyplot as plt
+
+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 = 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 = 100
+DeltaF0 = VF0 * np.sqrt(3)/(np.pi*duration)
+DeltaF1 = VF1 * np.sqrt(180)/(np.pi*duration**2)
+DeltaAlpha = 1e-2
+DeltaDelta = 1e-2
+theta_prior = {'F0': {'type': 'unif', 'lower': F0-DeltaF0/2.,
+                      'upper': F0+DeltaF0/2},
+               'F1': {'type': 'unif', 'lower': F1-DeltaF1/2.,
+                      'upper': F1+DeltaF1/2},
+               'F2': F2,
+               'Alpha': Alpha,
+               'Delta': Delta
+               }
+
+ntemps = 2
+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, subtractions=[F0, F1], labelpad=0.01,
+    plot_det_stat=False, return_fig=True, context='paper', 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)
diff --git a/examples/semi_coherent_search_using_MCMC.py b/examples/semi_coherent_search_using_MCMC.py
index ac10782ab9de4df64920d4d061f6f4507b24a1d6..1049a77da96807af5eec32ac528aac1c9c7e2073 100644
--- a/examples/semi_coherent_search_using_MCMC.py
+++ b/examples/semi_coherent_search_using_MCMC.py
@@ -1,18 +1,47 @@
 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 = 362750407.0
+Alpha = np.radians(83.6292)
+Delta = np.radians(22.0144)
+tref = .5*(tstart+tend)
 
-tstart = 1000000000
-duration = 100*86400
-tend = tstart + duration
+depth = 10
+h0 = sqrtSX / depth
+label = 'semi_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.
+print '\nV={:1.2e}, VF0={:1.2e}, VF1={:1.2e}\n'.format(VF0*VF1, VF0, VF1)
 
-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)},
+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
@@ -21,13 +50,14 @@ theta_prior = {'F0': {'type': 'unif', 'lower': F0*(1-1e-6), 'upper': F0*(1+1e-5)
 ntemps = 1
 log10beta_min = -1
 nwalkers = 100
-nsteps = [500, 500, 500]
+nsteps = [300, 300]
 
 mcmc = pyfstat.MCMCSemiCoherentSearch(
-    label='semi_coherent_search_using_MCMC', outdir='data', nsegs=20,
-    sftfilepattern='data/*basic*sft', theta_prior=theta_prior, tref=tref,
-    minStartTime=tstart, maxStartTime=tend, nsteps=nsteps, nwalkers=nwalkers,
-    ntemps=ntemps, log10beta_min=log10beta_min)
+    label=label, outdir=outdir, nsegs=3,
+    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.run()
 mcmc.plot_corner(add_prior=True)
 mcmc.print_summary()
diff --git a/examples/using_initialisation.py b/examples/using_initialisation.py
new file mode 100644
index 0000000000000000000000000000000000000000..5e7ea5e019f0d8b3abf155c41d861d6b7a0ad825
--- /dev/null
+++ b/examples/using_initialisation.py
@@ -0,0 +1,64 @@
+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 = .5*(tstart+tend)
+
+depth = 10
+h0 = sqrtSX / depth
+label = 'using_initialisation'
+outdir = 'data'
+
+data = pyfstat.Writer(
+    label=label, outdir=outdir, tref=tref,
+    tstart=tstart, F0=F0, F1=F1, F2=F2, duration=duration, Alpha=Alpha,
+    Delta=Delta, h0=h0, sqrtSX=sqrtSX)
+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
+log10beta_min = -1
+nwalkers = 100
+nsteps = [100, 100]
+
+mcmc = pyfstat.MCMCSearch(
+    label=label, outdir=outdir,
+    sftfilepattern='{}/*{}*sft'.format(outdir, label),
+    theta_prior=theta_prior, tref=tref, minStartTime=tstart, maxStartTime=tend,
+    nsteps=nsteps, nwalkers=nwalkers, ntemps=ntemps,
+    log10beta_min=log10beta_min)
+mcmc.setup_initialisation(100, scatter_val=1e-10)
+mcmc.run(subtractions=[F0, F1])
+mcmc.plot_corner(add_prior=True)
+mcmc.print_summary()
diff --git a/examples/weak_signal_follow_up.py b/examples/weak_signal_follow_up.py
deleted file mode 100644
index f4c87966945922f83aba9d9d721465785c28706c..0000000000000000000000000000000000000000
--- a/examples/weak_signal_follow_up.py
+++ /dev/null
@@ -1,69 +0,0 @@
-import pyfstat
-import numpy as np
-import matplotlib.pyplot as plt
-
-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 = 'weak_signal_follow_up_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()
-
-# 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
-log10beta_min = -0.5
-nwalkers = 100
-scatter_val = 1e-10
-nsteps = [100, 100]
-
-mcmc = pyfstat.MCMCFollowUpSearch(
-    label='weak_signal_follow_up', outdir='data',
-    sftfilepattern='data/*'+data_label+'*sft', theta_prior=theta_prior, tref=tref,
-    minStartTime=tstart, maxStartTime=tend, nwalkers=nwalkers, nsteps=nsteps,
-    ntemps=ntemps, log10beta_min=log10beta_min,
-    scatter_val=scatter_val)
-
-fig, axes = plt.subplots(nrows=2, ncols=2)
-mcmc.run(
-    R=10, Nsegs0=100, subtractions=[F0, F1, Alpha, Delta], context='paper',
-    fig=fig, axes=axes, plot_det_stat=False, return_fig=True)
-
-mcmc.plot_corner(add_prior=True)
-mcmc.print_summary()