Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • toolbox/WindEnergyToolbox
  • tlbl/WindEnergyToolbox
  • cpav/WindEnergyToolbox
  • frza/WindEnergyToolbox
  • borg/WindEnergyToolbox
  • mmpe/WindEnergyToolbox
  • ozgo/WindEnergyToolbox
  • dave/WindEnergyToolbox
  • mmir/WindEnergyToolbox
  • wluo/WindEnergyToolbox
  • welad/WindEnergyToolbox
  • chpav/WindEnergyToolbox
  • rink/WindEnergyToolbox
  • shfe/WindEnergyToolbox
  • shfe1/WindEnergyToolbox
  • acdi/WindEnergyToolbox
  • angl/WindEnergyToolbox
  • wliang/WindEnergyToolbox
  • mimc/WindEnergyToolbox
  • wtlib/WindEnergyToolbox
  • cmos/WindEnergyToolbox
  • fabpi/WindEnergyToolbox
22 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1770 additions and 482 deletions
wetb/signal/fix/_rotor_position.py
View file @ 62f0f2e4
......@@ -7,9 +7,10 @@ import numpy as np
from wetb.signal.fit._spline_fit import spline_fit
from wetb.signal.filters._differentiation import differentiation
def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt=None):
"""Rotor position fitted with spline
Parameters
----------
rotor_position : array_like
......@@ -25,61 +26,63 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt
Note that a significant speed up is achievable by specifying the parameter
plt : PyPlot or None
If PyPlot a visual interpretation is plotted
Returns
-------
y : nd_array
Fitted rotor position
"""
from wetb.signal.subset_mean import revolution_trigger
t = np.arange(len(rotor_position))
indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=4)
rp = rotor_position[:].copy()
for i in indexes:
rp[i:] += 360
if fix_dt is None:
fix_dt = find_fix_dt(rotor_position, sample_frq, rotor_speed)
N = int(np.round(fix_dt*sample_frq))
N2 = N//2
N = int(np.round(fix_dt * sample_frq))
N2 = N // 2
if plt:
a = (rp.max()-rp.min())/t.max()
plt.plot(t/sample_frq,rp-t*a,label='Continus rotor position (detrended)')
a = (rp.max() - rp.min()) / t.max()
plt.plot(t / sample_frq, rp - t * a, label='Continuous rotor position (detrended)')
points = []
for j, i in enumerate(range(0,len(rp), N)):
#indexes for subsets for overlapping a and b polynomials
i1 = max(0,i-N2)
i2 = min(len(rp),i+N2)
#fit a polynomial
if i1<len(rp):
poly_coef = np.polyfit(t[i1:i2]-t[i1], rp[i1:i2], 1)
points.append((t[i],np.poly1d(poly_coef)(t[i]-t[i1])))
for j, i in enumerate(range(0, len(rp), N)):
# indexes for subsets for overlapping a and b polynomials
i1 = max(0, i - N2)
i2 = min(len(rp), i + N2)
# fit a polynomial
if i1 < len(rp):
poly_coef = np.polyfit(t[i1:i2] - t[i1], rp[i1:i2], 1)
points.append((t[i], np.poly1d(poly_coef)(t[i] - t[i1])))
if plt:
plt.plot(t[i1:i2]/sample_frq, np.poly1d(poly_coef)(t[i1:i2]-t[i1])- t[i1:i2]*a, 'mc'[j%2], label=('', "Line fit for points (detrended)")[j<2])
x,y = np.array(points).T
plt.plot(t[i1:i2] / sample_frq, np.poly1d(poly_coef)(t[i1:i2] - t[i1]) - t[i1:i2]
* a, 'mc'[j % 2], label=('', "Line fit for points (detrended)")[j < 2])
x, y = np.array(points).T
if plt:
plt.plot(x/sample_frq,y-x*a,'.', label='Fit points (detrended)')
plt.plot(t/sample_frq, spline_fit(x,y)(t)-t*a, label='Spline (detrended)')
plt.plot(x / sample_frq, y - x * a, '.', label='Fit points (detrended)')
plt.plot(t / sample_frq, spline_fit(x, y)(t) - t * a, label='Spline (detrended)')
plt.legend()
plt.show()
fit = spline_fit(x,y)(t)
fit[differentiation(fit, "left")<0]=np.nan
return fit%360
fit = spline_fit(x, y)(t)
fit[differentiation(fit, "left") < 0] = np.nan
return fit % 360
def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
"""Find the optimal fix_dt parameter for fix_rotor_position (function above).
Optimal is defined as the value that minimizes the sum of squared differences
between differentiated rotor position and rotor speed
Parameters
----------
rotor_position : array_like
......@@ -90,35 +93,37 @@ def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
Rotor speed [RPM]
plt : pyplot or None
If pyplot, a visual interpretation is plotted
Returns
-------
y : int
Optimal value for the fix_dt parameter for fix_rotor_position
"""
from wetb.signal.filters import differentiation
def err(i):
drp = differentiation(fix_rotor_position(rotor_position, sample_frq, rotor_speed, i))
rpm_pos = drp%180 / 360 * sample_frq * 60
rpm_pos = drp % 180 / 360 * sample_frq * 60
return np.sum((rpm_pos - rotor_speed)**2)
best = 27
for step in [9,3,1]:
x_lst = np.arange(-2,3)*step + best
for step in [9, 3, 1]:
x_lst = np.arange(-2, 3) * step + best
res = [err(x) for x in x_lst]
if plt is not None:
plt.plot(x_lst, res,'.-')
plt.plot(x_lst, res, '.-')
best = x_lst[np.argmin(res)]
if plt is not None:
plt.show()
return best
return best
def check_rotor_position(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None):
"""Rotor position fitted with spline
Parameters
----------
rotor_position : array_like
......@@ -134,50 +139,53 @@ def check_rotor_position(rotor_position, sample_frq, rotor_speed, max_rev_diff=1
Note that a significant speed up is achievable by specifying the parameter
plt : PyPlot or None
If PyPlot a visual interpretation is plotted
Returns
-------
y : nd_array
Fitted rotor position
"""
from wetb.signal.subset_mean import revolution_trigger
t = np.arange(len(rotor_position))/sample_frq
t = np.arange(len(rotor_position)) / sample_frq
indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=max_rev_diff)
rotor_position = rotor_position[:]%360
rotor_position = rotor_position[:] % 360
rp_fit = fix_rotor_position(rotor_position, sample_frq, rotor_speed, None)
if plt:
#a = (rp.max()-rp.min())/t.max()
#plt.plot(t/sample_frq,rp-t*a,label='Continus rotor position (detrended)')
f, axarr = plt.subplots(2)
print (rp_fit)
axarr[0].plot(t, rotor_position)
axarr[0].plot(indexes/sample_frq, rotor_position[indexes],'.')
axarr[0].plot(t, rp_fit)
axarr[1].plot(t, rotor_speed)
axarr[1].plot(indexes/sample_frq,60/( differentiation(indexes)/sample_frq))
print (t[170:200])
print(rp_fit)
axarr[0].plot(t, rotor_position, label='rotor position')
axarr[0].plot(indexes / sample_frq, rotor_position[indexes], '.', label='Revolution trigger')
axarr[0].plot(t, rp_fit, label='Fitted rotor position')
axarr[0].legend()
axarr[1].plot(t, rotor_speed, label='Rotor speed')
axarr[1].plot(indexes / sample_frq, 60 / (differentiation(indexes) /
sample_frq), label='Rotor speed (revolution average)')
print(t[170:200])
drp = differentiation(rp_fit)
#drp[(drp<0)]= 0
axarr[1].plot(t, drp%180 / 360 * sample_frq * 60, label='fix')
plt.legend()
i1,i2 = indexes[0], indexes[-1]
print ("Rev from rotor position", np.sum(np.diff(rotor_position[i1:i2])%180)/360)
print ("Rev from rotor speed", np.mean(rotor_speed[i1:i2])*(i2-i1)/60/sample_frq)
print ("Rev from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360)
print ("Mean RPM from rotor speed", np.mean(rotor_speed))
print ("Mean RPM from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360 / ((i2-i1)/60/sample_frq))
spr1 = (np.diff(indexes)/sample_frq)
axarr[1].plot(t, drp % 180 / 360 * sample_frq * 60, label='Fix(calc from fitted rotor position)')
axarr[1].legend()
i1, i2 = indexes[0], indexes[-1]
print("Rev from rotor position", np.sum(np.diff(rotor_position[i1:i2]) % 180) / 360)
print("Rev from rotor speed", np.mean(rotor_speed[i1:i2]) * (i2 - i1) / 60 / sample_frq)
print("Rev from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2]) % 180) / 360)
print("Mean RPM from rotor speed", np.mean(rotor_speed))
print("Mean RPM from fitted rotor position", np.sum(
np.diff(rp_fit[i1:i2]) % 180) / 360 / ((i2 - i1) / 60 / sample_frq))
spr1 = (np.diff(indexes) / sample_frq)
#rs1 = 60/( np.diff(indexes)/sample_frq)
spr2 = np.array([60/rotor_speed[i1:i2].mean() for i1,i2 in zip(indexes[:-1], indexes[1:])])
err = spr1-spr2
print (err.max())
\ No newline at end of file
spr2 = np.array([60 / rotor_speed[i1:i2].mean() for i1, i2 in zip(indexes[:-1], indexes[1:])])
err = spr1 - spr2
print(err.max())
wetb/signal/interpolation.py
View file @ 62f0f2e4
......@@ -43,8 +43,8 @@ def interpolate(x, xp, yp, max_xp_step=None, max_dydxp=None, cyclic_range=None,
[359,nan,0,175,350]
"""
xp = np.array(xp, dtype=np.float)
yp = np.array(yp, dtype=np.float)
xp = np.array(xp, dtype=float)
yp = np.array(yp, dtype=float)
assert xp.shape[0] == yp.shape[0], "xp and yp must have same length (%d!=%d)" % (xp.shape[0], yp.shape[0])
non_nan = ~(np.isnan(xp) & np.isnan(yp))
yp = yp[non_nan]
......
wetb/signal/subset_mean.py
View file @ 62f0f2e4
......@@ -5,14 +5,12 @@ Created on 24/06/2016
'''
import numpy as np
import unittest
from wetb.utils.geometry import rpm2rads
from _collections import deque
from tables.tests.test_index_backcompat import Indexes2_0TestCase
from wetb.signal.filters._differentiation import differentiation
def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1) :
def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1):
"""Calculate the density normalized mean power, taking acceleration of the rotor into account
Parameters
......@@ -46,21 +44,24 @@ def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225,
rs1 = rotor_speed[trigger_indexes[:-1]]
rs2 = rotor_speed[trigger_indexes[1:] - 1]
else:
rs = np.array([rotor_speed[max(i - rotor_speed_mean_samples, 0):i - 1 + rotor_speed_mean_samples].mean() for i in trigger_indexes])
rs = np.array([rotor_speed[max(i - rotor_speed_mean_samples, 0):i - 1 +
rotor_speed_mean_samples].mean() for i in trigger_indexes])
rs1 = rs[:-1]
rs2 = rs[1:]
power = np.array([np.nanmean(power[i1:i2], 0) for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
power = np.array([np.nanmean(power[i1:i2], 0)
for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
if isinstance(air_density, (int, float)):
if air_density != 1.225:
power = power / air_density * 1.225
else:
air_density = np.array([np.nanmean(air_density[i1:i2], 0) for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
air_density = np.array([np.nanmean(air_density[i1:i2], 0)
for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
power = power / air_density * 1.225
return power + 1 / 2 * I * (rs2 ** 2 - rs1 ** 2) / (time[trigger_indexes[1:] - 1] - time[trigger_indexes[:-1]])
def power_mean_func_kW(I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1) :
def power_mean_func_kW(I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1):
"""Return a power mean function [kW] used to Calculate the density normalized mean power, taking acceleration of the rotor into account
Parameters
......@@ -87,14 +88,14 @@ def power_mean_func_kW(I, rotor_speed, time, air_density=1.225, rotor_speed_mean
wsp_mean, power_mean = subset_mean([wsp, power],trigger_indexes,mean_func={1:turbine_power_mean})
"""
def mean_power(power, trigger_indexes):
return power_mean(power * 1000, trigger_indexes, I, rotor_speed, time , air_density, rotor_speed_mean_samples) / 1000
return power_mean(power * 1000, trigger_indexes, I, rotor_speed, time, air_density, rotor_speed_mean_samples) / 1000
return mean_power
def subset_mean(data, trigger_indexes, mean_func={}):
if isinstance(data, list):
data = np.array(data).T
if len(data.shape)==1:
if len(data.shape) == 1:
no_sensors = 1
else:
no_sensors = data.shape[1]
......@@ -103,7 +104,8 @@ def subset_mean(data, trigger_indexes, mean_func={}):
steps = np.diff(triggers[:, 0])
lengths = np.diff(triggers)[:, 0]
if np.all(steps == steps[0]) and np.all(lengths == lengths[0]):
subset_mean = np.mean(np.r_[data.reshape(data.shape[0],no_sensors),np.empty((steps[0],no_sensors))+np.nan][triggers[0][0]:triggers.shape[0] * steps[0] + triggers[0][0]].reshape(triggers.shape[0], steps[0], no_sensors)[:, :lengths[0]], 1)
subset_mean = np.mean(np.r_[data.reshape(data.shape[0], no_sensors), np.empty((steps[0], no_sensors)) + np.nan][triggers[0]
[0]:triggers.shape[0] * steps[0] + triggers[0][0]].reshape(triggers.shape[0], steps[0], no_sensors)[:, :lengths[0]], 1)
else:
subset_mean = np.array([np.mean(data[i1:i2], 0) for i1, i2 in trigger_indexes])
for index, func in mean_func.items():
......@@ -113,15 +115,17 @@ def subset_mean(data, trigger_indexes, mean_func={}):
steps = np.diff(trigger_indexes)
if np.all(steps == steps[0]):
#equal distance
subset_mean = np.mean(data[trigger_indexes[0]:trigger_indexes[-1]].reshape([ len(trigger_indexes) - 1, steps[0], data.shape[1]]), 1)
# equal distance
subset_mean = np.mean(data[trigger_indexes[0]:trigger_indexes[-1]
].reshape([len(trigger_indexes) - 1, steps[0], data.shape[1]]), 1)
else:
subset_mean = np.array([np.mean(data[i1:i2], 0) for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
subset_mean = np.array([np.mean(data[i1:i2], 0) for i1, i2 in zip(
trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
for index, func in mean_func.items():
att = data[:, index]
subset_mean[:, index] = func(att, trigger_indexes)
if len(data.shape)==1 and len(subset_mean.shape)==2:
return subset_mean[:,0]
if len(data.shape) == 1 and len(subset_mean.shape) == 2:
return subset_mean[:, 0]
else:
return subset_mean
......@@ -137,9 +141,10 @@ def cycle_trigger(values, trigger_value=None, step=1, ascending=True, tolerance=
else:
return np.where((values[1:] < trigger_value - tolerance) & (values[:-1] >= trigger_value + tolerance))[0][::step]
def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None):
"""Returns one index per revolution (minimum rotor position)
Parameters
----------
rotor_position : array_like
......@@ -148,84 +153,75 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1,
Sample frequency [Hz]
rotor_speed : array_like
Rotor speed [RPM]
Returns
-------
nd_array : Array of indexes
"""
if isinstance(rotor_speed, (float, int)):
rotor_speed = np.ones_like(rotor_position)*rotor_speed
deg_per_sample = rotor_speed*360/60/sample_frq
thresshold = deg_per_sample.max()*3
drp = (np.diff(rotor_position)+thresshold)%360-thresshold
rotor_speed = np.ones_like(rotor_position) * rotor_speed
deg_per_sample = rotor_speed * 360 / 60 / sample_frq
thresshold = deg_per_sample.max() * 3
drp = (np.diff(rotor_position) + thresshold) % 360 - thresshold
rp = rotor_position
rp = np.array(rotor_position).copy() % 360
# filter degree increase > thresshold
rp[np.r_[False, np.diff(rp) > thresshold]] = 180
upper_indexes = np.where((rp[:-1] > (360 - thresshold)) & (rp[1:] < (360 - thresshold)))[0]
lower_indexes = np.where((rp[:-1] > thresshold) & (rp[1:] < thresshold))[0] + 1
rp = np.array(rotor_position).copy()%360
#filter degree increase > thresshold
rp[np.r_[False, np.diff(rp)>thresshold]] = 180
upper_indexes = np.where((rp[:-1]>(360-thresshold))&(rp[1:]<(360-thresshold)))[0]
lower_indexes = np.where((rp[:-1]>thresshold)&(rp[1:]<thresshold))[0] +1
if plt:
plt.plot(rp)
plt.plot(lower_indexes, rp[lower_indexes],'.')
plt.plot(upper_indexes, rp[upper_indexes],'.')
plt.plot(lower_indexes, rp[lower_indexes], '.')
plt.plot(upper_indexes, rp[upper_indexes], '.')
# Best lower is the first lower after upper
best_lower = lower_indexes[np.searchsorted(lower_indexes, upper_indexes)]
upper2lower = best_lower - upper_indexes
trigger_indexes = best_lower[upper2lower<upper2lower.mean()*2].tolist()
trigger_indexes = best_lower[upper2lower < upper2lower.mean() * 2].tolist()
if len(trigger_indexes) > 1:
if len(trigger_indexes)>1:
rpm_rs = np.array([rev.mean() for rev in np.split(rotor_speed, trigger_indexes)[1:-1]])
rpm_i = 1/np.diff(trigger_indexes)*60*sample_frq
spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
rpm_i = 1 / np.diff(trigger_indexes) * 60 * sample_frq
spr_rs = np.array([rev.mean() for rev in np.split(1 / rotor_speed * 60 * sample_frq, trigger_indexes)[1:-1]])
spr_i = np.diff(trigger_indexes)
while np.any(spr_rs>spr_i*1.9):
i = np.where(spr_rs>spr_i*1.9)[0][0]
if np.abs(spr_i[i-1] + spr_i[i] - spr_rs[i])< np.abs(spr_i[i] + spr_i[i+1] - spr_rs[i]):
while np.any(spr_rs > spr_i * 1.9):
i = np.where(spr_rs > spr_i * 1.9)[0][0]
if np.abs(spr_i[i - 1] + spr_i[i] - spr_rs[i]) < np.abs(spr_i[i] + spr_i[i + 1] - spr_rs[i]):
del trigger_indexes[i]
else:
del trigger_indexes[i+1]
del trigger_indexes[i + 1]
spr_i = np.diff(trigger_indexes)
spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
spr_rs = np.array([rev.mean() for rev in np.split(
1 / rotor_speed * 60 * sample_frq, trigger_indexes)[1:-1]])
# if a revolution is too long add trigger
if np.any(rpm_rs>rpm_i*2.1):
if np.any(rpm_rs > rpm_i * 2.1):
# >1 missing triggers
raise NotImplementedError
trigger_indexes.extend([np.mean(trigger_indexes[i:i+2]) for i in np.where(rpm_rs>rpm_i*1.9)[0]])
trigger_indexes = np.sort(trigger_indexes).astype(np.int)
i1,i2 = trigger_indexes[0], trigger_indexes[-1]
nround_rotor_speed = np.nansum(rotor_speed[i1:i2]/60/sample_frq)
#mod = [v for v in [5,10,30,60,90] if v>thresshold][0]
nround_rotor_position = len(trigger_indexes)-1 #np.nansum(np.diff(rotor_position[i1:i2])%mod)/360
trigger_indexes.extend([np.mean(trigger_indexes[i:i + 2]) for i in np.where(rpm_rs > rpm_i * 1.9)[0]])
trigger_indexes = np.sort(trigger_indexes).astype(int)
i1, i2 = trigger_indexes[0], trigger_indexes[-1]
nround_rotor_speed = np.nansum(rotor_speed[i1:i2] / 60 / sample_frq)
#mod = [v for v in [5,10,30,60,90] if v>thresshold][0]
nround_rotor_position = len(trigger_indexes) - 1 # np.nansum(np.diff(rotor_position[i1:i2])%mod)/360
if max_rev_diff is not None:
diff_pct = abs(nround_rotor_position-nround_rotor_speed)/nround_rotor_position*100
assert diff_pct<max_rev_diff, "No of revolution mismatch: rotor_position (%d), rotor_speed (%.1f), diff %.1f%%"%(nround_rotor_position, nround_rotor_speed, diff_pct)
diff_pct = abs(nround_rotor_position - nround_rotor_speed) / nround_rotor_position * 100
assert diff_pct < max_rev_diff, "No of revolution mismatch: rotor_position (%d), rotor_speed (%.1f), diff %.1f%%" % (
nround_rotor_position, nround_rotor_speed, diff_pct)
return trigger_indexes
def revolution_trigger_old(values, rpm_dt=None, dmin=5, dmax=10, ):
"""Return indexes where values are > max(values)-dmin and decreases more than dmax
If RPM and time step is provided, triggers steps < time of 1rpm is removed
Parameters
---------
values : array_like
......@@ -237,40 +233,39 @@ def revolution_trigger_old(values, rpm_dt=None, dmin=5, dmax=10, ):
Minimum normal position increase between samples
dmax : float or int, optional
Maximum normal position increase between samples
Returns
-------
trigger indexes
[i1,i2,...,in] if rpm_dt is not provided
[(start1,stop1),(start2,stop2),...,(startn, stopn)] if rpm_dt is provided
"""
values = np.array(values)
indexes = np.where((np.diff(values)<-dmax)&(values[:-1]>values.max()-dmax))[0]
indexes = np.where((np.diff(values) < -dmax) & (values[:-1] > values.max() - dmax))[0]
if rpm_dt is None:
return indexes
else:
index_pairs = []
rpm, dt = rpm_dt
d_deg = rpm *360/60*dt
d_deg = rpm * 360 / 60 * dt
cum_d_deg = np.cumsum(d_deg)
lbound, ubound = values.max()-dmax, values.max()+dmax
index_pairs = [(i1,i2) for i1,i2, deg in zip(indexes[:-1], indexes[1:], cum_d_deg[indexes[1:]-1]-cum_d_deg[indexes[:-1]])
if deg > lbound and deg<ubound]
lbound, ubound = values.max() - dmax, values.max() + dmax
index_pairs = [(i1, i2) for i1, i2, deg in zip(indexes[:-1], indexes[1:], cum_d_deg[indexes[1:] - 1] - cum_d_deg[indexes[:-1]])
if deg > lbound and deg < ubound]
return index_pairs
def time_trigger(time, step, start=None, stop=None):
if start is None:
start = time[0]
decimals = int(np.ceil(np.log10(1 / np.nanmin(np.diff(time)))))
time = np.round(time - start, decimals)
steps = np.round(np.diff(time), decimals)
if np.sum(steps == steps[0])/len(time)>.99: #np.all(steps == steps[0]):
if np.sum(steps == steps[0]) / len(time) > .99: # np.all(steps == steps[0]):
# equal time step
time = np.r_[time, time[-1] + max(set(steps), key=list(steps).count)]
if stop is None:
......@@ -286,10 +281,9 @@ def non_nan_index_trigger(sensor, step):
i = 0
nan_indexes = deque(np.where(np.isnan(sensor))[0].tolist() + [len(sensor)])
while i + step <= sensor.shape[0]:
if i+step<=nan_indexes[0]:
trigger.append((i,i+step))
i+=step
if i + step <= nan_indexes[0]:
trigger.append((i, i + step))
i += step
else:
i = nan_indexes.popleft()+1
i = nan_indexes.popleft() + 1
return trigger
wetb/signal/tests/test_first_order.py
View file @ 62f0f2e4
......@@ -14,9 +14,8 @@ from wetb.signal.fit._fourier_fit import F2x, x2F, rx2F
class Test_first_order_filters(unittest.TestCase):
def test_low_pass(self):
a = np.random.randint(0,100,100).astype(np.float)
a = np.random.randint(0, 100, 100000).astype(float)
b = first_order.low_pass(a, 1, 1)
self.assertLess(b.std(), a.std())
if 0:
......@@ -24,30 +23,29 @@ class Test_first_order_filters(unittest.TestCase):
plt.plot(a)
plt.plot(b)
plt.show()
def test_low_pass2(self):
t = np.linspace(0, 1.0, 2001)
xlow = np.sin(2 * np.pi * 5 * t)
xhigh = np.sin(2 * np.pi * 250 * t)
x = xlow + xhigh
b, a = signal.butter(8, 0.125)
y = signal.filtfilt(b, a, x, padlen=150)
self.assertAlmostEqual(np.abs(y - xlow).max(),0,4)
self.assertAlmostEqual(np.abs(y - xlow).max(), 0, 4)
if 0:
import matplotlib.pyplot as plt
plt.plot(x)
plt.plot(y)
plt.show()
# def test_low_pass3(self):
# t = np.linspace(0, 1.0, 2001)
# # xlow = np.sin(2 * np.pi * 5 * t)
# # xhigh = np.sin(2 * np.pi * 250 * t)
# # x = xlow + xhigh
# x = np.sum([np.sin(t*x*2*np.pi) for x in range(1,200)],0)
# x = np.sum([np.sin(t*x*2*np.pi) for x in range(1,200)],0)
# cutoff = .2
# b, a = signal.butter(8,cutoff)
# w, h = signal.freqs(b, a)
......@@ -55,19 +53,19 @@ class Test_first_order_filters(unittest.TestCase):
# F = rx2F(x, max_nfft=len(t))
# tF = np.linspace(0, len(F)/t[-1],len(F))
# Fy = rx2F(y, max_nfft=len(t))
#
#
# if 1:
# import matplotlib.pyplot as plt
# # plt.plot(x)
# # plt.plot(y)
# # plt.show()
#
#
# plt.plot(tF, np.abs(F))
# plt.plot(tF, np.abs(Fy))
# plt.xlim([0,260])
# plt.show()
# print (b,a)
#
#
# plt.plot(w, 20 * np.log10(abs(h)))
# plt.xscale('log')
# plt.title('Butterworth filter frequency response')
......@@ -77,7 +75,7 @@ class Test_first_order_filters(unittest.TestCase):
# plt.grid(which='both', axis='both')
# plt.axvline(cutoff, color='green') # cutoff frequency
# plt.show()
#
#
# def test_low_pass2(self):
# t = np.arange(100000)/10000
# dt = t[1]-t[0]
......@@ -89,11 +87,11 @@ class Test_first_order_filters(unittest.TestCase):
# print (sig.shape)
# F = rx2F(sig, max_nfft=len(t))
# tF = np.linspace(0, len(F)/t[-1],len(F))
#
#
#
#
# sig_lp1 = first_order.low_pass(sig, dt, 30*dt)
# F_lp1 = rx2F(sig_lp1, max_nfft=len(t))
#
#
# bb, ba = signal.butter(1,5, 'low', analog=True)
# print (bb, ba)
# w1, h1 = signal.freqs(bb, ba)
......@@ -104,24 +102,24 @@ class Test_first_order_filters(unittest.TestCase):
# # sig_lp3 = signal.lfilter(bb,ba, sig)
# # F_lp3 = rx2F(sig_lp3, max_nfft=len(t))
# # w2, h2 = signal.freqs(bb, ba)
# # #
# # #
# # self.assertLess(b.std(), a.std())
# if 1:
# import matplotlib.pyplot as plt
# plt.plot(t, sig)
# # plt.plot(t, np.sum([1/x*np.sin(t*x*2*np.pi) for x in [1]],0))
# #
# #
# plt.plot(t, sig_lp1)
# plt.plot(t, sig_lp2)
# #plt.plot(t, sig_lp3)
#
#
# plt.show()
#
#
# # plt.plot(tF, np.abs(F))
# # plt.plot(tF, np.abs(F_lp1))
# # # plt.plot(tF, np.abs(F_lp3))
# # # plt.plot(tF, np.abs(F_lp2))
# # plt.plot([0,1000],[0.708,.708])
# # plt.plot([0,1000],[0.708,.708])
# # plt.xlim([0,205])
# # plt.ylim([0,1])
# # plt.show()
......@@ -135,21 +133,20 @@ class Test_first_order_filters(unittest.TestCase):
# # plt.grid(which='both', axis='both')
# # plt.axvline(100, color='green') # cutoff frequency
# # plt.show()
# # def test_low_pass2(self):
# # F = [0,1,0,0,0,0,0,0,0,0,0,0,0,0,.1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,.1]
# # F +=[0]*(len(F)-1)
# # x = F2x(F)
# # a = np.tile(x, 3)
# # print (x.shape)
# # # a = np.random.randint(0,100,100).astype(np.float)
# # # a = np.random.randint(0,100,100).astype(float)
# # b = first_order.low_pass(a, .1, 1)
# # bb, ba = signal.butter(10,100, 'low', analog=True)
# # #c = signal.lfilter(bb,ba, a)
# # w, h = signal.freqs(bb, ba)
# #
# #
# # # self.assertLess(b.std(), a.std())
# # if 1:
# # import matplotlib.pyplot as plt
......@@ -168,12 +165,11 @@ class Test_first_order_filters(unittest.TestCase):
# # plt.grid(which='both', axis='both')
# # plt.axvline(100, color='green') # cutoff frequency
# # plt.show()
# #
# #
# #
# #
def test_high_pass(self):
a = np.random.randint(0,100,100).astype(np.float)
a = np.random.randint(0, 100, 100).astype(float)
b = first_order.high_pass(a, 1, 1)
self.assertLess(b.mean(), a.mean())
if 0:
......@@ -181,6 +177,8 @@ class Test_first_order_filters(unittest.TestCase):
plt.plot(a)
plt.plot(b)
plt.show()
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
unittest.main()
\ No newline at end of file
unittest.main()
wetb/signal/tests/test_interpolation.py
View file @ 62f0f2e4
......@@ -27,10 +27,10 @@ class TestInterpolation(unittest.TestCase):
def test_interpolate2(self):
x = np.arange(0, 100, 5, dtype=np.float)
xp = np.arange(0, 100, 10, dtype=np.float)
x = np.arange(0, 100, 5, dtype=float)
xp = np.arange(0, 100, 10, dtype=float)
xp = np.r_[xp[:3], xp[5:]]
yp = np.arange(10, dtype=np.float)
yp = np.arange(10, dtype=float)
yp[7:8] = np.nan
yp = np.r_[yp[:3], yp[5:]]
......@@ -39,10 +39,10 @@ class TestInterpolation(unittest.TestCase):
#yp = [0.0, 1.0, 2.0, 5.0, 6.0, nan, 8.0, 9.0]
y = signal.interpolate(x, xp, yp)
np.testing.assert_array_equal(y[~np.isnan(y)], [0., 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 8. , 8.5, 9. ])
np.testing.assert_array_equal(y[~np.isnan(y)], [0., 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6., 8. , 8.5, 9. ])
y = signal.interpolate(x, xp, yp, 10)
np.testing.assert_array_equal(y[~np.isnan(y)], [ 0. , 0.5, 1. , 1.5, 2. , 5. , 5.5, 8. , 8.5, 9. ])
np.testing.assert_array_equal(y[~np.isnan(y)], [ 0. , 0.5, 1. , 1.5, 2. , 5. , 5.5, 6., 8. , 8.5, 9. ])
......
wetb/signal/tests/test_subset_mean.py
View file @ 62f0f2e4
......@@ -39,7 +39,7 @@ class TestSubsetMean(unittest.TestCase):
self.assertEqual(p[1], data[2501:5000, 0].mean())
def test_non_nan_index_trigger(self):
sensor = np.arange(18).astype(np.float)
sensor = np.arange(18).astype(float)
sensor[[5, 11]] = np.nan
triggers = non_nan_index_trigger(sensor, 3)
for i1, i2 in triggers:
......@@ -48,7 +48,7 @@ class TestSubsetMean(unittest.TestCase):
def test_subset_mean_trigger_tuple(self):
sensor = np.arange(18).astype(np.float)
sensor = np.arange(18).astype(float)
triggers = non_nan_index_trigger(sensor, 3)
np.testing.assert_array_equal(subset_mean(sensor, triggers), [ 1., 4., 7., 10., 13., 16])
......@@ -58,7 +58,7 @@ class TestSubsetMean(unittest.TestCase):
np.testing.assert_array_equal(subset_mean(sensor, triggers), [ 2, 5, 8, 11, 14])
#nan in the middle, noneq step and len
sensor = np.arange(18).astype(np.float)
sensor = np.arange(18).astype(float)
sensor[[5, 11]] = np.nan
triggers = non_nan_index_trigger(sensor, 3)
......
wetb/soil/KlinkvortModel.py 0 → 100644
View file @ 62f0f2e4
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 18 13:34:12 2018
@author: shfe
Note:
This script is used for Klinkvort model used for cyclic accumulation calculation
"""
import numpy as np
import matplotlib.pyplot as plt
plt.close('all')
class KlinkvortModel(object):
def __init__(self, xi_b, xi_c, N):
"""
Initialize the input parameters
"""
self.xi_b = xi_b
self.xi_c = xi_c
self.N = N
def alpha(self):
"""
Calculate the parameter alpha
"""
# Calculate the Tb, Tc
if not 0 <= self.xi_b <= 1:
raise ValueError('The maximum load ratio should between 0 to 1')
if not -1 <= self.xi_c <= 1:
raise ValueError('The cyclic load ratio should between -1 to 1')
if self.xi_b <= 0.022:
Tb = 0
else:
Tb = 0.61*self.xi_b - 0.013
Tc = (self.xi_c + 0.63)*(self.xi_c - 1)*(self.xi_c - 1.64)
alpha = Tb * Tc
return Tb, Tc, alpha
def evoluation(self, y1):
"""
Evoluation model for cyclic accumulated rotation
"""
_, _, alpha = self.alpha()
yn = y1 * self.N**alpha
return yn
def equivalentnum(self, y1, yn, alpha):
"""
Equivalent number for cylces
"""
num = (yn/y1)**(1/alpha)
return num
def superposition(xi_b, xi_c, N_cases, y1_cases, s0 = 0):
# Check the array equal number of items
if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
raise ValueError('Size for input should be identical')
print(s0)
yn = np.zeros(np.size(xi_b))
for i, (b, c, N, y1) in enumerate(zip(xi_b, xi_c, N_cases, y1_cases)):
if i==0:
em = KlinkvortModel(b, c, N)
_, _, alpha = em.alpha()
if s0 == 0:
yn[i] = em.evoluation(y1)
else:
N_eq = (s0/y1)**(1/alpha)
if N_eq < 1e8:
# N_eq = (s0/y1)**(1/alpha)
yn[i] = y1*(N+N_eq)**alpha
else:
yn[i] = s0
else:
em = KlinkvortModel(b, c, N)
_, _, alpha = em.alpha()
N_eq = (yn[i-1]/y1)**(1/alpha)
if N_eq < 1e8:
# N_eq = (yn[i-1]/y1)**(1/alpha)
yn[i] = y1*(N_eq + N)**alpha
else:
yn[i] = yn[i-1]
return yn
def alpha_super(xi_b, xi_c, N_cases):
# Check the alpha value
if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
raise ValueError('Size for input should be identical')
alpha_list = np.zeros(np.size(xi_b))
for i, (b, c, N) in enumerate(zip(xi_b, xi_c, N_cases)):
em = KlinkvortModel(b, c, N)
_, _, alpha = em.alpha()
alpha_list[i] = alpha
return alpha_list
\ No newline at end of file
wetb/soil/LeblancModel.py 0 → 100644
View file @ 62f0f2e4
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 25 16:15:25 2018
@author: shfe@dtu.dk
The model proposed by Leblanc
"""
import matplotlib.pyplot as plt
import numpy as np
plt.close('all')
class LeblancModel(object):
def __init__(self, xi_b, xi_c, N):
"""
Initialize the input parameters
"""
self.xi_b = xi_b
self.xi_c = xi_c
self.N = N
self.alpha = 0.31
def coeff(self, rd = 0.38):
"""
Calculate the parameter alpha
"""
# Calculate the Tb, Tc
if not 0 <= self.xi_b <= 1:
raise ValueError('The maximum load ratio should between 0 to 1')
if not -1 <= self.xi_c <= 1:
raise ValueError('The cyclic load ratio should between -1 to 1')
if rd == 0.38:
if self.xi_b <= 0.051:
Tb = 0
else:
Tb = 0.4238*self.xi_b - 0.0217
elif rd == 0.04:
if self.xi_b <= 0.1461:
Tb = 0
else:
Tb = 0.3087*self.xi_b - 0.0451
else:
raise ValueError('The relative density value is not validated yet')
if -1 <= self.xi_c < -0.65:
Tc = 13.71*self.xi_c + 13.71
elif -0.65 <= self.xi_c < 0:
Tc = -5.54*self.xi_c + 1.2
else:
Tc = -1.2*self.xi_c + 1.2
return Tb, Tc
def evoluation(self, y1):
"""
Evoluation model for cyclic accumulated rotation
"""
Tb, Tc = self.coeff()
yn = y1 * Tb * Tc * self.N**0.31
return yn
def equivalentnum(self, y1, yn):
"""
Equivalent number for cylces
"""
Tb, Tc = self.coeff()
num = (yn/(y1*Tb*Tc))**(1/self.alpha)
return num
def superposition(xi_b, xi_c, N_cases, y1_cases, s0 = 0, y0 = 0):
# Check the array equal number of items
if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
raise ValueError('Size for input should be identical')
# print(s0)
dyn = np.zeros(np.size(xi_b))
yn = np.zeros(np.size(xi_b))
y_static = np.zeros(np.size(xi_b))
for i, (b, c, N, y1) in enumerate(zip(xi_b, xi_c, N_cases, y1_cases)):
if i==0:
em = LeblancModel(b, c, N)
Tb, Tc = em.coeff()
if s0 == 0:
dyn[i] = em.evoluation(y1)
else:
N_eq = (s0/(y1*Tb*Tc))**(1/em.alpha)
if N_eq < 1e8:
dyn[i] = (y1*Tb*Tc)*(N+N_eq)**em.alpha
else:
dyn[i] = s0
y_static[i] = max(y1, y0)
else:
em = LeblancModel(b, c, N)
Tb, Tc = em.coeff()
N_eq = (dyn[i-1]/(y1*Tb*Tc))**(1/em.alpha)
if N_eq < 1e8:
dyn[i] = (y1*Tb*Tc)*(N_eq + N)**em.alpha
else:
dyn[i] = dyn[i-1]
y_static[i] = max(y_static[i-1], y1)
yn[i] = dyn[i] + y_static[i]
return dyn, y_static, yn
def coeff_super(xi_b, xi_c, N_cases):
# Check the alpha value
if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
raise ValueError('Size for input should be identical')
Tb_list = np.zeros(np.size(xi_b))
Tc_list = np.zeros(np.size(xi_c))
for i, (b, c, N) in enumerate(zip(xi_b, xi_c, N_cases)):
em = LeblancModel(b, c, N)
Tb, Tc = em.coeff()
Tb_list[i] = Tb
Tc_list[i] = Tc
return Tb_list, Tc_list
def moment_ult(L, D, gamma):
"""
Ultimate static moment
Parameters
----------------
L -- Pile length
D -- Pile diameter
gamma -- Submerged unit weight
Output
----------------
m_ult -- Ultimate static moment
"""
m_ult = L**3*D*gamma*0.6
return m_ult
def rot_moment(m, L, D, gamma, pa = 101, rd = 0.04):
"""
Ultimate static moment
Parameters
----------------
L -- Pile length
D -- Pile diameter
gamma -- Submerged unit weight
m -- Moment
Output
----------------
theta - Rotational angle in radians
"""
if rd == 0.04:
m_ult = L**3*D*gamma*0.6/1e3 # unit: MN
theta_norm = 0.038 * (m/m_ult)**2.33
theta = theta_norm*np.sqrt((L*gamma)/pa)
else:
m_ult = L**3*D*gamma*1.2/1e3 # unit: MN
theta_norm = 0.042 * (m/m_ult)**1.92
theta = theta_norm*np.sqrt((L*gamma)/pa)
return theta
wetb/soil/__init__.py 0 → 100644
View file @ 62f0f2e4
wetb/soil/hyperbolic_model.py 0 → 100644
View file @ 62f0f2e4
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 22 17:06:06 2018
@author: shfe
Hyperbolic model used for soil load displacement behaviour
"""
import numpy as np
def hyperbolic_load(y, pu, a, b):
"""
Static hyperbolic model developed by Lee
Parameters
----------------
y -- displacement
pu -- ultimate load
a, b -- coefficients
Output
----------------
p -- soil resistence pressure
"""
p = y/(a + b*y) * pu
return p
def hyperbolic_disp(p, pu, a, b):
"""
Static hyperbolic model developed by Lee
Parameters
----------------
p -- load
pu -- ultimate load
a, b -- coefficients
Output
----------------
y -- displacement
"""
y = p * a/(pu - p*b)
return y
\ No newline at end of file
wetb/soil/py_model.py 0 → 100644
View file @ 62f0f2e4
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 25 16:15:25 2018
@author: shfe@dtu.dk
Static p-y soil model and cyclic degradation model
"""
import numpy as np
def py_static(y, pu, yc):
"""
Static p-y model developed by Matlock
Parameters
----------------
y -- displacement
pu -- ultimate stress
yc -- ultimate displacement
Output
----------------
p -- soil resistence pressure
"""
if y/yc <= 8:
p = 0.5*pu*(y/yc)**(1/3)
else:
p = pu
return p
def strain_static(p, pu, yc):
"""
Static strain using p-y model developed by Matlock
Parameters
----------------
p -- stress
pu -- ultimate stress
yc -- ultimate displacement
Output
----------------
y -- soil strain
"""
if p <= pu:
y = (p/(0.5*pu))**3*yc
else:
raise ValueError("The applied stress is LARGER than ultimate strength")
return y
def load_unload_py(p, pu, yc):
"""
load/unload module using p-y model developed by Matlock
Parameters
----------------
p -- stress
pu -- ultimate stress
yc -- ultimate displacement
Output
----------------
y -- soil strain
y_cum -- cumulative soil displacement
"""
y = strain_static(p, pu, yc)
if p <= 0.5*pu: # elastic
y_cum = 0
else:
y_cum = y - p/(0.5*pu/yc)
return y, y_cum
def degra_coeff(N, y1, b):
"""
Degradation factor for p-y model
Parameters
----------------
N -- number of cycle
y1 -- displacement predicted by static p-y curve
b -- pile diameter
Output
----------------
lambda -- Degradation factor
"""
# check
lambda_N = np.abs(y1/(0.2*b)*np.log(N))
try:
lambda_N <= 1
except:
print("soil must be degraded!!!")
return lambda_N
def py_cyclic(y, pu, yc, N, y1, b):
"""
Degradation p-y model with cyclic load
Parameters
----------------
y -- displacement
pu -- ultimate stress
yc -- ultimate displacement
N -- number of cycle
y1 -- displacement predicted by static p-y curve
b -- pile diameter
Output
----------------
p -- degraded soil resistence pressure
"""
lambda_N = degra_coeff(N, y1, b)
p = py_static(y, pu, yc) * (1 - lambda_N)
return p
def py_parcel(load_parcel, cycle_parcel, pu, yc, b):
"""
Degradation p-y model after a parcel of cyclic load
Parameters
----------------
load_parcel -- load parcel
cycle_parcel -- number of cycle of load parcel
yc -- ultimate displacement
b -- pile diameter
Output
----------------
p -- degraded soil resistence pressure
"""
pu_N = np.zeros(np.size(load_parcel)+1)
pu_N[0] = pu
lambda_N = np.zeros(np.size(load_parcel))
y_N = np.zeros(np.size(load_parcel))
for i in np.arange(np.size(load_parcel)):
load_i = load_parcel[i]
cycle_i = cycle_parcel[i]
if np.isnan(load_i):
y_i = 0
else:
y_i = strain_static(load_i, pu_N[i], yc)
y_N[i] = y_i
lambda_i = degra_coeff(cycle_i, y_i, b)
lambda_N[i] = lambda_i
pu_N[i+1] = pu_N[i] * (1-lambda_i)
return y_N, lambda_N, pu_N
wetb/standard_models/power_curve.py 0 → 100644
View file @ 62f0f2e4
import numpy as np
def standard_power_curve(power_norm, diameter, turbulence_intensity=.1, shear=(0, 100),
rho=1.225, max_cp=.49,
gear_loss_const=.01, gear_loss_var=.014, generator_loss=0.03, converter_loss=.03):
"""Generate standard power curve
The method is extracted from an Excel spreadsheet made by Kenneth Thomsen, DTU Windenergy and
ported into python by Mads M. Pedersen, DTU Windenergy.
Parameters
----------
power_norm : int or float
Nominal power [kW]
diameter : int or float
Rotor diameter [m]
turbulence_intensity : float
Turbulence intensity [%]
shear : (power shear coefficient, hub height)
Power shear arguments\n
- Power shear coeficient, alpha\n
- Hub height [m]
rho : float optional
Density of air [kg/m^3], defualt is 1.225
max_cp : float
Maximum power coefficient
gear_loss_const : float
Constant gear loss [%]
gear_loss_var : float
Variable gear loss [%]
generator_loss : float
Generator loss [%]
converter_loss : float
Examples
--------
wsp, power = standard_power_curve(10000, 178.3)
plot(wsp, power)
show()
"""
area = (diameter / 2) ** 2 * np.pi
wsp_lst = np.arange(0.5, 25, .5)
sigma_lst = wsp_lst * turbulence_intensity
def norm_dist(x, my, sigma):
if turbulence_intensity > 0:
return 1 / np.sqrt(2 * np.pi * sigma ** 2) * np.exp(-(x - my) ** 2 / (2 * sigma ** 2))
else:
return x == my
p_aero = .5 * rho * area * wsp_lst ** 3 * max_cp / 1000
# calc power - gear, generator and conv loss
# gear_loss = gear_loss_const * power_norm + gear_loss_var * p_aero
# p_gear = p_aero - gear_loss
# p_gear[p_gear < 0] = 0
# gen_loss = generator_loss * p_gear
# p_gen = p_gear - gen_loss
# converter_loss = converter_loss * p_gen
# p_raw = p_gen - converter_loss
# p_raw[p_raw > power_norm] = power_norm
p_raw = p_aero - power_loss(p_aero, power_norm, gear_loss_const, gear_loss_var, generator_loss, converter_loss)
powers = []
shear_weighted_wsp = []
alpha, hub_height = shear
r = diameter / 2
z = np.linspace(hub_height - r, hub_height + r, 100, endpoint=True)
shear_factors = (z / hub_height) ** alpha
rotor_width = 2 * np.sqrt(r ** 2 - (hub_height - z) ** 2)
for wsp, sigma in zip(wsp_lst, sigma_lst):
shear_weighted_wsp.append(wsp * (np.trapz(shear_factors ** 3 * rotor_width, z) / (area)) ** (1 / 3))
ndist = norm_dist(wsp_lst, wsp, sigma)
powers.append((ndist * p_raw).sum() / ndist.sum())
return wsp_lst, lambda wsp : np.interp(wsp, wsp_lst, powers)
return wsp_lst, np.interp(shear_weighted_wsp, wsp_lst, powers)
def power_loss(power_aero, power_norm, gear_loss_const=.01, gear_loss_var=.014, generator_loss=0.03, converter_loss=.03):
gear_loss = gear_loss_const * power_norm + gear_loss_var * power_aero
p_gear = power_aero - gear_loss
p_gear[p_gear < 0] = 0
gen_loss = generator_loss * p_gear
p_gen = p_gear - gen_loss
converter_loss = converter_loss * p_gen
p_electric = p_gen - converter_loss
p_electric[p_electric > power_norm] = power_norm
return power_aero - p_electric
if __name__ == '__main__':
from matplotlib.pyplot import plot, show
plot(*standard_power_curve(10000, 178.3, 0., (0, 119), gear_loss_const=.0, gear_loss_var=.0, generator_loss=0.0, converter_loss=.0))
plot(*standard_power_curve(10000, 178.3, 0.03, (0, 119), gear_loss_const=.0, gear_loss_var=.0, generator_loss=0.0, converter_loss=.0))
show()
wetb/utils/caching.py
View file @ 62f0f2e4
......@@ -3,15 +3,9 @@ Created on 07/02/2014
@author: MMPE
'''
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
import sys
from collections import OrderedDict
import os
standard_library.install_aliases()
import inspect
import numpy as np
def set_cache_property(obj, name, get_func, set_func=None):
......
wetb/utils/cluster_tools/cluster_resource.py
View file @ 62f0f2e4
......@@ -34,15 +34,16 @@ def unix_path(path, cwd=None, fail_on_missing=False):
if r:
path = r[0]
elif fail_on_missing:
raise FileExistsError("File or folder matching '%s' not found"%path)
raise FileExistsError("File or folder matching '%s' not found" % path)
if cwd:
path = os.path.relpath(path, cwd)
if os.path.isdir(path):
path+="/"
return path.replace("\\","/")
path += "/"
return path.replace("\\", "/")
def unix_path_old(filename):
filename = os.path.realpath(filename.replace("\\", "/")).replace("\\", "/")
filename = os.path.realpath(str(filename)).replace("\\", "/")
ufn, rest = os.path.splitdrive(filename)
ufn += "/"
for f in rest[1:].split("/"):
......@@ -51,22 +52,20 @@ def unix_path_old(filename):
f_lst = [f_ for f_ in f_lst if f_ == f]
elif len(f_lst) == 0:
raise IOError("'%s' not found in '%s'" % (f, ufn))
else: # one match found
else: # one match found
ufn = os.path.join(ufn, f_lst[0])
return ufn.replace("\\", "/")
class Resource(object):
def __init__(self, min_cpu, min_free):
self.min_cpu = min_cpu
self.min_free = min_free
self.cpu_free=0
self.cpu_free = 0
self.acquired = 0
self.no_cpu="?"
self.used_by_user=0
self.no_cpu = "?"
self.used_by_user = 0
self.resource_lock = threading.Lock()
def ok2submit(self):
......@@ -75,7 +74,7 @@ class Resource(object):
#print ("ok2submit")
total, free, user = self.check_resources()
user = max(user, self.acquired)
user_max = max((self.min_cpu-user), self.cpu_free - self.min_free)
user_max = max((self.min_cpu - user), self.cpu_free - self.min_free)
#print ("ok2submit", total, free, user, user_max)
return user_max
except:
......@@ -95,7 +94,6 @@ class Resource(object):
with self.resource_lock:
self.acquired -= 1
def update_resource_status(self):
try:
self.no_cpu, self.cpu_free, self.used_by_user = self.check_resources()
......@@ -107,28 +105,27 @@ class SSHPBSClusterResource(Resource):
finished = []
loglines = {}
is_executing = []
def __init__(self, sshclient, min_cpu, min_free):
Resource.__init__(self, min_cpu, min_free)
self.ssh = sshclient
self.resource_lock = threading.Lock()
def glob(self, filepattern, cwd="", recursive=False):
return self.ssh.glob(filepattern, cwd, recursive)
@property
def host(self):
return self.ssh.host
@property
def username(self):
return self.ssh.username
def new_ssh_connection(self):
from wetb.utils.cluster_tools.ssh_client import SSHClient
return SSHClient(self.host, self.ssh.username, self.ssh.password, self.ssh.port)
#return self.ssh
# return self.ssh
def check_resources(self):
with self.resource_lock:
......@@ -153,22 +150,20 @@ class SSHPBSClusterResource(Resource):
except Exception as e:
raise EnvironmentError(str(e))
def jobids(self, jobname_prefix):
_, output, _ = self.ssh.execute('qstat -u %s' % self.username)
return [l.split()[0].split(".")[0] for l in output.split("\n")[5:] if l.strip() != "" and l.split()[3].startswith("h2l")]
_, output, _ = self.ssh.execute('qstat -u %s' % self.username)
return [l.split()[0].split(".")[0] for l in output.split("\n")[5:] if l.strip() != "" and l.split()[3].startswith("h2l")]
def stop_pbsjobs(self, jobids):
if not hasattr(jobids, "len"):
jobids = list(jobids)
self.ssh.execute("qdel %s" % (" ".join(jobids)))
def setup_wine(self):
def setup_wine(self):
self.ssh.execute("""rm -f ./config-wine-hawc2.sh &&
wget https://gitlab.windenergy.dtu.dk/toolbox/pbsutils/raw/master/config-wine-hawc2.sh &&
chmod 777 config-wine-hawc2.sh &&
./config-wine-hawc2.sh""")
class LocalResource(Resource):
......@@ -178,8 +173,6 @@ class LocalResource(Resource):
#self.process_name = process_name
self.host = 'Localhost'
def check_resources(self):
import psutil
no_cpu = multiprocessing.cpu_count()
......
wetb/utils/cluster_tools/os_path.py 0 → 100644
View file @ 62f0f2e4
import os
def repl(path):
return path.replace("\\", "/")
def abspath(path):
return repl(os.path.abspath(path))
def relpath(path, start=None):
return repl(os.path.relpath(path, start))
def realpath(path):
return repl(os.path.realpath(path))
def pjoin(*path):
return repl(os.path.join(*path))
def fixcase(path):
path = realpath(str(path)).replace("\\", "/")
p, rest = os.path.splitdrive(path)
p += "/"
for f in rest[1:].split("/"):
f_lst = [f_ for f_ in os.listdir(p) if f_.lower() == f.lower()]
if len(f_lst) > 1:
# use the case sensitive match
f_lst = [f_ for f_ in f_lst if f_ == f]
if len(f_lst) == 0:
raise IOError("'%s' not found in '%s'" % (f, p))
# Use matched folder
p = pjoin(p, f_lst[0])
return p
def normpath(path):
return repl(os.path.normpath(path))
def cluster_path(path):
try:
import win32wnet
drive, folder = os.path.splitdrive(abspath(path))
path = abspath(win32wnet.WNetGetUniversalName(drive, 1) + folder)
except Exception:
path = repl(path)
path = path.replace("//jess.dtu.dk", "/home")
path = path.replace("//mimer.risoe.dk/aiolos", "/mnt/aiolos")
path = path.replace("//mimer.risoe.dk", "/mnt/mimer")
return path
wetb/utils/cluster_tools/pbsfile.py 0 → 100644
View file @ 62f0f2e4
import os
import inspect
from wetb.utils.cluster_tools.os_path import cluster_path, pjoin, normpath
class Template():
def __init__(self, template):
self.template = template
def __call__(self, **kwargs):
s = self.template
found = True
while found:
found = False
for k, v in dict(kwargs).items():
if "[%s]" % k in s:
found = True
s = s.replace("[%s]" % k, str(v))
return s
pbs_template = Template('''### Jobid
#PBS -N [jobname]
### Standard Output
#PBS -o [stdout_filename]
[stderr]
#PBS -W umask=0003
### Maximum wallclock time format HOURS:MINUTES:SECONDS
#PBS -l walltime=[walltime]
#PBS -l nodes=[nodes]:ppn=[ppn]
### Queue name
#PBS -q [queue]
cd "[workdir]"
mkdir -p "stdout"
if [ -z "$PBS_JOBID" ]; then echo "Run using qsub"; exit ; fi
pwd
[commands]
exit
''')
class PBSFile():
_walltime = "00:30:00"
def __init__(self, workdir, jobname, commands, queue='workq', walltime='00:10:00', nodes=1, ppn=1, merge_std=True):
"""Description
Parameters
----------
walltime : int, str
wall time as string ("hh:mm:ss") or second (integer)
"""
self.workdir = workdir
self.jobname = jobname
self.commands = commands
self.queue = queue
self.walltime = walltime
self.nodes = nodes
self.ppn = ppn
self.merge_std = merge_std
self.stdout_filename = normpath(pjoin(workdir, './stdout/%s.out' % jobname))
self.filename = "pbs_in/%s.in" % self.jobname
@property
def walltime(self):
return self._walltime
@walltime.setter
def walltime(self, walltime):
if isinstance(walltime, (float, int)):
from math import ceil
h = walltime // 3600
m = (walltime % 3600) // 60
s = ceil(walltime % 60)
self._walltime = "%02d:%02d:%02d" % (h, m, s)
else:
self._walltime = walltime
def __str__(self):
if self.merge_std:
stderr = "### merge stderr into stdout\n#PBS -j oe"
else:
stderr = '### Standard Error\n#PBS -e "./err/[jobname].err"'
if callable(self.commands):
commands = self.commands()
else:
commands = self.commands
return pbs_template(workdir=cluster_path(self.workdir),
stdout_filename=cluster_path(self.stdout_filename),
stderr=stderr,
jobname=self.jobname,
queue=self.queue,
walltime=self.walltime,
nodes=self.nodes,
ppn=self.ppn,
commands=commands)
def save(self, modelpath=None, filename=None):
modelpath = modelpath or self.workdir
self.filename = filename or self.filename
filename = os.path.join(modelpath, self.filename)
os.makedirs(os.path.dirname(filename), exist_ok=True)
with open(filename, 'w', newline='\n') as fid:
fid.write(str(self))
os.chmod(filename, 0o774)
multirunner_template = Template("""echo "[make_dict]
" | python
for node in `cat $PBS_NODEFILE | sort | uniq`
do
ssh -T $node << EOF &
cd "[workdir]"
python -c "[start_jobs]
"
EOF
done
wait
""")
class PBSMultiRunner(PBSFile):
def __init__(self, workdir, queue='workq', walltime='01:00:00', nodes=1, ppn=1, merge_std=True, pbsfiles=None,
jobname='pbs_multirunner'):
if pbsfiles:
def fmt(pbsfile):
if isinstance(pbsfile, PBSFile):
return pbsfile.filename
return pbsfile
self.PBS_FILES = "['" + "',\n '".join([fmt(pbsfile) for pbsfile in pbsfiles]) + "']"
else:
self.PBS_FILES = """[os.path.join(root, f) for root, folders, f_lst in os.walk('.') for f in f_lst if f.endswith('.in')]"""
commands = multirunner_template(make_dict=self.get_src(self.make_dict),
start_jobs=self.get_src(self.start_jobs),
workdir=cluster_path(workdir)).replace("self.ppn", str(ppn))
PBSFile.__init__(self, workdir, jobname, commands, queue, walltime=walltime,
nodes=nodes, ppn=ppn, merge_std=merge_std)
self.filename = "%s.%s" % (self.jobname, ("lst", "all")[pbsfiles is None])
self.pbsfiles = pbsfiles
def make_dict(self):
import os
import glob
import numpy as np
import re
# find available nodes
with open(os.environ['PBS_NODEFILE']) as fid:
nodes = set([f.strip() for f in fid.readlines() if f.strip() != ''])
pbs_files = eval(self.PBS_FILES)
# Make a list of [(pbs_in_filename, stdout_filename, walltime),...]
pat = re.compile(r'[\s\S]*#\s*PBS\s+-o\s+(.*)[\s\S]*(\d\d:\d\d:\d\d)[\s\S]*')
def get_info(f):
try:
with open(f) as fid:
return (f,) + pat.match(fid.read()).groups()
except Exception:
return (f, f.replace('.in', '.out'), '00:30:00')
pbs_info_lst = map(get_info, pbs_files)
# sort wrt walltime
pbs_info_lst = sorted(pbs_info_lst, key=lambda fow: tuple(map(int, fow[2].split(':'))))[::-1]
# make dict {node1: pbs_info_lst1, ...} and save
d = dict([(f, pbs_info_lst[i::len(nodes)]) for i, f in enumerate(nodes)])
with open('pbs.dict', 'w') as fid:
fid.write(str(d))
def start_jobs(self):
import os
import multiprocessing
import platform
import time
with open('pbs.dict') as fid:
pbs_info_lst = eval(fid.read())[platform.node()]
arg_lst = ['echo starting %s && mkdir -p "%s" && env PBS_JOBID=$PBS_JOBID "%s" &> "%s" && echo finished %s' %
(f, os.path.dirname(o), f, o, f) for f, o, _ in pbs_info_lst]
print(arg_lst[0])
print('Starting %d jobs on %s' % (len(arg_lst), platform.node()))
pool = multiprocessing.Pool(int('$PBS_NUM_PPN'))
res = pool.map_async(os.system, arg_lst)
t = time.time()
for (f, _, _), r in zip(pbs_info_lst, res.get()):
print('%-50s\t%s' % (f, ('Errorcode %d' % r, 'Done')[r == 0]))
print('Done %d jobs on %s in %ds' % (len(arg_lst), platform.node(), time.time() - t))
def get_src(self, func):
src_lines = inspect.getsource(func).split("\n")[1:]
indent = len(src_lines[0]) - len(src_lines[0].lstrip())
src = "\n".join([l[indent:] for l in src_lines])
if func == self.make_dict:
src = src.replace("eval(self.PBS_FILES)", self.PBS_FILES)
return src
if __name__ == '__main__':
pbsmr = PBSMultiRunner(workdir="W:/simple1", queue='workq', nodes=2, ppn=10, pbsfiles=['/mnt/t1', "/mnt/t2"])
print(pbsmr.get_src(pbsmr.make_dict))
# pbsmr.save("c:/tmp/")
wetb/utils/cluster_tools/pbsfile_hawc2.py 0 → 100644
View file @ 62f0f2e4
from wetb.hawc2.htc_file import HTCFile
import os
from os import path
pbs_template = '''### Standard Output
#PBS -N [jobname]
#PBS -o [pbsoutdir]/[jobname].out
### Standard Error
#PBS -e [pbsoutdir]/[jobname].err
#PBS -W umask=0003
### Maximum wallclock time format HOURS:MINUTES:SECONDS
#PBS -l walltime=[walltime]
#PBS -l nodes=1:ppn=1
### Queue name
#PBS -q workq
# ==============================================================================
# single PBS mode: one case per PBS job
# evaluates to true if LAUNCH_PBS_MODE is NOT set
if [ -z ${LAUNCH_PBS_MODE+x} ] ; then
### Create scratch directory and copy data to it
cd $PBS_O_WORKDIR
echo "current working dir (pwd):"
pwd
cp -R ./[modelzip] /scratch/$USER/$PBS_JOBID
fi
# ==============================================================================
# ==============================================================================
# single PBS mode: one case per PBS job
# evaluates to true if LAUNCH_PBS_MODE is NOT set
if [ -z ${LAUNCH_PBS_MODE+x} ] ; then
echo
echo 'Execute commands on scratch nodes'
cd /scratch/$USER/$PBS_JOBID
# create unique dir for each CPU
mkdir "1"; cd "1"
pwd
/usr/bin/unzip ../[modelzip]
mkdir -p [htcdir]
mkdir -p [resdir]
mkdir -p [logdir]
mkdir -p [turbdir]
cp -R $PBS_O_WORKDIR/[htcdir]/[jobname].htc ./[htcdir]
cp -R $PBS_O_WORKDIR/[turbdir][turbfileroot]*.bin [turbdir]
_HOSTNAME_=`hostname`
if [[ ${_HOSTNAME_:0:1} == "j" ]] ; then
WINEARCH=win64 WINEPREFIX=~/.wine winefix
fi
# ==============================================================================
# ------------------------------------------------------------------------------
# find+xargs mode: 1 PBS job, multiple cases
else
# with find+xargs we first browse to CPU folder
cd "$CPU_NR"
fi
# ------------------------------------------------------------------------------
echo ""
# ==============================================================================
# single PBS mode: one case per PBS job
# evaluates to true if LAUNCH_PBS_MODE is NOT set
if [ -z ${LAUNCH_PBS_MODE+x} ] ; then
echo "execute HAWC2, fork to background"
time WINEARCH=win64 WINEPREFIX=~/.wine wine hawc2-latest ./[htcdir]/[jobname].htc &
wait
# ==============================================================================
# ------------------------------------------------------------------------------
# find+xargs mode: 1 PBS job, multiple cases
else
echo "execute HAWC2, do not fork and wait"
(time WINEARCH=win64 WINEPREFIX=~/.wine numactl --physcpubind=$CPU_NR wine hawc2-latest ./[htcdir]/[jobname].htc) 2>&1 | tee [pbsoutdir]/[jobname].err.out
fi
# ------------------------------------------------------------------------------
### Epilogue
# ==============================================================================
# single PBS mode: one case per PBS job
# evaluates to true if LAUNCH_PBS_MODE is NOT set
if [ -z ${LAUNCH_PBS_MODE+x} ] ; then
### wait for jobs to finish
wait
echo ""
echo "Copy back from scratch directory"
mkdir -p $PBS_O_WORKDIR/[resdir]
mkdir -p $PBS_O_WORKDIR/[logdir]
mkdir -p $PBS_O_WORKDIR/animation/
mkdir -p $PBS_O_WORKDIR/[turbdir]
cp -R [resdir]. $PBS_O_WORKDIR/[resdir].
cp -R [logdir]. $PBS_O_WORKDIR/[logdir].
cp -R animation/. $PBS_O_WORKDIR/animation/.
echo ""
echo "COPY BACK TURB IF APPLICABLE"
cd turb/
for i in `ls *.bin`; do if [ -e $PBS_O_WORKDIR/[turbdir]$i ]; then echo "$i exists no copyback"; else echo "$i copyback"; cp $i $PBS_O_WORKDIR/[turbdir]; fi; done
cd /scratch/$USER/$PBS_JOBID/1/
echo "END COPY BACK TURB"
echo ""
echo "COPYBACK [copyback_files]/[copyback_frename]"
echo "END COPYBACK"
echo ""
echo ""
echo "following files are on node/cpu 1 (find .):"
find .
# ==============================================================================
# ------------------------------------------------------------------------------
# find+xargs mode: 1 PBS job, multiple cases
else
cd /scratch/$USER/$PBS_JOBID/$CPU_NR/
rsync -a --remove-source-files [resdir]. ../HAWC2SIM/[resdir].
rsync -a --remove-source-files [logdir]. ../HAWC2SIM/[logdir].
rsync -a --remove-source-files [pbsoutdir]/. ../HAWC2SIM/[pbsoutdir]/.
rsync -a --remove-source-files animation/. ../HAWC2SIM/animation/.
echo ""
echo "COPY BACK TURB IF APPLICABLE"
cd turb/
for i in `ls *.bin`; do if [ -e $PBS_O_WORKDIR/[turbdir]$i ]; then echo "$i exists no copyback"; else echo "$i copyback"; cp $i $PBS_O_WORKDIR/[turbdir]; fi; done
cd /scratch/$USER/$PBS_JOBID/$CPU_NR/
echo "END COPY BACK TURB"
echo ""
echo "COPYBACK [copyback_files]/[copyback_frename]"
echo "END COPYBACK"
echo ""
# ------------------------------------------------------------------------------
fi
exit
'''
def htc2pbs(htc_fn, walltime='00:40:00', zipfile=None):
"""
Creates a PBS launch file (.p) based on a HAWC2 .htc file.
- Assumes htc files are within a htc/[casename]/ directory relative to current directory.
- Assumes there is a .zip file in the current directory which contains the turbine model.
If there is none, the zip file is set to 'model.zip' by default
- Will place a .p fine in pbs_in/[casename]/ directory relative to current directory.
-
Parameters
----------
htc_fn : str
The file name and path to the .htc file relative to current directory.
walltime: str (default='00:40:00')
A string indicating the walltime of the job of the form 'HH:MM:SS'
zipfile: str (default=None)
The filename of the zipfile containing the wind turbine model files and
HAWC2 executable. if zipfile=None, searches the current directory for a
zip file. If none is found, sets zipfile to 'model.zip'
Returns
-------
str
The filename and path to the output .p file
Raises
------
FileNotFoundError: If the file structure is not correct.
"""
basename = path.relpath(path.dirname(htc_fn), 'htc')
jobname = path.splitext(path.basename(htc_fn))[0]
pbs_in_dir = path.join('pbs_in', basename)
if basename == '.':
raise FileNotFoundError('File structure is incorrect.')
if not zipfile:
try:
zipfile = [x for x in os.listdir() if x.lower().endswith('.zip')][0]
except:
print('No .zip file found in current directory. Set model zip to \'model.zip\'')
zipfile = 'model.zip'
# get the required parameters for the pbs file from the htc file
htc = HTCFile(htc_fn) #modelpath='../..')
p = {
'walltime' : walltime,
'modelzip' : zipfile,
'jobname' : jobname,
'htcdir' : 'htc/' + basename,
'logdir' : path.dirname(htc.simulation.logfile.str_values())[2:] + '/',
'resdir' : path.dirname(htc.output.filename.str_values())[2:] + '/',
'turbdir' : path.dirname(htc.wind.mann.filename_u.str_values()) + '/',
'turbfileroot' : path.basename(htc.wind.mann.filename_u.str_values()).split('u.')[0],
'pbsoutdir' : 'pbs_out/' + basename
}
#Write pbs file based on template file and tags
if not os.path.exists(pbs_in_dir):
os.makedirs(pbs_in_dir)
template = str(pbs_template)
for key, value in p.items():
template = template.replace('[' + key + ']', value)
with open(os.path.join(pbs_in_dir, jobname + '.p'), 'w') as f:
f.write(template)
wetb/utils/cluster_tools/pbsjob.py
View file @ 62f0f2e4
......@@ -4,6 +4,9 @@ Created on 04/12/2015
@author: mmpe
'''
import os
import io
from wetb.utils.cluster_tools.pbsfile import PBSFile
from wetb.utils.cluster_tools.os_path import relpath
NOT_SUBMITTED = "Job not submitted"
PENDING = "Pending"
......@@ -11,26 +14,35 @@ RUNNING = "Running"
DONE = "Done"
def pjoin(*args):
return os.path.join(*args).replace('\\', '/')
class SSHPBSJob(object):
_status = NOT_SUBMITTED
nodeid = None
jobid = None
def __init__(self, sshClient):
self.ssh = sshClient
def submit(self, job, cwd, pbs_out_file):
def submit(self, pbsfile, cwd=None, pbs_out_file=None):
self.cwd = cwd
self.pbs_out_file = os.path.relpath(cwd + pbs_out_file).replace("\\", "/")
self.nodeid = None
#self.execute()
if isinstance(pbsfile, PBSFile):
f = io.StringIO(str(pbsfile))
f.seek(0)
pbs_filename = pjoin(cwd, pbsfile.filename)
self.ssh.upload(f, pbs_filename)
self.pbs_out_file = pjoin(cwd, pbsfile.stdout_filename)
cwd = pbsfile.workdir
pbsfile = pbsfile.filename
else:
self.pbs_out_file = os.path.relpath(cwd + pbs_out_file).replace("\\", "/")
cmds = ['rm -f %s' % self.pbs_out_file]
if cwd != "":
cmds.append("cd %s" % cwd)
cmds.append("qsub %s" % job)
cmds.append("qsub %s" % pbsfile)
_, out, _ = self.ssh.execute(";".join(cmds))
self.jobid = out.split(".")[0]
self._status = PENDING
......@@ -65,7 +77,6 @@ class SSHPBSJob(object):
return
raise e
def is_executing(self):
try:
self.ssh.execute("qstat %s" % self.jobid)
......
wetb/utils/cluster_tools/ssh_client.py
View file @ 62f0f2e4
......@@ -7,7 +7,7 @@ Created on 27/11/2015
from io import StringIO
import sys
import paramiko
import os
import threading
from _collections import deque
......@@ -17,39 +17,40 @@ import zipfile
import glob
from sshtunnel import SSHTunnelForwarder, SSH_CONFIG_FILE
from wetb.utils.ui import UI
from contextlib import contextmanager
import io
import tempfile
class SSHInteractiveAuthTunnelForwarder(SSHTunnelForwarder):
def __init__(
self,
interactive_auth_handler,
ssh_address_or_host=None,
ssh_config_file=SSH_CONFIG_FILE,
ssh_host_key=None,
ssh_password=None,
ssh_pkey=None,
ssh_private_key_password=None,
ssh_proxy=None,
ssh_proxy_enabled=True,
ssh_username=None,
local_bind_address=None,
local_bind_addresses=None,
logger=None,
mute_exceptions=False,
remote_bind_address=None,
remote_bind_addresses=None,
set_keepalive=0.0,
threaded=True,
compression=None,
allow_agent=True, *
args, **
kwargs):
self,
interactive_auth_handler,
ssh_address_or_host=None,
ssh_config_file=SSH_CONFIG_FILE,
ssh_host_key=None,
ssh_password=None,
ssh_pkey=None,
ssh_private_key_password=None,
ssh_proxy=None,
ssh_proxy_enabled=True,
ssh_username=None,
local_bind_address=None,
local_bind_addresses=None,
logger=None,
mute_exceptions=False,
remote_bind_address=None,
remote_bind_addresses=None,
set_keepalive=0.0,
threaded=True,
compression=None,
allow_agent=True, *
args, **
kwargs):
self.interactive_auth_handler = interactive_auth_handler
SSHTunnelForwarder.__init__(self, ssh_address_or_host=ssh_address_or_host, ssh_config_file=ssh_config_file, ssh_host_key=ssh_host_key, ssh_password=ssh_password, ssh_pkey=ssh_pkey, ssh_private_key_password=ssh_private_key_password, ssh_proxy=ssh_proxy, ssh_proxy_enabled=ssh_proxy_enabled, ssh_username=ssh_username, local_bind_address=local_bind_address, local_bind_addresses=local_bind_addresses, logger=logger, mute_exceptions=mute_exceptions, remote_bind_address=remote_bind_address, remote_bind_addresses=remote_bind_addresses, set_keepalive=set_keepalive, threaded=threaded, compression=compression, allow_agent=allow_agent, *args, **kwargs)
SSHTunnelForwarder.__init__(self, ssh_address_or_host=ssh_address_or_host, ssh_config_file=ssh_config_file, ssh_host_key=ssh_host_key, ssh_password=ssh_password, ssh_pkey=ssh_pkey, ssh_private_key_password=ssh_private_key_password, ssh_proxy=ssh_proxy, ssh_proxy_enabled=ssh_proxy_enabled, ssh_username=ssh_username,
local_bind_address=local_bind_address, local_bind_addresses=local_bind_addresses, logger=logger, mute_exceptions=mute_exceptions, remote_bind_address=remote_bind_address, remote_bind_addresses=remote_bind_addresses, set_keepalive=set_keepalive, threaded=threaded, compression=compression, allow_agent=allow_agent, *args, **kwargs)
def _connect_to_gateway(self):
"""
Open connection to SSH gateway
......@@ -66,9 +67,9 @@ class SSHInteractiveAuthTunnelForwarder(SSHTunnelForwarder):
except paramiko.AuthenticationException:
self.logger.debug('Authentication error')
self._stop_transport()
self.logger.error('Could not open connection to gateway')
def _connect_to_gateway_old(self):
"""
Open connection to SSH gateway
......@@ -88,38 +89,49 @@ class SSHInteractiveAuthTunnelForwarder(SSHTunnelForwarder):
self._transport.connect(hostkey=self.ssh_host_key,
username=self.ssh_username,
password=self.ssh_password)
if self._transport.is_alive:
return
except paramiko.AuthenticationException:
self.logger.debug('Authentication error')
self._stop_transport()
self.logger.error('Could not open connection to gateway')
def unix_path(path):
drive, p = os.path.splitdrive(str(path))
return os.path.join(drive, os.path.splitdrive(os.path.realpath(p))[1]).replace("\\", "/")
class SSHClient(object):
"A wrapper of paramiko.SSHClient"
TIMEOUT = 4
def __init__(self, host, username, password=None, port=22, key=None, passphrase=None, interactive_auth_handler=None, gateway=None, ui=UI()):
def __init__(self, host, username, password=None, port=22, key=None,
passphrase=None, interactive_auth_handler=None, gateway=None, ui=UI()):
self.host = host
self.username = username
self.password = password
if password is None and key is None:
from os.path import expanduser
home = expanduser("~")
key = home + '/.ssh/id_rsa'
self.port = port
self.key = key
self.gateway=gateway
self.gateway = gateway
self.interactive_auth_handler = interactive_auth_handler
self.ui = ui
self.disconnect = 0
self.client = None
self.ssh_lock = threading.RLock()
#self.sftp = None
self._sftp = None
self.transport = None
self.counter_lock = threading.RLock()
self.counter=0
self.counter = 0
if key is not None:
self.key = paramiko.RSAKey.from_private_key(StringIO(key), password=passphrase)
with open(key) as fid:
self.key = paramiko.RSAKey.from_private_key(fid, password=passphrase)
def info(self):
return self.host, self.username, self.password, self.port
......@@ -147,50 +159,51 @@ class SSHClient(object):
ssh_password=self.gateway.password,
remote_bind_address=(self.host, self.port),
local_bind_address=('0.0.0.0', 10022)
)
)
else:
self.tunnel = SSHTunnelForwarder((self.gateway.host, self.gateway.port),
ssh_username=self.gateway.username,
ssh_password=self.gateway.password,
remote_bind_address=(self.host, self.port),
local_bind_address=('0.0.0.0', 10022)
)
print ("start tunnel")
)
print("start tunnel")
self.tunnel.start()
print ("self.client = paramiko.SSHClient()")
print("self.client = paramiko.SSHClient()")
self.client = paramiko.SSHClient()
self.client.set_missing_host_key_policy(paramiko.AutoAddPolicy())
print ('self.client.connect("127.0.0.1", 10022, username=self.username, password=self.password)')
print('self.client.connect("127.0.0.1", 10022, username=self.username, password=self.password)')
self.client.connect("127.0.0.1", 10022, username=self.username, password=self.password)
print ("done")
print("done")
elif self.password is None or self.password == "":
raise IOError("Password not set for %s"%self.host)
elif self.key is None and (self.password is None or self.password == ""):
raise IOError("Password not set for %s" % self.host)
else:
self.client = paramiko.SSHClient()
self.client.set_missing_host_key_policy(paramiko.AutoAddPolicy())
try:
self.client.connect(self.host, self.port, username=self.username, password=self.password, pkey=self.key, timeout=self.TIMEOUT)
self.client.connect(self.host, self.port, username=self.username,
password=self.password, pkey=self.key, timeout=self.TIMEOUT)
except paramiko.ssh_exception.SSHException as e:
transport = self.client.get_transport()
transport.auth_interactive(self.username, self.interactive_auth_handler)
assert self.client is not None
#self.sftp = paramiko.SFTPClient.from_transport(self.client._transport)
return self
def __del__(self):
self.close()
@property
def sftp(self):
return paramiko.SFTPClient.from_transport(self.client._transport)
if self._sftp is None:
self._sftp = paramiko.SFTPClient.from_transport(self.client._transport)
return self._sftp
# @sftp.setter
# def sftp(self, values):
# pass
......@@ -200,17 +213,16 @@ class SSHClient(object):
if self.disconnect == 0:
self.close()
def download(self, remotefilepath, localfile, verbose=False, retry=1, callback=None):
if verbose:
ret = None
print ("Download %s > %s" % (remotefilepath, str(localfile)))
print("Download %s > %s" % (remotefilepath, str(localfile)))
if callback is None:
callback = self.ui.progress_callback()
remotefilepath = remotefilepath.replace("\\", "/")
for i in range(retry):
if i>0:
print ("Retry download %s, #%d"%(remotefilepath, i))
if i > 0:
print("Retry download %s, #%d" % (remotefilepath, i))
try:
SSHClient.__enter__(self)
......@@ -219,112 +231,126 @@ class SSHClient(object):
elif hasattr(localfile, 'write'):
ret = self.sftp.putfo(remotefilepath, localfile, callback=callback)
break
except:
except Exception:
pass
finally:
SSHClient.__exit__(self)
print ("Download %s failed from %s"%(remotefilepath, self.host))
if verbose:
print (ret)
def upload(self, localfile, filepath, verbose=False, callback=None):
print("Download %s failed from %s" % (remotefilepath, self.host))
if verbose and ret:
print(ret)
def upload(self, localfile, filepath, chmod="770", verbose=False, callback=None):
if verbose:
print ("Upload %s > %s" % (localfile, filepath))
print("Upload %s > %s" % (localfile, filepath))
if callback is None:
callback = self.ui.progress_callback()
try:
SSHClient.__enter__(self)
self.execute('mkdir -p "%s"' % (os.path.dirname(filepath)))
sftp = self.sftp
if isinstance(localfile, (str, bytes, int)):
ret = self.sftp.put(localfile, filepath, callback=callback)
ret = sftp.put(localfile, filepath, callback=callback)
elif hasattr(localfile, 'read'):
size = len(localfile.read())
localfile.seek(0)
ret = self.sftp.putfo(localfile, filepath, file_size=size, callback=callback)
ret = sftp.putfo(localfile, filepath, file_size=size, callback=callback)
self.execute('chmod %s "%s"' % (chmod, filepath))
except Exception as e:
print ("upload failed ", str(e))
print("upload failed ", str(e))
raise e
finally:
SSHClient.__exit__(self)
if verbose:
print (ret)
if verbose and ret:
print(ret)
def upload_files(self, localpath, remotepath, file_lst=["."], compression_level=1, callback=None):
assert os.path.isdir(localpath)
if not isinstance(file_lst, (tuple, list)):
file_lst = [file_lst]
files = ([os.path.join(root, f) for fp in file_lst for root,_,files in os.walk(os.path.join(localpath, fp )) for f in files] +
[f for fp in file_lst for f in glob.glob(os.path.join(localpath, fp)) ])
files = set([os.path.abspath(f) for f in files])
files = ([os.path.join(root, f) for fp in file_lst for root, _, files in os.walk(os.path.join(localpath, fp)) for f in files] +
[f for fp in file_lst for f in glob.glob(os.path.join(localpath, fp))])
files = set([os.path.abspath(f) for f in files if os.path.isfile(f)])
compression_levels = {0:zipfile.ZIP_STORED, 1:zipfile.ZIP_DEFLATED, 2:zipfile.ZIP_BZIP2, 3:zipfile.ZIP_LZMA}
compression_levels = {0: zipfile.ZIP_STORED, 1: zipfile.ZIP_DEFLATED, 2: zipfile.ZIP_BZIP2, 3: zipfile.ZIP_LZMA}
with self.counter_lock:
self.counter+=1
zn = 'tmp_%s_%04d.zip'%(id(self),self.counter)
self.counter += 1
zn = 'tmp_%s_%04d.zip' % (id(self), self.counter)
zipf = zipfile.ZipFile(zn, 'w', compression_levels[compression_level])
try:
for f in files:
zipf.write(f, os.path.relpath(f, localpath))
zipf.close()
remote_zn = os.path.join(remotepath, zn).replace("\\","/")
remote_zn = os.path.join(remotepath, zn).replace("\\", "/")
with self:
self.execute("mkdir -p %s"%(remotepath))
self.execute('mkdir -p "%s"' % (remotepath))
self.upload(zn, remote_zn, callback=callback)
self.execute("unzip %s -d %s && rm %s"%(remote_zn, remotepath, remote_zn))
except:
print ("upload files failed", )
self.execute('unzip -o "%s" -d "%s" && rm "%s"' % (remote_zn, remotepath, remote_zn))
except Exception:
print("upload files failed", )
traceback.print_exc()
raise
finally:
os.remove(zn)
def download_files(self, remote_path, localpath, file_lst=["."], compression_level=1, callback=None):
if not isinstance(file_lst, (tuple, list)):
file_lst = [file_lst]
file_lst = [f.replace("\\","/") for f in file_lst]
file_lst = [f.replace("\\", "/") for f in file_lst]
with self.counter_lock:
self.counter+=1
zn = 'tmp_%s_%04d.zip'%(id(self),self.counter)
remote_zip = os.path.join(remote_path, zn).replace("\\","/")
self.execute("cd %s && zip -r %s %s"%(remote_path, zn, " ".join(file_lst)))
self.counter += 1
zn = 'tmp_%s_%04d.zip' % (id(self), self.counter)
remote_zip = os.path.join(remote_path, zn).replace("\\", "/")
self.execute('cd "%s" && zip -r "%s" "%s"' % (remote_path, zn, " ".join(file_lst)))
local_zip = os.path.join(localpath, zn)
if not os.path.isdir(localpath):
os.makedirs(localpath)
self.download(remote_zip, local_zip, callback=callback)
self.execute("rm -f %s" % remote_zip)
self.execute('rm -f "%s"' % remote_zip)
with zipfile.ZipFile(local_zip, "r") as z:
z.extractall(localpath)
os.remove(local_zip)
def close(self):
for x in ["client", 'tunnel' ]:
for x in ["_sftp", "client", "tunnel"]:
try:
getattr(self, x).close()
setattr(self, x, None)
except:
except Exception:
pass
self.disconnect = False
def file_exists(self, filename):
_, out, _ = (self.execute('[ -f %s ] && echo "File exists" || echo "File does not exists"' % filename.replace("\\", "/")))
_, out, _ = (self.execute(
'[ -f "%s" ] && echo "File exists" || echo "File does not exists"' % unix_path(filename)))
return out.strip() == "File exists"
def execute(self, command, sudo=False, verbose=False):
def isfile(self, filename):
return self.file_exists(filename)
def folder_exists(self, folder):
_, out, _ = (self.execute(
'[ -d "%s" ] && echo "Folder exists" || echo "Folder does not exists"' % unix_path(folder)))
return out.strip() == "Folder exists"
def isdir(self, folder):
return self.folder_exists(folder)
def execute(self, command, cwd='.', sudo=False, verbose=False):
feed_password = False
if sudo and self.username != "root":
command = "sudo -S -p '' %s" % command
feed_password = self.password is not None and len(self.password) > 0
if isinstance(command, (list, tuple)):
command = "\n".join(command)
#cwd = unix_path(cwd)
if verbose:
print (">>> " + command)
print("[%s]$ %s" % (cwd, command))
command = 'cd "%s" && %s' % (cwd, command)
with self as ssh:
if ssh is None:
exc_info = sys.exc_info()
......@@ -337,9 +363,8 @@ class SSHClient(object):
v, out, err = stdout.channel.recv_exit_status(), stdout.read().decode(), stderr.read().decode()
if v:
raise Warning ("out:\n%s\n----------\nerr:\n%s" % (out, err))
raise Warning("out:\n%s\n----------\nerr:\n%s" % (out, err))
elif verbose:
if out:
sys.stdout.write(out)
......@@ -348,7 +373,8 @@ class SSHClient(object):
return v, out, err
def append_wine_path(self, path):
ret = self.execute('wine regedit /E tmp.reg "HKEY_LOCAL_MACHINE\System\CurrentControlSet\Control\Session Manager\Environment"')
ret = self.execute(
r'wine regedit /E tmp.reg "HKEY_LOCAL_MACHINE\System\CurrentControlSet\Control\Session Manager\Environment"')
self.download('tmp.reg', 'tmp.reg')
with open('tmp.reg') as fid:
lines = fid.readlines()
......@@ -372,67 +398,79 @@ class SSHClient(object):
cwd = os.path.join(cwd, os.path.split(filepattern)[0]).replace("\\", "/")
filepattern = os.path.split(filepattern)[1]
if recursive:
_, out, _ = self.execute(r'find %s -type f -name "%s"' % (cwd, filepattern))
_, out, _ = self.execute(r'find "%s" -type f -name "%s"' % (cwd, filepattern))
else:
_, out, _ = self.execute(r'find %s -maxdepth 1 -type f -name "%s"' % (cwd, filepattern))
_, out, _ = self.execute(r'find "%s" -maxdepth 1 -type f -name "%s"' % (cwd, filepattern))
return [file for file in out.strip().split("\n") if file != ""]
def listdir(self, folder):
_, out, _ = self.execute('ls -p', cwd=folder)
return out.split()
@contextmanager
def open(self, filename, mode='r', encoding=None, newline=None):
with tempfile.NamedTemporaryFile(mode='x', delete=False, suffix=os.path.splitext(filename)[1]) as tmp:
tmp_name = tmp.name
if mode in 'ra+':
# try:
self.download(filename, tmp_name)
try:
fid = open(tmp_name, mode=mode, encoding=encoding, newline=newline)
yield fid
finally:
fid.close()
if mode in 'wa+':
self.upload(tmp_name, filename)
os.remove(tmp_name)
class SharedSSHClient(SSHClient):
def __init__(self, host, username, password=None, port=22, key=None, passphrase=None, interactive_auth_handler=None, gateway=None):
SSHClient.__init__(self, host, username, password=password, port=port, key=key, passphrase=passphrase, interactive_auth_handler=interactive_auth_handler, gateway=gateway)
def __init__(self, host, username, password=None, port=22, key=None,
passphrase=None, interactive_auth_handler=None, gateway=None):
SSHClient.__init__(self, host, username, password=password, port=port, key=key,
passphrase=passphrase, interactive_auth_handler=interactive_auth_handler, gateway=gateway)
self.shared_ssh_queue = deque()
self.next = None
def execute(self, command, sudo=False, verbose=False):
res = SSHClient.execute(self, command, sudo=sudo, verbose=verbose)
return res
def __enter__(self):
with self.ssh_lock:
SSHClient.__enter__(self)
#print ("request SSH", threading.currentThread())
# if len(self.shared_ssh_queue)>0 and self.shared_ssh_queue[0] == threading.get_ident():
# # SSH already allocated to this thread ( multiple use-statements in "with ssh:" block
# # SSH already allocated to this thread ( multiple use-statements in "with ssh:" block
# self.shared_ssh_queue.appendleft(threading.get_ident())
# else:
# self.shared_ssh_queue.append(threading.get_ident())
if len(self.shared_ssh_queue)>0 and self.shared_ssh_queue[0] == threading.get_ident():
# SSH already allocated to this thread ( multiple use-statements in "with ssh:" block
if len(self.shared_ssh_queue) > 0 and self.shared_ssh_queue[0] == threading.get_ident():
# SSH already allocated to this thread ( multiple use-statements in "with ssh:" block
self.shared_ssh_queue.popleft()
self.shared_ssh_queue.append(threading.get_ident())
while self.shared_ssh_queue[0] != threading.get_ident():
time.sleep(2)
return self.client
def __exit__(self, *args):
with self.ssh_lock:
if len(self.shared_ssh_queue)>0 and self.shared_ssh_queue[0] == threading.get_ident():
if len(self.shared_ssh_queue) > 0 and self.shared_ssh_queue[0] == threading.get_ident():
self.shared_ssh_queue.popleft()
if __name__ == "__main__":
from mmpe.ui.qt_ui import QtInputUI
q = QtInputUI(None)
x = None
username, password = "mmpe", x.password #q.get_login("mmpe")
client = SSHClient(host='gorm', port=22, username=username, password=password)
print (client.glob("*.*", ".hawc2launcher/medium1__1__"))
# ssh.upload('../News.txt', 'news.txt')
try:
import x
username, password = 'mmpe', x.password
client = SSHClient(host='jess.dtu.dk', port=22, username=username, password=password)
print(client.execute("echo hello $USER from $HOSTNAME")[1])
except ImportError:
x = None
wetb/utils/envelope.py
View file @ 62f0f2e4
......@@ -6,25 +6,16 @@ Created on Tue Feb 13 12:58:25 2018
@author: dave
"""
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import dict
from io import open
from builtins import zip
from builtins import range
from builtins import str
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import object
import numpy as np
import numpy.typing as npt
from typing import List, Tuple, Union
import scipy
from wetb.utils.rotation import projection_2d
def compute_env_of_env(envelope, dlc_list, Nx=300, Nsectors=12, Ntheta=181):
def compute_env_of_env(
envelope, dlc_list: list, Nx: int = 300, Nsectors: int = 12, Ntheta: int = 181
) -> npt.NDArray[np.float64]:
"""
The function computes load envelopes for given channels and a groups of
load cases starting from the envelopes computed for single simulations.
......@@ -64,20 +55,21 @@ def compute_env_of_env(envelope, dlc_list, Nx=300, Nsectors=12, Ntheta=181):
"""
# Group all the single DLCs
cloud = np.zeros(((Nx+1)*len(envelope),6))
cloud = np.zeros(((Nx + 1) * len(envelope), 6))
for i in range(len(envelope)):
cloud[(Nx+1)*i:(Nx+1)*(i+1),:] = envelope[dlc_list[i]]
cloud[(Nx + 1) * i : (Nx + 1) * (i + 1), :] = envelope[dlc_list[i]]
# Compute total Hull of all the envelopes
hull = scipy.spatial.ConvexHull(cloud[:,:2])
cc = np.append(cloud[hull.vertices,:2],
cloud[hull.vertices[0],:2].reshape(1,2),axis=0)
hull = scipy.spatial.ConvexHull(cloud[:, :2])
cc = np.append(
cloud[hull.vertices, :2], cloud[hull.vertices[0], :2].reshape(1, 2), axis=0
)
# Interpolate full envelope
cc_x,cc_up,cc_low,cc_int= int_envelope(cc[:,0], cc[:,1], Nx=Nx)
cc_x, cc_up, cc_low, cc_int = int_envelope(cc[:, 0], cc[:, 1], Nx=Nx)
# Project full envelope on given direction
cc_proj = proj_envelope(cc_x, cc_up, cc_low, cc_int, Nx, Nsectors, Ntheta)
env_proj = np.zeros([len(cc_proj),6])
env_proj[:,:2] = cc_proj
env_proj = np.zeros([len(cc_proj), 6])
env_proj[:, :2] = cc_proj
# Based on Mx and My, gather the remaining cross-sectional forces and
# moments
......@@ -87,19 +79,27 @@ def compute_env_of_env(envelope, dlc_list, Nx=300, Nsectors=12, Ntheta=181):
s0 = np.append(s0, s1, axis=0)
cc = np.append(cc, s0, axis=1)
_,_,_,extra_sensor = int_envelope(cc[:,0],cc[:,ich],Nx)
es = np.atleast_2d(np.array(extra_sensor[:,1])).T
cc_int = np.append(cc_int,es,axis=1)
_, _, _, extra_sensor = int_envelope(cc[:, 0], cc[:, ich], Nx)
es = np.atleast_2d(np.array(extra_sensor[:, 1])).T
cc_int = np.append(cc_int, es, axis=1)
for isec in range(Nsectors):
ids = (np.abs(cc_int[:,0]-cc_proj[isec,0])).argmin()
env_proj[isec,ich] = (cc_int[ids-1,ich]+cc_int[ids,ich]+\
cc_int[ids+1,ich])/3
ids = (np.abs(cc_int[:, 0] - cc_proj[isec, 0])).argmin()
env_proj[isec, ich] = (
cc_int[ids - 1, ich] + cc_int[ids, ich] + cc_int[ids + 1, ich]
) / 3
return env_proj
def int_envelope(ch1, ch2, Nx):
def int_envelope(
ch1, ch2, Nx: int
) -> Tuple[
npt.NDArray[np.float64],
npt.NDArray[np.float64],
npt.NDArray[np.float64],
npt.NDArray[np.float64],
]:
"""
Function to interpolate envelopes and output arrays of same length
......@@ -113,31 +113,47 @@ def int_envelope(ch1, ch2, Nx):
indmax = np.argmax(ch1)
indmin = np.argmin(ch1)
if indmax > indmin:
lower = np.array([ch1[indmin:indmax+1],ch2[indmin:indmax+1]]).T
upper = np.concatenate((np.array([ch1[indmax:],ch2[indmax:]]).T,
np.array([ch1[:indmin+1],ch2[:indmin+1]]).T),
axis=0)
lower = np.array([ch1[indmin : indmax + 1], ch2[indmin : indmax + 1]]).T
upper = np.concatenate(
(
np.array([ch1[indmax:], ch2[indmax:]]).T,
np.array([ch1[: indmin + 1], ch2[: indmin + 1]]).T,
),
axis=0,
)
else:
upper = np.array([ch1[indmax:indmin+1],ch2[indmax:indmin+1]]).T
lower = np.concatenate((np.array([ch1[indmin:],ch2[indmin:]]).T,
np.array([ch1[:indmax+1],ch2[:indmax+1]]).T),
axis=0)
int_1 = np.linspace(min(upper[:,0].min(),lower[:,0].min()),
max(upper[:,0].max(),lower[:,0].max()),Nx/2+1)
upper = np.array([ch1[indmax : indmin + 1], ch2[indmax : indmin + 1]]).T
lower = np.concatenate(
(
np.array([ch1[indmin:], ch2[indmin:]]).T,
np.array([ch1[: indmax + 1], ch2[: indmax + 1]]).T,
),
axis=0,
)
int_1 = np.linspace(
min(upper[:, 0].min(), lower[:, 0].min()),
max(upper[:, 0].max(), lower[:, 0].max()),
Nx / 2 + 1,
)
upper = np.flipud(upper)
int_2_up = np.interp(int_1,np.array(upper[:,0]),np.array(upper[:,1]))
int_2_low = np.interp(int_1,np.array(lower[:,0]),np.array(lower[:,1]))
int_2_up = np.interp(int_1, np.array(upper[:, 0]), np.array(upper[:, 1]))
int_2_low = np.interp(int_1, np.array(lower[:, 0]), np.array(lower[:, 1]))
int_env = np.concatenate((np.array([int_1[:-1],int_2_up[:-1]]).T,
np.array([int_1[::-1],int_2_low[::-1]]).T),
axis=0)
int_env = np.concatenate(
(
np.array([int_1[:-1], int_2_up[:-1]]).T,
np.array([int_1[::-1], int_2_low[::-1]]).T,
),
axis=0,
)
return int_1, int_2_up, int_2_low, int_env
def proj_envelope(env_x, env_up, env_low, env, Nx, Nsectors, Ntheta):
def proj_envelope(
env_x, env_up, env_low, env, Nx, Nsectors, Ntheta
) -> npt.NDArray[np.float64]:
"""
Function to project envelope on given angles
......@@ -146,72 +162,76 @@ def proj_envelope(env_x, env_up, env_low, env, Nx, Nsectors, Ntheta):
cartesian
"""
theta_int = np.linspace(-np.pi,np.pi,Ntheta)
sectors = np.linspace(-np.pi,np.pi,Nsectors+1)
proj = np.zeros([Nsectors,2])
theta_int = np.linspace(-np.pi, np.pi, Ntheta)
sectors = np.linspace(-np.pi, np.pi, Nsectors + 1)
proj = np.zeros([Nsectors, 2])
R_up = np.sqrt(env_x**2+env_up**2)
theta_up = np.arctan2(env_up,env_x)
R_up = np.sqrt(env_x**2 + env_up**2)
theta_up = np.arctan2(env_up, env_x)
R_low = np.sqrt(env_x**2+env_low**2)
theta_low = np.arctan2(env_low,env_x)
R_low = np.sqrt(env_x**2 + env_low**2)
theta_low = np.arctan2(env_low, env_x)
R = np.concatenate((R_up,R_low))
theta = np.concatenate((theta_up,theta_low))
R = np.concatenate((R_up, R_low))
theta = np.concatenate((theta_up, theta_low))
R = R[np.argsort(theta)]
theta = np.sort(theta)
R_int = np.interp(theta_int,theta,R,period=2*np.pi)
R_int = np.interp(theta_int, theta, R, period=2 * np.pi)
for i in range(Nsectors):
if sectors[i]>=-np.pi and sectors[i+1]<-np.pi/2:
indices = np.where(np.logical_and(theta_int >= sectors[i],
theta_int <= sectors[i+1]))
if sectors[i] >= -np.pi and sectors[i + 1] < -np.pi / 2:
indices = np.where(
np.logical_and(theta_int >= sectors[i], theta_int <= sectors[i + 1])
)
maxR = R_int[indices].max()
proj[i+1,0] = maxR*np.cos(sectors[i+1])
proj[i+1,1] = maxR*np.sin(sectors[i+1])
elif sectors[i]==-np.pi/2:
proj[i + 1, 0] = maxR * np.cos(sectors[i + 1])
proj[i + 1, 1] = maxR * np.sin(sectors[i + 1])
elif sectors[i] == -np.pi / 2:
continue
elif sectors[i]>-np.pi/2 and sectors[i+1]<=0:
indices = np.where(np.logical_and(theta_int >= sectors[i],
theta_int <= sectors[i+1]))
elif sectors[i] > -np.pi / 2 and sectors[i + 1] <= 0:
indices = np.where(
np.logical_and(theta_int >= sectors[i], theta_int <= sectors[i + 1])
)
maxR = R_int[indices].max()
proj[i,0] = maxR*np.cos(sectors[i])
proj[i,1] = maxR*np.sin(sectors[i])
elif sectors[i]>=0 and sectors[i+1]<np.pi/2:
indices = np.where(np.logical_and(theta_int >= sectors[i],
theta_int <= sectors[i+1]))
proj[i, 0] = maxR * np.cos(sectors[i])
proj[i, 1] = maxR * np.sin(sectors[i])
elif sectors[i] >= 0 and sectors[i + 1] < np.pi / 2:
indices = np.where(
np.logical_and(theta_int >= sectors[i], theta_int <= sectors[i + 1])
)
maxR = R_int[indices].max()
proj[i+1,0] = maxR*np.cos(sectors[i+1])
proj[i+1,1] = maxR*np.sin(sectors[i+1])
elif sectors[i]==np.pi/2:
proj[i + 1, 0] = maxR * np.cos(sectors[i + 1])
proj[i + 1, 1] = maxR * np.sin(sectors[i + 1])
elif sectors[i] == np.pi / 2:
continue
elif sectors[i]>np.pi/2 and sectors[i+1]<=np.pi:
indices = np.where(np.logical_and(theta_int >= sectors[i],
theta_int <= sectors[i+1]))
elif sectors[i] > np.pi / 2 and sectors[i + 1] <= np.pi:
indices = np.where(
np.logical_and(theta_int >= sectors[i], theta_int <= sectors[i + 1])
)
maxR = R_int[indices].max()
proj[i,0] = maxR*np.cos(sectors[i])
proj[i,1] = maxR*np.sin(sectors[i])
proj[i, 0] = maxR * np.cos(sectors[i])
proj[i, 1] = maxR * np.sin(sectors[i])
ind = np.where(sectors==0)
proj[ind,0] = env[:,0].max()
ind = np.where(sectors == 0)
proj[ind, 0] = env[:, 0].max()
ind = np.where(sectors==np.pi/2)
proj[ind,1] = env[:,1].max()
ind = np.where(sectors == np.pi / 2)
proj[ind, 1] = env[:, 1].max()
ind = np.where(sectors==-np.pi)
proj[ind,0] = env[:,0].min()
ind = np.where(sectors == -np.pi)
proj[ind, 0] = env[:, 0].min()
ind = np.where(sectors==-np.pi/2)
proj[ind,1] = env[:,1].min()
ind = np.where(sectors == -np.pi / 2)
proj[ind, 1] = env[:, 1].min()
return proj
def closed_contour(cloud):
"""Return indices of the vertices of the closed convex contour that
contain all points of the (x,y) cloud of pairs.
def closed_contour(
cloud: npt.NDArray[np.float64],
) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
"""Returns a tuple of the vertices of the closed convex contour, and the indices of the vertices in the input array.
Parameters
----------
......@@ -220,36 +240,38 @@ def closed_contour(cloud):
Coordinates of points to construct a convex hull from. Ndim should be
at least 2 or higher.
cloud_extra : np.ndarray(npoints, nr_extra_channels)
Returns
-------
ivertices : ndarray(nvertices)
Indices of the coordinates in cloud that make out the vertices of the closed convex contour.
vertices : ndarray(nvertices, 2+nr_extra_channels)
The coordinates of the vertices that make out the closed convex contour of cloud.
"""
if not cloud.shape[1] >= 2:
raise IndexError('Cloud dimension should be 2 or greater')
raise IndexError("Cloud dimension should be 2 or greater")
hull = scipy.spatial.ConvexHull(cloud[:,0:2])
hull = scipy.spatial.ConvexHull(cloud[:, 0:2])
# indices to the vertices containing the cloud of the first 2 dimensions
ivertices = np.ndarray((len(hull.vertices)+1,), dtype=np.int32)
ivertices = np.ndarray((len(hull.vertices) + 1,), dtype=np.int32)
ivertices[0:-1] = hull.vertices
ivertices[-1] = hull.vertices[0]
# actual vertices of all dimensions in the cloud, based on the first two
# dimensions
vertices = np.ndarray((len(ivertices), cloud.shape[1]))
vertices[0:-1,:] = cloud[ivertices[0:-1],:]
vertices[-1,:] = cloud[ivertices[-1],:]
vertices[0:-1, :] = cloud[ivertices[0:-1], :]
vertices[-1, :] = cloud[ivertices[-1], :]
return vertices, ivertices
def compute_envelope(cloud, int_env=False, Nx=300):
def compute_envelope(
cloud: npt.NDArray[np.float64], int_env: bool = False, Nx: int = 300
) -> npt.NDArray[np.float64]:
"""
The function computes load envelopes for given signals and a single
load case. Starting from Mx and My moments, the other cross-sectional
......@@ -258,7 +280,7 @@ def compute_envelope(cloud, int_env=False, Nx=300):
Parameters
----------
x, y : np.ndarray(n)
x, y : np.ndarray(int, 2)
2 components of the same time signal, e.g x and y.
int_env : boolean, default=False
......@@ -271,10 +293,13 @@ def compute_envelope(cloud, int_env=False, Nx=300):
Returns
-------
envelope : dictionary,
The dictionary has entries refered to the channels selected.
Inside the dictonary under each entry there is a matrix with 6
columns, each for the sectional forces and moments
vertices : np.ndarray(int, 2),
Returns an array of the vertices of the closed contour. The
# envelope : dictionary,
# The dictionary has entries refered to the channels selected.
# Inside the dictonary under each entry there is a matrix with 6
# columns, each for the sectional forces and moments
"""
......@@ -282,12 +307,142 @@ def compute_envelope(cloud, int_env=False, Nx=300):
# interpolate to a fixed location of equally spaced vertices
if int_env:
vert_int = np.ndarray((Nx+1, cloud.shape[1]))
_,_,_,vert_int[:,0:2] = int_envelope(vertices[:,0], vertices[:,1], Nx)
vert_int = np.ndarray((Nx + 1, cloud.shape[1]))
_, _, _, vert_int[:, 0:2] = int_envelope(vertices[:, 0], vertices[:, 1], Nx)
for i in range(2, cloud.shape[1]):
_,_,_,extra = int_envelope(vertices[:,0], vertices[:,i], Nx)
vert_int[:,i] = extra[:,1]
_, _, _, extra = int_envelope(vertices[:, 0], vertices[:, i], Nx)
vert_int[:, i] = extra[:, 1]
vertices = vert_int
return vertices
def projected_extremes(
signal: npt.NDArray,
angles: npt.NDArray = np.linspace(-150, 180, 12),
sweep_angle: Union[float, None] = None,
degrees: bool = True
) -> npt.NDArray:
"""_summary_
Parameters
----------
signal : npt.NDArray
2-Dimensional array, corresponding to a time-series or a load envelope for two given channels. The first column is treated as the x-coordinate, and the second column as the y-coordinate.
angles : npt.NDArray, optional
List of angles to project signals onto, before calculating the maximum value, by default numpy.linspace(-150,180, 12). Angles are given in degrees. The angles will be returned in the range (-180:180]
sweep_angle : float | None, optional
Sweep angle, the allowed deviation from the projection angle for the maximum loads, by default None. If None, no sweep angle is applied to the calculation. If no loads exist within the sweep angle, a load of 0, and the index nan is returned.
Returns
-------
npt.NDArray
_description_
"""
if degrees:
angles = np.deg2rad(angles)
if sweep_angle:
sweep_angle = np.deg2rad(sweep_angle)
# Initialize output variables
extremes = np.zeros(shape=(len(angles), 3))
# rearrange angles to (-pi; pi]
angles[angles > np.pi] -= 2*np.pi
angles[angles <= -np.pi] += 2*np.pi
# Calculate the angle of the signal in the time-series
signal_angles = np.angle(signal[:, 0]+1j*signal[:, 1])
for index, angle in enumerate(angles):
# Project signal into the desired angle, saving only the first component of the 2D load for extremes analysis
projected_signal = signal @ projection_2d(angle, degrees=False)
if sweep_angle:
# Remove the loads not covered by the main angle +/- the sweep angle
lower_bound = angle - sweep_angle
if lower_bound<-np.pi: lower_bound +=2*np.pi
upper_bound = angle + sweep_angle
if upper_bound>np.pi: upper_bound -=2*np.pi
if upper_bound > lower_bound:
# If lower bound is below upper bound, and the range is thus monotone, use the intersection of angles above lower- and below upper bound.
projected_signal = (
projected_signal
* (signal_angles > lower_bound)
* (signal_angles < upper_bound)
)
else:
# If lower bound is above upper bound, use the union of angles above lower- and below upper bounds.
projected_signal = (
projected_signal
* (signal_angles > lower_bound)
+ projected_signal
* (signal_angles < upper_bound)
)
if (projected_signal.any()) == True:
# Save the larges projected load that satisfies the angle+sweep criteria
idx = np.argmax(projected_signal)
val = projected_signal[idx]
else:
# If no loads exist within the swept area, set idx to nan, and value to 0
idx = np.nan
val = 0
else:
# If no sweep angle is defined, simply get the maximum value of the projected vector
idx = np.argmax(projected_signal)
val = projected_signal[idx]
if degrees:
angle = np.rad2deg(angle)
extremes[index, :] = (angle, val, idx)
return extremes
def compute_ensemble_2d_envelope(
channels: List[int], result_files: List[str]
) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
"""Calculate the ensemble envelope of an arbitrary number of 2D signals. The envelope is calculated as a convex hull.
Parameters
----------
channels : list[int]
List containing the two channels for calculation of the ensemble signal envelope across the provided result files.
result_files : list[str]
List of paths to HAWC2 result files. If the gtsdf format was not used as output format, the path to the *.sel files must be specified.
Returns
-------
ensemble_envelope : npt.NDArray[np.float64]
Coordinates of ensemble envelope. The shape of the array is (N, 2), where N is the number of points required to define the ensemble envelope.
individual_envelopes : list[npt.NDArray[np.float64]]
List of the arrays containing the coordinated of the individual signal envelopes calculated from the list of files provided.
Raises
------
ValueError
If the number of queried channels is not exaclty 2
ValueError
If the result files do not contain the queried channels
"""
# Assert that the number of channels is exactly 2.
try:
assert len(channels) == 2
except AssertionError:
raise ValueError("The list of channels must contain exactly two integers")
individual_envelopes = []
# Open result files in a loop, using the same variable to store the input to avoid excessive memory use
from wetb.hawc2.Hawc2io import ReadHawc2
for result_file in result_files:
res = ReadHawc2(f"{result_file}")
# Make sure the queried channels are available, otherwise raise error.
try:
assert res.NrCh >= (max(channels) - 1)
except AssertionError:
raise ValueError(
f"Result file {result_file} only has {res.NrCh} channels available. channels={channels} was provided."
)
res = res.ReadAll()[:, channels]
individual_envelopes.append(compute_envelope(res))
ensemble_env = compute_envelope(np.vstack(individual_envelopes))
return ensemble_env, individual_envelopes