interpretation speed up update

services/seq-interp/src/interfaces/rf_exporter.py

@@ -4,52 +4,19 @@ import numpy as np
 class RFExporter:
     """
     Экспорт RF в формате test1_full/srv_interp.py
-    """
-
-    @staticmethod
-    def _radio_ampl_convertation(rf_ampl, t_rf, rf_raster_local):
-        """
-        Resample the RF waveform to a dense raster grid.
-
-        pypulseq represents block (hard) pulses as just 2 samples with equal
-        complex amplitude at t_start and t_end.  The old loop-based approach
-        treated every large time-gap as a zero region, which caused block pulses
-        to disappear entirely.
-
-        np.interp handles all pulse shapes correctly:
-          - Shaped pulses  : many samples → linear interpolation (good approx)
-          - Block pulses   : 2 equal samples → constant amplitude throughout
-          - Inter-pulse gaps (amplitude ≈ 0) : np.interp fills with ~0
-        """
-        rf_ampl = np.asarray(rf_ampl, dtype=complex)
-        t_rf    = np.asarray(t_rf,    dtype=float)
-
-        if len(rf_ampl) == 0:
-            return np.empty(0)
-
-        rf_ampl_maximum = float(np.max(np.abs(rf_ampl)))
-        if rf_ampl_maximum == 0:
-            return np.empty(0)
-
-        proportional_cf_rf = 127.0 / rf_ampl_maximum
 
