Skip to content

Commit 1c288e1

Browse files
nasa-ddallenasa-ddalle
committed
Add getidx_ordered_stats() to apply confidnece level to ordered stats
1 parent 5f69aa1 commit 1c288e1

File tree

3 files changed

+100
-87
lines changed

3 files changed

+100
-87
lines changed

‎cape/cfdx/cntl.py‎

Lines changed: 0 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -4498,34 +4498,6 @@ def abspath(self, fname: str) -> str:
44984498

44994499
# Copy files
45004500
def copy_files(self, i: int):
4501-
r"""Copy files from *Mesh* section
4502-
4503-
This applies to both *CopyFiles* and *CopyAsFiles* in the
4504-
*Mesh* section. The former will copy a given file into the run
4505-
folder for case *i* using the base name of the original (source)
4506-
file. Using
4507-
4508-
.. code-block:: javascript
4509-
4510-
"Mesh": {
4511-
"CopyAsFiles": {
4512-
"inputs/mesh-config02.ugrid": "mesh.ugrid"
4513-
}
4514-
}
4515-
4516-
will copy the file ``inputs/mesh-config02.ugrid`` into the run
4517-
folder but name it ``mesh.ugrid`` there.
4518-
4519-
:Call:
4520-
>>> cntl.copy_files(i)
4521-
:Inputs:
4522-
*cntl*: :class:`cape.cfdx.cntl.Cntl`
4523-
Overall CAPE control instance
4524-
*i*: :class:`int`
4525-
Case index
4526-
:Versions:
4527-
* 2025-09-19 ``@ddalle``: v1.0
4528-
"""
45294501
# Ensure case index is set
45304502
self.opts.setx_i(i)
45314503
# Create case folder
@@ -4584,35 +4556,6 @@ def _copy_as_files(self, i: int):
45844556

45854557
# Link files
45864558
def link_files(self, i: int):
4587-
r"""Link files from *Mesh* section
4588-
4589-
This applies to both *LinkFiles* and *LinkAsFiles* in the
4590-
*Mesh* section. The former will copy a given file into the run
4591-
folder for case *i* using the base name of the original (source)
4592-
file. Using
4593-
4594-
.. code-block:: javascript
4595-
4596-
"Mesh": {
4597-
"LinkAsFiles": {
4598-
"inputs/mesh-config02.ugrid": "mesh.ugrid"
4599-
}
4600-
}
4601-
4602-
will create a link (using the absolute path) from
4603-
``inputs/mesh-config02.ugrid`` to ``mesh.ugrid`` in the case run
4604-
folder.
4605-
4606-
:Call:
4607-
>>> cntl.link_files(i)
4608-
:Inputs:
4609-
*cntl*: :class:`cape.cfdx.cntl.Cntl`
4610-
Overall CAPE control instance
4611-
*i*: :class:`int`
4612-
Case index
4613-
:Versions:
4614-
* 2025-09-19 ``@ddalle``: v1.0
4615-
"""
46164559
# Ensure case index is set
46174560
self.opts.setx_i(i)
46184561
# Create case folder

‎cape/cfdx/cntlbase.py‎

