From 967ae03268c52134250ab37350ec15bb4e04fab9 Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 14:57:27 -0400
Subject: [PATCH 01/10] nexus: introduce interval dist

---
 nexus/nexus/statistics.py | 294 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 294 insertions(+)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index e1e05fc3e1..1ff9b0a446 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -701,3 +701,297 @@ def series_stats(x,t_auto=None):
     x_stderr = np.std(x)/np.sqrt(N_eff)
     return x_mean,x_stderr,t_auto
 #end def series_stats
+
+
+############################################################################
+#                                                                          #
+#              Line-crossing and interval-distribution analysis            #
+#              ------------------------------------------------            #
+#                                                                          #
+# These functions represent a time series as intervals between neighboring #
+# values and count their overlap along the value axis.  The resulting      #
+# interval distribution is a line-crossing density: locations with many    #
+# overlapping segments identify values persistently traversed by the       #
+# series.                                                                  #
+#                                                                          #
+# The distribution and its peak provide robust center estimates that       #
+# emphasize locally stable, equilibrium-like portions of a fluctuating     #
+# series.  Rolling versions track this center over time, while related     #
+# utilities support broader interval-distribution analysis.                #
+#                                                                          #
+############################################################################
+
+
+def _int_dist_input(x1,x2=None):
+    # process input types
+    x1 = np.asarray(x1)
+    if x2 is not None:
+        x2 = np.asarray(x2)
+        if x1.ndim>1:
+            assert x1.size==max(x1.shape)
+            x1 = x1.ravel()
+        if x2.ndim>1:
+            assert x2.size==max(x2.shape)
+            x2 = x2.ravel()
+        xi = np.vstack((x1,x2)).T
+    else:
+        xi = x1
+    # xi is array of N intervals
+    assert xi.ndim==2
+    assert xi.shape[1]==2
+    # check endpoint ordering
+    assert (xi[:,1]-xi[:,0]).min()>=0
+    si = np.empty(xi.shape,dtype=int)
+    si[:,0] =  1
+    si[:,1] = -1
+    return xi,si
+#end def _int_dist_input
+
+
+
+def interval_distribution(x1,x2=None):
+    xi,si = _int_dist_input(x1,x2)
+    xi = xi.ravel()
+    si = si.ravel()
+    # organize by edge order
+    order = xi.argsort()
+    xi = xi[order]
+    si = si[order]
+    # make the interval counting distribution
+    cd = {}
+    n=0
+    for s,xv in zip(si,xi):
+        n += s
+        cd[xv] = n
+    # first organize into sorted point/edge arrays
+    x = []
+    c = []
+    for xv in sorted(cd.keys()):
+        x.append(xv)
+        c.append(cd[xv])
+    x = np.array(x)
+    c = np.array(c)
+    # next into interval array
+    xi = np.zeros((len(x)-1,2),dtype=x.dtype)
+    xi[:,0] = x[:-1]
+    xi[:,1] = x[1:]
+    ci = c[:-1].copy()
+    return xi,ci
+#end def interval_distribution
+
+
+
+def plot_interval_dist(xi,ci,style='b.-'):
+    import matplotlib.pyplot as plt
+    xif = xi.ravel()
+    cif = np.zeros(xif.shape)
+    cif[::2]  = ci
+    cif[1::2] = ci
+    plt.axhline(0,color='k')
+    plt.plot(xif,cif,style)
+#end def plot_interval_dist
+
+
+
+def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
+    if method=='interval_mid':
+        cm = ci.max()
+        xm = xi[ci==cm].mean()
+    elif method=='interval_rand':
+        cm    = ci.max()
+        xi    = xi[ci==cm]
+        u     = np.random.uniform(size=len(xi))
+        x1,x2 = xi.T
+        dx    = x2-x1
+        xmid  = (x2+x1)/2
+        x     = xmid + (u-0.5)*dx/2
+        xm    = x.mean()
+    elif method=='quad_peak':
+        # first find the overall max
+        imax = ci.argmax()
+        cf = peak_frac*ci[imax]
+        # flatten for search and fit
+        xif = xi.ravel()
+        cif = np.zeros(xif.shape)
+        cif[::2]  = ci
+        cif[1::2] = ci
+        N = len(cif)
+        # move left until below peak frac
+        i1 = imax
+        for n in range(N):
+            if i1==0 or cif[i1]<cf:
+                break
+            i1 -= 1
+        # move right until below peak frac
+        i2 = imax
+        for n in range(N):
+            if i2==N-1 or cif[i2]<cf:
+                break
+            i2 += 1
+        # fit the peak
+        xp = xif[i1:i2+1]
+        cp = cif[i1:i2+1]
+        p = np.polyfit(xp,cp,2)
+        # find the max
+        xm = -p[1]/(2*p[0])
+        cm = np.polyval(p,xm)
+    else:
+        raise ValueError(f'unrecognized int. dist. max method: "{method}"')
+    if not height:
+        return xm
+    else:
+        return xm,cm
+#end def interval_dist_peak
+
+
+
+def rolling_interval_dist_peak(
+        x1,
+        x2          = None,
+        window      = 10,
+        step        = 5,
+        method      = 'interval_mid',
+        ret_height  = False,
+        ret_windows = False,
+        ):
+    assert step<=window
+    interval_mid  = method=='interval_mid'
+    interval_rand = method=='interval_rand'
+    # map inputs to intervals
+    xia,sia = _int_dist_input(x1,x2)
+    N = len(xia)
+    assert N>=window
+    # find window segments
+    windows = []
+    i1 = 0
+    for n in range(N):
+        i2 = i1 + window
+        if i2<N:
+            windows.append((i1,i2))
+        elif i2==N:
+            windows.append((i1,i2))
+            break
+        else: #i2>N
+            windows.append((N-window,N))
+            break
+        i1 += step
+    assert n+1<N
+    assert len(windows)>0
+    assert windows[-1][1]==N
+    # find interval dist peaks in each window
+    xp = []
+    cp = []
+    for i1,i2 in windows:
+        xi = xia[i1:i2]
+        si = sia[i1:i2]
+        assert len(xi)==window
+        xi = xi.ravel()
+        si = si.ravel()
+        order = xi.argsort()
+        xi = xi[order]
+        si = si[order]
+        cd = {}
+        n=0
+        for s,xv in zip(si,xi):
+            n += s
+            cd[xv] = n
+        x = []
+        c = []
+        for xv in sorted(cd.keys()):
+            x.append(xv)
+            c.append(cd[xv])
+        x = np.array(x)
+        c = np.array(c)[:-1]
+        cm = c.max()
+        if interval_mid:
+            xmid = (x[:-1]+x[1:])/2
+            xm   = xmid[c==cm].mean()
+        elif interval_rand:
+            x1 = x[:-1]
+            x2 = x[1:]
+            xmid = (x1+x2)/2
+            dx   = x2-x1
+            u    = np.random.uniform(size=len(xmid))
+            xc   = xmid + (u-0.5)*dx/2
+            xm   = xc[c==cm].mean()
+        else:
+            raise ValueError(f'method "{method}" is unrecognized')
+        xp.append(xm)
+        cp.append(cm)
+    xp = np.array(xp)
+    ret = [xp]
+    if ret_height:
+        cp = np.array(cp)
+        ret.append(cp)
+    if ret_windows:
+        ret.append(windows)
+    if len(ret)==0:
+        return ret[0]
+    else:
+        return tuple(ret)
+#end def rolling_interval_dist_peak
+
+
+
+
+
+
+def time_series_intervals(x,t=None):
+    xi = np.empty((len(x)-1,2),dtype=x.dtype)
+    for n in range(len(x)-1):
+        xi[n,0] = x[n]
+        xi[n,1] = x[n+1]
+    xi = np.sort(xi,axis=1)
+    if t is None:
+        return xi,None
+    else:
+        ti = (t[:-1]+t[1:])/2
+        return xi,ti
+#end def time_series_intervals
+
+
+
+
+def line_crossing_distribution(x,nperm=0):
+    if nperm>1:
+        xperm = []
+        for n in range(nperm):
+            xp = x.copy()
+            np.random.shuffle(xp)
+            xperm.append(xp)
+        x = np.hstack(xperm)
+    xi,_ = time_series_intervals(x,t=None)
+    xi,ci = interval_distribution(xi)
+    return xi,ci
+#end def line_crossing_distribution
+
+
+
+def lcd_peak(x,method='interval_mid',peak_frac=0.5):
+    xi,ci = line_crossing_distribution(x)
+    xp = interval_dist_peak(xi,ci,method=method,peak_frac=peak_frac)
+    return xp
+#end def lcd_peak
+
+
+
+def lcd_smooth(x,
+               t      = None,
+               window = 10,
+               step   = 5,
+               method = 'interval_rand',
+               ):
+    xi,ti = time_series_intervals(x,t)
+    xp,windows = rolling_interval_dist_peak(
+        xi,
+        window      = window,
+        step        = step,
+        method      = method,
+        ret_windows = True,
+        )
+    if t is None:
+        return xp
+    else:
+        tp = np.array([ti[i1:i2].mean() for i1,i2 in windows])
+        return xp,tp
+#end def lcd_smooth

From db34898ae3a1bb9455911abefe1030a632242691 Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 15:00:30 -0400
Subject: [PATCH 02/10] nexus: add tests

---
 nexus/nexus/statistics.py            |  34 ++++----
 nexus/nexus/tests/test_statistics.py | 119 +++++++++++++++++++++++++++
 2 files changed, 134 insertions(+), 19 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index 1ff9b0a446..a3ddb71a8b 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -722,6 +722,21 @@ def series_stats(x,t_auto=None):
 ############################################################################
 
 
+def time_series_intervals(x,t=None):
+    xi = np.empty((len(x)-1,2),dtype=x.dtype)
+    for n in range(len(x)-1):
+        xi[n,0] = x[n]
+        xi[n,1] = x[n+1]
+    xi = np.sort(xi,axis=1)
+    if t is None:
+        return xi,None
+    else:
+        ti = (t[:-1]+t[1:])/2
+        return xi,ti
+#end def time_series_intervals
+
+
+
 def _int_dist_input(x1,x2=None):
     # process input types
     x1 = np.asarray(x1)
@@ -933,25 +948,6 @@ def rolling_interval_dist_peak(
 
 
 
-
-
-
-def time_series_intervals(x,t=None):
-    xi = np.empty((len(x)-1,2),dtype=x.dtype)
-    for n in range(len(x)-1):
-        xi[n,0] = x[n]
-        xi[n,1] = x[n+1]
-    xi = np.sort(xi,axis=1)
-    if t is None:
-        return xi,None
-    else:
-        ti = (t[:-1]+t[1:])/2
-        return xi,ti
-#end def time_series_intervals
-
-
-
-
 def line_crossing_distribution(x,nperm=0):
     if nperm>1:
         xperm = []
diff --git a/nexus/nexus/tests/test_statistics.py b/nexus/nexus/tests/test_statistics.py
index fd1f14a590..ab3234c893 100644
--- a/nexus/nexus/tests/test_statistics.py
+++ b/nexus/nexus/tests/test_statistics.py
@@ -368,3 +368,122 @@ def fake_autocorr_time(x_arg):
             ):
             statistics.series_stats(x,t_auto=t_auto_invalid)
 #end def test_series_stats
+
+
+
+def test_time_series_intervals():
+    """Check adjacent-value intervals and their associated midpoint times."""
+    x = np.array([3.,1.,2.])
+    t = np.array([0.,2.,5.])
+
+    intervals,times = statistics.time_series_intervals(x,t)
+    np.testing.assert_array_equal(intervals,[[1.,3.],[1.,2.]])
+    np.testing.assert_array_equal(times,[1.,3.5])
+
+    intervals,no_times = statistics.time_series_intervals(x)
+    np.testing.assert_array_equal(intervals,[[1.,3.],[1.,2.]])
+    assert(no_times is None)
+#end def test_time_series_intervals
+
+
+
+def test_interval_distribution_and_peak(monkeypatch):
+    """Check interval-overlap counts and the supported peak selections."""
+    endpoints = np.array([1.,2.,3.])
+    upper     = np.array([4.,5.,6.])
+    intervals,counts = statistics.interval_distribution(endpoints,upper)
+    expected_intervals = np.array(
+        [[1.,2.],[2.,3.],[3.,4.],[4.,5.],[5.,6.]]
+        )
+    np.testing.assert_array_equal(intervals,expected_intervals)
+    np.testing.assert_array_equal(counts,[1,2,3,2,1])
+
+    matrix_intervals,matrix_counts = statistics.interval_distribution(
+        np.column_stack((endpoints,upper))
+        )
+    np.testing.assert_array_equal(matrix_intervals,expected_intervals)
+    np.testing.assert_array_equal(matrix_counts,counts)
+
+    peak_intervals = np.array([[0.,1.],[1.,2.],[2.,3.]])
+    peak_counts    = np.array([1,3,3])
+    peak,height = statistics.interval_dist_peak(
+        peak_intervals,
+        peak_counts,
+        height=True,
+        )
+    assert(peak==pytest.approx(2.))
+    assert(height==3)
+
+    monkeypatch.setattr(
+        statistics.np.random,
+        'uniform',
+        lambda size: np.full(size,.5),
+        )
+    assert(
+        statistics.interval_dist_peak(
+            peak_intervals,
+            peak_counts,
+            method='interval_rand',
+            )
+        ==pytest.approx(2.)
+        )
+
+    with pytest.raises(ValueError,match=r'unrecognized int. dist. max method'):
+        statistics.interval_dist_peak(peak_intervals,peak_counts,method='invalid')
+#end def test_interval_distribution_and_peak
+
+
+
+def test_rolling_interval_dist_peak_and_lcd_smooth():
+    """Check rolling peak locations, heights, window bounds, and times."""
+    intervals = np.array([[0.,2.],[1.,3.],[2.,4.],[3.,5.]])
+    peaks,heights,windows = statistics.rolling_interval_dist_peak(
+        intervals,
+        window=2,
+        step=2,
+        ret_height=True,
+        ret_windows=True,
+        )
+    np.testing.assert_allclose(peaks,[1.5,3.5])
+    np.testing.assert_array_equal(heights,[2,2])
+    assert(windows==[(0,2),(2,4)])
+
+    x = np.array([0.,2.,1.,3.])
+    t = np.array([0.,1.,3.,6.])
+    smooth,times = statistics.lcd_smooth(x,t,window=2,step=1,method='interval_mid')
+    np.testing.assert_allclose(smooth,[1.5,1.5])
+    np.testing.assert_allclose(times,[1.25,3.25])
+    np.testing.assert_allclose(
+        statistics.lcd_smooth(x,window=2,step=1,method='interval_mid'),
+        smooth,
+        )
+#end def test_rolling_interval_dist_peak_and_lcd_smooth
+
+
+
+def test_line_crossing_distribution_and_lcd_peak():
+    """Check LCD interval counts and their peak-derived center."""
+    x = np.array([0.,2.,1.])
+    intervals,counts = statistics.line_crossing_distribution(x)
+    np.testing.assert_array_equal(intervals,[[0.,1.],[1.,2.]])
+    np.testing.assert_array_equal(counts,[1,2])
+    assert(statistics.lcd_peak(x)==pytest.approx(1.5))
+#end def test_line_crossing_distribution_and_lcd_peak
+
+
+
+def test_plot_interval_dist():
+    """Check that interval-distribution plotting adds the expected lines."""
+    matplotlib = pytest.importorskip('matplotlib')
+    matplotlib.use('Agg')
+    import matplotlib.pyplot as plt
+
+    figure,axis = plt.subplots()
+    plt.sca(axis)
+    statistics.plot_interval_dist(
+        np.array([[0.,1.],[1.,2.]]),
+        np.array([1,2]),
+        )
+    assert(len(axis.lines)==2)
+    plt.close(figure)
+#end def test_plot_interval_dist

From c52dedeaae03e6358cf5eddcccd583c6eb12b43e Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 15:13:05 -0400
Subject: [PATCH 03/10] nexus: guard inputs

---
 nexus/nexus/statistics.py            | 175 +++++++++++++++++----------
 nexus/nexus/tests/test_statistics.py |  58 +++++++++
 2 files changed, 170 insertions(+), 63 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index a3ddb71a8b..53a25f7ff6 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -73,6 +73,29 @@ def _paired_real_arrays(x,y):
 #end def _paired_real_arrays
 
 
