Expose numeric trend segment formulas from native reconstruction models

python/set_devices/plot_formula.py

@@ -0,0 +1,107 @@
+"""Local polynomial representation of the curves actually drawn on a plot."""
+from dataclasses import dataclass
+import math
+
+
+@dataclass(frozen=True)
+class FormulaPiece:
+    left: float
+    right: float
+    coefficients: tuple  # ascending powers of u=(x-left)/(right-left)
+
+    def value(self, x):
+        u = (x - self.left) / (self.right - self.left)
+        value = 0.
+        for coefficient in reversed(self.coefficients):
+            value = value * u + coefficient
+        return value
+
+    def text(self, origin=0., unit=''):
+        def number(value):
+            return format(value, '.6g')
+        terms = []
+        for power, coefficient in enumerate(self.coefficients):
+            if coefficient == 0 and power:
+                continue
+            term = number(abs(coefficient))
+            if power:
+                term += '·u' + ('²' if power == 2 else '³' if power == 3 else f'^{power}' if power > 1 else '')
+            terms.append(('−' if coefficient < 0 else '+' if terms else '') + term)
+        left = self.left - origin
+        shift = ('−' + number(left)) if left >= 0 else ('+' + number(-left))
+        return ('y = ' + ' '.join(terms) + '\n'
+                + f'u = (x{shift})/{number(self.right - self.left)}; '
+                + f'{number(left)} ≤ x ≤ {number(self.right - origin)} {unit}').rstrip()
+
+
+def linear_piece(timestamps, values, x):
+    """A source polyline has a linear equation, not an inferred analytic model."""
+    previous = None
+    for stamp, value in zip(timestamps, values):
+        if not math.isfinite(stamp) or not math.isfinite(value):
+            previous = None
+            continue
+        if previous is not None:
+            t0, y0 = previous
+            if t0 <= x <= stamp and stamp > t0:
+                return FormulaPiece(t0, stamp, (y0, value - y0))
+        previous = stamp, value
+    return None
+
+
+def polynomial_coefficients(nodes, values):
+    """Convert evaluations of an existing polynomial to power coefficients."""
+    divided = list(values)
+    for order in range(1, len(nodes)):
+        for i in range(len(nodes) - 1, order - 1, -1):
+            divided[i] = (divided[i] - divided[i - 1]) / (nodes[i] - nodes[i - order])
+    result = [divided[-1]]
+    for i in range(len(nodes) - 2, -1, -1):
+        product = [0.] * (len(result) + 1)
+        for power, coefficient in enumerate(result):
+            product[power] -= nodes[i] * coefficient
+            product[power + 1] += coefficient
+        product[0] += divided[i]
+        result = product
+    return tuple(result)
+
+
+def reconstruction_pieces(work, count, unique_count, method, degree, output_y, endpoint):
+    """Read the native reconstruction workspace; no second interpolation fit.
+
+    set_signal.c stores sorted original pairs, normalized x/y and derivatives
+    in the first five count-sized blocks. Derivatives are with respect to the
+    normalized domain, and spline derivatives are second derivatives.
+    """
+    origin, end = work[0], work[2 * (count - 1)]
+    span = end - origin
+    scale = 1.
+    i = 0
+    while i < count:
+        j = i + 1
+        while j < count and work[2 * j] == work[2 * i]:
+            j += 1
+        scale = max(scale, abs(sum(work[2 * k + 1] / (j - i) for k in range(i, j))))
+        i = j
+    if method == 'polynomial':
+        indices = [round(i * (len(output_y) - 1) / degree) for i in range(degree + 1)]
+        nodes = [i / (len(output_y) - (1 if endpoint else 0)) for i in indices]
+        coefficients = polynomial_coefficients(nodes, [output_y[i] for i in indices])
+        return (FormulaPiece(origin, end, coefficients),)
+    pieces = []
+    for i in range(unique_count - 1):
+        x0, x1 = work[2 * count + i], work[2 * count + i + 1]
+        y0, y1 = work[3 * count + i], work[3 * count + i + 1]
+        d0, d1 = work[4 * count + i], work[4 * count + i + 1]
+        h = x1 - x0
+        if method == 'linear':
+            coefficients = (y0, y1 - y0)
+        elif method == 'pchip':
+            coefficients = (y0, h * d0, 3 * (y1 - y0) - h * (2 * d0 + d1),
+                            2 * (y0 - y1) + h * (d0 + d1))
+        else:
+            coefficients = (y0, y1 - y0 - h * h * (2 * d0 + d1) / 6,
+                            h * h * d0 / 2, h * h * (d1 - d0) / 6)
+        pieces.append(FormulaPiece(origin + span * x0, origin + span * x1,
+                                   tuple(c * scale for c in coefficients)))
+    return tuple(pieces)