-        # Dense time grid aligned to raster
-        t_start   = float(t_rf[0])
-        t_end     = float(t_rf[-1])
-        n_samples = max(1, round((t_end - t_start) / rf_raster_local))
-        t_dense   = np.linspace(t_start, t_end, n_samples)
-
-        # Linear interpolation — exact for block pulses (constant), good approx
-        # for shaped pulses.  Extrapolation is clamped to boundary values.
-        real_dense = np.interp(t_dense, t_rf, rf_ampl.real)
-        imag_dense = np.interp(t_dense, t_rf, rf_ampl.imag)
+    Файл .bin — это непрерывный поток int8-пар [re, im] на плотной растровой
+    сетке от первого до последнего RF-сэмпла последовательности. Для длинных
+    последовательностей (минуты при растре ~1 мкс) сетка достигает сотен
+    миллионов отсчётов, поэтому экспорт потоковый: файл пишется кусками
+    ограниченного размера, а межимпульсные паузы, гарантированно дающие
+    нулевые байты, эмитятся из нулевого буфера без интерполяции.
+    """
 
-        # Interleaved [re, im, re, im, ...]; np.round is round-half-even,
-        # same as builtin round() used previously.
-        out = np.empty(2 * n_samples)
-        out[0::2] = np.round(real_dense * proportional_cf_rf)
-        out[1::2] = np.round(imag_dense * proportional_cf_rf)
-        return out
+    # Отсчётов плотной сетки на один интерполируемый кусок (~24 МБ временных
+    # массивов) и размер нулевого буфера в байтах.
+    _CHUNK = 1 << 20
+    _ZERO_BUF = bytes(1 << 22)
 
     def export(self, waveforms: dict, params: dict, output_dir: str):
         rf_raster_local = params["rf_raster_time"]
@@ -63,25 +30,134 @@ class RFExporter:
         else:
             empty_block_time_delay = 0
 
-        head = np.zeros(int(2 * (empty_block_time_delay // rf_raster_local)))
-        body = self._radio_ampl_convertation(
-            waveforms["rf"],
-            waveforms["t_rf"],
-            rf_raster_local,
-        )
-
+        n_head = int(2 * (empty_block_time_delay // rf_raster_local))
         scale_rf = params.get("scale_rf", 1.0)
-        rf_out = np.round(np.concatenate((head, body)) * scale_rf)
-
-        if len(rf_out) and (rf_out.min() < -128 or rf_out.max() > 127):
-            raise ValueError(
-                f"RF amplitude out of int8 range after scale_rf={scale_rf}: "
-                f"[{rf_out.min():.0f}, {rf_out.max():.0f}]"
-            )
 
         file_path = f"{output_dir}/rf_{rf_raster_local}_raster.bin"
         with open(file_path, "wb") as file_rf:
-            file_rf.write(rf_out.astype(np.int8).tobytes())
+            # Ведущие нули (задержка пустого блока); нули не зависят от scale_rf.
+            self._write_zeros(file_rf, n_head)
+            self._stream_body(
+                file_rf,
+                waveforms["rf"],
+                waveforms["t_rf"],
+                rf_raster_local,
+                scale_rf,
+            )
 
         np.savetxt(f"{output_dir}/rf_time.txt", np.transpose(waveforms["t_rf"]))
         np.savetxt(f"{output_dir}/rf_ampl.txt", np.transpose(waveforms["rf"]))
+
+    # ------------------------------------------------------------------
+    # Streaming implementation
+    # ------------------------------------------------------------------
+
+    @classmethod
+    def _write_zeros(cls, f, n_bytes: int):
+        buf = cls._ZERO_BUF
+        while n_bytes > 0:
+            n = min(n_bytes, len(buf))
+            f.write(buf[:n])
+            n_bytes -= n
+
+    @classmethod
+    def _stream_body(cls, f, rf_ampl, t_rf, rf_raster_local, scale_rf):
+        """
+        Пишет тело файла: поток int8-пар [re, im], эквивалентный
+        поэлементно np.interp по np.linspace(t_start, t_end, n_samples)
+        с двухступенчатым округлением (как в прежней невекторизованной
+        и векторизованной версиях).
+
+        pypulseq представляет блочные (hard) импульсы двумя сэмплами равной
+        амплитуды, поэтому линейная интерполяция корректна для всех форм:
+        shaped-импульсы — хорошее приближение, block-импульсы — константа,
+        межимпульсные паузы — линейный спад ~0.
+        """
+        rf_ampl = np.asarray(rf_ampl, dtype=complex)
+        t_rf    = np.asarray(t_rf,    dtype=float)
+
+        if len(rf_ampl) == 0:
+            return
+
+        rf_ampl_maximum = float(np.max(np.abs(rf_ampl)))
+        if rf_ampl_maximum == 0:
+            return
+
+        cf = 127.0 / rf_ampl_maximum
+
+        t_start   = float(t_rf[0])
+        t_end     = float(t_rf[-1])
+        n_samples = max(1, round((t_end - t_start) / rf_raster_local))
+        # Шаг сетки, битово повторяющей np.linspace(t_start, t_end,
+        # n_samples): t_k = k * step + t_start, t_[n-1] = t_end.
+        step = (t_end - t_start) / max(n_samples - 1, 1)
+
+        for k0, k1, is_zero in cls._segments(t_rf, rf_ampl, cf, step, n_samples):
+            if is_zero:
+                cls._write_zeros(f, 2 * (k1 - k0))
+            else:
+                for c0 in range(k0, k1, cls._CHUNK):
+                    c1 = min(c0 + cls._CHUNK, k1)
+                    cls._write_interp_chunk(
+                        f, c0, c1, t_rf, rf_ampl, cf, scale_rf,
+                        t_start, t_end, step, n_samples,
+                    )
+
+    @staticmethod
+    def _segments(t_rf, rf_ampl, cf, step, n_samples):
+        """
+        Разбивает плотную сетку [0, n_samples) на сегменты (k0, k1, is_zero).
+
+        Zero-сегмент — интервал строго внутри паузы между соседними
+        RF-сэмплами, оба из которых квантуются в 0 (|re|·cf и |im|·cf ≤ 0.5).
+        Линейная интерполяция между такими точками по модулю не превышает
+        максимума концов, поэтому все плотные отсчёты внутри дают 0 байт —
+        их можно эмитить без вычислений. Границы берутся с запасом в 2
+        отсчёта: краевые точки уходят в обычную интерполяцию, так что
+        погрешности плавающей точки на результат не влияют.
+        """
+        re_zero = np.abs(np.round(rf_ampl.real * cf)) == 0
+        im_zero = np.abs(np.round(rf_ampl.imag * cf)) == 0
+        pt_zero = re_zero & im_zero
+
+        gaps = np.diff(t_rf)
+        wide = gaps > 8 * step
+        zero_gap = wide & pt_zero[:-1] & pt_zero[1:]
+
+        t0 = t_rf[0]
+        segments = []
+        cursor = 0
+        for i in np.flatnonzero(zero_gap):
+            k_lo = int(np.ceil((t_rf[i] - t0) / step)) + 2
+            k_hi = int(np.floor((t_rf[i + 1] - t0) / step)) - 2
+            k_lo = max(k_lo, cursor)
+            k_hi = min(k_hi, n_samples)
+            if k_hi - k_lo < 16:
+                continue
+            if k_lo > cursor:
+                segments.append((cursor, k_lo, False))
+            segments.append((k_lo, k_hi, True))
+            cursor = k_hi
+        if cursor < n_samples:
+            segments.append((cursor, n_samples, False))
+        return segments
+
+    @staticmethod
+    def _write_interp_chunk(f, k0, k1, t_rf, rf_ampl, cf, scale_rf,
+                            t_start, t_end, step, n_samples):
+        k = np.arange(k0, k1, dtype=np.float64)
+        t = k * step + t_start
+        if k1 == n_samples and n_samples > 1:
+            t[-1] = t_end
+
+        out = np.empty(2 * (k1 - k0))
+        out[0::2] = np.round(np.interp(t, t_rf, rf_ampl.real) * cf)
+        out[1::2] = np.round(np.interp(t, t_rf, rf_ampl.imag) * cf)
+        out = np.round(out * scale_rf)
+
+        if out.min() < -128 or out.max() > 127:
+            raise ValueError(
+                f"RF amplitude out of int8 range after scale_rf={scale_rf}: "
+                f"[{out.min():.0f}, {out.max():.0f}]"
+            )
+        f.write(out.astype(np.int8).tobytes())

services/seq-interp/src/interfaces/xml_generator.py

@@ -1,11 +1,15 @@
 from itertools import accumulate
 
-from yattag import Doc, indent
-
 
 class XMLGenerator:
     """
     Генератор XML синхронизации в формате, совместимом с test1_full/srv_interp.py
