From 5bf698d1f4aec01c4dfacd4d25814dc2c0770dc9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carlos=20A=2E=20Michel=C3=A9n=20Str=C3=B6fer?= Date: Wed, 22 Nov 2023 11:05:29 -0700 Subject: [PATCH] reset --- wecopttool/core.py | 261 +++++++++++++++++++++------------------------ 1 file changed, 119 insertions(+), 142 deletions(-) diff --git a/wecopttool/core.py b/wecopttool/core.py index 86302314..55b49616 100644 --- a/wecopttool/core.py +++ b/wecopttool/core.py @@ -24,7 +24,7 @@ "derivative2_mat", "mimo_transfer_mat", "vec_to_dofmat", - "dofmat_to_vec", ## + "dofmat_to_vec", "real_to_complex", "complex_to_real", "fd_to_td", @@ -43,14 +43,14 @@ "add_linear_friction", "wave_excitation", "hydrodynamic_impedance", - "atleast_2d", ## + "atleast_2d", "degrees_to_radians", "subset_close", "scale_dofs", "decompose_state", "frequency_parameters", "time_results", - "set_fb_centers", # + "set_fb_centers", ] @@ -363,8 +363,8 @@ def from_bem( inertia_matrix = hydro_data['inertia_matrix'].values # frequency array - f1, nfreq, nfreq_start = frequency_parameters( - hydro_data.omega.values/(2*np.pi)) + f1, nfreq = frequency_parameters( + hydro_data.omega.values/(2*np.pi), False) # check real part of damping diagonal > 0 if min_damping is not None: @@ -532,7 +532,7 @@ def from_impedance( If :python:`impedance` does not have the correct size: :python:`(ndof, ndof, nfreq)`. """ - f1, nfreq, nfreq_start = frequency_parameters(freqs) + f1, nfreq = frequency_parameters(freqs, False) # impedance matrix shape shape = impedance.shape @@ -1269,61 +1269,43 @@ def td_to_fd( def ncomponents( nfreq : int, - nfreq_start: Optional[int] = 1, zero_freq: Optional[bool] = True, - total: Optional[bool] = False, ) -> int: - """Number of Fourier components for each DOF. + """Number of Fourier components (:python:`2*nfreq`) for each + DOF. The sine component of the highest frequency (the 2-point wave) + is excluded as it will always evaluate to zero. - For most frequencies this is two compenents per frequency. - However, the sine component of the highest frequency - (the 2-point wave) is excluded as it will always evaluate to zero. - Similarly, the zero-frequency (DC) component, if included via - :python:`zero_freq=True` (default) , only has a real part. - - With the DC component included the total number of components is - :python:`2*nfreq`, corresponding to :math:`x = [X0, Re(X1), Im(x1), ..., Re(Xn)]` - were :math:`X0` corresponds to frequency :python:`0`, - :math:`X1` corresponds to frequency :python:`nfreq_start*f1` and - :math:`Xn` corresponds to frequency :python:`(nfreq_start+nfreq-1)*f1`. - - Without the DC component the totral number of components is one - less, or :python:`2*nfreq-1` + If :python:`zero_freq = False` (not default), the mean (DC) component + :python:`X0` is excluded, and the number of components is reduced by 1. Parameters ---------- nfreq Number of frequencies. - nfreq_start - Frequency index at which frequency vector starts. - total - Whether to return the total number of frequencies as if it - started at :python:`nfreq_start=1`. zero_freq - Whether to include the zero (DC) component. + Whether to include the zero-frequency. """ - nfreq = nfreq if (not total) else (nfreq+nfreq_start-1) - ncomp = 2*nfreq if zero_freq else 2*nfreq-1 + ncomp = 2*nfreq - 1 + if zero_freq: + ncomp = ncomp + 1 return ncomp def frequency( f1: float, nfreq: int, - nfreq_start: Optional[int] = 1, zero_freq: Optional[bool] = True, ) -> ndarray: """Construct equally spaced frequency array. - Returns the frequency array - :python:`[0, f1*nfreq_start, f1*(nfreq_start+1), ..., f1*(nfreq_start+nfreq-1)]`. - :python:`f1` is fundamental frequency (1st harmonic) and - :python:`n_start` allows the vector starting at a frequency other - than :python:`f1`. - Larger values of :python:`nfreq_start` can be useful in reducing the - computational cost when there is little or no energy in the lower - frequency components, specially when using a very fine - discretization (small :python:`f1`). + The array includes :python:`0` and has length of :python:`nfreq+1`. + :python:`f1` is fundamental frequency (1st harmonic). + + Returns the frequency array, e.g., + :python:`freqs = [0, f1, 2*f1, ..., nfreq*f1]`. + + If :python:`zero_freq = False` (not default), the mean (DC) component + :python:`0` is excluded, and the vector length is reduced by 1. Parameters ---------- @@ -1331,8 +1313,6 @@ def frequency( Fundamental frequency :python:`f1` [:math:`Hz`]. nfreq Number of frequencies. - nfreq_start - Frequency index at which frequency vector starts. zero_freq Whether to include the zero-frequency. """ @@ -1344,16 +1324,15 @@ def frequency( def time( f1: float, nfreq: int, - nfreq_start: Optional[int] = 1, nsubsteps: Optional[int] = 1, ) -> ndarray: """Assemble the time vector with :python:`nsubsteps` subdivisions. Returns the 1D time vector, in seconds, starting at time :python:`0`, and not containing the end time :python:`tf=1/f1`. - The time vector has length :python:`2*(nfreq+nfreq_start-1)*nsubsteps`. + The time vector has length :python:`(2*nfreq)*nsubsteps`. The timestep length is :python:`dt = dt_default * 1/nsubsteps`, - where :python:`dt_default=tf/(2*(nfreq+nfreq_start-1))`. + where :python:`dt_default=tf/(2*nfreq)`. Parameters ---------- @@ -1361,21 +1340,19 @@ def time( Fundamental frequency :python:`f1` [:math:`Hz`]. nfreq Number of frequencies. - nfreq_start - Frequency index at which frequency vector starts. nsubsteps Number of steps between the default (implied) time steps. A value of :python:`1` corresponds to the default step length. """ - ncomp = ncomponents(nfreq, nfreq_start, total=True, zero_freq=True) - nsteps = nsubsteps * ncomp + if nsubsteps < 1: + raise ValueError("'nsubsteps' must be 1 or greater") + nsteps = nsubsteps * ncomponents(nfreq) return np.linspace(0, 1/f1, nsteps, endpoint=False) def time_mat( f1: float, nfreq: int, - nfreq_start: Optional[int] = 1, nsubsteps: Optional[int] = 1, zero_freq: Optional[bool] = True, ) -> ndarray: @@ -1383,16 +1360,16 @@ def time_mat( time-series. For a state :math:`x` consisting of the mean (DC) component - (if :python:`nfreq_start=0`) followed by the real and imaginary - components of the Fourier coefficients (excluding the imaginary - component of the 2-point wave) as + followed by the real and imaginary components of the Fourier + coefficients (excluding the imaginary component of the 2-point wave) as :math:`x=[X0, Re(X1), Im(X1), ..., Re(Xn)]`, the response vector in the time-domain (:math:`x(t)`) is given as :math:`Mx`, where :math:`M` is the time matrix. - If :python:`nfreq_start>0` the state vector starts at components for - frequency :python:`nfreq_start*f1` instead of :python:`f1`. - If :python:`zero_freq=False` the DC component is assumed missing - from the state vector. + + The time matrix has size :python:`(nfreq*2, nfreq*2)`. + + If :python:`zero_freq = False` (not default), the mean (DC) component + :python:`X0` is excluded, and the matrix/vector length is reduced by 1. Parameters --------- @@ -1400,33 +1377,28 @@ def time_mat( Fundamental frequency :python:`f1` [:math:`Hz`]. nfreq Number of frequencies. - nfreq_start - Frequency index at which the frequency vector starts. nsubsteps Number of steps between the default (implied) time steps. A value of :python:`1` corresponds to the default step length. zero_freq - Whether the state vector includes the zero (DC) component. + Whether the first frequency should be zero. """ - t = time(f1, nfreq, nfreq_start, nsubsteps) - omega = frequency(f1, nfreq, nfreq_start, zero_freq=False) * 2*np.pi - ncomp = ncomponents(nfreq, nfreq_start, zero_freq) - time_mat = np.empty((nsubsteps*len(t), ncomp)) - wt = np.outer(t, omega) - if zero_freq: - time_mat[:, 0] = 1.0 - c0, s0 = 1,2 - else: - c0,s0 = 0,1 - time_mat[:, c0::2] = np.cos(wt) - time_mat[:, s0::2] = -np.sin(wt[:, :-1]) # remove 2pt wave sine component + t = time(f1, nfreq, nsubsteps) + omega = frequency(f1, nfreq) * 2*np.pi + wt = np.outer(t, omega[1:]) + ncomp = ncomponents(nfreq) + time_mat = np.empty((nsubsteps*ncomp, ncomp)) + time_mat[:, 0] = 1.0 + time_mat[:, 1::2] = np.cos(wt) + time_mat[:, 2::2] = -np.sin(wt[:, :-1]) # remove 2pt wave sine component + if not zero_freq: + time_mat = time_mat[:, 1:] return time_mat def derivative_mat( f1: float, nfreq: int, - nfreq_start: Optional[int] = 1, zero_freq: Optional[bool] = True, ) -> ndarray: """Assemble the derivative matrix that converts the state vector of @@ -1434,13 +1406,15 @@ def derivative_mat( For a state :math:`x` consisting of the mean (DC) component followed by the real and imaginary components of the Fourier - coefficients (excluding the imaginary component of the 2-point wave) - as :math:`x=[X0, Re(X1), Im(X1), ..., Re(Xn)]`, the state of its - derivative is given as :math:`Dx`, where :math:`D` is the derivative - matrix. - If :python:`zero_freq=False` the state vector does not include the - DC component and starts at components for frequency - :python:`nfreq_start*f1`. + coefficients (excluding the imaginary component of the 2-point wave) as + :math:`x=[X0, Re(X1), Im(X1), ..., Re(Xn)]`, + the state of its derivative is given as :math:`Dx`, where + :math:`D` is the derivative matrix. + + The time matrix has size :python:`(nfreq*2, nfreq*2)`. + + If :python:`zero_freq = False` (not default), the mean (DC) component + :python:`X0` is excluded, and the matrix/vector length is reduced by 1. Parameters --------- @@ -1448,37 +1422,36 @@ def derivative_mat( Fundamental frequency :python:`f1` [:math:`Hz`]. nfreq Number of frequencies. - nfreq_start - Frequency index at which frequency vector starts. zero_freq - Whether the state vector includes the zero (DC) component. + Whether the first frequency should be zero. """ def block(n): return np.array([[0, -1], [1, 0]]) * n*f1 * 2*np.pi - blocks = [block(n) for n in range(nfreq_start, nfreq+nfreq_start)] + blocks = [block(n+1) for n in range(nfreq)] if zero_freq: blocks = [0.0] + blocks - deriv_mat = block_diag(*blocks)[:-1, :-1] # remove 2pt wave sine component - return deriv_mat + deriv_mat = block_diag(*blocks) + return deriv_mat[:-1, :-1] # remove 2pt wave sine component def derivative2_mat( f1: float, nfreq: int, - nfreq_start: Optional[int] = 0, zero_freq: Optional[bool] = True, ) -> ndarray: - """Assemble the second derivative matrix that converts the state - vector of a response to the state vector of its second derivative. + """Assemble the second derivative matrix that converts the state vector of + a response to the state vector of its second derivative. For a state :math:`x` consisting of the mean (DC) component followed by the real and imaginary components of the Fourier - coefficients (excluding the imaginary component of the 2-point wave) - as :math:`x=[X0, Re(X1), Im(X1), ..., Re(Xn)]`, the state of its - second derivative is given as :math:`(DD)x`, where :math:`DD` is the - second derivative matrix. - If :python:`zero_freq=False` the state vector does not include the - DC component and starts at components for frequency - :python:`nfreq_start*f1`. + coefficients (excluding the imaginary component of the 2-point wave) as + :math:`x=[X0, Re(X1), Im(X1), ..., Re(Xn)]`, + the state of its second derivative is given as :math:`(DD)x`, where + :math:`DD` is the second derivative matrix. + + The time matrix has size :python:`(nfreq*2, nfreq*2)`. + + If :python:`zero_freq = False` (not default), the mean (DC) component + :python:`X0` is excluded, and the matrix/vector length is reduced by 1. Parameters --------- @@ -1486,13 +1459,11 @@ def derivative2_mat( Fundamental frequency :python:`f1` [:math:`Hz`]. nfreq Number of frequencies. - nfreq_start - Frequency index at which frequency vector starts. zero_freq - Whether the state vector includes the zero (DC) component. + Whether the first frequency should be zero. """ - vals = [((n)*f1 * 2*np.pi)**2 for n in range(nfreq_start, nfreq+nfreq_start)] - diagonal = np.squeeze(np.repeat(-np.ones(nfreq-1) * vals, 2))[:-1] # remove 2pt wave sine + vals = [((n+1)*f1 * 2*np.pi)**2 for n in range(nfreq)] + diagonal = np.repeat(-np.ones(nfreq) * vals, 2)[:-1] # remove 2pt wave sine if zero_freq: diagonal = np.concatenate(([0.0], diagonal)) return np.diag(diagonal) @@ -1512,23 +1483,23 @@ def mimo_transfer_mat( :python`(nfreq, ndof, ndof)`. Returns the 2D real matrix that transform the state representation - of the input variables to the state representation of the output - variables. + of the input variable variable to the state representation of the + output variable. Here, a state representation :python:`x` consists of the mean (DC) - component (if included) followed by the real and imaginary - components of the Fourier coefficients (excluding the imaginary - component of the 2-point wave) as + component followed by the real and imaginary components of the + Fourier coefficients (excluding the imaginary component of the + 2-point wave) as :python:`x=[X0, Re(X1), Im(X1), ..., Re(Xn)]`. - For multiple degrees of freedom, the state vector is stacked, as - :python:`[[x1], [x2], ...]`. + + If :python:`zero_freq = False` (not default), the mean (DC) component + :python:`X0` is excluded, and the matrix/vector length is reduced by 1. Parameters ---------- transfer_mat Complex transfer matrix. zero_freq - Whether the first elements :python:`tranfser_mat[:,:,0]` - correspond to frequency zero. + Whether the first frequency should be zero. """ ndof = transfer_mat.shape[1] assert transfer_mat.shape[2] == ndof @@ -1547,9 +1518,7 @@ def block(re, im): return np.array([[re, -im], [im, re]]) re = np.real(Zp) im = np.imag(Zp) blocks = [block(ire, iim) for (ire, iim) in zip(re[:-1], im[:-1])] - blocks = blocks + [re[-1]] - if zero_freq: - blocks = [Zp0] + blocks + blocks = [Zp0] + blocks + [re[-1]] elem[idof][jdof] = block_diag(*blocks) return np.block(elem) @@ -2431,7 +2400,6 @@ def decompose_state( state: ndarray, ndof: int, nfreq: int, - nfreq_start: Optional[int] = 0, ) -> tuple[ndarray, ndarray]: """Split the state vector into the WEC dynamics state and the optimization (control) state. @@ -2449,8 +2417,6 @@ def decompose_state( Number of degrees of freedom for the WEC dynamics. nfreq Number of frequencies. - nfreq_start - Frequency index at which frequency vector starts. Returns ------- @@ -2459,55 +2425,66 @@ def decompose_state( state_opt Optimization (control) state. """ - nstate_wec = ndof * ncomponents(nfreq, nfreq_start) + nstate_wec = ndof * ncomponents(nfreq) return state[:nstate_wec], state[nstate_wec:] def frequency_parameters( freqs: ArrayLike, - precision: Optional[int] = 10, + zero_freq: bool = True, ) -> tuple[float, int]: - """Return the fundamental frequency, the number of frequencies, and - first frequency in a frequency array. + """Return the fundamental frequency and the number of frequencies + in a frequency array. This function can be used as a check for inputs to other functions since it raises an error if the frequency vector does not have - the correct format (equally spaced). + the correct format :python:`freqs = [0, f1, 2*f1, ..., nfreq*f1]` + (or :python:`freqs = [f1, 2*f1, ..., nfreq*f1]` if + :python:`zero_freq = False`). Parameters ---------- freqs - The frequency array with equal spacing. - precision - Controls rounding of fundamental frequency. + The frequency array, starting at zero and having equal spacing. + zero_freq + Whether the first frequency should be zero. Returns ------- f1 - Fundamental frequency :python:`f1` [:math:`Hz`]. - This is also the spacing. + Fundamental frequency :python:`f1` [:math:`Hz`] nfreq - Number of frequencies. - nfreq_start - Frequency (index) at which vector starts. + Number of frequencies (not including zero frequency), + i.e., :python:`freqs = [0, f1, 2*f1, ..., nfreq*f1]`. Raises ------ ValueError If the frequency vector is not evenly spaced. + ValueError + If the zero-frequency was expected but not included or not + expected but included. """ - f1 = freqs[1] - freqs[0] - f1 = f1 if precision is None else round(f1, precision) - nfreq_start = round(f reqs[0]/f1) - assert np.isclose(nfreq_start, freqs[0]/f1) - nfreq = len(freqs) - f_check = np.arange(nfreq_start, nfreq_start+nfreq)*f1 - if not np.allclose(f_check, freqs): - print(f_check) - print(freqs) - raise ValueError("Frequency array must be evenly spaced by " + - "the fundamental frequency ") - return f1, nfreq, nfreq_start + if np.isclose(freqs[0], 0.0): + if zero_freq: + freqs0 = freqs[:] + else: + raise ValueError('Zero frequency was included.') + else: + if zero_freq: + raise ValueError( + 'Frequency array must start with the zero frequency.') + else: + freqs0 = np.concatenate([[0.0,], freqs]) + + f1 = freqs0[1] + nfreq = len(freqs0) - 1 + f_check = np.arange(0, f1*(nfreq+0.5), f1) + if not np.allclose(f_check, freqs0): + raise ValueError("Frequency array 'omega' must be evenly spaced by" + + "the fundamental frequency " + + "(i.e., 'omega = [0, f1, 2*f1, ..., nfreq*f1]')") + return f1, nfreq def time_results(fd: DataArray, time: DataArray) -> ndarray: