From 5d115632f89045f0462909fa465da449d0d89703 Mon Sep 17 00:00:00 2001
From: "gregory.ashton" <gregory.ashton@ligo.org>
Date: Thu, 2 Feb 2017 14:02:45 +0100
Subject: [PATCH] Adds test to see if samples are railing

Also adds thetas_at_ref data to the Writer object
---
 pyfstat/core.py                |  1 +
 pyfstat/mcmc_based_searches.py | 20 ++++++++++++++++++++
 2 files changed, 21 insertions(+)

diff --git a/pyfstat/core.py b/pyfstat/core.py
index 0fd5517..130039d 100755
--- a/pyfstat/core.py
+++ b/pyfstat/core.py
@@ -116,6 +116,7 @@ class BaseSearchClass(object):
                 post_theta_at_ith_glitch = pre_theta_at_ith_glitch + dt
                 thetas.append(self.shift_coefficients(
                     post_theta_at_ith_glitch, self.tref - tbounds[i+1]))
+        self.thetas_at_tref = thetas
         return thetas
 
     def generate_loudest(self):
diff --git a/pyfstat/mcmc_based_searches.py b/pyfstat/mcmc_based_searches.py
index 8b3ddab..a32046a 100644
--- a/pyfstat/mcmc_based_searches.py
+++ b/pyfstat/mcmc_based_searches.py
@@ -851,6 +851,26 @@ class MCMCSearch(BaseSearchClass):
             d[k+'_std'] = np.std(s)
         return d
 
+    def check_if_samples_are_railing(self, threshold=0.01):
+        return_flag = False
+        for s, k in zip(self.samples.T, self.theta_keys):
+            prior = self.theta_prior[k]
+            if prior['type'] == 'unif':
+                prior_range = prior['upper'] - prior['lower']
+                edges = []
+                fracs = []
+                for l in ['lower', 'upper']:
+                    bools = np.abs(s - prior[l])/prior_range < threshold
+                    if np.any(bools):
+                        edges.append(l)
+                        fracs.append(str(100*float(np.sum(bools))/len(bools)))
+                if len(edges) > 0:
+                    logging.warning(
+                        '{}% of the {} posterior is railing on the {} edges'
+                        .format('% & '.join(fracs), k, ' & '.join(edges)))
+                    return_flag = True
+        return return_flag
+
     def write_par(self, method='med'):
         """ Writes a .par of the best-fit params with an estimated std """
         logging.info('Writing {}/{}.par using the {} method'.format(
-- 
GitLab