+def _real_vector(x,name):
+    """Return a finite real vector, flattening vector-shaped arrays."""
+    x = np.asarray(x)
+    if np.iscomplexobj(x):
+        msg = f'{name} must be real-valued'
+        raise ValueError(msg)
+    if x.ndim>1 and np.max(x.shape)==x.size:
+        x = x.ravel()
+    if x.ndim!=1:
+        msg = f'{name} must be one-dimensional'
+        raise ValueError(msg)
+    try:
+        x = np.asarray(x,dtype=float)
+    except (TypeError,ValueError):
+        msg = f'{name} must be numeric'
+        raise ValueError(msg) from None
+    if not np.all(np.isfinite(x)):
+        msg = f'{name} must contain only finite values'
+        raise ValueError(msg)
+    return x
+#end def _real_vector
+
+
 def theil_sen(x,y):
     """Return the Theil--Sen slope and intercept for paired observations.
 
@@ -723,6 +746,15 @@ def series_stats(x,t_auto=None):
 
 
 def time_series_intervals(x,t=None):
+    x = _real_vector(x,'data array')
+    if len(x)<2:
+        msg = 'data array must contain at least two values'
+        raise ValueError(msg)
+    if t is not None:
+        t = _real_vector(t,'time array')
+        if len(t)!=len(x):
+            msg = 'time array must have the same length as data array'
+            raise ValueError(msg)
     xi = np.empty((len(x)-1,2),dtype=x.dtype)
     for n in range(len(x)-1):
         xi[n,0] = x[n]
@@ -738,24 +770,34 @@ def time_series_intervals(x,t=None):
 
 
 def _int_dist_input(x1,x2=None):
-    # process input types
-    x1 = np.asarray(x1)
     if x2 is not None:
-        x2 = np.asarray(x2)
-        if x1.ndim>1:
-            assert x1.size==max(x1.shape)
-            x1 = x1.ravel()
-        if x2.ndim>1:
-            assert x2.size==max(x2.shape)
-            x2 = x2.ravel()
+        x1 = _real_vector(x1,'lower endpoints')
+        x2 = _real_vector(x2,'upper endpoints')
+        if len(x1)!=len(x2):
+            msg = 'interval endpoint arrays must have equal lengths'
+            raise ValueError(msg)
         xi = np.vstack((x1,x2)).T
     else:
-        xi = x1
+        xi = np.asarray(x1)
+        if np.iscomplexobj(xi):
+            msg = 'interval array must be real-valued'
+            raise ValueError(msg)
     # xi is array of N intervals
-    assert xi.ndim==2
-    assert xi.shape[1]==2
+    if xi.ndim!=2 or xi.shape[1]!=2:
+        msg = 'interval array must have shape (n,2)'
+        raise ValueError(msg)
+    if len(xi)==0:
+        msg = 'interval array must not be empty'
+        raise ValueError(msg)
+    try:
+        xi = np.asarray(xi,dtype=float)
+    except (TypeError,ValueError):
+        msg = 'interval array must be numeric'
+        raise ValueError(msg) from None
     # check endpoint ordering
-    assert (xi[:,1]-xi[:,0]).min()>=0
+    if np.any(xi[:,1]<xi[:,0]):
+        msg = 'interval upper endpoints must not be less than lower endpoints'
+        raise ValueError(msg)
     si = np.empty(xi.shape,dtype=int)
     si[:,0] =  1
     si[:,1] = -1
@@ -798,6 +840,11 @@ def interval_distribution(x1,x2=None):
 
 def plot_interval_dist(xi,ci,style='b.-'):
     import matplotlib.pyplot as plt
+    xi,_ = _int_dist_input(xi)
+    ci = _real_vector(ci,'interval counts')
+    if len(ci)!=len(xi):
+        msg = 'interval counts must have the same length as intervals'
+        raise ValueError(msg)
     xif = xi.ravel()
     cif = np.zeros(xif.shape)
     cif[::2]  = ci
@@ -809,6 +856,22 @@ def plot_interval_dist(xi,ci,style='b.-'):
 
 
 def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
+    xi,_ = _int_dist_input(xi)
+    ci = _real_vector(ci,'interval counts')
+    if len(ci)!=len(xi):
+        msg = 'interval counts must have the same length as intervals'
+        raise ValueError(msg)
+    if not isinstance(method,str):
+        msg = 'peak method must be a string'
+        raise ValueError(msg)
+    try:
+        peak_frac = float(peak_frac)
+    except (TypeError,ValueError):
+        msg = 'peak fraction must be a finite number'
+        raise ValueError(msg) from None
+    if not np.isfinite(peak_frac):
+        msg = 'peak fraction must be a finite number'
+        raise ValueError(msg)
     if method=='interval_mid':
         cm = ci.max()
         xm = xi[ci==cm].mean()
@@ -869,68 +932,45 @@ def rolling_interval_dist_peak(
         ret_height  = False,
         ret_windows = False,
         ):
-    assert step<=window
+    for value,name in ((window,'window'),(step,'step')):
+        if isinstance(value,(bool,np.bool_)) or not isinstance(
+            value,(int,np.integer)
+            ) or value<1:
+            msg = f'{name} must be a positive integer'
+            raise ValueError(msg)
+    if step>window:
+        msg = 'step must not exceed window'
+        raise ValueError(msg)
     interval_mid  = method=='interval_mid'
     interval_rand = method=='interval_rand'
+    if not interval_mid and not interval_rand:
+        msg = f'method "{method}" is unrecognized'
+        raise ValueError(msg)
     # map inputs to intervals
     xia,sia = _int_dist_input(x1,x2)
     N = len(xia)
-    assert N>=window
+    if N<window:
+        msg = 'window must not exceed the number of intervals'
+        raise ValueError(msg)
     # find window segments
-    windows = []
-    i1 = 0
-    for n in range(N):
-        i2 = i1 + window
-        if i2<N:
-            windows.append((i1,i2))
-        elif i2==N:
-            windows.append((i1,i2))
-            break
-        else: #i2>N
-            windows.append((N-window,N))
-            break
-        i1 += step
-    assert n+1<N
-    assert len(windows)>0
-    assert windows[-1][1]==N
+    starts = list(range(0,N-window+1,step))
+    if starts[-1]!=N-window:
+        starts.append(N-window)
+    windows = [(i1,i1+window) for i1 in starts]
     # find interval dist peaks in each window
     xp = []
     cp = []
     for i1,i2 in windows:
-        xi = xia[i1:i2]
-        si = sia[i1:i2]
-        assert len(xi)==window
-        xi = xi.ravel()
-        si = si.ravel()
-        order = xi.argsort()
-        xi = xi[order]
-        si = si[order]
-        cd = {}
-        n=0
-        for s,xv in zip(si,xi):
-            n += s
-            cd[xv] = n
-        x = []
-        c = []
-        for xv in sorted(cd.keys()):
-            x.append(xv)
-            c.append(cd[xv])
-        x = np.array(x)
-        c = np.array(c)[:-1]
-        cm = c.max()
+        xi,ci = interval_distribution(xia[i1:i2])
+        if len(ci)==0:
+            msg = 'each rolling window must span a nonzero interval'
+            raise ValueError(msg)
         if interval_mid:
-            xmid = (x[:-1]+x[1:])/2
-            xm   = xmid[c==cm].mean()
-        elif interval_rand:
-            x1 = x[:-1]
-            x2 = x[1:]
-            xmid = (x1+x2)/2
-            dx   = x2-x1
-            u    = np.random.uniform(size=len(xmid))
-            xc   = xmid + (u-0.5)*dx/2
-            xm   = xc[c==cm].mean()
+            xm,cm = interval_dist_peak(xi,ci,height=True)
         else:
-            raise ValueError(f'method "{method}" is unrecognized')
+            xm,cm = interval_dist_peak(
+                xi,ci,method='interval_rand',height=True
+                )
         xp.append(xm)
         cp.append(cm)
     xp = np.array(xp)
@@ -949,6 +989,15 @@ def rolling_interval_dist_peak(
 
 
 def line_crossing_distribution(x,nperm=0):
+    x = _real_vector(x,'data array')
+    if len(x)<2:
+        msg = 'data array must contain at least two values'
+        raise ValueError(msg)
+    if isinstance(nperm,(bool,np.bool_)) or not isinstance(
+        nperm,(int,np.integer)
+        ) or nperm<0:
+        msg = 'number of permutations must be a nonnegative integer'
+        raise ValueError(msg)
     if nperm>1:
         xperm = []
         for n in range(nperm):
diff --git a/nexus/nexus/tests/test_statistics.py b/nexus/nexus/tests/test_statistics.py
index ab3234c893..8865a404cd 100644
--- a/nexus/nexus/tests/test_statistics.py
+++ b/nexus/nexus/tests/test_statistics.py
@@ -387,6 +387,64 @@ def test_time_series_intervals():
 
 
 
+def test_interval_and_lcd_input_validation():
+    """Check diagnostics for malformed interval and LCD inputs."""
+    with pytest.raises(ValueError,match=r'data array must contain at least two values'):
+        statistics.time_series_intervals([1.])
+    with pytest.raises(ValueError,match=r'data array must contain only finite values'):
+        statistics.time_series_intervals([1.,np.nan])
+    with pytest.raises(ValueError,match=r'time array must have the same length'):
+        statistics.time_series_intervals([1.,2.],[0.])
+    with pytest.raises(ValueError,match=r'time array must be real-valued'):
+        statistics.time_series_intervals([1.,2.],[0.+1.j,1.+1.j])
+
+    invalid_intervals = [
+        (np.array([1.,2.]),None,r'interval array must have shape'),
+        (np.empty((0,2)),None,r'interval array must not be empty'),
+        (np.array([[2.,1.]]),None,r'upper endpoints must not be less'),
+        (np.array([[1.,np.nan]]),None,r'interval array must contain only finite values'),
+        (np.array([1.,2.]),np.array([3.]),r'endpoint arrays must have equal lengths'),
+        ]
+    for lower,upper,message in invalid_intervals:
+        with pytest.raises(ValueError,match=message):
+            statistics.interval_distribution(lower,upper)
+
+    intervals = np.array([[0.,1.],[1.,2.]])
+    for counts,message in [
+        ([1.],r'counts must have the same length'),
+        ([1.,np.nan],r'counts must contain only finite values'),
+        ]:
+        with pytest.raises(ValueError,match=message):
+            statistics.interval_dist_peak(intervals,counts)
+    with pytest.raises(ValueError,match=r'peak method must be a string'):
+        statistics.interval_dist_peak(intervals,[1.,2.],method=1)
+    with pytest.raises(ValueError,match=r'peak fraction must be a finite number'):
+        statistics.interval_dist_peak(intervals,[1.,2.],peak_frac=np.nan)
+
+    for window,step,message in [
+        (0,1,r'window must be a positive integer'),
+        (1,0,r'step must be a positive integer'),
+        (1,2,r'step must not exceed window'),
+        (3,1,r'window must not exceed the number of intervals'),
+        ]:
+        with pytest.raises(ValueError,match=message):
+            statistics.rolling_interval_dist_peak(intervals,window=window,step=step)
+    with pytest.raises(ValueError,match=r'method "invalid" is unrecognized'):
+        statistics.rolling_interval_dist_peak(
+            intervals,window=1,step=1,method='invalid'
+            )
+
+    for x,nperm,message in [
+        ([1.],0,r'data array must contain at least two values'),
+        ([1.,2.],-1,r'number of permutations must be a nonnegative integer'),
+        ([1.,2.],True,r'number of permutations must be a nonnegative integer'),
+        ]:
+        with pytest.raises(ValueError,match=message):
+            statistics.line_crossing_distribution(x,nperm=nperm)
+#end def test_interval_and_lcd_input_validation
+
+
+
 def test_interval_distribution_and_peak(monkeypatch):
     """Check interval-overlap counts and the supported peak selections."""
     endpoints = np.array([1.,2.,3.])

From a22d88919ebb2594056b4af6ecb50d71b0a4d136 Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 15:55:31 -0400
Subject: [PATCH 04/10] nexus: fixes

---
 nexus/nexus/statistics.py            | 65 ++++++++++++++--------------
 nexus/nexus/tests/test_statistics.py | 57 ++++++++++++++++++++++--
 2 files changed, 87 insertions(+), 35 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index 53a25f7ff6..0ae2908348 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -885,34 +885,30 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
         x     = xmid + (u-0.5)*dx/2
         xm    = x.mean()
     elif method=='quad_peak':
-        # first find the overall max
         imax = ci.argmax()
         cf = peak_frac*ci[imax]
-        # flatten for search and fit
-        xif = xi.ravel()
-        cif = np.zeros(xif.shape)
-        cif[::2]  = ci
-        cif[1::2] = ci
-        N = len(cif)
-        # move left until below peak frac
+        # Locate the contiguous high-count region in interval coordinates.
+        # ``imax`` indexes ``ci``, not the flattened endpoint array.
         i1 = imax
-        for n in range(N):
-            if i1==0 or cif[i1]<cf:
-                break
+        while i1>0 and ci[i1-1]>=cf:
             i1 -= 1
-        # move right until below peak frac
         i2 = imax
-        for n in range(N):
-            if i2==N-1 or cif[i2]<cf:
-                break
+        while i2+1<len(ci) and ci[i2+1]>=cf:
             i2 += 1
-        # fit the peak
-        xp = xif[i1:i2+1]
-        cp = cif[i1:i2+1]
-        p = np.polyfit(xp,cp,2)
-        # find the max
-        xm = -p[1]/(2*p[0])
-        cm = np.polyval(p,xm)
+        xp = xi[i1:i2+1].ravel()
+        cp = np.repeat(ci[i1:i2+1],2)
+        if len(np.unique(xp))<3:
+            xm = xi[ci==ci[imax]].mean()
+            cm = ci[imax]
+        else:
+            p = np.polyfit(xp,cp,2)
+            if not np.isfinite(p[0]) or p[0]>=0.:
+                xm = xi[ci==ci[imax]].mean()
+                cm = ci[imax]
+            else:
+                xm = -p[1]/(2*p[0])
+                xm = np.clip(xm,xp.min(),xp.max())
+                cm = np.polyval(p,xm)
     else:
         raise ValueError(f'unrecognized int. dist. max method: "{method}"')
     if not height:
@@ -998,16 +994,21 @@ def line_crossing_distribution(x,nperm=0):
         ) or nperm<0:
         msg = 'number of permutations must be a nonnegative integer'
         raise ValueError(msg)
-    if nperm>1:
-        xperm = []
-        for n in range(nperm):
-            xp = x.copy()
-            np.random.shuffle(xp)
-            xperm.append(xp)
-        x = np.hstack(xperm)
-    xi,_ = time_series_intervals(x,t=None)
-    xi,ci = interval_distribution(xi)
-    return xi,ci
+
+    # permutation-free (typical) case
+    if nperm==0:
+        xi,_ = time_series_intervals(x,t=None)
+        return interval_distribution(xi)
+
+    # use permutation shuffling
+    permutation_intervals = []
+    for n in range(nperm):
+        xp = x.copy()
+        np.random.shuffle(xp)
+        xi,_ = time_series_intervals(xp,t=None)
+        permutation_intervals.append(xi)
+    xi,ci = interval_distribution(np.vstack(permutation_intervals))
+    return xi,ci/nperm
 #end def line_crossing_distribution
 
 
diff --git a/nexus/nexus/tests/test_statistics.py b/nexus/nexus/tests/test_statistics.py
index 8865a404cd..6f23b62683 100644
--- a/nexus/nexus/tests/test_statistics.py
+++ b/nexus/nexus/tests/test_statistics.py
@@ -383,6 +383,13 @@ def test_time_series_intervals():
     intervals,no_times = statistics.time_series_intervals(x)
     np.testing.assert_array_equal(intervals,[[1.,3.],[1.,2.]])
     assert(no_times is None)
+
+    column_intervals,column_times = statistics.time_series_intervals(
+        x.reshape(-1,1),
+        t.reshape(-1,1),
+        )
+    np.testing.assert_array_equal(column_intervals,intervals)
+    np.testing.assert_array_equal(column_times,times)
 #end def test_time_series_intervals
 
 
@@ -402,7 +409,6 @@ def test_interval_and_lcd_input_validation():
         (np.array([1.,2.]),None,r'interval array must have shape'),
         (np.empty((0,2)),None,r'interval array must not be empty'),
         (np.array([[2.,1.]]),None,r'upper endpoints must not be less'),
-        (np.array([[1.,np.nan]]),None,r'interval array must contain only finite values'),
         (np.array([1.,2.]),np.array([3.]),r'endpoint arrays must have equal lengths'),
         ]
     for lower,upper,message in invalid_intervals:
@@ -441,6 +447,9 @@ def test_interval_and_lcd_input_validation():
         ]:
         with pytest.raises(ValueError,match=message):
             statistics.line_crossing_distribution(x,nperm=nperm)
+
+    with pytest.raises(ValueError,match=r'window must not exceed the number of intervals'):
+        statistics.lcd_smooth([1.,2.],window=2,step=1)
 #end def test_interval_and_lcd_input_validation
 
 
@@ -462,6 +471,13 @@ def test_interval_distribution_and_peak(monkeypatch):
     np.testing.assert_array_equal(matrix_intervals,expected_intervals)
     np.testing.assert_array_equal(matrix_counts,counts)
 
+    column_intervals,column_counts = statistics.interval_distribution(
+        endpoints.reshape(-1,1),
+        upper.reshape(1,-1),
+        )
+    np.testing.assert_array_equal(column_intervals,expected_intervals)
+    np.testing.assert_array_equal(column_counts,counts)
+
     peak_intervals = np.array([[0.,1.],[1.,2.],[2.,3.]])
     peak_counts    = np.array([1,3,3])
     peak,height = statistics.interval_dist_peak(
@@ -472,6 +488,26 @@ def test_interval_distribution_and_peak(monkeypatch):
     assert(peak==pytest.approx(2.))
     assert(height==3)
 
+    column_peak = statistics.interval_dist_peak(
+        peak_intervals,
+        peak_counts.reshape(-1,1),
+        )
+    assert(column_peak==pytest.approx(peak))
+
+    quadratic_intervals = np.array(
+        [[0.,1.],[1.,2.],[2.,3.],[3.,4.],[4.,5.]]
+        )
+    quadratic_counts = np.array([1.,2.,3.,2.,1.])
+    quadratic_peak,quadratic_height = statistics.interval_dist_peak(
+        quadratic_intervals,
+        quadratic_counts,
+        method='quad_peak',
+        height=True,
+    )
+    assert(quadratic_peak==pytest.approx(2.5))
+    assert(np.isfinite(quadratic_height))
+    assert(quadratic_height>0.)
+
     monkeypatch.setattr(
         statistics.np.random,
         'uniform',
@@ -519,13 +555,28 @@ def test_rolling_interval_dist_peak_and_lcd_smooth():
 
 
 
-def test_line_crossing_distribution_and_lcd_peak():
-    """Check LCD interval counts and their peak-derived center."""
+def test_line_crossing_distribution_and_lcd_peak(monkeypatch):
+    """Check LCD counts, peaks, and independently accumulated permutations."""
     x = np.array([0.,2.,1.])
+
+    def fail_shuffle(values):
+        pytest.fail('the nperm=0 path must not shuffle data')
+    #end def fail_shuffle
+
+    monkeypatch.setattr(statistics.np.random,'shuffle',fail_shuffle)
     intervals,counts = statistics.line_crossing_distribution(x)
     np.testing.assert_array_equal(intervals,[[0.,1.],[1.,2.]])
     np.testing.assert_array_equal(counts,[1,2])
     assert(statistics.lcd_peak(x)==pytest.approx(1.5))
+
+    def reverse(values):
+        values[:] = values[::-1]
+    #end def reverse
+
+    monkeypatch.setattr(statistics.np.random,'shuffle',reverse)
+    intervals,counts = statistics.line_crossing_distribution(x,nperm=2)
+    np.testing.assert_array_equal(intervals,[[0.,1.],[1.,2.]])
+    np.testing.assert_array_equal(counts,[1.,2.])
 #end def test_line_crossing_distribution_and_lcd_peak
 
 

From 489a5c320f252272e6983f2c99aca6e70fa9c9ae Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 17:04:12 -0400
Subject: [PATCH 05/10] nexus: docstrings/quad_peak

---
 nexus/nexus/statistics.py            | 65 ++++++++++++++++++++--------
 nexus/nexus/tests/test_statistics.py | 28 +++++++++---
 2 files changed, 71 insertions(+), 22 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index 0ae2908348..7971cc76b9 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -74,7 +74,7 @@ def _paired_real_arrays(x,y):
 
 
 def _real_vector(x,name):
-    """Return a finite real vector, flattening vector-shaped arrays."""
+    """Return a real vector, flattening vector-shaped arrays."""
     x = np.asarray(x)
     if np.iscomplexobj(x):
         msg = f'{name} must be real-valued'
@@ -89,9 +89,6 @@ def _real_vector(x,name):
     except (TypeError,ValueError):
         msg = f'{name} must be numeric'
         raise ValueError(msg) from None
-    if not np.all(np.isfinite(x)):
-        msg = f'{name} must contain only finite values'
-        raise ValueError(msg)
     return x
 #end def _real_vector
 
@@ -746,6 +743,10 @@ def series_stats(x,t_auto=None):
 
 
 def time_series_intervals(x,t=None):
+    """Return ordered intervals between adjacent time-series values.
+
+    If times are supplied, return their adjacent-pair midpoints as well.
+    """
     x = _real_vector(x,'data array')
     if len(x)<2:
         msg = 'data array must contain at least two values'
@@ -770,6 +771,7 @@ def time_series_intervals(x,t=None):
 
 
 def _int_dist_input(x1,x2=None):
+    """Normalize one interval matrix or paired lower and upper endpoints."""
     if x2 is not None:
         x1 = _real_vector(x1,'lower endpoints')
         x2 = _real_vector(x2,'upper endpoints')
@@ -807,6 +809,7 @@ def _int_dist_input(x1,x2=None):
 
 
 def interval_distribution(x1,x2=None):
+    """Return spans between interval edges and their overlap counts."""
     xi,si = _int_dist_input(x1,x2)
     xi = xi.ravel()
     si = si.ravel()
@@ -839,6 +842,7 @@ def interval_distribution(x1,x2=None):
 
 
 def plot_interval_dist(xi,ci,style='b.-'):
+    """Plot an interval distribution as a piecewise-constant curve."""
     import matplotlib.pyplot as plt
     xi,_ = _int_dist_input(xi)
     ci = _real_vector(ci,'interval counts')
@@ -856,6 +860,11 @@ def plot_interval_dist(xi,ci,style='b.-'):
 
 
 def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
+    """Return a representative location at the peak of an interval distribution.
+
+    ``method`` selects a peak-interval midpoint, random interior samples, or
+    a quadratic fit over the region at least ``peak_frac`` of the maximum.
+    """
     xi,_ = _int_dist_input(xi)
     ci = _real_vector(ci,'interval counts')
     if len(ci)!=len(xi):
@@ -869,8 +878,8 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
     except (TypeError,ValueError):
         msg = 'peak fraction must be a finite number'
         raise ValueError(msg) from None
-    if not np.isfinite(peak_frac):
-        msg = 'peak fraction must be a finite number'
+    if not np.isfinite(peak_frac) or not 0.<peak_frac<=1.:
+        msg = 'peak fraction must be in the interval (0,1]'
         raise ValueError(msg)
     if method=='interval_mid':
         cm = ci.max()
@@ -887,6 +896,7 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
     elif method=='quad_peak':
         imax = ci.argmax()
         cf = peak_frac*ci[imax]
+        peak_mean = xi[ci==ci[imax]].mean()
         # Locate the contiguous high-count region in interval coordinates.
         # ``imax`` indexes ``ci``, not the flattened endpoint array.
         i1 = imax
@@ -898,12 +908,14 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
         xp = xi[i1:i2+1].ravel()
         cp = np.repeat(ci[i1:i2+1],2)
         if len(np.unique(xp))<3:
-            xm = xi[ci==ci[imax]].mean()
+            # A single usable span cannot determine a quadratic peak.
+            xm = peak_mean
             cm = ci[imax]
         else:
             p = np.polyfit(xp,cp,2)
             if not np.isfinite(p[0]) or p[0]>=0.:
-                xm = xi[ci==ci[imax]].mean()
+                # A non-concave fit has no interior maximum.
+                xm = peak_mean
                 cm = ci[imax]
             else:
                 xm = -p[1]/(2*p[0])
@@ -925,9 +937,15 @@ def rolling_interval_dist_peak(
         window      = 10,
         step        = 5,
         method      = 'interval_mid',
+        peak_frac  = 0.5,
         ret_height  = False,
         ret_windows = False,
         ):
+    """Return interval-distribution peaks for overlapping input windows.
+
+    ``method`` accepts ``'interval_mid'``, ``'interval_rand'``, and
+    ``'quad_peak'``.  ``peak_frac`` is used by the quadratic method.
+    """
     for value,name in ((window,'window'),(step,'step')):
         if isinstance(value,(bool,np.bool_)) or not isinstance(
             value,(int,np.integer)
@@ -939,7 +957,8 @@ def rolling_interval_dist_peak(
         raise ValueError(msg)
     interval_mid  = method=='interval_mid'
     interval_rand = method=='interval_rand'
-    if not interval_mid and not interval_rand:
+    interval_quad = method=='quad_peak'
+    if not interval_mid and not interval_rand and not interval_quad:
         msg = f'method "{method}" is unrecognized'
         raise ValueError(msg)
     # map inputs to intervals
@@ -961,12 +980,9 @@ def rolling_interval_dist_peak(
         if len(ci)==0:
             msg = 'each rolling window must span a nonzero interval'
             raise ValueError(msg)
-        if interval_mid:
-            xm,cm = interval_dist_peak(xi,ci,height=True)
-        else:
-            xm,cm = interval_dist_peak(
-                xi,ci,method='interval_rand',height=True
-                )
+        xm,cm = interval_dist_peak(
+            xi,ci,method=method,peak_frac=peak_frac,height=True
+            )
         xp.append(xm)
         cp.append(cm)
     xp = np.array(xp)
@@ -985,6 +1001,11 @@ def rolling_interval_dist_peak(
 
 
 def line_crossing_distribution(x,nperm=0):
+    """Return the line-crossing distribution of a series or its permutations.
+
+    Positive ``nperm`` values average independently shuffled line-crossing
+    distributions without connecting successive permutations.
+    """
     x = _real_vector(x,'data array')
     if len(x)<2:
         msg = 'data array must contain at least two values'
@@ -1013,8 +1034,12 @@ def line_crossing_distribution(x,nperm=0):
 
 
 
-def lcd_peak(x,method='interval_mid',peak_frac=0.5):
-    xi,ci = line_crossing_distribution(x)
+def lcd_peak(x,method='interval_mid',peak_frac=0.5,nperm=0):
+    """Return a peak of a series line-crossing distribution.
+
+    ``nperm`` is passed to :func:`line_crossing_distribution`.
+    """
+    xi,ci = line_crossing_distribution(x,nperm=nperm)
     xp = interval_dist_peak(xi,ci,method=method,peak_frac=peak_frac)
     return xp
 #end def lcd_peak
@@ -1026,13 +1051,19 @@ def lcd_smooth(x,
                window = 10,
                step   = 5,
                method = 'interval_rand',
+               peak_frac = 0.5,
                ):
+    """Return rolling line-crossing-distribution peaks for a time series.
+
+    ``peak_frac`` is forwarded to the rolling quadratic-peak method.
+    """
     xi,ti = time_series_intervals(x,t)
     xp,windows = rolling_interval_dist_peak(
         xi,
         window      = window,
         step        = step,
         method      = method,
+        peak_frac   = peak_frac,
         ret_windows = True,
         )
     if t is None:
diff --git a/nexus/nexus/tests/test_statistics.py b/nexus/nexus/tests/test_statistics.py
index 6f23b62683..c6d605ffe6 100644
--- a/nexus/nexus/tests/test_statistics.py
+++ b/nexus/nexus/tests/test_statistics.py
@@ -398,8 +398,6 @@ def test_interval_and_lcd_input_validation():
     """Check diagnostics for malformed interval and LCD inputs."""
     with pytest.raises(ValueError,match=r'data array must contain at least two values'):
         statistics.time_series_intervals([1.])
-    with pytest.raises(ValueError,match=r'data array must contain only finite values'):
-        statistics.time_series_intervals([1.,np.nan])
     with pytest.raises(ValueError,match=r'time array must have the same length'):
         statistics.time_series_intervals([1.,2.],[0.])
     with pytest.raises(ValueError,match=r'time array must be real-valued'):
@@ -418,14 +416,14 @@ def test_interval_and_lcd_input_validation():
     intervals = np.array([[0.,1.],[1.,2.]])
     for counts,message in [
         ([1.],r'counts must have the same length'),
-        ([1.,np.nan],r'counts must contain only finite values'),
         ]:
         with pytest.raises(ValueError,match=message):
             statistics.interval_dist_peak(intervals,counts)
     with pytest.raises(ValueError,match=r'peak method must be a string'):
         statistics.interval_dist_peak(intervals,[1.,2.],method=1)
-    with pytest.raises(ValueError,match=r'peak fraction must be a finite number'):
-        statistics.interval_dist_peak(intervals,[1.,2.],peak_frac=np.nan)
+    for peak_frac in (np.nan,0.,-1.,1.1):
+        with pytest.raises(ValueError,match=r'peak fraction must be in the interval'):
+            statistics.interval_dist_peak(intervals,[1.,2.],peak_frac=peak_frac)
 
     for window,step,message in [
         (0,1,r'window must be a positive integer'),
@@ -508,6 +506,15 @@ def test_interval_distribution_and_peak(monkeypatch):
     assert(np.isfinite(quadratic_height))
     assert(quadratic_height>0.)
 
+    fallback_peak,fallback_height = statistics.interval_dist_peak(
+        np.array([[0.,2.]]),
+        np.array([4.]),
+        method='quad_peak',
+        height=True,
+        )
+    assert(fallback_peak==pytest.approx(1.))
+    assert(fallback_height==4.)
+
     monkeypatch.setattr(
         statistics.np.random,
         'uniform',
@@ -542,6 +549,16 @@ def test_rolling_interval_dist_peak_and_lcd_smooth():
     np.testing.assert_array_equal(heights,[2,2])
     assert(windows==[(0,2),(2,4)])
 
+    quadratic_peaks,quadratic_heights = statistics.rolling_interval_dist_peak(
+        np.array([[1.,4.],[2.,5.],[3.,6.]]),
+        window=3,
+        step=1,
+        method='quad_peak',
+        ret_height=True,
+        )
+    np.testing.assert_allclose(quadratic_peaks,[3.5])
+    assert(quadratic_heights[0]>0.)
+
     x = np.array([0.,2.,1.,3.])
     t = np.array([0.,1.,3.,6.])
     smooth,times = statistics.lcd_smooth(x,t,window=2,step=1,method='interval_mid')
@@ -577,6 +594,7 @@ def reverse(values):
     intervals,counts = statistics.line_crossing_distribution(x,nperm=2)
     np.testing.assert_array_equal(intervals,[[0.,1.],[1.,2.]])
     np.testing.assert_array_equal(counts,[1.,2.])
+    assert(statistics.lcd_peak(x,nperm=2)==pytest.approx(1.5))
 #end def test_line_crossing_distribution_and_lcd_peak
 
 

From 3a0adf8b851bad7f9502431a866a0a44986d16af Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 17:08:43 -0400
Subject: [PATCH 06/10] nexus: better docs, expand tests

---
 nexus/nexus/statistics.py            | 227 ++++++++++++++++++++++-----
 nexus/nexus/tests/test_statistics.py |  60 +++++++
 2 files changed, 249 insertions(+), 38 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index 7971cc76b9..22dee4306a 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -745,7 +745,22 @@ def series_stats(x,t_auto=None):
 def time_series_intervals(x,t=None):
     """Return ordered intervals between adjacent time-series values.
 
