Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix: Restore changes from PR #127 overwritten by PR #138 #154

Merged
merged 1 commit into from
Dec 19, 2024
Merged

Fix: Restore changes from PR #127 overwritten by PR #138 #154

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 101 additions & 14 deletions mhkit/tests/Loads_TestExtreme.m
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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


Loading