Lines changed: 41 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -2525,34 +2525,69 @@ def abspath(self, fname: str) -> str:
25252525
# Copy files
25262526
@abstractmethod
25272527
def copy_files(self, i: int):
2528-
r"""Copy specified files to case *i* run folder
2528+
r"""Copy files from *Mesh* section
2529+
2530+
This applies to both *CopyFiles* and *CopyAsFiles* in the
2531+
*Mesh* section. The former will copy a given file into the run
2532+
folder for case *i* using the base name of the original (source)
2533+
file. Using
2534+
2535+
.. code-block:: javascript
2536+
2537+
"Mesh": {
2538+
"CopyAsFiles": {
2539+
"inputs/mesh-config02.ugrid": "mesh.ugrid"
2540+
}
2541+
}
2542+
2543+
will copy the file ``inputs/mesh-config02.ugrid`` into the run
2544+
folder but name it ``mesh.ugrid`` there.
25292545
25302546
:Call:
25312547
>>> cntl.copy_files(i)
25322548
:Inputs:
2533-
*cntl*: :class:`Cntl`
2534-
CAPE main control instance
2549+
*cntl*: :class:`cape.cfdx.cntl.Cntl`
2550+
Overall CAPE control instance
25352551
*i*: :class:`int`
25362552
Case index
25372553
:Versions:
25382554
* 2025-03-26 ``@ddalle``: v1.0
2555+
* 2025-09-19 ``@ddalle``: v1.1; *CopyAsFiles*
25392556
"""
25402557
pass
25412558

25422559
# Link files
25432560
@abstractmethod
25442561
def link_files(self, i: int):
2545-
r"""Link specified files to case *i* run folder
2562+
r"""Link files from *Mesh* section
2563+
2564+
This applies to both *LinkFiles* and *LinkAsFiles* in the
2565+
*Mesh* section. The former will copy a given file into the run
2566+
folder for case *i* using the base name of the original (source)
2567+
file. Using
2568+
2569+
.. code-block:: javascript
2570+
2571+
"Mesh": {
2572+
"LinkAsFiles": {
2573+
"inputs/mesh-config02.ugrid": "mesh.ugrid"
2574+
}
2575+
}
2576+
2577+
will create a link (using the absolute path) from
2578+
``inputs/mesh-config02.ugrid`` to ``mesh.ugrid`` in the case run
2579+
folder.
25462580
25472581
:Call:
25482582
>>> cntl.link_files(i)
25492583
:Inputs:
2550-
*cntl*: :class:`Cntl`
2551-
CAPE main control instance
2584+
*cntl*: :class:`cape.cfdx.cntl.Cntl`
2585+
Overall CAPE control instance
25522586
*i*: :class:`int`
25532587
Case index
25542588
:Versions:
25552589
* 2025-03-26 ``@ddalle``: v1.0
2590+
* 2025-09-19 ``@ddalle``: v1.1; *LinkAsFiles*
25562591
"""
25572592
pass
25582593

‎cape/statutils.py‎

Lines changed: 59 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -28,11 +28,46 @@
2828

2929
# Statistics modules from SciPy
3030
try:
31+
from scipy.stats import binom
3132
from scipy.stats import norm
3233
from scipy.stats import t as student
3334
except ImportError:
3435
pass
3536