-    If times are supplied, return their adjacent-pair midpoints as well.
+    Parameters
+    ----------
+    x : array_like
+        Real one-dimensional series with at least two values. Vector-shaped
+        arrays are flattened.
+
+    t : array_like, optional
+        Real times paired with ``x``. If supplied, must have the same length.
+
+    Returns
+    -------
+    xi : ndarray
+        ``(len(x)-1, 2)`` array of ordered adjacent-value intervals.
+
+    ti : ndarray or None
+        Adjacent-pair time midpoints, or ``None`` when ``t`` is omitted.
     """
     x = _real_vector(x,'data array')
     if len(x)<2:
@@ -771,7 +786,10 @@ def time_series_intervals(x,t=None):
 
 
 def _int_dist_input(x1,x2=None):
-    """Normalize one interval matrix or paired lower and upper endpoints."""
+    """Normalize one interval matrix or paired lower and upper endpoints.
+
+    Returns ordered endpoint pairs and corresponding ``+1/-1`` edge signs.
+    """
     if x2 is not None:
         x1 = _real_vector(x1,'lower endpoints')
         x2 = _real_vector(x2,'upper endpoints')
@@ -809,7 +827,25 @@ def _int_dist_input(x1,x2=None):
 
 
 def interval_distribution(x1,x2=None):
-    """Return spans between interval edges and their overlap counts."""
+    """Return spans between interval edges and their overlap counts.
+
+    Parameters
+    ----------
+    x1 : array_like
+        ``(n,2)`` ordered interval array, or lower endpoints when ``x2`` is
+        supplied.
+
+    x2 : array_like, optional
+        Upper endpoints paired with ``x1``.
+
+    Returns
+    -------
+    xi : ndarray
+        Consecutive spans between sorted unique interval endpoints.
+
+    ci : ndarray
+        Number of input intervals overlapping each span in ``xi``.
+    """
     xi,si = _int_dist_input(x1,x2)
     xi = xi.ravel()
     si = si.ravel()
@@ -842,7 +878,19 @@ def interval_distribution(x1,x2=None):
 
 
 def plot_interval_dist(xi,ci,style='b.-'):
-    """Plot an interval distribution as a piecewise-constant curve."""
+    """Plot an interval distribution as a piecewise-constant curve.
+
+    Parameters
+    ----------
+    xi : array_like
+        ``(n,2)`` interval-distribution spans.
+
+    ci : array_like
+        Counts paired with ``xi``.
+
+    style : str, optional
+        Matplotlib style specification for the distribution line.
+    """
     import matplotlib.pyplot as plt
     xi,_ = _int_dist_input(xi)
     ci = _real_vector(ci,'interval counts')
@@ -862,8 +910,31 @@ def plot_interval_dist(xi,ci,style='b.-'):
 def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
     """Return a representative location at the peak of an interval distribution.
 
