Commit ce971a07 authored by Frederik Zahle's avatar Frederik Zahle
Browse files

added ks function for smooth min/max

parent 0340b1b5
Pipeline #1573 passed with stage
......@@ -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__':
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment