From bd776d5a72988e5a0c3fd31258cc797ed2ec6b25 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Fri, 18 Mar 2016 14:01:10 +0000 Subject: [PATCH] Added compactness fitness function. --- user_libs/optimisation_taguchi_fitness.py | 153 +++++++++++++++++++--- 1 file changed, 134 insertions(+), 19 deletions(-) diff --git a/user_libs/optimisation_taguchi_fitness.py b/user_libs/optimisation_taguchi_fitness.py index 50f23525..cecaa806 100644 --- a/user_libs/optimisation_taguchi_fitness.py +++ b/user_libs/optimisation_taguchi_fitness.py @@ -5,6 +5,7 @@ # # Please use the attribution at http://dx.doi.org/10.1190/1.3548506 +import sys import h5py import numpy as np np.seterr(divide='ignore') @@ -33,10 +34,10 @@ def minvalue(filename, args): nrx = f.attrs['nrx'] for rx in range(1, nrx + 1): - tmp = f['/rxs/rx' + str(rx) + '/'] - if tmp.attrs['Name'] in args['outputs']: - fieldname = list(tmp.keys())[0] - minvalue = np.amin(tmp[fieldname]) + output = f['/rxs/rx' + str(rx) + '/'] + if output.attrs['Name'] in args['outputs']: + outputname = list(output.keys())[0] + minvalue = np.amin(output[outputname]) return np.abs(minvalue) @@ -56,10 +57,10 @@ def maxvalue(filename, args): nrx = f.attrs['nrx'] for rx in range(1, nrx + 1): - tmp = f['/rxs/rx' + str(rx) + '/'] - if tmp.attrs['Name'] in args['outputs']: - fieldname = list(tmp.keys())[0] - maxvalue = np.amax(tmp[fieldname]) + output = f['/rxs/rx' + str(rx) + '/'] + if output.attrs['Name'] in args['outputs']: + outputname = list(output.keys())[0] + maxvalue = np.amax(output[outputname]) return maxvalue @@ -88,16 +89,16 @@ def xcorr(filename, args): modeltime = np.arange(0, f.attrs['dt'] * f.attrs['Iterations'], f.attrs['dt']) for rx in range(1, nrx + 1): - tmp = f['/rxs/rx' + str(rx) + '/'] - if tmp.attrs['Name'] in args['outputs']: - fieldname = list(tmp.keys())[0] - modelresp = tmp[fieldname] + output = f['/rxs/rx' + str(rx) + '/'] + if output.attrs['Name'] in args['outputs']: + outputname = list(output.keys())[0] + modelresp = output[outputname] # Convert field value (V/m) to voltage - if fieldname == 'Ex': + if outputname == 'Ex': modelresp *= -1 * f.attrs['dx, dy, dz'][0] - elif fieldname == 'Ey': + elif outputname == 'Ey': modelresp *= -1 * f.attrs['dx, dy, dz'][1] - if fieldname == 'Ez': + if outputname == 'Ez': modelresp *= -1 * f.attrs['dx, dy, dz'][2] # Normalise respose from output file @@ -160,10 +161,10 @@ def min_sum_diffs(filename, args): diffdB = 0 outputs = 0 for rx in range(1, nrx + 1): - tmp = f['/rxs/rx' + str(rx) + '/'] - if tmp.attrs['Name'] in args['outputs']: - fieldname = list(tmp.keys())[0] - modelresp = np.array(tmp[fieldname]) + output = f['/rxs/rx' + str(rx) + '/'] + if output.attrs['Name'] in args['outputs']: + outputname = list(output.keys())[0] + modelresp = np.array(output[outputname]) # Calculate sum of differences tmp = 20 * np.log10(np.abs(modelresp - refresp) / np.amax(np.abs(refresp))) tmp = np.abs(np.sum(tmp[-np.isneginf(tmp)])) / len(tmp[-np.isneginf(tmp)]) @@ -173,6 +174,56 @@ def min_sum_diffs(filename, args): return diffdB / outputs +def compactness(filename, args): + """A measure of the compactness of a time domain signal. + + Args: + filename (str): Name of output file + args (dict): 'outputs' key with a list of names (IDs) of outputs (rxs) from input file + + Returns: + compactness (float): Compactness value of signal + """ + + f = h5py.File(filename, 'r') + nrx = f.attrs['nrx'] + dt = f.attrs['dt'] + iterations = f.attrs['Iterations'] + time = np.linspace(0, 1, iterations) + time *= (iterations * dt) + + for rx in range(1, nrx + 1): + output = f['/rxs/rx' + str(rx) + '/'] + if output.attrs['Name'] in args['outputs']: + outputname = list(output.keys())[0] + outputdata = output[outputname] + + # Get absolute maximum value in signal + peak = np.amax([np.amax(outputdata), np.abs(np.amin(outputdata))]) + + # Get peaks and troughs in signal + delta = peak / 50 # Considered a peak/trough if it has the max/min value, and was preceded (to the left) by a value lower by delta + maxtab, mintab = peakdet(outputdata, delta) + peaks = maxtab + mintab + peaks.sort() + + # Remove any peaks smaller than a threshold + thresholdpeak = 1e-3 + peaks = [peak for peak in peaks if np.abs(outputdata[peak]) > thresholdpeak] + + # Percentage of maximum value to measure compactness of signal + durationthreshold = 2 + # Check if there is a peak/trough smaller than threshold + durationthresholdexist = np.where(np.abs(outputdata[peaks]) < (peak * (durationthreshold / 100)))[0] + if durationthresholdexist.size == 0: + compactness = time[peaks[-1]] + else: + time2threshold = time[peaks[durationthresholdexist[0]]] + compactness = time2threshold - time[min(peaks)] + + return (1 / compactness) + + ###################################### # Helper methods for signal analyses # @@ -213,6 +264,70 @@ def zero_crossings(x): npos = ~pos indexzeros = ((pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:])).nonzero()[0] return indexzeros + + +def peakdet(v, delta, x = None): + """Detect peaks and troughs in a vector (adapted from MATLAB script at http://billauer.co.il/peakdet.html). + A point is considered a maximum peak if it has the maximal value, and was preceded (to the left) by a value lower by delta. + + Eli Billauer, 3.4.05 (Explicitly not copyrighted). + This function is released to the public domain; Any use is allowed. + + Args: + v (float): 1D array + delta (float): threshold for determining peaks/troughs + + Returns: + maxtab, mintab (list): Indices of peak/trough locations + """ + + maxtab = [] + mintab = [] + + if x is None: + x = np.arange(len(v)) + + v = np.asarray(v) + + if len(v) != len(x): + sys.exit('Input vectors v and x must have same length') + + if not np.isscalar(delta): + sys.exit('Input argument delta must be a scalar') + + if delta <= 0: + sys.exit('Input argument delta must be positive') + + mn, mx = np.Inf, -np.Inf + mnpos, mxpos = np.NaN, np.NaN + + lookformax = True + + for i in np.arange(len(v)): + this = v[i] + if this > mx: + mx = this + mxpos = x[i] + if this < mn: + mn = this + mnpos = x[i] + + if lookformax: + if this < mx-delta: + if int(mxpos) != 0: + maxtab.append(int(mxpos)) + mn = this + mnpos = x[i] + lookformax = False + else: + if this > mn+delta: + if int(mnpos) != 0: + mintab.append(int(mnpos)) + mx = this + mxpos = x[i] + lookformax = True + + return maxtab, mintab