-    ``method`` selects a peak-interval midpoint, random interior samples, or
-    a quadratic fit over the region at least ``peak_frac`` of the maximum.
+    Parameters
+    ----------
+    xi : array_like
+        ``(n,2)`` interval-distribution spans.
+
+    ci : array_like
+        Counts paired with ``xi``.
+
+    method : {'interval_mid', 'interval_rand', 'quad_peak'}, optional
+        Peak estimator. The first averages all maximum-count intervals, the
+        second averages random interior samples of those intervals, and the
+        third fits each separated high-count peak region quadratically.
+
+    peak_frac : float, optional
+        Fraction of the maximum count retained for each quadratic-fit region.
+        Must lie in ``(0,1]``.
+
+    height : bool, optional
+        If true, return the peak location and its estimated height.
+
+    Returns
+    -------
+    peak : float or (float, float)
+        Peak location, optionally followed by peak height. Separated
+        equal-height modes are averaged.
     """
     xi,_ = _int_dist_input(xi)
     ci = _real_vector(ci,'interval counts')
@@ -894,33 +965,38 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
         x     = xmid + (u-0.5)*dx/2
         xm    = x.mean()
     elif method=='quad_peak':
-        imax = ci.argmax()
-        cf = peak_frac*ci[imax]
-        peak_mean = xi[ci==ci[imax]].mean()
-        # Locate the contiguous high-count region in interval coordinates.
-        # ``imax`` indexes ``ci``, not the flattened endpoint array.
-        i1 = imax
-        while i1>0 and ci[i1-1]>=cf:
-            i1 -= 1
-        i2 = imax
-        while i2+1<len(ci) and ci[i2+1]>=cf:
-            i2 += 1
-        xp = xi[i1:i2+1].ravel()
-        cp = np.repeat(ci[i1:i2+1],2)
-        if len(np.unique(xp))<3:
-            # A single usable span cannot determine a quadratic peak.
-            xm = peak_mean
-            cm = ci[imax]
-        else:
-            p = np.polyfit(xp,cp,2)
-            if not np.isfinite(p[0]) or p[0]>=0.:
-                # A non-concave fit has no interior maximum.
-                xm = peak_mean
-                cm = ci[imax]
+        cm = ci.max()
+        cf = peak_frac*cm
+        high = ci>=cf
+        edges = np.flatnonzero(np.diff(np.r_[False,high,False]))
+        xpeaks = []
+        cpeaks = []
+        for i1,i2 in zip(edges[::2],edges[1::2]-1):
+            ci_region = ci[i1:i2+1]
+            if ci_region.max()!=cm:
+                continue
+            xi_region = xi[i1:i2+1]
+            peak_mean = xi_region[ci_region==cm].mean()
+            xp = xi_region.ravel()
+            cp = np.repeat(ci_region,2)
+            if len(np.unique(xp))<3:
+                # A single usable span cannot determine a quadratic peak.
+                xp = peak_mean
+                cp = cm
             else:
-                xm = -p[1]/(2*p[0])
-                xm = np.clip(xm,xp.min(),xp.max())
-                cm = np.polyval(p,xm)
+                p = np.polyfit(xp,cp,2)
+                if not np.isfinite(p[0]) or p[0]>=0.:
+                    # A non-concave fit has no interior maximum.
+                    xp = peak_mean
+                    cp = cm
+                else:
+                    xp = -p[1]/(2*p[0])
+                    xp = np.clip(xp,xi_region.min(),xi_region.max())
+                    cp = np.polyval(p,xp)
+            xpeaks.append(xp)
+            cpeaks.append(cp)
+        xm = np.mean(xpeaks)
+        cm = np.mean(cpeaks)
     else:
         raise ValueError(f'unrecognized int. dist. max method: "{method}"')
     if not height:
@@ -943,8 +1019,41 @@ def rolling_interval_dist_peak(
         ):
     """Return interval-distribution peaks for overlapping input windows.
 
-    ``method`` accepts ``'interval_mid'``, ``'interval_rand'``, and
-    ``'quad_peak'``.  ``peak_frac`` is used by the quadratic method.
+    Parameters
+    ----------
+    x1 : array_like
+        ``(n,2)`` interval array, or lower endpoints when ``x2`` is given.
+
+    x2 : array_like, optional
+        Upper endpoints paired with ``x1``.
+
+    window : int, optional
+        Number of input intervals in each rolling distribution.
+
+    step : int, optional
+        Number of intervals between successive window starts. It must not
+        exceed ``window``.
+
+    method : {'interval_mid', 'interval_rand', 'quad_peak'}, optional
+        Peak estimator passed to :func:`interval_dist_peak`.
+
+    peak_frac : float, optional
+        Quadratic peak-region threshold passed to :func:`interval_dist_peak`.
+
+    ret_height : bool, optional
+        Include a peak-height array in the returned tuple.
+
+    ret_windows : bool, optional
+        Include ``(start, stop)`` bounds for each returned window.
+
+    Returns
+    -------
+    result : tuple
+        A tuple beginning with the peak-location array. When requested, it
+        then contains the peak-height array and/or window-bound list, in that
+        order. Windows advance by ``step``; a final window ending at the last
+        input interval is appended when the regular sequence does not reach
+        it exactly.
     """
     for value,name in ((window,'window'),(step,'step')):
         if isinstance(value,(bool,np.bool_)) or not isinstance(
@@ -1003,8 +1112,23 @@ def rolling_interval_dist_peak(
 def line_crossing_distribution(x,nperm=0):
     """Return the line-crossing distribution of a series or its permutations.
 
