diff --git a/mhkit/tests/Loads_TestExtreme.m b/mhkit/tests/Loads_TestExtreme.m index ce2575da..588eeceb 100644 --- a/mhkit/tests/Loads_TestExtreme.m +++ b/mhkit/tests/Loads_TestExtreme.m @@ -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) @@ -33,11 +31,11 @@ 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); @@ -45,15 +43,15 @@ function test_mler_wave_amp_normalize(testCase) 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); @@ -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 \ No newline at end of file diff --git a/mhkit/wave/graphics/plot_environmental_contours.m b/mhkit/wave/graphics/plot_environmental_contours.m index 21445258..b07014e9 100644 --- a/mhkit/wave/graphics/plot_environmental_contours.m +++ b/mhkit/wave/graphics/plot_environmental_contours.m @@ -3,7 +3,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Plots an overlay of the x1 and x2 variables to the calculated % environmental contours. -% +% % Parameters % ------------ % x1: vector @@ -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) % @@ -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 @@ -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 \ No newline at end of file diff --git a/mhkit/wave/resource/environmental_contours.m b/mhkit/wave/resource/environmental_contours.m index 779c3118..209fc1d7 100644 --- a/mhkit/wave/resource/environmental_contours.m +++ b/mhkit/wave/resource/environmental_contours.m @@ -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. @@ -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 @@ -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); @@ -106,3 +106,4 @@ environmental_contour.fit = struct(data_struct.(varfit)); end +