NDBC functions convertsion

mhkit/loads/extreme/short_term_extreme.m

@@ -1,80 +1,183 @@
 function ste = short_term_extreme(t, data, t_st, type, x, method, options)
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% SHORT_TERM_EXTREME Estimate short-term extreme value statistics
 %
-%     Estimate the peaks distribution by fitting a Weibull
-%     distribution to the peaks of the response.
+%   ste = short_term_extreme(t, data, t_st, type, x, method)
 %
-%     To learn more about the methods that can be used in this function,
-%     refer to the scipy.stats.rv_continuous methods documentation!
+%   Estimates extreme value distributions from a timeseries of the response
+%   using either block maxima or peak-over-threshold methods.
 %
-%     Parameters
-%     ----------
-%         t : array
-%             Time array
-%         data : array
-%             Response timeseries
-%         t_st : double or int
-%             Short time period
-%         type : string
-%             Method for estimating the short-term extreme distribution
-%         x : array or double
-%             Input for the statistical function/method
-%         method : str
-%             Statistical method to apply to resulting data. Options to
-%             choose from are: "pdf", "cdf", "ppf", or "sf"
-%         output_py : bool (optional)
-%             Select if you want to return the native python result for use
-%             in long term extreme calculations.
+% Parameters
+% ----------
+%   t       : time vector
+%   data    : response signal
+%   t_st    : short-term period (scalar)
+%   type    : method string:
+%             'peaks_weibull', 'peaks_weibull_tail_fit',
+%             'peaks_over_threshold', 'block_maxima_gev', 'block_maxima_gumbel'
+%   x       : vector of values to evaluate the distribution at
+%   method  : 'pdf', 'cdf', 'ppf' (inverse CDF), or 'sf' (1 - CDF)
+%   options : struct with optional fields:
+%             .threshold (required for peaks_over_threshold)
 %
-%      Returns
-%      -------
-%         ste : array or double
-%             Probability distribution of the peaks
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Returns
+% -------
+%   ste     : values of the estimated distribution at x
 
 arguments
-    t
-    data
-    t_st
-    type
-    x
-    method
-    options.output_py=false;
+    t (:,1) double
+    data (:,1) double
+    t_st (1,1) double
+    type (1,:) char
+    x (:,1) double
+    method (1,:) char
+    options.threshold double = NaN
+end
+
+% Error checking
+if length(t) ~= length(data)
+    error('Time and data must be the same length');
 end
 
+% Sampling interval and frequency
+dt = mean(diff(t));
+fs = 1/dt;
+
+switch lower(type)
+    case 'peaks_weibull'
+        peaks = find_peaks(data);
+        nst = round(length(peaks) * t_st / duration);
+        pd_st = fit_weibull(peaks, nst, x);
+
+    case 'peaks_weibull_tail_fit'
+        peaks = find_peaks(data);
+        pd_st = fit_weibull_tail(peaks, x);
+
+    case 'peaks_over_threshold'
+        if isnan(options.threshold)
+            error('You must provide options.threshold for this method.');
+        end
+        peaks = data(data > options.threshold);
+        pd_st = fit_pareto(peaks, x);
 
-if ~isa(t,'numeric')
-    error('ERROR: t must be a double array')
+    case 'block_maxima_gev'
+        block_size = round(t_st * fs);
+        blocks = reshape_truncated(data, block_size);
+        maxima = max(blocks, [], 1);
+        pd_st = fit_gumbel(maxima, x);
+
+    case 'block_maxima_gumbel'
+        block_size = round(t_st * fs);
+        blocks = reshape_truncated(data, block_size);
+        maxima = max(blocks, [], 1);
+        pd_st = fit_gev(maxima, x);
+
+    otherwise
+        error('Unknown method: %s', type);
 end