-    Positive ``nperm`` values average independently shuffled line-crossing
-    distributions without connecting successive permutations.
+    Parameters
+    ----------
+    x : array_like
+        Real one-dimensional series with at least two values.
+
+    nperm : int, optional
+        Number of independently shuffled series to average. Zero evaluates
+        the input series directly.
+
+    Returns
+    -------
+    xi : ndarray
+        Line-crossing distribution spans.
+
+    ci : ndarray
+        Crossing counts, averaged over permutations when ``nperm`` is
+        positive. Permutations are never connected to one another.
     """
     x = _real_vector(x,'data array')
     if len(x)<2:
@@ -1037,7 +1161,19 @@ def line_crossing_distribution(x,nperm=0):
 def lcd_peak(x,method='interval_mid',peak_frac=0.5,nperm=0):
     """Return a peak of a series line-crossing distribution.
 
-    ``nperm`` is passed to :func:`line_crossing_distribution`.
+    Parameters
+    ----------
+    x : array_like
+        Time series supplied to :func:`line_crossing_distribution`.
+
+    method, peak_frac, nperm
+        Options forwarded to :func:`interval_dist_peak` and
+        :func:`line_crossing_distribution`.
+
+    Returns
+    -------
+    float
+        Estimated line-crossing-distribution peak.
     """
     xi,ci = line_crossing_distribution(x,nperm=nperm)
     xp = interval_dist_peak(xi,ci,method=method,peak_frac=peak_frac)
@@ -1055,7 +1191,22 @@ def lcd_smooth(x,
                ):
     """Return rolling line-crossing-distribution peaks for a time series.
 
-    ``peak_frac`` is forwarded to the rolling quadratic-peak method.
+    Parameters
+    ----------
+    x : array_like
+        Time series to smooth with rolling line-crossing peaks.
+
+    t : array_like, optional
+        Times paired with ``x``.
+
+    window, step, method, peak_frac
+        Options forwarded to :func:`rolling_interval_dist_peak`.
+
+    Returns
+    -------
+    peaks : ndarray or (ndarray, ndarray)
+        Rolling peak locations, optionally paired with their mean window
+        times when ``t`` is supplied.
     """
     xi,ti = time_series_intervals(x,t)
     xp,windows = rolling_interval_dist_peak(
diff --git a/nexus/nexus/tests/test_statistics.py b/nexus/nexus/tests/test_statistics.py
index c6d605ffe6..2ae0e64a87 100644
--- a/nexus/nexus/tests/test_statistics.py
+++ b/nexus/nexus/tests/test_statistics.py
@@ -476,6 +476,30 @@ def test_interval_distribution_and_peak(monkeypatch):
     np.testing.assert_array_equal(column_intervals,expected_intervals)
     np.testing.assert_array_equal(column_counts,counts)
 
+    touching_intervals,touching_counts = statistics.interval_distribution(
+        np.array([[1.,2.],[2.,3.],[3.,4.]])
+        )
+    np.testing.assert_array_equal(
+        touching_intervals,
+        [[1.,2.],[2.,3.],[3.,4.]],
+        )
+    np.testing.assert_array_equal(touching_counts,[1,1,1])
+
+    repeated_intervals,repeated_counts = statistics.interval_distribution(
+        np.array([[1.,1.],[1.,2.],[2.,2.]])
+        )
+    np.testing.assert_array_equal(repeated_intervals,[[1.,2.]])
+    np.testing.assert_array_equal(repeated_counts,[1])
+
+    irregular_intervals,irregular_counts = statistics.interval_distribution(
+        np.array([[0.,10.],[.5,.75]])
+        )
+    np.testing.assert_array_equal(
+        irregular_intervals,
+        [[0.,.5],[.5,.75],[.75,10.]],
+        )
+    np.testing.assert_array_equal(irregular_counts,[1,2,1])
+
     peak_intervals = np.array([[0.,1.],[1.,2.],[2.,3.]])
     peak_counts    = np.array([1,3,3])
     peak,height = statistics.interval_dist_peak(
@@ -515,6 +539,27 @@ def test_interval_distribution_and_peak(monkeypatch):
     assert(fallback_peak==pytest.approx(1.))
     assert(fallback_height==4.)
 
+    multimodal_intervals = np.array([[0.,1.],[1.,2.],[2.,3.],[3.,4.]])
+    multimodal_counts    = np.array([4.,1.,1.,4.])
+    assert(
+        statistics.interval_dist_peak(
+            multimodal_intervals,multimodal_counts,method='quad_peak'
+            )
+        ==pytest.approx(2.)
+        )
+
+    monkeypatch.setattr(
+        statistics.np,
+        'polyfit',
+        lambda x,y,degree: np.array([1.,0.,0.]),
+        )
+    assert(
+        statistics.interval_dist_peak(
+            quadratic_intervals,quadratic_counts,method='quad_peak'
+            )
+        ==pytest.approx(2.5)
+        )
+
     monkeypatch.setattr(
         statistics.np.random,
         'uniform',
@@ -595,6 +640,21 @@ def reverse(values):
     np.testing.assert_array_equal(intervals,[[0.,1.],[1.,2.]])
     np.testing.assert_array_equal(counts,[1.,2.])
     assert(statistics.lcd_peak(x,nperm=2)==pytest.approx(1.5))
+
+    permutations = [
+        np.array([0.,1.,3.,6.]),
+        np.array([0.,3.,1.,6.]),
+        ]
+    def set_permutation(values):
+        values[:] = permutations.pop(0)
+    #end def set_permutation
+
+    monkeypatch.setattr(statistics.np.random,'shuffle',set_permutation)
+    intervals,counts = statistics.line_crossing_distribution(
+        np.array([0.,1.,3.,6.]),nperm=2
+        )
+    np.testing.assert_array_equal(intervals,[[0.,1.],[1.,3.],[3.,6.]])
+    np.testing.assert_array_equal(counts,[1.,2.,1.])
 #end def test_line_crossing_distribution_and_lcd_peak
 
 

From ba2f124db70dcc4fe4e7ce288bb0b53ab4586c9f Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 17:12:27 -0400
Subject: [PATCH 07/10] nexus: width-weighting

---
 nexus/nexus/statistics.py            | 55 ++++++++++++++++++++++++----
 nexus/nexus/tests/test_statistics.py | 37 +++++++++++++++++++
 2 files changed, 85 insertions(+), 7 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index 22dee4306a..2580157ba9 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -845,6 +845,13 @@ def interval_distribution(x1,x2=None):
 
     ci : ndarray
         Number of input intervals overlapping each span in ``xi``.
+
+    Notes
+    -----
+    Each row of ``xi`` denotes the open span between consecutive unique
+    endpoints. Endpoint membership is not counted separately: touching
+    intervals occupy adjacent spans, and zero-width intervals contribute no
+    span to the returned distribution.
     """
     xi,si = _int_dist_input(x1,x2)
     xi = xi.ravel()
@@ -907,7 +914,8 @@ def plot_interval_dist(xi,ci,style='b.-'):
 
 
 
-def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
+def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False,
+                       quad_weighting='endpoint'):
     """Return a representative location at the peak of an interval distribution.
 
     Parameters
@@ -930,6 +938,12 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
     height : bool, optional
         If true, return the peak location and its estimated height.
 
+    quad_weighting : {'endpoint', 'width'}, optional
+        Weighting used only by ``'quad_peak'``. ``'endpoint'`` gives every
+        duplicated interval endpoint equal fit weight. ``'width'`` weights
+        each endpoint by the square root of its interval width, making the
+        least-squares objective proportional to interval width.
+
     Returns
     -------
     peak : float or (float, float)
@@ -952,6 +966,11 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
     if not np.isfinite(peak_frac) or not 0.<peak_frac<=1.:
         msg = 'peak fraction must be in the interval (0,1]'
         raise ValueError(msg)
+    if not isinstance(quad_weighting,str) or quad_weighting not in (
+        'endpoint','width'
+        ):
+        msg = 'quadratic weighting must be "endpoint" or "width"'
+        raise ValueError(msg)
     if method=='interval_mid':
         cm = ci.max()
         xm = xi[ci==cm].mean()
@@ -979,12 +998,23 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False):
             peak_mean = xi_region[ci_region==cm].mean()
             xp = xi_region.ravel()
             cp = np.repeat(ci_region,2)
+            weights = None
+            if quad_weighting=='width':
+                widths = xi_region[:,1]-xi_region[:,0]
+                weights = np.repeat(np.sqrt(widths),2)
+                nonzero = weights>0.
+                xp = xp[nonzero]
+                cp = cp[nonzero]
+                weights = weights[nonzero]
             if len(np.unique(xp))<3:
                 # A single usable span cannot determine a quadratic peak.
                 xp = peak_mean
                 cp = cm
             else:
-                p = np.polyfit(xp,cp,2)
+                if weights is None:
+                    p = np.polyfit(xp,cp,2)
+                else:
+                    p = np.polyfit(xp,cp,2,w=weights)
                 if not np.isfinite(p[0]) or p[0]>=0.:
                     # A non-concave fit has no interior maximum.
                     xp = peak_mean
@@ -1014,6 +1044,7 @@ def rolling_interval_dist_peak(
         step        = 5,
         method      = 'interval_mid',
         peak_frac  = 0.5,
+        quad_weighting = 'endpoint',
         ret_height  = False,
         ret_windows = False,
         ):
@@ -1040,6 +1071,9 @@ def rolling_interval_dist_peak(
     peak_frac : float, optional
         Quadratic peak-region threshold passed to :func:`interval_dist_peak`.
 
+    quad_weighting : {'endpoint', 'width'}, optional
+        Quadratic-fit weighting passed to :func:`interval_dist_peak`.
+
     ret_height : bool, optional
         Include a peak-height array in the returned tuple.
 
@@ -1090,7 +1124,8 @@ def rolling_interval_dist_peak(
             msg = 'each rolling window must span a nonzero interval'
             raise ValueError(msg)
         xm,cm = interval_dist_peak(
-            xi,ci,method=method,peak_frac=peak_frac,height=True
+            xi,ci,method=method,peak_frac=peak_frac,height=True,
+            quad_weighting=quad_weighting,
             )
         xp.append(xm)
         cp.append(cm)
@@ -1158,7 +1193,8 @@ def line_crossing_distribution(x,nperm=0):
 
 
 
-def lcd_peak(x,method='interval_mid',peak_frac=0.5,nperm=0):
+def lcd_peak(x,method='interval_mid',peak_frac=0.5,nperm=0,
+             quad_weighting='endpoint'):
     """Return a peak of a series line-crossing distribution.
 
     Parameters
@@ -1166,7 +1202,7 @@ def lcd_peak(x,method='interval_mid',peak_frac=0.5,nperm=0):
     x : array_like
         Time series supplied to :func:`line_crossing_distribution`.
 
-    method, peak_frac, nperm
+    method, peak_frac, quad_weighting, nperm
         Options forwarded to :func:`interval_dist_peak` and
         :func:`line_crossing_distribution`.
 