37+
38+
# Get orderered stats with confidence level
39+
def getidx_ordered_stats(p: float, cl: float, n: int) -> float:
40+
r"""Get number of cases to exclude for ordered statistics w/ CL
41+
42+
:Call:
43+
>>> m = getidx_ordered_stats(p, cl, n)
44+
:Inputs:
45+
*p*: :class:`float` | :class:`np.ndarray`\ [:class:`float`]
46+
Coverage fraction, probability
47+
*cl*: :class:`float` | :class:`np.ndarray`\ [:class:`float`]
48+
Confidence interval, or array thereof
49+
*n*: :class:`int`
50+
Number of samples (must be scalar)
51+
:Versions:
52+
* 2025-09-24 ``@ddalle``: v1.0
53+
"""
54+
# Check for lower bounds
55+
if isinstance(p, (float, np.floating)) and (p < 0.5):
56+
return n - getidx_ordered_stats(1-p, cl, n)
57+
# Array of failure counts
58+
m_excl = np.arange(n + 1)
59+
m_keep = m_excl[::-1]
60+
# Use binomial distribution to get p(exactly k failures)
61+
p_excl = binom.pmf(m_keep, n, p)
62+
# Cumulative probability
63+
total_p_excl = np.cumsum(p_excl)
64+
total_confidence = 1 - total_p_excl
65+
# Get unique values to avoid issues with confidence < machine prec.
66+
tmp, idx = np.unique(total_confidence[::-1], return_index=True)
67+
# Interpolate to requested confidnece level
68+
return np.interp(cl, tmp, m_excl[idx])
69+
70+
3671
# Get ordered stats
3772
def get_ordered_stats(V, cov=None, onesided=False, **kw):
3873
r"""Calculate coverage using ordered statistics
@@ -61,7 +96,7 @@ def get_ordered_stats(V, cov=None, onesided=False, **kw):
6196
*vlim*: :class:`float`
6297
Upper limit of one-sided coverage interval
6398
:Versions:
64-
* 2021-09-30 ``@ddalle``: Version 1.0
99+
* 2021-09-30 ``@ddalle``: v1.0
65100
"""
66101
# Get standard deviation counts
67102
ksig = kw.get("ksig")
@@ -119,7 +154,7 @@ def get_ordered_lower(V, cov):
119154
than or equal to *v*; may be interpolated between sorted
120155
values of *V*
121156
:Versions:
122-
* 2021-09-30 ``@ddalle``: Version 1.0
157+
* 2021-09-30 ``@ddalle``: v1.0
123158
"""
124159
# Get size
125160
n = len(V)
@@ -164,7 +199,7 @@ def get_ordered_upper(V, cov):
164199
or equal to *v*; may be interpolated between sorted values
165200
of *V*
166201
:Versions:
167-
* 2021-09-30 ``@ddalle``: Version 1.0
202+
* 2021-09-30 ``@ddalle``: v1.0
168203
"""
169204
# Get size
170205
n = len(V)
@@ -190,16 +225,16 @@ def get_ordered_upper(V, cov):
190225
cova = ia / float(n)
191226
# Interpolate ... n = 1 / (covb-cova)
192227
return va + n*(cov-cova)*(vb-va)
193-
228+
194229