-if ~isa(data,'numeric')
-    error('ERROR: data must be a double array')
+
+switch lower(method)
+    case 'pdf'
+        ste = pd_st.pdf;
+    case 'cdf'
+        ste = pd_st.cdf;
+    case {'ppf', 'inv'}
+        ste = pd_st.ppf;
+    case 'sf'
+        ste = 1 - pd_st.cdf;
+    case 'expect'
+        % Numerical expectation: integral of x * pdf(x)
+        dx = mean(diff(x));
+        ste = sum(x .* pd_st.pdf) * dx;
+    otherwise
+        error('Invalid method: %s. Choose from "pdf", "cdf", "ppf", "sf", "expect".', method);
 end
-if ~isa(t_st,'numeric')
-    error('ERROR: t_st must be a double')
 end
-if ~isa(type,'string')
-    error('ERROR: type must be a string')
+
+
+function peaks = find_peaks(data)
+    locs = find(data(2:end-1) > data(1:end-2) & data(2:end-1) > data(3:end)) + 1;
+    peaks = data(locs);
 end
-if ~isa(x,'numeric')
-    error('ERROR: x must be a double array')
+
+function blocks = reshape_truncated(data, block_size)
+    N = floor(length(data)/block_size) * block_size;
+    blocks = reshape(data(1:N), block_size, []);
 end
-if ~isa(method,'string')
-    error('ERROR: method must be a string')
+
+function pd = fit_weibull_tail(data, x)
+    % Manual peak filtering
+    locs = find(data(2:end-1) > data(1:end-2) & data(2:end-1) > data(3:end)) + 1;
+    peaks = data(locs);
+
+    % Tail selection
+    threshold = prctile(peaks, 95);
+    tail_peaks = peaks(peaks > threshold);
+
+    % Grid search over Weibull parameters
+    scale_vals = linspace(mean(tail_peaks)*0.5, mean(tail_peaks)*2, 50);
+    shape_vals = linspace(0.5, 5, 50);
+
+    min_nll = inf;
+
+    for a = scale_vals
+        for b = shape_vals
+            pdf_vals = weibull_pdf(tail_peaks, a, b);
+            if all(pdf_vals > 0 & isfinite(pdf_vals))
+                nll = -sum(log(pdf_vals));
+                if nll < min_nll
+                    min_nll = nll;
+                    best_scale = a;
+                    best_shape = b;
+                end
+            end
+        end
+    end
+
+    pd.type = 'Weibull';
+    pd.scale = best_scale;
+    pd.shape = best_shape;
+    pd.pdf = weibull_pdf(x, best_scale, best_shape);
+    pd.cdf = weibull_cdf(x, best_scale, best_shape);
+    pd.ppf = weibull_ppf(x, best_scale, best_shape);
 end
 
-py.importlib.import_module('mhkit');
+function pd = fit_weibull(data, nst, x)
+    % Same fitting approach as the tail, just on all peaks
+    scale_vals = linspace(mean(data)*0.5, mean(data)*2, 50);
+    shape_vals = linspace(0.5, 5, 50);
 
-t = py.numpy.array(t);
-data = py.numpy.array(data);
+    min_nll = inf;
 
-result = py.mhkit.loads.extreme.short_term_extreme(t, data, t_st, type);
+    for a = scale_vals
+        for b = shape_vals
+            pdf_vals = weibull_pdf(data, a, b);
+            if all(pdf_vals > 0 & isfinite(pdf_vals))
+                nll = -sum(log(pdf_vals));
+                if nll < min_nll
+                    min_nll = nll;
+                    best_scale = a;
+                    best_shape = b;
+                end
+            end
+        end
+    end
 
-if options.output_py
-    ste = result;
-else
-    stat = py.mhkit_python_utils.scipy_stats.convert_to_array(result, method, x);
-    ste = double(stat);
+    pd.type = 'Weibull';
+    pd.scale = best_scale;
+    pd.shape = best_shape;
+    pd.pdf = weibull_pdf(x, best_scale, best_shape);
+    pd.cdf = weibull_cdf(x, best_scale, best_shape).^nst;  % scaled for ST realization
+    pd.ppf = weibull_ppf(x, best_scale, best_shape);
 end
 
-end
+function p = weibull_pdf(x, scale, shape)
+    p = (shape./scale) .* (x./scale).^(shape - 1) .* exp(-(x./scale).^shape);
+end
+
+function c = weibull_cdf(x, scale, shape)
+    c = 1 - exp(-(x./scale).^shape);
+end
+
+function x_ppf = weibull_ppf(p, scale, shape)
+    x_ppf = scale .* (-log(1 - p)).^(1/shape);
+end

mhkit/wave/resource/energy_period.m

@@ -1,49 +1,74 @@
-function Te=energy_period(S,varargin)
-
+function Te = energy_period(S, varargin)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
+% Calculates wave energy period Te from wave spectra
 %
 % Parameters
 % ------------
-%   S: Spectral Density (m^2/Hz)
-%       Pandas data frame
-%           To make a pandas data frame from user supplied frequency and spectra
-%           use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(frequency,spectra)
-%
-%       OR
-%
-%       structure of form:
-%           S.spectrum: Spectral Density (m^2/Hz)
-%
-%           S.type: String of the spectra type, i.e. Bretschneider,
-%           time series, date stamp etc.
-%
-%           S.frequency: frequency (Hz)
-%
-%     frequency_bins: vector (optional)
-%       Bin widths for frequency of S. Required for unevenly sized bins
+%   S: structure or numeric array
+%       If structure:
+%           S.spectrum   - Spectral Density (m^2/Hz)
+%           S.frequency  - Frequency (Hz)
+%       If numeric:
+%           S is spectral density array (vector or matrix)
+%           varargin{1} must contain frequency vector
 %
+%   frequency_bins: vector (optional)
+%       Frequency bin widths [Hz]. Required for unevenly spaced bins.
 %
 % Returns
 % ---------
-%    Te: float
-%        Wave energy Period (s)
-%
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%   Te: double
+%       Wave energy period [s]
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-% assign feq_bin
-if nargin == 2
-    freq_bins = py.numpy.array(varargin{1});
-elseif nargin == 1
-    freq_bins = py.None;
-else
-    ME = MException('MATLAB:energy_period','Incorrect number of input arguments');
-        throw(ME);
-end
+    % Extract spectrum and frequency
+    if isstruct(S)
+        spectrum = S.spectrum;
+        frequency = S.frequency;
+    elseif isnumeric(S)
+        if nargin < 2
+            error('When S is numeric, frequency vector must be provided as second argument');
+        end
+        spectrum = S;
+        frequency = varargin{1};
+        varargin(1) = [];
+    else
+        error('Input S must be a struct or numeric array');
+    end
+
+    % Determine frequency bin widths
+    if ~isempty(varargin)
+        freq_bins = varargin{1};
+    else
+        df = diff(frequency);
+        freq_bins = mean(df);
+    end
 
-S_py = typecast_spectra_to_mhkit_python(S);
+    % Ensure proper shapes
+    frequency = frequency(:);
+    if isvector(spectrum)
+        spectrum = spectrum(:);
+    end
 
-Te = py.mhkit.wave.resource.energy_period(S_py, pyargs('frequency_bins',freq_bins));
+    % Check size consistency
+    if length(frequency) ~= size(spectrum,1)
+        error('Length of frequency vector must match number of rows in spectrum');
+    end
 
-Te = typecast_from_mhkit_python(Te).data;
+    % Calculate moments: m0 and m-1
+    if isscalar(freq_bins)
+        m0 = sum(spectrum .* freq_bins, 1);
+        m_neg1 = sum((spectrum ./ frequency) .* freq_bins, 1);
+    else
+        freq_bins = freq_bins(:);
+        if length(freq_bins) ~= length(frequency)
+            error('Length of freq_bins must match frequency vector');
+        end
+        m0 = sum(spectrum .* freq_bins, 1);
+        m_neg1 = sum((spectrum ./ frequency) .* freq_bins, 1);
+    end
+
+    % Calculate energy period
+    Te = m_neg1 ./ m0;
+
+end