Skip to content

Commit 66ad3a6

Browse files
heptaflarheptaflar
committed
Update: OCX model functions and tests
1 parent c18d672 commit 66ad3a6

2 files changed

Lines changed: 227 additions & 110 deletions

File tree

‎test/oceancolour/test_ocx.py‎

Lines changed: 120 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,9 @@
77
import numpy as np
88
import pandas as pd
99

10+
from uncertaintyx.oceancolour.ocx import CI
1011
from uncertaintyx.oceancolour.ocx import OC4
12+
from uncertaintyx.oceancolour.ocx import OCI
1113

1214

1315
def read_test_data(
@@ -34,7 +36,57 @@ def read_test_data(
3436

3537

3638
class OCxTest(unittest.TestCase):
37-
"""Tests OC4 and OCI model functions on optical water type classes."""
39+
"""Tests OCX model functions on optical water type classes."""
40+
41+
def test_ci(self):
42+
w, R, u, M, m = read_test_data( # noqa : N806
43+
"test.resources", "owt.csv"
44+
)
45+
W = np.broadcast_to(w, R.shape) # noqa : N806
46+
47+
f = CI()
48+
x = np.stack([W[:, [1, 4, 5]], R[:, [1, 4, 5]]], axis=1)
49+
u = np.stack([np.broadcast_to(0.5, (M, 3)), u[:, [1, 4, 5]]], axis=1)
50+
p = f.estimate()
51+
y = f.eval(p, x)
52+
53+
self.assertEqual((M,), y.shape)
54+
self.assertAlmostEqual(0.022, y[0], delta=0.001)
55+
self.assertAlmostEqual(0.033, y[1], delta=0.001)
56+
self.assertAlmostEqual(0.048, y[2], delta=0.001)
57+
self.assertAlmostEqual(0.074, y[3], delta=0.001)
58+
self.assertAlmostEqual(0.111, y[4], delta=0.001)
59+
self.assertAlmostEqual(0.152, y[5], delta=0.001)
60+
self.assertAlmostEqual(0.205, y[6], delta=0.001)
61+
self.assertAlmostEqual(0.259, y[7], delta=0.001)
62+
self.assertAlmostEqual(0.350, y[8], delta=0.001)
63+
self.assertAlmostEqual(0.426, y[9], delta=0.001)
64+
self.assertAlmostEqual(0.538, y[10], delta=0.001)
65+
self.assertAlmostEqual(3.897, y[11], delta=0.001)
66+
# 12 and 13 are out of domain
67+
self.assertTrue(np.isfinite(y[12]))
68+
self.assertTrue(np.isfinite(y[13]))
69+
70+
U = np.square(u) # noqa : N806
71+
V = f.lpu_x(p, x, U) # noqa : N806
72+
self.assertEqual((M,), V.shape)
73+
74+
v = np.sqrt(V)
75+
self.assertAlmostEqual(0.026, v[0], delta=0.001)
76+
self.assertAlmostEqual(0.024, v[1], delta=0.001)
77+
self.assertAlmostEqual(0.042, v[2], delta=0.001)
78+
self.assertAlmostEqual(0.068, v[3], delta=0.001)
79+
self.assertAlmostEqual(0.070, v[4], delta=0.001)
80+
self.assertAlmostEqual(0.074, v[5], delta=0.001)
81+
self.assertAlmostEqual(0.107, v[6], delta=0.001)
82+
self.assertAlmostEqual(0.100, v[7], delta=0.001)
83+
self.assertAlmostEqual(0.174, v[8], delta=0.001)
84+
self.assertAlmostEqual(0.204, v[9], delta=0.001)
85+
self.assertAlmostEqual(0.200, v[10], delta=0.001)
86+
self.assertAlmostEqual(5.461, v[11], delta=0.001)
87+
# 12 and 13 are out of domain
88+
self.assertTrue(np.isfinite(v[12]))
89+
self.assertTrue(np.isfinite(y[13]))
3890

3991
def test_oc4(self):
4092
_, R, u, M, m = read_test_data( # noqa : N806
@@ -43,6 +95,7 @@ def test_oc4(self):
4395

4496
f = OC4()
4597
x = R[:, 1:-1]
98+
u = u[:, 1:-1]
4699
p = f.estimate()
47100
y = f.eval(p, x)
48101

@@ -62,10 +115,74 @@ def test_oc4(self):
62115
self.assertAlmostEqual(3.295, y[12], delta=0.001)
63116
self.assertAlmostEqual(3.427, y[13], delta=0.001)
64117

65-
U = np.square(u[:, 1:-1]) # noqa : N806
66-
V = f.propagate_x_diag(p, x, U)
118+
U = np.square(u) # noqa : N806
119+
V = f.lpu_x(p, x, U) # noqa : N806
67120
self.assertEqual((M,), V.shape)
68121

122+
v = np.sqrt(V)
123+
self.assertAlmostEqual(0.116, v[0], delta=0.001)
124+
self.assertAlmostEqual(0.093, v[1], delta=0.001)
125+
self.assertAlmostEqual(0.136, v[2], delta=0.001)
126+
self.assertAlmostEqual(0.171, v[3], delta=0.001)
127+
self.assertAlmostEqual(0.135, v[4], delta=0.001)
128+
self.assertAlmostEqual(0.120, v[5], delta=0.001)
129+
self.assertAlmostEqual(0.159, v[6], delta=0.001)
130+
self.assertAlmostEqual(0.159, v[7], delta=0.001)
131+
self.assertAlmostEqual(0.360, v[8], delta=0.001)
132+
self.assertAlmostEqual(0.559, v[9], delta=0.001)
133+
self.assertAlmostEqual(1.158, v[10], delta=0.001)
134+
self.assertAlmostEqual(6.411, v[11], delta=0.001)
135+
self.assertAlmostEqual(4.170, v[12], delta=0.001)
136+
self.assertAlmostEqual(4.362, v[13], delta=0.001)
137+
138+
def test_oci(self):
139+
w, R, u, M, m = read_test_data( # noqa : N806
140+
"test.resources", "owt.csv"
141+
)
142+
W = np.broadcast_to(w, R.shape) # noqa : N806
143+
144+
f = OCI()
145+
x = np.stack([W[:, 1:], R[:, 1:]], axis=1)
146+
u = np.stack([np.broadcast_to(0.5, (M, 5)), u[:, 1:]], axis=1)
147+
p = f.estimate()
148+
y = f.eval(p, x)
149+
150+
self.assertEqual((M,), y.shape)
151+
self.assertAlmostEqual(0.022, y[0], delta=0.001)
152+
self.assertAlmostEqual(0.033, y[1], delta=0.001)
153+
self.assertAlmostEqual(0.048, y[2], delta=0.001)
154+
self.assertAlmostEqual(0.074, y[3], delta=0.001)
155+
self.assertAlmostEqual(0.111, y[4], delta=0.001)
156+
self.assertAlmostEqual(0.152, y[5], delta=0.001)
157+
self.assertAlmostEqual(0.205, y[6], delta=0.001)
158+
self.assertAlmostEqual(0.260, y[7], delta=0.001) # blend
159+
self.assertAlmostEqual(0.426, y[8], delta=0.001)
160+
self.assertAlmostEqual(0.572, y[9], delta=0.001)
161+
self.assertAlmostEqual(1.022, y[10], delta=0.001)
162+
self.assertAlmostEqual(4.174, y[11], delta=0.001)
163+
self.assertAlmostEqual(3.295, y[12], delta=0.001)
164+
self.assertAlmostEqual(3.427, y[13], delta=0.001)
165+
166+
U = np.square(u) # noqa : N806
167+
V = f.lpu_x(p, x, U) # noqa : N806
168+
self.assertEqual((M,), V.shape)
169+
170+
v = np.sqrt(V)
171+
self.assertAlmostEqual(0.026, v[0], delta=0.001)
172+
self.assertAlmostEqual(0.024, v[1], delta=0.001)
173+
self.assertAlmostEqual(0.042, v[2], delta=0.001)
174+
self.assertAlmostEqual(0.068, v[3], delta=0.001)
175+
self.assertAlmostEqual(0.070, v[4], delta=0.001)
176+
self.assertAlmostEqual(0.074, v[5], delta=0.001)
177+
self.assertAlmostEqual(0.107, v[6], delta=0.001)
178+
self.assertAlmostEqual(0.117, v[7], delta=0.001) # blend
179+
self.assertAlmostEqual(0.360, v[8], delta=0.001)
180+
self.assertAlmostEqual(0.559, v[9], delta=0.001)
181+
self.assertAlmostEqual(1.158, v[10], delta=0.001)
182+
self.assertAlmostEqual(6.411, v[11], delta=0.001)
183+
self.assertAlmostEqual(4.170, v[12], delta=0.001)
184+
self.assertAlmostEqual(4.362, v[13], delta=0.001)
185+
69186

70187
if __name__ == "__main__":
71188
unittest.main()

‎uncertaintyx/oceancolour/ocx.py‎

Lines changed: 107 additions & 107 deletions
Original file line numberDiff line numberDiff line change
@@ -28,112 +28,9 @@
2828
from ..m.jax import ToM
2929

3030

31-
class Blended(ToM):
31+
class CI(ToM):
3232
"""
33-
The blended OC4/OCI model function.
34-
35-
The two or three nearest wavebands to 412, 443, 490 and 510 nm
36-
are used for the blue, while the nearest waveband to 555 nm is
37-
used for the green, and the nearest waveband to 670 nm is used
38-
for the red.
39-
40-
The blending occurs when :class:`OC4` is between, e.g.,
41-
0.25-0.35 mg m-3, creating a smooth handover between :class:`OCI`
42-
(low chlorophyll specialist) and :class:`OC4` (baseline). The
43-
blending uses perfectly normalized weights. Low :class:`OCI`
44-
values favour :class:`OCI` while high :class:`OCI` values
45-
favour :class:`OC4` through self-weighting. A quadratic term
46-
creates curvature that eliminates boundary discontinuities while
47-
:class:`OC4` acts as regime detector.
48-
"""
49-
50-
def __init__(self, b: Literal[0, 1] = 0):
51-
"""
52-
Creates a new model function instance.
53-
54-
:param b: The index of the blue waveband near 443 nm.
55-
"""
56-
self._oci = OCI()
57-
self._ocx = OC4()
58-
59-
def f(p, x):
60-
r"""
61-
The blended OCI/OCX model function.
62-
63-
The blending is self-adjusting.
64-
65-
Let :math:`p = \in \mathbb{R}^{k}, k = 2 + 2 + 5` be the
66-
model parameters and let :math:`\lambda = (\lambda_1,\dots
67-
\lambda_\mathrm{m}) \in \mathbb{R}^{m}` denote the central
68-
wavelengths of :math:`2 \le m-2 \le 3` wavebands in the blue,
69-
and a single waveband :math:`\lambda_{m-1}` in the green, and
70-
a single waveband :math:`\lambda_{m}` in the red. Further let
71-
72-
.. math::
73-
x = (\lambda, R_\mathrm{rs}(\lambda))
74-
\in \mathbb{R}^{2 \times m}
75-
76-
denote the function inputs. Then:
77-
78-
:param p: The parameters :math:`p \in \mathbb{R}^{k}`.
79-
:param x: The inputs :math:`x \in \mathbb{R}^{2 \times m}`.
80-
:returns: The chlorophyll concentration (mg m-3).
81-
"""
82-
83-
def blend(p, ci, cx):
84-
"""Returns the blended chlorophyll concentration."""
85-
ti = p[1] - ci
86-
tx = ci - p[0]
87-
return (ci * ti + cx * tx) / (p[1] - p[0])
88-
89-
ci = self._oci.f(p[2:4], x[:, [b, -2, -1]])
90-
cx = self._ocx.f(p[4:9], x[1, 0:-1])
91-
92-
return jnp.where(
93-
cx < p[0],
94-
ci,
95-
jnp.where(cx > p[1], cx, blend(p, ci, cx)),
96-
)
97-
98-
super().__init__(f)
99-
100-
def estimate(
101-
self,
102-
x: np.ndarray | None = None,
103-
y: np.ndarray | None = None,
104-
preset: Literal[
105-
"CZCS",
106-
"GOCI",
107-
"HAWKEYE",
108-
"MERIS",
109-
"MODIS",
110-
"OCTS",
111-
"OLCI",
112-
"PACE",
113-
"SEAWIFS",
114-
"VIIRS20",
115-
"VIIRS21",
116-
]
117-
| None = None,
118-
) -> np.ndarray:
119-
"""
120-
Returns the blended OCI/OCX default parameter values.
121-
122-
Elements ``[0:2]`` refer to the blending, elements ``[2:4]``
123-
refer to OCI, and elements ``[4:9]`` refer to OCX.
124-
"""
125-
return np.concatenate(
126-
(
127-
np.array([0.25, 0.35]),
128-
self._oci.estimate(x, y),
129-
self._ocx.estimate(x, y, preset),
130-
)
131-
)
132-
133-
134-
class OCI(ToM):
135-
"""
136-
NASA's ocean colour chlorophyll index (OCI) model function.
33+
NASA's ocean colour chlorophyll index (CI) model function.
13734
13835
The nearest wavebands to 443, 555, and 670 nm are used for the blue,
13936
green, and red, respectively, for all sensors.
@@ -228,9 +125,9 @@ def estimate(
228125

229126
class OC4(ToM):
230127
"""
231-
NASA's OC3/OC4 chlorophyll model function.
128+
NASA's OC4 (and OC3) chlorophyll model function.
232129
233-
The two or three nearest wavebands to 412, 443, 490 and 510 nm
130+
The three (or two) nearest wavebands to 412, 443, 490 and 510 nm
234131
are used for the blue, while the nearest waveband to 555 nm is
235132
used for the green, for all sensors.
236133
"""
@@ -314,3 +211,106 @@ def estimate(
314211
case _:
315212
pass
316213
return np.array(params)
214+
215+
216+
class OCI(ToM):
217+
"""
218+
The blended OC4/CI model function.
219+
220+
The two or three nearest wavebands to 412, 443, 490 and 510 nm
221+
are used for the blue, while the nearest waveband to 555 nm is
222+
used for the green, and the nearest waveband to 670 nm is used
223+
for the red.
224+
225+
The blending occurs when :class:`OC4` is between, e.g.,
226+
0.25-0.35 mg m-3, creating a smooth handover between :class:`OCI`
227+
(low chlorophyll specialist) and :class:`OC4` (baseline). The
228+
blending uses perfectly normalized weights. Low :class:`CI`
229+
values favour :class:`CI` while high :class:`CI` values
230+
favour :class:`OC4` through self-weighting. A quadratic term
231+
creates curvature that eliminates boundary discontinuities while
232+
:class:`OC4` acts as regime detector.
233+
"""
234+
235+
def __init__(self, b: Literal[0, 1] = 0):
236+
"""
237+
Creates a new model function instance.
238+
239+
:param b: The index of the blue waveband near 443 nm.
240+
"""
241+
self.m_ci: ToM = CI()
242+
self.m_oc: ToM = OC4()
243+
244+
def f(p, x):
245+
r"""
246+
The blended OCI/OCX model function.
247+
248+
The blending is self-adjusting.
249+
250+
Let :math:`p = \in \mathbb{R}^{k}, k = 2 + 2 + 5` be the
251+
model parameters and let :math:`\lambda = (\lambda_1,\dots
252+
\lambda_\mathrm{m}) \in \mathbb{R}^{m}` denote the central
253+
wavelengths of :math:`2 \le m-2 \le 3` wavebands in the blue,
254+
and a single waveband :math:`\lambda_{m-1}` in the green, and
255+
a single waveband :math:`\lambda_{m}` in the red. Further let
256+
257+
.. math::
258+
x = (\lambda, R_\mathrm{rs}(\lambda))
259+
\in \mathbb{R}^{2 \times m}
260+
261+
denote the function inputs. Then:
262+
263+
:param p: The parameters :math:`p \in \mathbb{R}^{k}`.
264+
:param x: The inputs :math:`x \in \mathbb{R}^{2 \times m}`.
265+
:returns: The chlorophyll concentration (mg m-3).
266+
"""
267+
268+
def blend(p, ci, cx):
269+
"""Returns the blended chlorophyll concentration."""
270+
ti = p[1] - ci
271+
tx = ci - p[0]
272+
return (ci * ti + cx * tx) / (p[1] - p[0])
273+
274+
ci = self.m_ci.f(p[2:4], x[:, [b, -2, -1]])
275+
oc = self.m_oc.f(p[4:9], x[1, 0:-1])
276+
277+
return jnp.where(
278+
oc < p[0],
279+
ci,
280+
jnp.where(oc > p[1], oc, blend(p, ci, oc)),
281+
)
282+
283+
super().__init__(f)
284+
285+
def estimate(
286+
self,
287+
x: np.ndarray | None = None,
288+
y: np.ndarray | None = None,
289+
preset: Literal[
290+
"CZCS",
291+
"GOCI",
292+
"HAWKEYE",
293+
"MERIS",
294+
"MODIS",
295+
"OCTS",
296+
"OLCI",
297+
"PACE",
298+
"SEAWIFS",
299+
"VIIRS20",
300+
"VIIRS21",
301+
]
302+
| None = None,
303+
) -> np.ndarray:
304+
"""
305+
Returns the blended OCI/OCX default parameter values.
306+
307+
Elements ``[0:2]`` refer to the blending, elements ``[2:4]``
308+
refer to CI, and elements ``[4:9]`` refer to OC4.
309+
"""
310+
return np.concatenate(
311+
(
312+
np.array([0.25, 0.35]),
313+
self.m_ci.estimate(x, y),
314+
self.m_oc.estimate(x, y, preset),
315+
)
316+
)

0 commit comments

Comments
 (0)