Serialized magnets: mixed-function families, strengths/hardwares arrays and partial series writes

pyaml/arrays/magnet_array.py

@@ -67,18 +67,16 @@ def set(self, value: np.array):
         nvalue = np.ones(self.__nb) * value if isinstance(value, float) else value
         if not self.__aggregator:
             # Magnets of a series share one power supply: it is written once, from the strengths of the
-            # magnets of the series present in the array
+            # magnets of the series present in the array, the others following the resulting hardware value
             series: dict[int, tuple[Magnet, dict[int, float]]] = {}
             for idx, m in enumerate(self.__magnets):
                 if m.model is not None and m.model.shares_power_supply():
                     series.setdefault(id(m.model), (m, {}))[1][m.strength.index()] = nvalue[idx]
                 else:
                     m.strength.set(nvalue[idx])
             for m, requested in series.values():
-                if len(requested) == 1:
-                    m.strength.set(next(iter(requested.values())))
-                else:
-                    m.hardware.set(m.model.compute_hardware_values_partial(None, requested)[0])
+                current_hardware = np.array([m.hardware.get()])
+                m.hardware.set(m.model.compute_hardware_values_partial(current_hardware, requested)[0])
         else:
             self.__aggregator.set(nvalue)
 

pyaml/control/abstract_impl.py

@@ -431,7 +431,8 @@ def add_magnet(self, magnet: Magnet, devs: list[DeviceAccess]):
             Hardware devices associated with the magnet model.
         """
         # Magnets of a serialized group share the serialized model (one power supply for N strengths):
-        # the power-supply current is computed by the model from the strengths of the magnets present in the array.
+        # the power-supply current is computed by the model from the strengths of the magnets present in the array,
+        # the others following the resulting current.
         is_indexed = isinstance(magnet.strength, (abstract.RWMapper, RWSerializedStrengthScalar))
         strengthIndex = magnet.strength.index() if is_indexed else 0
         if magnet.model not in self.__models:

pyaml/lattice/simulator.py

@@ -271,12 +271,12 @@ def _fill_serialized_magnets(self, magnets: SerializedMagnets) -> None:
         strengths = []
         for index, magnet in enumerate(magnets.get_magnets()):
             current = (
-                RWHardwareScalar(self.get_at_elems(magnet), magnets.polynom, magnets.model.get_sub_model(index))
+                RWHardwareScalar(self.get_at_elems(magnet), magnets.polynoms[index], magnets.model.get_sub_model(index))
                 if magnets.model.has_hardware()
                 else None
             )
             strength = (
-                RWStrengthScalar(self.get_at_elems(magnet), magnets.polynom, magnets.model.get_sub_model(index))
+                RWStrengthScalar(self.get_at_elems(magnet), magnets.polynoms[index], magnets.model.get_sub_model(index))
                 if magnets.model.has_physics()
                 else None
             )

pyaml/magnet/linear_serialized_model.py

@@ -21,7 +21,7 @@
 
 logger = logging.getLogger(__name__)
 
-# Relative spread above which the currents requested by members of a series are considered inconsistent
+# Relative gap above which strengths requested for a series are considered unreachable with a single hardware value
 _SPREAD_RTOL = 1e-6
 
 
@@ -180,6 +180,8 @@ class LinearSerializedMagnetModel(MagnetModel, DynamicValidation):
         Return the average strength of the series for a hardware value.
     compute_hardware_from_mean_strength(strength)
         Return the hardware value giving an average strength of the series.
+    compute_hardware_from_strengths(strengths)
+        Return the hardware value for a requested strength of each magnet of the series.
     shares_power_supply()
         Return True, all the magnets of the series share one power supply.
     get_strength_units()
@@ -355,37 +357,37 @@ def compute_strengths(self, currents: np.array) -> np.array:
         return np.array([model.compute_strengths([current])[0] for model in self.__sub_models])
 
     def compute_hardware_values_partial(
-        self, current_hardware: np.ndarray | None, strengths_by_index: dict[int, float]
+        self, current_hardware: np.ndarray, strengths_by_index: dict[int, float]
     ) -> np.ndarray:
         """
         Convert the strengths of some magnets of the series to the shared hardware value.
 
-        Each requested strength implies a hardware value through its own submodel. The power supply gets the mean of
-        the hardware values implied by the requested magnets only: the magnets that are not requested do not take part.
-        A warning is logged when the requested strengths are not consistent with a single hardware value.
+        Only the requested magnets define the hardware value: it solves ``sum(f_i(I)) == sum(k_i)`` over them. The
+        other magnets of the series follow the resulting hardware value. A warning is logged when the requested
+        strengths cannot all be reached with one hardware value, or when they are outside the excitation curves.
 
         Parameters
         ----------
-        current_hardware : np.array | None
-            Current hardware value of the series. Not used: the requested strengths define the hardware value.
+        current_hardware : np.array
+            Current hardware value of the series, one value per power supply. Only its shape is used.
         strengths_by_index : dict[int, float]
             Requested strengths, indexed by the position of the magnet in the series.
 
         Returns
         -------
         np.array
-            Hardware value of the power supply.
+            Hardware value of the power supply, with the shape of ``current_hardware``.
+
+        Raises
+        ------
+        PyAMLException
+            If ``current_hardware`` does not hold one value per power supply of the series, or if the average strength
+            of the requested magnets is not monotonic in the hardware value.
         """
-        currents = np.array(
-            [self.__sub_models[idx].compute_hardware_values([strength])[0] for idx, strength in strengths_by_index.items()]
-        )
-        current = np.mean(currents)
-        if np.ptp(currents) > _SPREAD_RTOL * max(abs(current), 1.0):
-            logger.warning(
-                f"Inconsistent strengths requested for magnets {list(strengths_by_index.keys())} of a series on "
-                f"{self._powerconverter}: implied hardware values {currents}, the mean {current} is applied"
-            )
-        return np.array([current])
+        _check_len(current_hardware, "current_hardware", len(self.get_device_names()))
+        indices = list(strengths_by_index.keys())
+        strengths = np.array(list(strengths_by_index.values()), dtype=float)
+        return np.array([self._compute_hardware_from_strengths(strengths, indices)])
 
     def compute_mean_strength(self, current: float) -> float:
         """
@@ -408,7 +410,8 @@ def compute_hardware_from_mean_strength(self, strength: float) -> float:
         Return the hardware value giving an average strength of the series.
 
         The submodels are piecewise linear in the hardware value, so the average strength is piecewise linear on the
-        union of the curve breakpoints and is inverted exactly by interpolation.
+        union of the curve breakpoints and is inverted exactly by interpolation. A strength outside the range reachable
+        by the excitation curves is clamped to the end of the curves, and a warning is logged.
 
         Parameters
         ----------
@@ -425,16 +428,104 @@ def compute_hardware_from_mean_strength(self, strength: float) -> float:
         PyAMLException
             If the average strength is not monotonic in the hardware value.
         """
-        currents = np.unique(np.concatenate([curve.get_curve()[:, 0] for curve in self.__curves]))
-        strengths = np.array([self.compute_mean_strength(current) for current in currents])
+        return self._compute_hardware_from_mean_strength(strength, list(range(self.__nbMagnets)))[0]
+
+    def compute_hardware_from_strengths(self, strengths: np.ndarray) -> float:
+        """
+        Return the hardware value for a requested strength of each magnet of the series.
+
+        The series has a single hardware value, so it solves ``sum(f_i(I)) == sum(k_i)`` for ``I``. When the requested
+        strengths cannot all be reached with one hardware value, the result is a compromise and a warning is logged.
+
+        Parameters
+        ----------
+        strengths : np.ndarray
+            Requested strength ``k_i`` of each magnet of the series, in declaration order.
+
+        Returns
+        -------
+        float
+            Hardware value of the power supply.
+
+        Raises
+        ------
+        PyAMLException
+            If the number of strengths does not match the number of magnets, or if the average strength is not
+            monotonic in the hardware value.
+        """
+        strengths = np.asarray(strengths, dtype=float)
+        _check_len(strengths, "strengths", self.__nbMagnets)
+        return self._compute_hardware_from_strengths(strengths, list(range(self.__nbMagnets)))
+
+    def _compute_hardware_from_strengths(self, strengths: np.ndarray, indices: list[int]) -> float:
+        """
+        Return the hardware value solving ``sum(f_i(I)) == sum(k_i)`` over some magnets of the series.
+
+        Parameters
+        ----------
+        strengths : np.ndarray
+            Requested strength ``k_i`` of each magnet in ``indices``.
+        indices : list[int]
+            Positions of the requested magnets in the series.
+
+        Returns
+        -------
+        float
+            Hardware value of the power supply.
+        """
+        current, clamped = self._compute_hardware_from_mean_strength(float(np.mean(strengths)), indices)
+        if not clamped:
+            reached = np.array([self.__sub_models[idx].compute_strengths([current])[0] for idx in indices])
+            gaps = reached - strengths
+            if np.any(np.abs(gaps) > _SPREAD_RTOL * np.maximum(np.abs(strengths), 1.0)):
+                worst = int(np.argmax(np.abs(gaps)))
+                logger.warning(
+                    f"Inconsistent strengths requested for the series on {self._powerconverter}: the hardware value "
+                    f"{current} misses the requested strengths by up to {gaps[worst]:.3g} (magnet {indices[worst]})"
+                )
+        return current
+
+    def _compute_hardware_from_mean_strength(self, strength: float, indices: list[int]) -> tuple[float, bool]:
+        """
+        Return the hardware value giving an average strength of some magnets of the series.
+
+        Parameters
+        ----------
+        strength : float
+            Requested average strength of the magnets in ``indices``.
+        indices : list[int]
+            Positions of the magnets taken into account in the series.
+
+        Returns
+        -------
+        tuple[float, bool]
+            Hardware value, and whether it was clamped to the end of the excitation curves.
+
+        Raises
+        ------
+        PyAMLException
+            If the average strength of these magnets is not monotonic in the hardware value.
+        """
+        currents = np.unique(np.concatenate([self.__curves[idx].get_curve()[:, 0] for idx in indices]))
+        strengths = np.array(
+            [np.mean([self.__sub_models[idx].compute_strengths([current])[0] for idx in indices]) for current in currents]
+        )
         steps = np.diff(strengths)
         if np.all(steps <= 0):
             currents, strengths = currents[::-1], strengths[::-1]
         elif not np.all(steps >= 0):
             raise PyAMLException(
                 f"Average strength of the series on {self._powerconverter} is not monotonic, it cannot be inverted"
             )
-        return float(np.interp(strength, strengths, currents))
+        tolerance = _SPREAD_RTOL * max(abs(strengths[0]), abs(strengths[-1]), 1.0)
+        clamped = bool(strength < strengths[0] - tolerance or strength > strengths[-1] + tolerance)
+        if clamped:
+            logger.warning(
+                f"Strength {strength:.6g} requested for the series on {self._powerconverter} is outside the range "
+                f"[{strengths[0]:.6g}, {strengths[-1]:.6g}] reachable by the excitation curves, the hardware value is "
+                f"clamped to the end of the curves"
+            )
+        return float(np.interp(strength, strengths, currents)), clamped
 
     def shares_power_supply(self) -> bool:
         """Return True, all the magnets of the series share one power supply."""