@@ -1176,7 +1212,10 @@ def lcd_peak(x,method='interval_mid',peak_frac=0.5,nperm=0):
         Estimated line-crossing-distribution peak.
     """
     xi,ci = line_crossing_distribution(x,nperm=nperm)
-    xp = interval_dist_peak(xi,ci,method=method,peak_frac=peak_frac)
+    xp = interval_dist_peak(
+        xi,ci,method=method,peak_frac=peak_frac,
+        quad_weighting=quad_weighting,
+        )
     return xp
 #end def lcd_peak
 
@@ -1188,6 +1227,7 @@ def lcd_smooth(x,
                step   = 5,
                method = 'interval_rand',
                peak_frac = 0.5,
+               quad_weighting = 'endpoint',
                ):
     """Return rolling line-crossing-distribution peaks for a time series.
 
@@ -1199,7 +1239,7 @@ def lcd_smooth(x,
     t : array_like, optional
         Times paired with ``x``.
 
-    window, step, method, peak_frac
+    window, step, method, peak_frac, quad_weighting
         Options forwarded to :func:`rolling_interval_dist_peak`.
 
     Returns
@@ -1215,6 +1255,7 @@ def lcd_smooth(x,
         step        = step,
         method      = method,
         peak_frac   = peak_frac,
+        quad_weighting = quad_weighting,
         ret_windows = True,
         )
     if t is None:
diff --git a/nexus/nexus/tests/test_statistics.py b/nexus/nexus/tests/test_statistics.py
index 2ae0e64a87..203328d4b1 100644
--- a/nexus/nexus/tests/test_statistics.py
+++ b/nexus/nexus/tests/test_statistics.py
@@ -576,10 +576,47 @@ def test_interval_distribution_and_peak(monkeypatch):
 
     with pytest.raises(ValueError,match=r'unrecognized int. dist. max method'):
         statistics.interval_dist_peak(peak_intervals,peak_counts,method='invalid')
+    with pytest.raises(ValueError,match=r'quadratic weighting'):
+        statistics.interval_dist_peak(
+            peak_intervals,peak_counts,quad_weighting='invalid'
+            )
 #end def test_interval_distribution_and_peak
 
 
 
+def test_quad_peak_width_weighting(monkeypatch):
+    """Check width-based quadratic weights and convenience-API forwarding."""
+    intervals = np.array(
+        [[0.,1.],[1.,3.],[3.,6.],[6.,10.],[10.,15.]]
+        )
+    counts = np.array([1.,2.,3.,2.,1.])
+    polyfit = statistics.np.polyfit
+    weights = []
+
+    def capture_polyfit(x,y,degree,**kwargs):
+        weights.append(kwargs.get('w'))
+        return polyfit(x,y,degree,**kwargs)
+    #end def capture_polyfit
+
+    monkeypatch.setattr(statistics.np,'polyfit',capture_polyfit)
+    peak = statistics.interval_dist_peak(
+        intervals,counts,method='quad_peak',quad_weighting='width'
+        )
+    assert(np.isfinite(peak))
+    np.testing.assert_allclose(weights[0],np.sqrt([2.,2.,3.,3.,4.,4.]))
+
+    rolling_peak, = statistics.rolling_interval_dist_peak(
+        intervals,
+        window=5,
+        step=1,
+        method='quad_peak',
+        quad_weighting='width',
+        )
+    assert(np.isfinite(rolling_peak[0]))
+#end def test_quad_peak_width_weighting
+
+
+
 def test_rolling_interval_dist_peak_and_lcd_smooth():
     """Check rolling peak locations, heights, window bounds, and times."""
     intervals = np.array([[0.,2.],[1.,3.],[2.,4.],[3.,5.]])

From 33d5074904146955f2f01a1625caaf6912f9656d Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 17:16:15 -0400
Subject: [PATCH 08/10] nexus: optimize speed

---
 nexus/nexus/statistics.py | 39 +++++++++++++++------------------------
 1 file changed, 15 insertions(+), 24 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index 2580157ba9..b9dc23fcf9 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -854,31 +854,22 @@ def interval_distribution(x1,x2=None):
     span to the returned distribution.
     """
     xi,si = _int_dist_input(x1,x2)
-    xi = xi.ravel()
-    si = si.ravel()
     # organize by edge order
-    order = xi.argsort()
-    xi = xi[order]
-    si = si[order]
-    # make the interval counting distribution
-    cd = {}
-    n=0
-    for s,xv in zip(si,xi):
-        n += s
-        cd[xv] = n
-    # first organize into sorted point/edge arrays
-    x = []
-    c = []
-    for xv in sorted(cd.keys()):
-        x.append(xv)
-        c.append(cd[xv])
-    x = np.array(x)
-    c = np.array(c)
-    # next into interval array
-    xi = np.zeros((len(x)-1,2),dtype=x.dtype)
-    xi[:,0] = x[:-1]
-    xi[:,1] = x[1:]
-    ci = c[:-1].copy()
+    edges = xi.ravel()
+    signs = si.ravel()
+    order = edges.argsort()
+    edges = edges[order]
+    signs = signs[order]
+
+    # Combine all coincident edges before accumulating their net change.
+    # This vectorized sweep avoids one Python dictionary entry and two list
+    # appends per edge while preserving the span counts between unique edges.
+    values,starts = np.unique(edges,return_index=True)
+    counts = np.cumsum(np.add.reduceat(signs,starts))
+    xi = np.empty((len(values)-1,2),dtype=values.dtype)
+    xi[:,0] = values[:-1]
+    xi[:,1] = values[1:]
+    ci = counts[:-1]
     return xi,ci
 #end def interval_distribution
 

From d5a86a8e057ed19636a88f9dcf69a388b998feee Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Fri, 11 Sep 2026 17:23:41 -0400
Subject: [PATCH 09/10] nexus: format

---
 nexus/nexus/statistics.py | 103 ++++++++++++++++++++++++++------------
 1 file changed, 70 insertions(+), 33 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index b9dc23fcf9..5fcbc86e87 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -248,7 +248,13 @@ def theil_sen_stoch_reblock(x,y):
 #end def theil_sen_stoch_reblock
 
 
-def reblocked_autocorr_time(x,min_blocks=10,plot=False,show=False):
+def reblocked_autocorr_time(
+        x,
+        min_blocks = 10,
+        *,
+        plot       = False,
+        show       = False,
+        ):
     """Estimate autocorrelation time from the growth of blocked errors.
 
     This estimator currently overestimates the autocorrelation times in a 
@@ -393,7 +399,7 @@ def reblocked_autocorr_time(x,min_blocks=10,plot=False,show=False):
 
 
 
-def acf_autocorr_time(x,reliability=False):
+def acf_autocorr_time(x,*,reliability=False):
     """Estimate autocorrelation time from a windowed sample ACF.
 
     Best for long chains.  Generally prefer the Geyer method.
@@ -505,7 +511,13 @@ def acf_autocorr_time(x,reliability=False):
 
 
 
-def geyer_ims_autocorr_time(x,c=5.0,reliability=False,acf_fallback=True):
+def geyer_ims_autocorr_time(
+        x,
+        c            = 5.0,
+        *,
+        reliability  = False,
+        acf_fallback = True,
+        ):
     """Estimate integrated autocorrelation time with Geyer's IMS method.
 
     This is the single best autocorrelation estimator.
@@ -637,7 +649,7 @@ def geyer_ims_autocorr_time(x,c=5.0,reliability=False,acf_fallback=True):
 
 
 
-def autocorr_time(x,reliability=False):
+def autocorr_time(x,*,reliability=False):
     """Conservatively combine autocorrelation-time estimates.
 
     The ACF and Geyer initial-monotone-sequence probe the correlation
@@ -875,7 +887,11 @@ def interval_distribution(x1,x2=None):
 
 
 
-def plot_interval_dist(xi,ci,style='b.-'):
+def plot_interval_dist(
+        xi,
+        ci,
+        style = 'b.-',
+        ):
     """Plot an interval distribution as a piecewise-constant curve.
 
     Parameters
@@ -905,8 +921,15 @@ def plot_interval_dist(xi,ci,style='b.-'):
 
 
 
-def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False,
-                       quad_weighting='endpoint'):
+def interval_dist_peak(
+        xi,
+        ci,
+        method         = 'interval_mid',
+        peak_frac      = 0.5,
+        *,
+        height         = False,
+        quad_weighting = 'endpoint',
+        ):
     """Return a representative location at the peak of an interval distribution.
 
     Parameters
@@ -1030,14 +1053,15 @@ def interval_dist_peak(xi,ci,method='interval_mid',peak_frac=0.5,height=False,
 
 def rolling_interval_dist_peak(
         x1,
-        x2          = None,
-        window      = 10,
-        step        = 5,
-        method      = 'interval_mid',
-        peak_frac  = 0.5,
+        x2             = None,
+        window         = 10,
+        step           = 5,
+        method         = 'interval_mid',
+        peak_frac      = 0.5,
+        *,
         quad_weighting = 'endpoint',
-        ret_height  = False,
-        ret_windows = False,
+        ret_height     = False,
+        ret_windows    = False,
         ):
     """Return interval-distribution peaks for overlapping input windows.
 
@@ -1115,8 +1139,12 @@ def rolling_interval_dist_peak(
             msg = 'each rolling window must span a nonzero interval'
             raise ValueError(msg)
         xm,cm = interval_dist_peak(
-            xi,ci,method=method,peak_frac=peak_frac,height=True,
-            quad_weighting=quad_weighting,
+            xi,
+            ci,
+            method         = method,
+            peak_frac      = peak_frac,
+            height         = True,
+            quad_weighting = quad_weighting,
             )
         xp.append(xm)
         cp.append(cm)
@@ -1184,8 +1212,13 @@ def line_crossing_distribution(x,nperm=0):
 
 
 
-def lcd_peak(x,method='interval_mid',peak_frac=0.5,nperm=0,
-             quad_weighting='endpoint'):
+def lcd_peak(
+        x,
+        method         = 'interval_mid',
+        peak_frac      = 0.5,
+        nperm          = 0,
+        quad_weighting = 'endpoint',
+        ):
     """Return a peak of a series line-crossing distribution.
 
     Parameters
@@ -1204,22 +1237,26 @@ def lcd_peak(x,method='interval_mid',peak_frac=0.5,nperm=0,
     """
     xi,ci = line_crossing_distribution(x,nperm=nperm)
     xp = interval_dist_peak(
-        xi,ci,method=method,peak_frac=peak_frac,
-        quad_weighting=quad_weighting,
+        xi,
+        ci,
+        method         = method,
+        peak_frac      = peak_frac,
+        quad_weighting = quad_weighting,
         )
     return xp
 #end def lcd_peak
 
 
 
-def lcd_smooth(x,
-               t      = None,
-               window = 10,
-               step   = 5,
-               method = 'interval_rand',
-               peak_frac = 0.5,
-               quad_weighting = 'endpoint',
-               ):
+def lcd_smooth(
+        x,
+        t              = None,
+        window         = 10,
+        step           = 5,
+        method         = 'interval_rand',
+        peak_frac      = 0.5,
+        quad_weighting = 'endpoint',
+        ):
     """Return rolling line-crossing-distribution peaks for a time series.
 
     Parameters
@@ -1242,12 +1279,12 @@ def lcd_smooth(x,
     xi,ti = time_series_intervals(x,t)
     xp,windows = rolling_interval_dist_peak(
         xi,
-        window      = window,
-        step        = step,
-        method      = method,
-        peak_frac   = peak_frac,
+        window         = window,
+        step           = step,
+        method         = method,
+        peak_frac      = peak_frac,
         quad_weighting = quad_weighting,
-        ret_windows = True,
+        ret_windows    = True,
         )
     if t is None:
         return xp

From cd9e2164e2da18e456ac5f66824a8c561a51aad7 Mon Sep 17 00:00:00 2001
From: Jaron Krogel <krogeljt@ornl.gov>
Date: Mon, 14 Sep 2026 09:50:35 -0400
Subject: [PATCH 10/10] nexus: docs and constant case

---
 nexus/nexus/statistics.py            | 88 ++++++++++++++++++++++++++--
 nexus/nexus/tests/test_statistics.py | 45 +++++++++++++-
 2 files changed, 126 insertions(+), 7 deletions(-)

diff --git a/nexus/nexus/statistics.py b/nexus/nexus/statistics.py
index bab2cb5241..ef8fd4e05e 100644
--- a/nexus/nexus/statistics.py
+++ b/nexus/nexus/statistics.py
@@ -775,6 +775,8 @@ def series_stats(x,t_auto=None):
 # emphasize locally stable, equilibrium-like portions of a fluctuating     #
 # series.  Rolling versions track this center over time, while related     #
 # utilities support broader interval-distribution analysis.                #
+############################################################################
+
 
 def time_series_intervals(x,t=None):
     """Return ordered intervals between adjacent time-series values.
@@ -860,7 +862,37 @@ def _int_dist_input(x1,x2=None):
 
 
 
-def interval_distribution(x1,x2=None):
+def _perturb_constant_intervals(xi,perturb_const):
+    """Expand constant intervals by a fixed number of floating-point steps."""
+    if isinstance(perturb_const,(bool,np.bool_)) or not isinstance(
+        perturb_const,(int,np.integer)
+        ) or perturb_const<0:
+        msg = 'constant perturbation must be a nonnegative integer'
+        raise ValueError(msg)
+    if perturb_const==0:
+        return xi
+    constant = xi[:,0]==xi[:,1]
+    if not constant.any():
+        return xi
+    xi = xi.copy()
+    lower = xi[constant,0]
+    upper = xi[constant,1]
+    for n in range(perturb_const):
+        lower = np.nextafter(lower,-np.inf)
+        upper = np.nextafter(upper,np.inf)
+    xi[constant,0] = lower
+    xi[constant,1] = upper
+    return xi
+#end def _perturb_constant_intervals
+
+
+
+def interval_distribution(
+        x1,
+        x2             = None,
+        *,
+        perturb_const  = 1,
+        ):
     """Return spans between interval edges and their overlap counts.
 
     Parameters
@@ -872,6 +904,12 @@ def interval_distribution(x1,x2=None):
     x2 : array_like, optional
         Upper endpoints paired with ``x1``.
 
+    perturb_const : int, optional
+        Number of floating-point steps used to expand each zero-width input
+        interval by equal step counts toward negative and positive infinity.
+        One gives a deterministic ULP-scale representation of constant
+        intervals; zero leaves them unexpanded.
+
     Returns
     -------
     xi : ndarray
@@ -884,10 +922,13 @@ def interval_distribution(x1,x2=None):
     -----
     Each row of ``xi`` denotes the open span between consecutive unique
     endpoints. Endpoint membership is not counted separately: touching
-    intervals occupy adjacent spans, and zero-width intervals contribute no
-    span to the returned distribution.
+    intervals occupy adjacent spans.  By default, zero-width intervals are
+    expanded by ``perturb_const`` representable floating-point values on each
+    side before the distribution is constructed.  This preserves a narrow,
+    deterministic LCD contribution for repeated adjacent time-series values.
     """
     xi,si = _int_dist_input(x1,x2)
+    xi = _perturb_constant_intervals(xi,perturb_const)
     # organize by edge order
     edges = xi.ravel()
     signs = si.ravel()
@@ -951,6 +992,7 @@ def interval_dist_peak(
         *,
         height         = False,
         quad_weighting = 'endpoint',
+        perturb_const  = 1,
         ):
     """Return a representative location at the peak of an interval distribution.
 
@@ -980,6 +1022,11 @@ def interval_dist_peak(
         each endpoint by the square root of its interval width, making the
         least-squares objective proportional to interval width.
 
+    perturb_const : int, optional
+        Number of floating-point steps used to expand any zero-width spans
+        before locating the peak.  One gives a deterministic ULP-scale
+        representation; zero leaves the spans unexpanded.
+
     Returns
     -------
     peak : float or (float, float)
@@ -987,6 +1034,7 @@ def interval_dist_peak(
         equal-height modes are averaged.
     """
     xi,_ = _int_dist_input(xi)
+    xi = _perturb_constant_intervals(xi,perturb_const)
     ci = _real_vector(ci,'interval counts')
     if len(ci)!=len(xi):
         msg = 'interval counts must have the same length as intervals'
@@ -1156,7 +1204,7 @@ def rolling_interval_dist_peak(
     xp = []
     cp = []
     for i1,i2 in windows:
-        xi,ci = interval_distribution(xia[i1:i2])
+        xi,ci = interval_distribution(xia[i1:i2],perturb_const=1)
         if len(ci)==0:
             msg = 'each rolling window must span a nonzero interval'
             raise ValueError(msg)
@@ -1167,6 +1215,7 @@ def rolling_interval_dist_peak(
             peak_frac      = peak_frac,
             height         = True,
             quad_weighting = quad_weighting,
+            perturb_const  = 1,
             )
         xp.append(xm)
         cp.append(cm)
@@ -1205,6 +1254,30 @@ def line_crossing_distribution(x,nperm=0):
     ci : ndarray
         Crossing counts, averaged over permutations when ``nperm`` is
         positive. Permutations are never connected to one another.
+
+    Notes
+    -----
+    Each adjacent pair defines an interval, and the distribution count at a
+    value is the number of such intervals that span it.  For a continuous
+    equilibrium series with independent samples ``X`` and ``Y`` drawn from
+    CDF ``F``, the corresponding crossing probability is
+
+    .. math::
+
+       L(z) = P(\min(X,Y) < z < \max(X,Y)) = 2F(z)[1-F(z)].
+
+    Thus, for a series with ``N`` samples, the expected count is
+    ``(N - 1) L(z)``.  The distribution is maximized at a median of the
+    sampled distribution, which motivates its use as a robust equilibrium
+    location estimator.  It is a crossing-rate curve rather than a normalized
+    probability density; when ``E[|X-Y|]`` is finite, its normalized form is
+    ``2 F(z) [1-F(z)] / E[|X-Y|]``.
+
+    In the ideal continuous i.i.d. case, a probability-scale LCD can be
+    inverted to obtain ``F(z) = (1 - sqrt(1 - 2 L(z))) / 2`` below a median
+    and ``F(z) = (1 + sqrt(1 - 2 L(z))) / 2`` above one, followed by
+    differentiation to obtain the density.  Empirical inversion is noisy,
+    and the LCD does not uniquely determine distributions with atoms or gaps.
     """
     x = _real_vector(x,'data array')
     if len(x)<2:
@@ -1219,7 +1292,7 @@ def line_crossing_distribution(x,nperm=0):
     # permutation-free (typical) case
     if nperm==0:
         xi,_ = time_series_intervals(x,t=None)
-        return interval_distribution(xi)
+        return interval_distribution(xi,perturb_const=1)
 
     # use permutation shuffling
     permutation_intervals = []
@@ -1228,7 +1301,9 @@ def line_crossing_distribution(x,nperm=0):
         np.random.shuffle(xp)
         xi,_ = time_series_intervals(xp,t=None)
         permutation_intervals.append(xi)
-    xi,ci = interval_distribution(np.vstack(permutation_intervals))
+    xi,ci = interval_distribution(
+        np.vstack(permutation_intervals),perturb_const=1
+        )
     return xi,ci/nperm
 #end def line_crossing_distribution
 
@@ -1264,6 +1339,7 @@ def lcd_peak(
         method         = method,
         peak_frac      = peak_frac,
         quad_weighting = quad_weighting,
+        perturb_const  = 1,
         )
     return xp
 #end def lcd_peak
diff --git a/nexus/nexus/tests/test_statistics.py b/nexus/nexus/tests/test_statistics.py
index 34861ebdd3..e84b612f6a 100644
--- a/nexus/nexus/tests/test_statistics.py
+++ b/nexus/nexus/tests/test_statistics.py
@@ -486,11 +486,28 @@ def test_interval_distribution_and_peak(monkeypatch):
     np.testing.assert_array_equal(touching_counts,[1,1,1])
 
     repeated_intervals,repeated_counts = statistics.interval_distribution(
-        np.array([[1.,1.],[1.,2.],[2.,2.]])
+        np.array([[1.,1.],[1.,2.],[2.,2.]]),perturb_const=0
         )
     np.testing.assert_array_equal(repeated_intervals,[[1.,2.]])
     np.testing.assert_array_equal(repeated_counts,[1])
 
+    constant_intervals = np.array([[1.,1.]])
+    perturbed_intervals,perturbed_counts = statistics.interval_distribution(
+        constant_intervals
+        )
+    np.testing.assert_array_equal(
+        perturbed_intervals,
+        [[np.nextafter(1.,-np.inf),np.nextafter(1.,np.inf)]],
+        )
+    np.testing.assert_array_equal(perturbed_counts,[1])
+    constant_peak,constant_height = statistics.interval_dist_peak(
+        constant_intervals,
+        [3.],
+        height=True,
+        )
+    assert(constant_peak==pytest.approx(1.))
+    assert(constant_height==3.)
+
     irregular_intervals,irregular_counts = statistics.interval_distribution(
         np.array([[0.,10.],[.5,.75]])
         )
@@ -580,6 +597,21 @@ def test_interval_distribution_and_peak(monkeypatch):
         statistics.interval_dist_peak(
             peak_intervals,peak_counts,quad_weighting='invalid'
             )
+    for perturb_const in (-1,1.5,True):
+        with pytest.raises(
+            ValueError,
+            match=r'constant perturbation must be a nonnegative integer',
+            ):
+            statistics.interval_distribution(
+                constant_intervals,perturb_const=perturb_const
+                )
+        with pytest.raises(
+            ValueError,
+            match=r'constant perturbation must be a nonnegative integer',
+            ):
+            statistics.interval_dist_peak(
+                constant_intervals,[1.],perturb_const=perturb_const
+                )
 #end def test_interval_distribution_and_peak
 
 
@@ -668,6 +700,17 @@ def fail_shuffle(values):
     np.testing.assert_array_equal(counts,[1,2])
     assert(statistics.lcd_peak(x)==pytest.approx(1.5))
 
+    constant = np.full(4,5.)
+    constant_intervals,constant_counts = statistics.line_crossing_distribution(
+        constant
+        )
+    np.testing.assert_array_equal(
+        constant_intervals,
+        [[np.nextafter(5.,-np.inf),np.nextafter(5.,np.inf)]],
+        )
+    np.testing.assert_array_equal(constant_counts,[3])
+    assert(statistics.lcd_peak(constant)==pytest.approx(5.))
+
     def reverse(values):
         values[:] = values[::-1]
     #end def reverse