Skip to content

Commit aa8d512

Browse files
ali-cetin-4ssali-cetin-4ss
authored
feat: add "equal-energy" discrete direction method for spreading function classes (#81)
1 parent 81bf68a commit aa8d512

File tree

2 files changed

+324
-3
lines changed

2 files changed

+324
-3
lines changed

‎src/waveresponse/_core.py‎

Lines changed: 59 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,9 @@
44
from numbers import Number
55

66
import numpy as np
7-
from scipy.integrate import trapezoid
7+
from scipy.integrate import quad, trapezoid
88
from scipy.interpolate import RegularGridInterpolator as RGI
9+
from scipy.optimize import root_scalar
910
from scipy.special import gamma
1011

1112
from ._utils import _robust_modulus, complex_to_polar, polar_to_complex
@@ -1538,7 +1539,7 @@ def from_spectrum1d(
15381539
1-D array of grid frequency coordinates. Positive and monotonically increasing.
15391540
dirs : array-like
15401541
1-D array of grid direction coordinates. Positive and monotonically increasing.
1541-
Must cover the directional range [0, 360) degrees (or [0, 2 * numpy.pi) radians).
1542+
Must cover the directional range [0, 360) degrees (or [0, 2 * pi) radians).
15421543
spectrum1d : array-like
15431544
1-D array of non-directional spectrum density values. These 1-D spectrum
15441545
values will be scaled according to the spreading function, and distributed
@@ -2418,6 +2419,62 @@ def _spread_fun(self, omega, theta):
24182419
"""
24192420
raise NotImplementedError()
24202421

2422+
def discrete_directions(self, n, direction_offset=0.0):
2423+
"""
2424+
Split the spreading function into discrete direction bins with
2425+
"equal energy", i.e. equal area under the curve. The direcitons
2426+
representing the bins are chosen to have equal area under the curve on
2427+
each side within the bin.
2428+
2429+
Parameters
2430+
----------
2431+
n : int
2432+
Number of discrete directions.
2433+
direction_offset : float, default
2434+
A offset to add to the discrete directions. Units should be
2435+
according to the `degrees` flag given during initialization.
2436+
2437+
Returns
2438+
-------
2439+
ndarray
2440+
A sequence of direction representing "equal energy" bins with range
2441+
wrapped to [0, 360) degrees or [0, 2 * pi) radians according
2442+
to the `degrees` flag given during initialization.
2443+
"""
2444+
if self._degrees:
2445+
x_lb = -180.0
2446+
x_ub = 180.0
2447+
periodicity = 360.0
2448+
else:
2449+
x_lb = -np.pi
2450+
x_ub = np.pi
2451+
periodicity = 2.0 * np.pi
2452+
2453+
total_area, _ = quad(
2454+
lambda theta: self(None, theta), x_lb, x_ub, epsabs=1.0e-6, epsrel=0.0
2455+
)
2456+
2457+
half_bin_edges = np.empty(2 * n - 1)
2458+
2459+
x_prev = x_lb
2460+
for i in range(1, 2 * n):
2461+
target_area = total_area * i / (2 * n)
2462+
res = root_scalar(
2463+
lambda x: quad(
2464+
lambda theta: self(None, theta), x_lb, x, epsabs=1.0e-6, epsrel=0.0
2465+
)[0]
2466+
- target_area,
2467+
bracket=[x_prev, x_ub],
2468+
)
2469+
2470+
if not res.converged:
2471+
raise RuntimeError(f"Failed find the directions: {res.flag}")
2472+
2473+
x_prev = res.root
2474+
half_bin_edges[i - 1] = x_prev
2475+
2476+
return _robust_modulus(half_bin_edges[::2] + direction_offset, periodicity)
2477+
24212478

24222479
class CosineHalfSpreading(BaseSpreading):
24232480
"""

‎tests/test_core.py‎

Lines changed: 265 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99
from scipy.interpolate import RegularGridInterpolator as RGI
1010

1111
import waveresponse as wr
12-
from waveresponse import ( # BinGrid,
12+
from waveresponse import (
1313
RAO,
1414
CosineFullSpreading,
1515
CosineHalfSpreading,
@@ -5970,6 +5970,138 @@ def test_independent_of_frequency(self):
59705970

59715971
assert len(np.unique(np.array(spread_out_list))) == 1
59725972

5973+
@pytest.mark.parametrize(
5974+
"n,dirs_expect",
5975+
[
5976+
[2, np.array([-26.4, 26.4])],
5977+
[3, np.array([-37.6, 0.0, 37.6])],
5978+
[4, np.array([-44.6, -12.5, 12.5, 44.6])],
5979+
[5, np.array([-49.5, -20.5, 0.0, 20.5, 49.5])],
5980+
],
5981+
)
5982+
def test_discrete_directions_no_offset_full(self, n, dirs_expect):
5983+
spreading = CosineFullSpreading(4, degrees=True)
5984+
np.testing.assert_allclose(
5985+
(spreading.discrete_directions(n) + 1e-8) % 360.0,
5986+
dirs_expect % 360.0,
5987+
rtol=0.0,
5988+
atol=0.1,
5989+
)
5990+
5991+
spreading = CosineFullSpreading(4, degrees=False)
5992+
np.testing.assert_allclose(
5993+
np.degrees(spreading.discrete_directions(n) + 1e-8) % 360.0,
5994+
dirs_expect % 360.0,
5995+
rtol=0.0,
5996+
atol=0.1,
5997+
)
5998+
5999+
@pytest.mark.parametrize(
6000+
"n,dirs_expect",
6001+
[
6002+
[2, np.array([-16.4, 36.4])],
6003+
[3, np.array([-27.6, 10.0, 47.6])],
6004+
[4, np.array([-34.6, -2.5, 22.5, 54.6])],
6005+
[5, np.array([-39.5, -10.5, 10.0, 30.5, 59.5])],
6006+
],
6007+
)
6008+
def test_discrete_directions_offset_full(self, n, dirs_expect):
6009+
spreading = CosineFullSpreading(4, degrees=True)
6010+
np.testing.assert_allclose(
6011+
(spreading.discrete_directions(n, 10.0) + 1e-8) % 360.0,
6012+
dirs_expect % 360.0,
6013+
rtol=0.0,
6014+
atol=0.1,
6015+
)
6016+
6017+
spreading = CosineFullSpreading(4, degrees=False)
6018+
np.testing.assert_allclose(
6019+
np.degrees(spreading.discrete_directions(n, np.radians(10)) + 1e-8) % 360.0,
6020+
dirs_expect % 360.0,
6021+
rtol=0.0,
6022+
atol=0.1,
6023+
)
6024+
6025+
@pytest.mark.parametrize(
6026+
"n,dirs_expect",
6027+
[
6028+
[2, np.array([-36.4, 16.4])],
6029+
[3, np.array([-47.6, -10.0, 27.6])],
6030+
[4, np.array([-54.6, -22.5, 2.5, 34.6])],
6031+
[5, np.array([-59.5, -30.5, -10.0, 10.5, 39.5])],
6032+
],
6033+
)
6034+
def test_discrete_directions_neg_offset_full(self, n, dirs_expect):
6035+
spreading = CosineFullSpreading(4, degrees=True)
6036+
np.testing.assert_allclose(
6037+
(spreading.discrete_directions(n, -10.0) + 1e-8) % 360.0,
6038+
dirs_expect % 360.0,
6039+
rtol=0.0,
6040+
atol=0.1,
6041+
)
6042+
6043+
spreading = CosineFullSpreading(4, degrees=False)
6044+
np.testing.assert_allclose(
6045+
np.degrees(spreading.discrete_directions(n, np.radians(-10)) + 1e-8)
6046+
% 360.0,
6047+
dirs_expect % 360.0,
6048+
rtol=0.0,
6049+
atol=0.1,
6050+
)
6051+
6052+
@pytest.mark.parametrize(
6053+
"n,dirs_expect",
6054+
[
6055+
[2, np.array([332.6, 385.4])],
6056+
[3, np.array([321.4, 359.0, 396.6])],
6057+
[4, np.array([314.4, 346.5, 371.5, 403.6])],
6058+
[5, np.array([309.5, 338.5, 359.0, 379.5, 408.5])],
6059+
],
6060+
)
6061+
def test_discrete_directions_large_offset_full(self, n, dirs_expect):
6062+
spreading = CosineFullSpreading(4, degrees=True)
6063+
np.testing.assert_allclose(
6064+
(spreading.discrete_directions(n, 359) + 1e-8) % 360.0,
6065+
dirs_expect % 360.0,
6066+
rtol=0.0,
6067+
atol=0.1,
6068+
)
6069+
6070+
spreading = CosineFullSpreading(4, degrees=False)
6071+
np.testing.assert_allclose(
6072+
np.degrees(spreading.discrete_directions(n, np.radians(359)) + 1e-8)
6073+
% 360.0,
6074+
dirs_expect % 360.0,
6075+
rtol=0.0,
6076+
atol=0.1,
6077+
)
6078+
6079+
@pytest.mark.parametrize(
6080+
"n,dirs_expect",
6081+
[
6082+
[2, np.array([-11.0, 31.0])],
6083+
[3, np.array([-20.0, 10.0, 40.0])],
6084+
[4, np.array([-25.6, 0.1, 19.9, 45.6])],
6085+
[5, np.array([-29.5, -6.3, 10.0, 26.3, 49.5])],
6086+
],
6087+
)
6088+
def test_discrete_directions_other_s_full(self, n, dirs_expect):
6089+
spreading = CosineFullSpreading(6.5, degrees=True)
6090+
np.testing.assert_allclose(
6091+
(spreading.discrete_directions(n, 10) + 1e-8) % 360.0,
6092+
dirs_expect % 360.0,
6093+
rtol=0.0,
6094+
atol=0.1,
6095+
)
6096+
6097+
spreading = CosineFullSpreading(6.5, degrees=False)
6098+
np.testing.assert_allclose(
6099+
np.degrees(spreading.discrete_directions(n, np.radians(10)) + 1e-8) % 360.0,
6100+
dirs_expect % 360.0,
6101+
rtol=0.0,
6102+
atol=0.1,
6103+
)
6104+
59736105

59746106
class Test_CosineHalfSpreading:
59756107
def test__init__(self):
@@ -6085,3 +6217,135 @@ def test__call__degrees(self, f, d, s, spread_expect):
60856217
def test__call__radians(self, f, d, s, spread_expect):
60866218
spreading = CosineHalfSpreading(s, degrees=False)
60876219
assert spreading(f, d) == pytest.approx(spread_expect)
6220+
6221+
@pytest.mark.parametrize(
6222+
"n,dirs_expect",
6223+
[
6224+
[2, np.array([-13.2, 13.2])],
6225+
[3, np.array([-18.8, 0.0, 18.8])],
6226+
[4, np.array([-22.3, -6.3, 6.3, 22.3])],
6227+
[5, np.array([-24.8, -10.3, 0.0, 10.3, 24.8])],
6228+
],
6229+
)
6230+
def test_discrete_directions_no_offset(self, n, dirs_expect):
6231+
spreading = CosineHalfSpreading(4, degrees=True)
6232+
np.testing.assert_allclose(
6233+
(spreading.discrete_directions(n) + 1e-8) % 360.0,
6234+
dirs_expect % 360.0,
6235+
rtol=0.0,
6236+
atol=0.1,
6237+
)
6238+
6239+
spreading = CosineHalfSpreading(4, degrees=False)
6240+
np.testing.assert_allclose(
6241+
np.degrees(spreading.discrete_directions(n) + 1e-8) % 360.0,
6242+
dirs_expect % 360.0,
6243+
rtol=0.0,
6244+
atol=0.1,
6245+
)
6246+
6247+
@pytest.mark.parametrize(
6248+
"n,dirs_expect",
6249+
[
6250+
[2, np.array([-3.2, 23.2])],
6251+
[3, np.array([-8.8, 10.0, 28.8])],
6252+
[4, np.array([-12.3, 3.7, 16.3, 32.3])],
6253+
[5, np.array([-14.8, -0.3, 10.0, 20.3, 34.8])],
6254+
],
6255+
)
6256+
def test_discrete_directions_offset(self, n, dirs_expect):
6257+
spreading = CosineHalfSpreading(4, degrees=True)
6258+
np.testing.assert_allclose(
6259+
(spreading.discrete_directions(n, 10.0) + 1e-8) % 360.0,
6260+
dirs_expect % 360.0,
6261+
rtol=0.0,
6262+
atol=0.1,
6263+
)
6264+
6265+
spreading = CosineHalfSpreading(4, degrees=False)
6266+
np.testing.assert_allclose(
6267+
np.degrees(spreading.discrete_directions(n, np.radians(10)) + 1e-8) % 360.0,
6268+
dirs_expect % 360.0,
6269+
rtol=0.0,
6270+
atol=0.1,
6271+
)
6272+
6273+
@pytest.mark.parametrize(
6274+
"n,dirs_expect",
6275+
[
6276+
[2, np.array([-23.2, 3.2])],
6277+
[3, np.array([-28.8, -10.0, 8.8])],
6278+
[4, np.array([-32.3, -16.3, -3.7, 12.3])],
6279+
[5, np.array([-34.8, -20.3, -10.0, 0.3, 14.8])],
6280+
],
6281+
)
6282+
def test_discrete_directions_neg_offset(self, n, dirs_expect):
6283+
spreading = CosineHalfSpreading(4, degrees=True)
6284+
np.testing.assert_allclose(
6285+
(spreading.discrete_directions(n, -10.0) + 1e-8) % 360.0,
6286+
dirs_expect % 360.0,
6287+
rtol=0.0,
6288+
atol=0.1,
6289+
)
6290+
6291+
spreading = CosineHalfSpreading(4, degrees=False)
6292+
np.testing.assert_allclose(
6293+
np.degrees(spreading.discrete_directions(n, np.radians(-10)) + 1e-8)
6294+
% 360.0,
6295+
dirs_expect % 360.0,
6296+
rtol=0.0,
6297+
atol=0.1,
6298+
)
6299+
6300+
@pytest.mark.parametrize(
6301+
"n,dirs_expect",
6302+
[
6303+
[2, np.array([345.8, 372.2])],
6304+
[3, np.array([340.2, 359, 377.8])],
6305+
[4, np.array([336.7, 352.7, 365.3, 381.3])],
6306+
[5, np.array([334.2, 348.7, 359.0, 369.3, 383.8])],
6307+
],
6308+
)
6309+
def test_discrete_directions_large_offset(self, n, dirs_expect):
6310+
spreading = CosineHalfSpreading(4, degrees=True)
6311+
np.testing.assert_allclose(
6312+
(spreading.discrete_directions(n, 359) + 1e-8) % 360.0,
6313+
dirs_expect % 360.0,
6314+
rtol=0.0,
6315+
atol=0.1,
6316+
)
6317+
6318+
spreading = CosineHalfSpreading(4, degrees=False)
6319+
np.testing.assert_allclose(
6320+
np.degrees(spreading.discrete_directions(n, np.radians(359)) + 1e-8)
6321+
% 360.0,
6322+
dirs_expect % 360.0,
6323+
rtol=0.0,
6324+
atol=0.1,
6325+
)
6326+
6327+
@pytest.mark.parametrize(
6328+
"n,dirs_expect",
6329+
[
6330+
[2, np.array([-0.5, 20.5])],
6331+
[3, np.array([-5.0, 10.0, 25.0])],
6332+
[4, np.array([-7.8, 5.0, 15.0, 27.8])],
6333+
[5, np.array([-9.8, 1.8, 10.0, 18.2, 29.8])],
6334+
],
6335+
)
6336+
def test_discrete_directions_other_s(self, n, dirs_expect):
6337+
spreading = CosineHalfSpreading(6.5, degrees=True)
6338+
np.testing.assert_allclose(
6339+
(spreading.discrete_directions(n, 10) + 1e-8) % 360.0,
6340+
dirs_expect % 360.0,
6341+
rtol=0.0,
6342+
atol=0.1,
6343+
)
6344+
6345+
spreading = CosineHalfSpreading(6.5, degrees=False)
6346+
np.testing.assert_allclose(
6347+
np.degrees(spreading.discrete_directions(n, np.radians(10)) + 1e-8) % 360.0,
6348+
dirs_expect % 360.0,
6349+
rtol=0.0,
6350+
atol=0.1,
6351+
)

0 commit comments

Comments
 (0)