195230
# Calculate range
196231
def get_range(R, cov=None, **kw):
197232
r"""Calculate Student's t-distribution confidence range
198-
233+
199234
If the nominal application of the Student's t-distribution fails to
200235
cover a high enough fraction of the data, the bounds are extended
201236
until the data is covered.
202-
237+
203238
:Call:
204239
>>> width = get_range(R, cov, **kw)
205240
:Inputs:
@@ -218,8 +253,8 @@ def get_range(R, cov=None, **kw):
218253
*width*: :class:`float`
219254
Half-width of confidence region
220255
:Versions:
221-
* 2018-09-28 ``@ddalle``: Version 1.0
222-
* 2021-09-20 ``@ddalle``: Version 1.1
256+
* 2018-09-28 ``@ddalle``: v1.0
257+
* 2021-09-20 ``@ddalle``: v1.1
223258
- use :func:`_parse_options`
224259
- allow 100% coverage
225260
- remove confusing *kcov* vs *ksig* scaling
@@ -273,11 +308,11 @@ def get_range(R, cov=None, **kw):
273308
# Calculate interval
274309
def get_coverage(dx, cov=None, **kw):
275310
r"""Calculate Student's *t*\ -distribution confidence range
276-
311+
277312
If the nominal application of the Student's t-distribution fails to
278313
cover a high enough fraction of the data, the bounds are extended
279314
until *cov* (user-defined fraction) of the data is covered.
280-
315+
281316
:Call:
282317
>>> width = get_coverage(dx, cov, **kw)
283318
:Inputs:
@@ -296,8 +331,8 @@ def get_coverage(dx, cov=None, **kw):
296331
*width*: :class:`float`
297332
Half-width of confidence region
298333
:Versions:
299-
* 2019-02-04 ``@ddalle``: Version 1.0
300-
* 2021-09-20 ``@ddalle``: Version 1.1
334+
* 2019-02-04 ``@ddalle``: v1.0
335+
* 2021-09-20 ``@ddalle``: v1.1
301336
- use :func:`_parse_options`
302337
- allow 100% coverage
303338
- remove confusing *kcov* vs *ksig* scaling
@@ -311,11 +346,11 @@ def get_coverage(dx, cov=None, **kw):
311346
# Calculate interval
312347
def get_cov_interval(dx, cov=None, **kw):
313348
r"""Calculate Student's *t*\ -distribution confidence range
314-
349+
315350
If the nominal application of the Student's t-distribution fails to
316351
cover a high enough fraction of the data, the bounds are extended
317352
until *cov* (user-defined fraction) of the data is covered.
318-
353+
319354
:Call:
320355
>>> a, b = get_cov_interval(dx, cov, **kw)
321356
:Inputs:
@@ -336,8 +371,8 @@ def get_cov_interval(dx, cov=None, **kw):
336371
*b*: :class:`float`
337372
Upper bound of coverage interval
338373
:Versions:
339-
* 2019-02-04 ``@ddalle``: Version 1.0
340-
* 2021-09-20 ``@ddalle``: Version 1.1
374+
* 2019-02-04 ``@ddalle``: v1.0
375+
* 2021-09-20 ``@ddalle``: v1.1
341376
- use :func:`_parse_options`
342377
- allow 100% coverage
343378
- remove confusing *kcov* vs *ksig* scaling
@@ -372,7 +407,7 @@ def get_cov_interval(dx, cov=None, **kw):
372407
b = vmu + width
373408
# --- Coverage Check ---
374409
# Filter cases that are outside bounds
375-
J = np.logical_and(a<=dx, dx<=b)
410+
J = np.logical_and(a <= dx, dx <= b)
376411
# Count cases outside the bounds
377412
ncov = np.count_nonzero(J)
378413
# Check coverage
@@ -396,7 +431,7 @@ def get_cov_interval(dx, cov=None, **kw):
396431
# Filter outliers
397432
def check_outliers_range(R, cov=None, **kw):
398433
r"""Find outliers in an array of ranges
399-
434+
400435
:Call:
401436
>>> I = check_outliers_range(R, cov, **kw)
402437
:Inputs:
@@ -415,7 +450,7 @@ def check_outliers_range(R, cov=None, **kw):
415450
*I*: :class:`np.ndarray`\ [:class:`bool`]
416451
Flags for non-outlier cases, ``False`` if case is an outlier
417452
:Versions:
418-
* 2021-02-20 ``@ddalle``: Version 1.0
453+
* 2021-02-20 ``@ddalle``: v1.0
419454
"""
420455
# --- Setup ---
421456
# Enforce array (copy to preserve original data)
@@ -452,8 +487,8 @@ def check_outliers_range(R, cov=None, **kw):
452487
# Update degrees of freedom
453488
df = N - n1
454489
# Recalculate statistics
455-
vmu = np.mean(dx[I])
456-
vstd = np.std(dx[I])
490+
vmu = np.mean(R[I])
491+
vstd = np.std(R[I])
457492
# Find outliers
458493
I = R / vstd <= osig
459494
J = np.logical_not(I)
@@ -466,7 +501,7 @@ def check_outliers_range(R, cov=None, **kw):
466501
# Filter outliers
467502
def check_outliers(dx, cov=None, **kw):
468503
r"""Find outliers in a data set
469-
504+
470505
:Call:
471506
>>> I = check_outliers(dx, cov, **kw)
472507
:Inputs:
@@ -485,8 +520,8 @@ def check_outliers(dx, cov=None, **kw):
485520
*I*: :class:`np.ndarray`\ [:class:`bool`]
486521
Flags for non-outlier cases, ``False`` if case is an outlier
487522
:Versions:
488-
* 2019-02-04 ``@ddalle``: Version 1.0
489-
* 2021-09-20 ``@ddalle``: Version 1.1
523+
* 2019-02-04 ``@ddalle``: v1.0
524+
* 2021-09-20 ``@ddalle``: v1.1
490525
- use :func:`_parse_options`
491526
- allow 100% coverage
492527
"""

0 commit comments

Comments
 (0)