Skip to content

Commit

Permalink
Fix: Restore changes from PR #127 overwritten by PR #138
Browse files Browse the repository at this point in the history 
PR #138's merge inadvertently overwrote WDRT module updates from PR#127
due to an incorrect merge conflict resolution here: simmsa@a0d07da
  • Loading branch information
@simmsa
simmsa committed Dec 19, 2024
1 parent caf8263 commit 3cf9f34
Show file tree
Hide file tree
Showing 3 changed files with 142 additions and 54 deletions.
115 changes: 101 additions & 14 deletions mhkit/tests/Loads_TestExtreme.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,13 +16,11 @@ function test_mler_coefficients(testCase)
fpath = '../../examples/data/loads/mler.csv';
validation = readtable(fpath);
wave_freq = linspace(0,1,500);
% js = jonswap_spectrum(wave_freq, 15.1, 9);
this_wave = pierson_moskowitz_spectrum(wave_freq, 15.1, 9);
js = jonswap_spectrum(wave_freq,15.1,9);
response_desired = 1;
RAO = validation.RAO;
% execute function
% mler = mler_coefficients(RAO, js, response_desired);
mler = mler_coefficients(RAO, this_wave, response_desired);
mler = mler_coefficients(RAO, js, response_desired);
% assertions
assertEqual(testCase, mler.conditioned_spectrum, validation.Res_Spec, 'RelTol',0.005)
assertEqual(testCase, mler.phase, validation.phase, 'RelTol',0.001)
Expand All @@ -33,27 +31,27 @@ function test_mler_simulation(testCase)
T = linspace(-150,150,301);
X = linspace(-300,300,601);
sim = mler_simulation();

assertEqual(testCase, sim.T, T)
assertEqual(testCase, sim.X, X)
end

function test_mler_wave_amp_normalize(testCase)
fpath = '../../examples/data/loads/mler.csv';
validation = readtable(fpath);
mler.conditioned_spectrum = validation.Res_Spec;
mler.phase = validation.phase;
wave_freq = linspace(0,1,500);
mler.frequency = wave_freq';

k = wave_number(wave_freq, 70);
k.values = fillmissing(k.values,'constant',0);
sim = mler_simulation();
mler_norm = mler_wave_amp_normalize(4.5*1.9, mler, sim, k.values);

assertEqual(testCase, mler_norm.conditioned_spectrum, validation.Norm_Spec, 'AbsTol', 0.002)
end

function test_mler_export_time_series(testCase)
fpath = '../../examples/data/loads/mler_ts.csv';
validation = readtable(fpath);
Expand All @@ -65,13 +63,102 @@ function test_mler_export_time_series(testCase)
mler.frequency = wave_freq';
RAO = normed.RAO;
k = wave_number(wave_freq, 70);
k.values = fillmissing(k.values,'constant',0);
k.values = fillmissing(k.values,'constant',0);
sim = mler_simulation();
mler_ts = mler_export_time_series(RAO, mler, sim, k.values);

assertEqual(testCase, mler_ts.linear_response, validation.LinearResponse, 'AbsTol', 0.00005)
end
end

function test_global_peaks(testCase)
fpath = '../../examples/data/loads/sine_wave.csv';
validation = readtable(fpath);
Fs = 200; % samples per second
dt = 1/Fs; % seconds per sample
StopTime = 0.25; % seconds
t = (0:dt:StopTime-dt)'; % seconds
Fc = 60; % hertz
x = cos(2*pi*Fc*t);
[t_peaks, peaks] = global_peaks(t, x);

end
assertEqual(testCase, t_peaks, validation.t_peaks', 'AbsTol', 0.00005)
assertEqual(testCase, peaks, validation.peaks', 'AbsTol', 0.00005)
end

function test_number_of_short_term_peaks(testCase)
n = 10;
t = 100;
t_st = 50;
n_st = number_of_short_term_peaks(n, t, t_st);

assertEqual(testCase, n_st, 5)
end

function test_block_maxima(testCase)
load('../../examples/data/loads/block_maxima.mat', 'block_max');
load('../../examples/data/loads/example_qoi.mat', 'qoi')
% time in seconds
time = linspace(0, 10800, 2*10801+1);
% split time series into 10-minute (600s) segments
bm = block_maxima(time, qoi', 600);

assertEqual(testCase, bm, block_max)
end

function test_ste_block_maxima_gev(testCase)
load('../../examples/data/loads/block_maxima.mat', 'block_max');
load('../../examples/data/loads/ste.mat', 'ste_gev_pdf')
load('../../examples/data/loads/ste.mat', 'ste_gev_cdf')
x = linspace(0,3,1000);
stegev_pdf = ste_block_maxima_gev(block_max, x, "pdf");
stegev_cdf = ste_block_maxima_gev(block_max, x, "cdf");

assertEqual(testCase, stegev_pdf, ste_gev_pdf, 'AbsTol', 0.00005)
assertEqual(testCase, stegev_cdf, ste_gev_cdf, 'AbsTol', 0.00005)
end

function test_short_term_extremes(testCase)
txtfile = readtable("../../examples/data/loads/time_series_for_extremes.txt");
txtfile = table2array(txtfile);
t_st = 1 * 60 * 60;
x = 1.6;
cdfs_1 = [
0.006750456316537166;
0.5921659393757381;
0.6156789503874247;
0.6075807789811315;
0.9033574618279865];
methods = [
"peaks_weibull";
"peaks_weibull_tail_fit";
"peaks_over_threshold";
"block_maxima_gev";
"block_maxima_gumbel"];
cdfs = zeros(1,length(methods))';
for i = 1:length(methods)
cdfs(i) = short_term_extreme(txtfile(:,1), txtfile(:,2), t_st, methods(i), x, "cdf");
end

assertEqual(testCase, cdfs, cdfs_1, 'AbsTol', 0.00005)
end

function test_long_term_extreme(testCase)
txtfile = readtable("../../examples/data/loads/time_series_for_extremes.txt");
txtfile = table2array(txtfile);
t_st = 1 * 60 * 60;
x = linspace(0,10,1000);
z = 0.8;
ste1 = short_term_extreme(txtfile(:,1), txtfile(:,2), t_st, "peaks_weibull", z, "cdf");
ste{1} = short_term_extreme(txtfile(:,1), txtfile(:,2), t_st, "peaks_weibull", x, "cdf", output_py=true);
ste2 = short_term_extreme(txtfile(:,1), txtfile(:,2), t_st, "peaks_weibull", z, "cdf");
ste{2} = short_term_extreme(txtfile(:,1), txtfile(:,2), t_st, "peaks_weibull", x, "cdf", output_py=true);
pyste = py.list(ste);
w = [0.25, 0.75];

lte_cdf = full_seastate_long_term_extreme(pyste, w, z, "cdf");

end
assertEqual(testCase, lte_cdf, w(1)*ste1 + w(2)*ste2, 'AbsTol', 0.00005)
end

end
end
20 changes: 10 additions & 10 deletions mhkit/wave/graphics/plot_environmental_contours.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots an overlay of the x1 and x2 variables to the calculated
% environmental contours.
%
%
% Parameters
% ------------
% x1: vector
Expand All @@ -27,12 +27,12 @@
% to call: plot_environmental_contours(x1,x2,x1_contour,x2_contour,"y_label",y_label)
%
% data_label: string (optional)
% Legend label for x1, x2 data (e.g. 'Buoy 46022').
% Legend label for x1, x2 data (e.g. 'Buoy 46022').
% Default None.
% to call: plot_environmental_contours(x1,x2,x1_contour,x2_contour,"data_label",data_label)
%
% contour_label: string or array of strings (optional)
% Legend label for x1_contour, x2_contour countor data
% Legend label for x1_contour, x2_contour countor data
% (e.g. '100-year contour'). Default None.
% to call: plot_environmental_contours(x1,x2,x1_contour,x2_contour,"contour_label",contour_label)
%
Expand All @@ -43,12 +43,12 @@
% savepath: string (optional)
% path and filename to save figure.
% to call: plot_environmental_contours(x1,x2,x1_contour,x2_contour,"savepath",savepath)
%
%
% Returns
% ---------
% figure: figure
% Envoronmental contour plot
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

arguments
Expand All @@ -71,18 +71,18 @@
end
grid on
xlabel(options.x_label)
ylabel(options.y_label)
ylabel(options.y_label)

title(options.title)
len = strlength(options.data_label);
len2 = strlength(options.contour_label);
if len > 1 || len2 > 1
label = [ {options.data_label} , options.contour_label];
legend(label);
label = [ {options.data_label} , options.contour_label];
legend(label,'Location','bestoutside');
end
len = strlength(options.savepath);
if len > 1
saveas(figure, options.savepath);
end
end

hold off
hold off
61 changes: 31 additions & 30 deletions mhkit/wave/resource/environmental_contours.m
View file
Original file line number Diff line number Diff line change
@@ -1,49 +1,49 @@
function environmental_contour=environmental_contours(x1, x2, dt, period, method, options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates environmental contours of extreme sea
% states using the improved joint probability distributions
% with the inverse first-order reliability method (IFORM)
% probability for the desired return period (period). Given the
% Calculates environmental contours of extreme sea
% states using the improved joint probability distributions
% with the inverse first-order reliability method (IFORM)
% probability for the desired return period (period). Given the
% period of interest a circle of iso-probability is created in the
% in the PCA joint probability (x1, x2) reference frame.
% Using the joint probability value the CDF of the marginal
% distribution is used to find the quantile of each component.
% in the PCA joint probability (x1, x2) reference frame.
% Using the joint probability value the CDF of the marginal
% distribution is used to find the quantile of each component.
% Finally, using the improved PCA methodology
% the component 2 contour lines are calculated from component 1 using
% the relationships defined in Exkert-Gallup et. al. 2016.
%
% Eckert-Gallup, A. C., Sallaberry, C. J., Dallman, A. R., &
% Neary, V. S. (2016). Application of principal component
% analysis (PCA) and improved joint probability distributions to
% the inverse first-order reliability method (I-FORM) for predicting
% extreme sea states. Ocean Engineering, 112, 307-319.
%
% the component 2 contour lines are calculated from component 1 using
% the relationships defined in Exkert-Gallup et. al. 2016.
%
% Eckert-Gallup, A. C., Sallaberry, C. J., Dallman, A. R., &
% Neary, V. S. (2016). Application of principal component
% analysis (PCA) and improved joint probability distributions to
% the inverse first-order reliability method (I-FORM) for predicting
% extreme sea states. Ocean Engineering, 112, 307-319.
%
% Parameters
% ------------
% x1 : vector
% component 1 data
%
% x2 : vector
% component 2 data
% component 2 data
%
% dt : double
% x1 and x2 sample rate (seconds)
%
%
% period : scalar or vector
% Desired return period (years) for calculation of environmental
% contour, can be a scalar or a vector.
%
%
% PCA: Structure (optional)
% principal component analysis dictionary from previous function
% call. When supplied the function will skip the PCA calculation
% for the passe x1, and x2.
% principal component analysis dictionary from previous function
% call. When supplied the function will skip the PCA calculation
% for the passe x1, and x2.
% to call: environmental_contour(x1,x2,dt,period,"PCA",PCA)
%
%
% bin_size : double (optional)
% Data points in each bin
% Data points in each bin
% to call: environmental_contour(x1,x2,dt,period,"bin_size",bin_size)
%
%
% nb_steps : int (optional)
% Discretization of the circle in the normal space used for
% IFORM calculation.
Expand All @@ -52,17 +52,17 @@
% return_PCA: boolean
% Default False, if True will retun the PCA dictionary
% to call: environmental_contour(x1,x2,dt,period,"return_PCA",return_PCA)
%
%
%
%
% Returns
% ---------
% environmental_contour: Structure
% environmental_contour: Structure
% Structure with fields contour1, contour2, and optionally PCA
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

arguments
arguments
x1
x2
dt
Expand All @@ -83,7 +83,7 @@
if isscalar(period)
period_py = period;
elseif isvector(period)
period_py = py.numpy.array(period);
period_py = py.numpy.array(period);
else
ME = MException('MATLAB:environmental_contour','period must be a vector or scalar');
throw(ME);
Expand All @@ -106,3 +106,4 @@
environmental_contour.fit = struct(data_struct.(varfit));
end


0 comments on commit 3cf9f34

Please sign in to comment.