In [3]:
import nbpresent
nbpresent.__version__
Out[3]:
In [4]:
import time
import numpy as np
import pandas as pd
import scipy.stats as ss
import matplotlib.pyplot as plt
from peak_finder import fetch_timetrace, global_peaks, weibull_fit, weib2gumb
%matplotlib inline
plt.rcParams['figure.figsize'] =(16.0, 6.0)
In [5]:
# autoreload modules after changes
%load_ext autoreload
%autoreload 2
In [17]:
# %% Load data
pth = r'I:\Stavanger\02 STVR BU ManCom and Dept\05 - Engineering\051 Analysis\Developments\RigidMethod\rayl_weib_gumb\01_sens\runs\Hs3.60_Tp08.00_WD75_seed1timetraces.txt'
ttraces = fetch_timetrace(pth)
# %% Also load the 30 seeds runs for the multi-seed Gumbel method
df = pd.read_table(r'I:\Stavanger\02 STVR BU ManCom and Dept\05 - Engineering\051 Analysis\Developments\RigidMethod\rayl_weib_gumb\01_sens\Results.txt')
In [18]:
def process(signal, times, extreme_peaks_gumbel):
# %% Detect peaks
sample, spl_times = global_peaks(signal, times)
sorted_sample = sorted(sample)
# %% Fit a Weibull distrubution to the peaks
# shape, loc, scale
wb_params_pwm = weibull_fit(sorted_sample)
wb = ss.frechet_r(*wb_params_pwm)
# Also calculate the Rayleigh MPM
MPM_Rayleigh = sum(sample)/len(sample) + np.array(sample).std() * (2*np.log(len(sample)))**0.5
# %% Calculate a Gumbel distribution for the extreme response
gb_params = weib2gumb(wb_params_pwm, len(sample))
gd = ss.gumbel_r(*gb_params)
# Fit a Gumbel distribution to the peaks from 30 seeds
params_gb30 = ss.gumbel_r.fit(extreme_peaks_gumbel)
gb30 = ss.gumbel_r(*params_gb30)
# %% Print statistics of peaks of 1 simulations
print('Weibull method - 1 simulation')
print('*'*30)
print('Signal mean', sum(signal)/len(signal))
print('Peaks', len(sample))
print('Min', min(sample))
print('Mean', sum(sample)/len(sample))
print('Max', max(sample))
print('Std', np.array(sample).std())
print('P10', np.percentile(sample, 10))
print('P50', np.percentile(sample, 50))
print('P90', np.percentile(sample, 90))
print('P99', np.percentile(sample, 99))
print('MPM Gumbel', gd.ppf(0.37))
print('MPM 1/N', wb.ppf(1-1/len(sample)))
print('MPM Rayleigh', MPM_Rayleigh)
print('Weibull parameters', wb_params_pwm)
print('Gumbel parameters', gb_params)
print('Weibull Percentile of Gumbel P37', 100*wb.cdf(gd.ppf(0.37)))
print()
print('Gumbel method - 30 simulations')
print('*'*30)
print('Min', min(extreme_peaks_gumbel))
print('Mean', sum(extreme_peaks_gumbel)/len(extreme_peaks_gumbel))
print('Max', max(extreme_peaks_gumbel))
print('Std', np.array(extreme_peaks_gumbel).std())
print('P10', np.percentile(extreme_peaks_gumbel, 10))
print('P50', np.percentile(extreme_peaks_gumbel, 50))
print('P90', np.percentile(extreme_peaks_gumbel, 90))
print('Gumbel P90', gb30.ppf(0.9))
print('MPM', gb30.ppf(0.37))
print('Parameters', params_gb30)
# Quality plots
# #############
# %% Gumbel extreme
xg = np.linspace(gd.ppf(0.0001), gd.ppf(0.999),100)
yg = gd.pdf(xg)
plt.subplot(221)
plt.plot(sorted_sample, wb.pdf(sorted_sample), label='Weibull of peaks')
plt.plot(xg, yg, label='Gumbel from Weibull')
plt.plot(xg, gb30.pdf(xg), label='Gumbel from 30 seeds')
plt.hist(sample, bins=20, normed=True, alpha=0.5)
plt.plot([wb.ppf(1-1/len(sample))]*2, [0, 0.045], label='Weibull MPM 1/N')
plt.plot([MPM_Rayleigh]*2, [0, 0.045], label='Rayleigh MPM')
plt.title('Peaks, Weibull and Extreme Gumbel')
plt.legend(loc='best')
plt.grid()
#plt.show()
#plt.close()
# %% Log-Log paperof Weibull
x = np.linspace(sorted_sample[0], sorted_sample[-1])
N = len(sample)
y = np.array([(i+1)/(N+1) for i in range(N)])
loc = wb_params_pwm[1]
plt.subplot(222)
plt.plot(np.log(sorted_sample-loc), np.log(np.log(1/(1-y))), '+')
plt.plot(np.log(x-loc), np.log(np.log(1/(1-wb.cdf(x)))))
plt.title('Weibull of Peaks - loglog x log scale')
plt.grid()
#plt.show()
#plt.close()
# %% Gumbel extremes
N = 30
extremes = sorted_sample[-N:]
# %% Timetrace, peaks, and extreme peaks
plt.subplot(223)
plt.plot(times, signal)
plt.plot(spl_times, sample, 'o')
signal_list = list(signal)
time_extremes = [times[signal_list.index(x)] for x in extremes]
plt.plot(time_extremes,extremes, 'sr')
plt.grid()
plt.title('Timetraces, peaks and extreme peaks')
plt.xlabel('Time')
plt.ylabel('Variate')
#plt.show()
#plt.close()
# Gumbel paper
x = np.linspace(extremes[0], extremes[-1])
# Not sure which empirical distribution to use ... here's 2 options:
y = np.array([(i-0.4+1)/(N+0.4) for i in range(N)])
# y = np.linspace(gd.cdf(extremes[0]), gd.cdf(extremes[-1]), N)
y30 = np.array([(i-0.4+1)/(30+0.4) for i in range(30)])
plt.subplot(224)
plt.plot(extremes, -np.log(-np.log(y)), '+', label='pks 1 sim')
plt.plot(x, -np.log(-gd.logcdf(x)), label='Wb->Gb')
plt.plot(sorted(extreme_peaks_gumbel), -np.log(-np.log(y30)), 'x', label='30 sims')
plt.plot(x, -np.log(-gb30.logcdf(x)), label='Gb 30')
plt.plot((extremes[0], extremes[-1]), -np.log(-np.log([0.37, 0.37])))
plt.grid(); plt.title('Gumble of extreme peaks'); plt.legend(loc='best');
plt.show()
#plt.close()
In [19]:
# Check what we've got
print('Time Traces')
for key in ttraces:
print(' ', key)
print('Results')
for cl in df.columns:
print(' ', cl)
In [20]:
list_vars = [
('Pipe DNV OS F101 Load Controlled EndA', 'Pipe Max DNV OS F101 Load Controlled EndA'),
('Pipe Effective Tension enda', 'Pipe Max Effective Tension EndA'),
('Pipe Bend Moment touchdown', 'Pipe Max Bend Moment TDP'),
]
def run(tt, evs):
print('\n\n', '*' * 80, '\n', tt, '\n', '*' * 80, '\n')
# Time trace for Weibull method
times = ttraces['Time']
signal = ttraces[tt]
# And extreme peaks from 30 seeds
extreme_peaks_gumbel = df[evs]
process(signal, times, extreme_peaks_gumbel)
In [21]:
run(*list_vars[0])
In [22]:
run(*list_vars[1])
In [23]:
run(*list_vars[2])
In [43]:
print(time.asctime())