Added flag to gmx_stats_get_point

diff --git a/include/gmx_statistics.h b/include/gmx_statistics.h
index 5d28891b98faa401154d8f801dd8a0eb4480e17c..afc124894a4bac075ab8f17ec85e3f7ac16bdea9 100644 (file)
--- a/include/gmx_statistics.h
+++ b/include/gmx_statistics.h
@@ -60,6 +60,8 @@ int gmx_stats_done(gmx_stats_t stats);
    sigma. Level needs to be larger than one obviously. */
 int gmx_stats_remove_outliers(gmx_stats_t stats,double level);
 
+
+
 int gmx_stats_add_point(gmx_stats_t stats,double x,double y,
                               double dx,double dy);
 
@@ -72,9 +74,11 @@ int gmx_stats_add_points(gmx_stats_t stats,int n,real *x,real *y,
    more points, and returns estatsNOPOINTS when the last point has
    been returned. Should be used in a while loop. Variables for either
    pointer may be NULL, in which case the routine can be used as an
-   expensive point counter. */
+   expensive point counter. 
+   If level > 0 then the outliers outside level*sigma are reported
+   only. */
 int gmx_stats_get_point(gmx_stats_t stats,real *x,real *y,
-                              real *dx,real *dy);
+                       real *dx,real *dy,real level);
 
 /* Fit the data to y = ax + b, possibly weighted, if uncertainties
    have been input. Returns slope in *a and intercept in b, *return

diff --git a/src/gmxlib/statistics/gmx_statistics.c b/src/gmxlib/statistics/gmx_statistics.c
index eb723bed61e13c1f96c3338abbd16fc9f5e2b7b1..74367c76dc7d24386b04eec20449d0c133c85c30 100644 (file)
--- a/src/gmxlib/statistics/gmx_statistics.c
+++ b/src/gmxlib/statistics/gmx_statistics.c
@@ -129,20 +129,34 @@ int gmx_stats_add_point(gmx_stats_t gstats,double x,double y,
 }
 
 int gmx_stats_get_point(gmx_stats_t gstats,real *x,real *y,
-                        real *dx,real *dy)
+                        real *dx,real *dy,real level)
 {
     gmx_stats *stats = (gmx_stats *) gstats;
-  
-    if (stats->np_c < stats->np) 
+    int  ok,outlier;
+    real rmsd;
+    
+    if ((ok = gmx_stats_get_rmsd(gstats,&rmsd)) != estatsOK)
+    {
+        return ok;
+    }
+    outlier = 0;
+    while ((outlier == 0) && (stats->np_c < stats->np))
     {
-        if (NULL != x)  *x  = stats->x[stats->np_c];
-        if (NULL != y)  *y  = stats->y[stats->np_c];
-        if (NULL != dx) *dx = stats->dx[stats->np_c];
-        if (NULL != dy) *dy = stats->dy[stats->np_c];
+        r = fabs(stats->x[stats->np_c] - stats->y[stats->np_c]);
+        outlier = (r > rmsd*level);
+        if (outlier)
+        {
+            if (NULL != x)  *x  = stats->x[stats->np_c];
+            if (NULL != y)  *y  = stats->y[stats->np_c];
+            if (NULL != dx) *dx = stats->dx[stats->np_c];
+            if (NULL != dy) *dy = stats->dy[stats->np_c];
+        }
         stats->np_c++;
-    
-        return estatsOK;
+        
+        if (outlier)    
+            return estatsOK;
     }
+    
     stats->np_c = 0;
   
     return estatsNO_POINTS;