diff --git a/c/set-protocol/include/set_signal.h b/c/set-protocol/include/set_signal.h index 0816e77..a27fcd2 100644 --- a/c/set-protocol/include/set_signal.h +++ b/c/set-protocol/include/set_signal.h @@ -10,6 +10,12 @@ enum { SET_SIGNAL_POLYNOMIAL, SET_SIGNAL_LINEAR, SET_SIGNAL_PCHIP, SET_SIGNAL_SP * Inputs may be unsorted; duplicate times are averaged. No extrapolation. * meta[0..2] = original count, unique count, RMSE at original measurements. * endpoint=1 includes the last time; 0 samples a cyclic period [first,last). + * On success, workspace blocks used by host formula annotations are: + * work[0..2*count): sorted original interleaved x,y pairs; + * work[2*count..3*count): unique x normalized to [0,1]; + * work[3*count..4*count): unique averaged y divided by max(1,max(abs(y))); + * work[4*count..5*count): PCHIP first / spline second derivatives in these + * normalized coordinates (only meta[1] entries in each unique-data block). * Return: 0 success, 1 bounds, 2 nonfinite/degenerate data, 3 rank deficient. */ PCAN_ABI_API int set_signal_reconstruct(const double *x, const double *y, size_t count, int method, unsigned degree, size_t output_count, int endpoint, diff --git a/python/set_devices/plot_formula.py b/python/set_devices/plot_formula.py new file mode 100644 index 0000000..ac4c8a2 --- /dev/null +++ b/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) diff --git a/python/set_devices/plot_processing.py b/python/set_devices/plot_processing.py index 1b626cc..afd5a78 100644 --- a/python/set_devices/plot_processing.py +++ b/python/set_devices/plot_processing.py @@ -76,6 +76,13 @@ class Curve: input_count: int unique_count: int rmse: float + pieces: tuple = () + + def piece_at(self, x): + for piece in self.pieces: + if piece.left <= x <= piece.right: + return piece + return None @property def label(self): @@ -100,8 +107,8 @@ def prepare(snapshot, key, method="pchip", output_count=1000, degree=2): def process(request): if request is None: raise ValueError("Нет доступного аналогового канала") - result = reconstruct(request.series.points, request.method, request.output_count, request.degree) - return Curve(request, tuple(result.points), result.input_count, result.unique_count, result.rmse) + result = reconstruct(request.series.points, request.method, request.output_count, request.degree, with_model=True) + return Curve(request, tuple(result.points), result.input_count, result.unique_count, result.rmse, result.pieces) def write_csv(curve, stream): diff --git a/python/set_devices/qt_ports/plot_annotations.py b/python/set_devices/qt_ports/plot_annotations.py new file mode 100644 index 0000000..44529b8 --- /dev/null +++ b/python/set_devices/qt_ports/plot_annotations.py @@ -0,0 +1,28 @@ +"""Small right-aligned mathematical annotations shared by Qt plots.""" +try: + from PySide6.QtCore import QRectF, Qt + from PySide6.QtGui import QColor +except ImportError: + from PySide2.QtCore import QRectF, Qt + from PySide2.QtGui import QColor + + +def paint_formulas(painter, rect, entries): + """Draw (color, text) entries in the right half, wrapping long formulas.""" + painter.save() + painter.setClipRect(rect) + metrics = painter.fontMetrics() + width = max(1., rect.width() * .48 - 12) + top = rect.top() + 5 + flags = int(Qt.AlignRight | Qt.AlignTop | Qt.TextWordWrap) + for color, text in entries: + box = QRectF(rect.right() - width - 6, top, width, max(1., rect.bottom() - top)) + height = painter.boundingRect(box, flags, text).height() + box.setHeight(height) + painter.fillRect(box.adjusted(-3, -1, 3, 1), QColor('#0b1119')) + painter.setPen(QColor(color)) + painter.drawText(box, flags, text) + top += height + metrics.height() * .4 + if top >= rect.bottom(): + break + painter.restore() diff --git a/python/set_devices/qt_ports/plot_processing.py b/python/set_devices/qt_ports/plot_processing.py index 25c8ac3..307a95d 100644 --- a/python/set_devices/qt_ports/plot_processing.py +++ b/python/set_devices/qt_ports/plot_processing.py @@ -332,6 +332,17 @@ class PlotProcessingAttachment(QObject): self._dirty = True self._refresh_timer.start(0) + def formula_entries(self, x, origin=0., unit=''): + entries = [] + offset = self._external_offset if self._external_curve is not None else 0. + for curve in self.curves: + piece = curve.piece_at(x + offset) + if piece is not None: + color_index = self.panel._result_ids[id(curve)] if self._external_curve is None else 0 + text = curve.label + '\n' + piece.text(origin + offset, unit) + entries.append((RESULT_COLORS[color_index % len(RESULT_COLORS)], text)) + return entries + def refresh(self): if self.panel is not None: self.panel.set_snapshot(self.snapshot()) diff --git a/python/set_devices/signal_reconstruction.py b/python/set_devices/signal_reconstruction.py index 776d6ee..3172526 100644 --- a/python/set_devices/signal_reconstruction.py +++ b/python/set_devices/signal_reconstruction.py @@ -14,6 +14,7 @@ class Reconstruction: input_count: int unique_count: int rmse: float + pieces: tuple = () @lru_cache(maxsize=1) def library(): @@ -29,7 +30,7 @@ def library(): except (AttributeError, OSError, RuntimeError) as error: raise ValueError("Пересоберите SETProtocol с set_signal.c и set_wavegen.c") from error -def reconstruct(points, method="pchip", output_count=1000, degree=2, *, endpoint=True): +def reconstruct(points, method="pchip", output_count=1000, degree=2, *, endpoint=True, with_model=False): if method not in METHODS or type(output_count) is not int or not 2 <= output_count <= 10000: raise ValueError("Неизвестный метод или число выходных точек вне 2…10000") count = len(points) @@ -44,7 +45,14 @@ def reconstruct(points, method="pchip", output_count=1000, degree=2, *, endpoint raise ValueError({1: "Недостаточно точек для выбранной степени или неверные параметры", 2: "Нужны конечные значения и минимум две различные временные точки", 3: "Неустойчивая аппроксимация: уменьшите степень"}.get(code, "Ошибка расчёта")) - return Reconstruction(list(zip(ox, oy)), int(meta[0]), int(meta[1]), meta[2]) + pieces = () + if with_model: + from .plot_formula import reconstruction_pieces + if method == 'polynomial' and output_count <= degree: + pieces = reconstruct(points, method, degree + 1, degree, with_model=True).pieces + else: + pieces = reconstruction_pieces(work, count, int(meta[1]), method, degree, oy, endpoint) + return Reconstruction(list(zip(ox, oy)), int(meta[0]), int(meta[1]), meta[2], pieces) def dac12(volts, vref=3.3): data = (C.c_double * len(volts))(*volts) diff --git a/python/tests/test_plot_formula.py b/python/tests/test_plot_formula.py new file mode 100644 index 0000000..e30915d --- /dev/null +++ b/python/tests/test_plot_formula.py @@ -0,0 +1,32 @@ +"""Formula coefficients must reproduce the native reconstruction, not a refit.""" +import unittest +from set_devices.signal_reconstruction import reconstruct +from set_devices.plot_formula import linear_piece + + +class PlotFormulaTests(unittest.TestCase): + def test_all_methods_match_native_values_on_irregular_and_duplicate_knots(self): + points = [(0, 2), (.3, 8), (.3, 4), (1.1, -3), (2, 7), (4, 1), (5, 2)] + for method in ('linear', 'pchip', 'spline', 'polynomial'): + result = reconstruct(points, method, 101, 3, with_model=True) + for x, y in result.points: + piece = next(piece for piece in result.pieces if piece.left <= x <= piece.right) + self.assertAlmostEqual(y, piece.value(x), places=9) + + def test_polynomial_model_with_fewer_output_samples_than_coefficients(self): + points = [(i, 1 + i * .2 - i ** 2 * .03 + i ** 5 * .0001) for i in range(9)] + result = reconstruct(points, 'polynomial', 2, 5, with_model=True) + for x, y in points: + self.assertAlmostEqual(y, result.pieces[0].value(x), places=9) + + def test_source_line_and_formula_use_displayed_origin(self): + piece = linear_piece([1000, 1100, 1200], [1, 3, 0], 1050) + self.assertAlmostEqual(2, piece.value(1050)) + text = piece.text(1000, 'мс') + self.assertIn('y = 1 +2·u', text) + self.assertIn('0 ≤ x ≤ 100 мс', text) + self.assertIsNone(linear_piece([0, 1, 2], [1, float('nan'), 3], .5)) + + +if __name__ == '__main__': + unittest.main()