+
+    XML собирается прямой конкатенацией строк (формат прежнего yattag +
+    indent(): вложенность по 4 пробела, переводы строк '\\r', листовые теги в
+    одну строку). yattag создавал контекст-менеджер на каждый тег, а indent()
+    повторно парсил весь документ — на десятках тысяч блоков это занимало
+    заметную долю всего пайплайна.
     """
 
     def generate(self, sync_data: dict, path: str, hw):
@@ -17,81 +21,81 @@ class XMLGenerator:
         synchro_block_timer = sync_data["synchro_block_timer"]
         min_block_time = sync_data["min_block_time"]
 
-        doc, tag, text = Doc().tagtext()
-
-        with tag("root"):
-            with tag("ParamCount"):
-                text(number_of_blocks + 1)
-
-            adc_times_values = []
-            adc_times_starts = []
-
-            # Префиксные суммы длительностей: duration_prefix[k] ==
-            # sum(blocks_duration[0:k]). Считаются один раз, иначе суммирование
-            # на каждом ADC-фронте давало O(N^2) на длинных последовательностях.
-            duration_prefix = list(accumulate(blocks_duration, initial=0))
-
-            # Массивы синхронизатора имеют ведущий seed-элемент на индексе 0
-            # (стартовый блок), значимые блоки идут по индексам 1..nb. Seed
-            # представлен заголовочным тегом (*1), поэтому в цикле эмитим
-            # элементы 1..nb: block_iter k -> индекс массива k + 1.
-            # (Совпадает с эталонной логикой srv_interp.synchronization.)
-
-            with tag("RF"):
-                with tag("RF1"):
-                    text(0)
-                for rf_iter in range(number_of_blocks):
-                    with tag("RF" + str(rf_iter + 2)):
-                        text(gate_rf[rf_iter + 1])
-
-            with tag("SW"):
-                with tag("SW1"):
-                    text(1)
-                for sw_iter in range(number_of_blocks):
-                    with tag("SW" + str(sw_iter + 2)):
-                        text(gate_tr_switch[sw_iter + 1])
-
-            with tag("ADC"):
-                with tag("ADC1"):
-                    text(0)
-                for adc_iter in range(number_of_blocks):
-                    idx = adc_iter + 1
-                    # Непрерывный участок HIGH (в т.ч. собранный из разбитых
-                    # блоков) считается одним событием ADC: засекаем фронт и
-                    # суммируем длительность всего участка для points-окна.
-                    is_rising_edge = (
-                        gate_adc[idx] == 1 and gate_adc[idx - 1] == 0
-                    )
-                    if is_rising_edge:
-                        adc_times_starts.append(duration_prefix[idx])
-                        run_duration = 0
-                        run_iter = idx
-                        while (run_iter < len(gate_adc)
-                               and gate_adc[run_iter] == 1):
-                            run_duration += blocks_duration[run_iter]
-                            run_iter += 1
-                        adc_times_values.append(run_duration)
-                    with tag("ADC" + str(adc_iter + 2)):
-                        text(gate_adc[idx])
-
-            with tag("GR"):
-                with tag("GR1"):
-                    text(1)
-                for gr_iter in range(number_of_blocks):
-                    # Последнее событие импульсной последовательности
-                    # переводим из LOW в HIGH (включение градиентной системы)
-                    gr_value = 1 if gr_iter == number_of_blocks - 1 else 0
-                    with tag("GR" + str(gr_iter + 2)):
-                        text(gr_value)
-
-            with tag("CL"):
-                with tag("CL1"):
-                    text(int(min_block_time / synchro_block_timer))
-                for cl_iter in range(number_of_blocks):
-                    with tag("CL" + str(cl_iter + 2)):
-                        text(int(blocks_duration[cl_iter + 1] / synchro_block_timer))
-
-        xml_string = indent(doc.getvalue(), indentation=" " * 4, newline="\r")
+        adc_times_values = []
+        adc_times_starts = []
+
+        # Префиксные суммы длительностей: duration_prefix[k] ==
+        # sum(blocks_duration[0:k]). Считаются один раз, иначе суммирование
+        # на каждом ADC-фронте давало O(N^2) на длинных последовательностях.
+        duration_prefix = list(accumulate(blocks_duration, initial=0))
+
+        # Массивы синхронизатора имеют ведущий seed-элемент на индексе 0
+        # (стартовый блок), значимые блоки идут по индексам 1..nb. Seed
+        # представлен заголовочным тегом (*1), поэтому эмитим элементы
+        # 1..nb: block_iter k -> индекс массива k + 1.
+        # (Совпадает с эталонной логикой srv_interp.synchronization.)
+
+        lines = ["<root>"]
+        lines.append(f"    <ParamCount>{number_of_blocks + 1}</ParamCount>")
+
+        lines.append("    <RF>")
+        lines.append("        <RF1>0</RF1>")
+        for rf_iter in range(number_of_blocks):
+            n = rf_iter + 2
+            lines.append(f"        <RF{n}>{gate_rf[rf_iter + 1]}</RF{n}>")
+        lines.append("    </RF>")
+
+        lines.append("    <SW>")
+        lines.append("        <SW1>1</SW1>")
+        for sw_iter in range(number_of_blocks):
+            n = sw_iter + 2
+            lines.append(f"        <SW{n}>{gate_tr_switch[sw_iter + 1]}</SW{n}>")
+        lines.append("    </SW>")
+
+        lines.append("    <ADC>")
+        lines.append("        <ADC1>0</ADC1>")
+        for adc_iter in range(number_of_blocks):
+            idx = adc_iter + 1
+            # Непрерывный участок HIGH (в т.ч. собранный из разбитых
+            # блоков) считается одним событием ADC: засекаем фронт и
+            # суммируем длительность всего участка для points-окна.
+            is_rising_edge = (
+                gate_adc[idx] == 1 and gate_adc[idx - 1] == 0
+            )
+            if is_rising_edge:
+                adc_times_starts.append(duration_prefix[idx])
+                run_duration = 0
+                run_iter = idx
+                while (run_iter < len(gate_adc)
+                       and gate_adc[run_iter] == 1):
+                    run_duration += blocks_duration[run_iter]
+                    run_iter += 1
+                adc_times_values.append(run_duration)
+            n = adc_iter + 2
+            lines.append(f"        <ADC{n}>{gate_adc[idx]}</ADC{n}>")
+        lines.append("    </ADC>")
+
+        lines.append("    <GR>")
+        lines.append("        <GR1>1</GR1>")
+        for gr_iter in range(number_of_blocks):
+            # Последнее событие импульсной последовательности
+            # переводим из LOW в HIGH (включение градиентной системы)
+            gr_value = 1 if gr_iter == number_of_blocks - 1 else 0
+            n = gr_iter + 2
+            lines.append(f"        <GR{n}>{gr_value}</GR{n}>")
+        lines.append("    </GR>")
+
+        lines.append("    <CL>")
+        lines.append(f"        <CL1>{int(min_block_time / synchro_block_timer)}</CL1>")
+        for cl_iter in range(number_of_blocks):
+            n = cl_iter + 2
+            cl_value = int(blocks_duration[cl_iter + 1] / synchro_block_timer)
+            lines.append(f"        <CL{n}>{cl_value}</CL{n}>")
+        lines.append("    </CL>")
+
+        lines.append("</root>")
+
+        xml_string = "\r".join(lines)
         with open(path, "w", encoding="utf-8") as f:
             f.write(xml_string)