diff --git a/wetb/prepost/windIO.py b/wetb/prepost/windIO.py
index 170e6e1d8f943ec3fa02b87556b3cb3774b65774..426bbab98a7a7a7f699a162800368752ecda9852 100755
--- a/wetb/prepost/windIO.py
+++ b/wetb/prepost/windIO.py
@@ -803,7 +803,7 @@ class LoadResults(ReadHawc2):
         return slice_, window, zoomtype, time_range
 
     # TODO: general signal method, this is not HAWC2 specific, move out
-    def calc_stats(self, sig, i0=0, i1=None):
+    def calc_stats(self, sig, i0=0, i1=None, ks=False, rho=50.):
 
         stats = {}
         # calculate the statistics values:
@@ -815,6 +815,14 @@ class LoadResults(ReadHawc2):
         stats['absmax'] = np.absolute(sig[i0:i1, :]).max(axis=0)
         stats['rms'] = np.sqrt(np.mean(sig[i0:i1, :]*sig[i0:i1, :], axis=0))
         stats['int'] = integrate.trapz(sig[i0:i1, :], x=sig[i0:i1, 0], axis=0)
+        if ks:
+            stats['max_ks'] = np.zeros(sig.shape[1])
+            stats['min_ks'] = np.zeros(sig.shape[1])
+            for i in range(sig.shape[1]):
+                p = sig[i0:i1, i] / np.abs(stats['max'][i])
+                stats['max_ks'][i] = ksfunc(p, side=1, rho=rho) * np.abs(stats['max'][i])
+                p = sig[i0:i1, i] / np.abs(stats['min'][i])
+                stats['min_ks'][i] = ksfunc(p, side=-1, rho=rho) * np.abs(stats['min'][i])
         return stats
 
     # TODO: general signal method, this is not HAWC2 specific, move out
@@ -1638,6 +1646,28 @@ class Bladed(object):
         df.index.name = lines[0][:-2]
         return df
 
+def ksfunc(p, rho=50., side=1.):
+    """
+    Kreisselmeier and Steinhauser constraint aggregation function
+
+    params
+    --------
+    p: 1-d array
+         vector of constraints
+    beta: float
+        aggregation parameter
+    side: float
+        side=1 computes the max, side=-1. computes the min
+    """
+    p = p.flatten()
+    side = np.float(side)
+    if side == 1:
+        pmax = p.max()
+    elif side == -1:
+        pmax = p.min()
+
+    return pmax + side * np.log(np.sum(np.exp(side * rho * (p - pmax)))) / rho
+
 
 if __name__ == '__main__':