diff --git a/pyaml/arrays/magnet_array.py b/pyaml/arrays/magnet_array.py index 7f05f35a..ec2669ea 100644 --- a/pyaml/arrays/magnet_array.py +++ b/pyaml/arrays/magnet_array.py @@ -67,7 +67,7 @@ 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(): @@ -75,10 +75,8 @@ def set(self, value: np.array): 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) diff --git a/pyaml/control/abstract_impl.py b/pyaml/control/abstract_impl.py index ffea3234..176697ab 100644 --- a/pyaml/control/abstract_impl.py +++ b/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: diff --git a/pyaml/lattice/simulator.py b/pyaml/lattice/simulator.py index e0922ae3..4e4f92e0 100644 --- a/pyaml/lattice/simulator.py +++ b/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 ) diff --git a/pyaml/magnet/linear_serialized_model.py b/pyaml/magnet/linear_serialized_model.py index 6915c862..c00858b9 100644 --- a/pyaml/magnet/linear_serialized_model.py +++ b/pyaml/magnet/linear_serialized_model.py @@ -21,8 +21,60 @@ 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 +# Normalized residual below which a breakpoint is considered to be exactly on the solution (rounding errors) +_ZERO_ATOL = 1e-12 + + +def _first_zero(currents: np.ndarray, residuals: np.ndarray) -> float: + """ + Return the lowest hardware value where an increasing piecewise-linear residual reaches zero. + + Parameters + ---------- + currents : np.ndarray + Breakpoints, in ascending order. + residuals : np.ndarray + Residual at each breakpoint, increasing. + + Returns + ------- + float + Hardware value, interpolated on the segment where the residual becomes non-negative. + """ + above = np.nonzero(residuals >= 0.0)[0] + if len(above) == 0: + return float(currents[-1]) + j = above[0] + if j == 0: + return float(currents[0]) + return float(currents[j - 1] - residuals[j - 1] * (currents[j] - currents[j - 1]) / (residuals[j] - residuals[j - 1])) + + +def _last_zero(currents: np.ndarray, residuals: np.ndarray) -> float: + """ + Return the highest hardware value where an increasing piecewise-linear residual is still zero. + + Parameters + ---------- + currents : np.ndarray + Breakpoints, in ascending order. + residuals : np.ndarray + Residual at each breakpoint, increasing. + + Returns + ------- + float + Hardware value, interpolated on the segment where the residual becomes positive. + """ + below = np.nonzero(residuals <= 0.0)[0] + if len(below) == 0: + return float(currents[0]) + j = below[-1] + if j == len(currents) - 1: + return float(currents[-1]) + return float(currents[j] - residuals[j] * (currents[j + 1] - currents[j]) / (residuals[j + 1] - residuals[j])) def _get_length(elem) -> int: @@ -142,7 +194,8 @@ class LinearSerializedMagnetModel(MagnetModel, DynamicValidation): The serialized magnets have a single degree of freedom, the power-supply current ``I``. Magnet ``i`` has the strength ``k_i = f_i(I)``, where ``f_i`` is its sub-model. The family strength is the average strength - ``sum(f_i(I)) / N``. + ``sum(p_i * f_i(I)) / N``, each strength taken with the polarity of the first + magnet (``p_i = -1`` for a magnet wired in reverse, ``+1`` otherwise). Parameters ---------- @@ -180,6 +233,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() @@ -332,10 +387,9 @@ def compute_hardware_values(self, strengths: np.array) -> np.array: Returns ------- np.array - Hardware value obtained from the submodels: the mean of the hardware values implied by each strength. + Hardware value of the power supply, computed by :meth:`compute_hardware_from_strengths`. """ - currents = [model.compute_hardware_values([s])[0] for s, model in zip(strengths, self.__sub_models, strict=True)] - return np.array([np.mean(currents)]) + return np.array([self.compute_hardware_from_strengths(strengths)]) def compute_strengths(self, currents: np.array) -> np.array: """ @@ -355,42 +409,45 @@ 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, solved as in :meth:`compute_hardware_from_strengths`. 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, if an excitation curve has + no strength range, or if an excitation curve is not monotonic. """ - 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: """ Return the average strength of the series for a hardware value. + Each strength is taken with the polarity of the first magnet, so that a magnet wired in reverse counts with + its magnitude instead of cancelling the others. + Parameters ---------- current : float @@ -399,16 +456,20 @@ def compute_mean_strength(self, current: float) -> float: Returns ------- float - Average strength ``sum(f_i(I)) / N``. + Average strength ``sum(p_i * f_i(I)) / N``, with ``p_i = +1`` when magnet ``i`` has the polarity of the first + magnet and ``-1`` otherwise. """ - return float(np.mean(self.compute_strengths([current]))) + return float(np.mean(self._polarities() * self.compute_strengths([current]))) 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. + The average strength is defined as in :meth:`compute_mean_strength`, so every term varies in the same direction + and the average is monotonic. The submodels are piecewise linear in the hardware value, so the average is + piecewise linear on the 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 ---------- @@ -418,15 +479,16 @@ def compute_hardware_from_mean_strength(self, strength: float) -> float: Returns ------- float - Hardware value ``I`` such that ``sum(f_i(I)) / N == strength``. + Hardware value ``I`` such that ``sum(p_i * f_i(I)) / N == strength``. Raises ------ PyAMLException If the average strength is not monotonic in the hardware value. """ + polarities = self._polarities() 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]) + strengths = np.array([np.mean(polarities * self.compute_strengths([current])) for current in currents]) steps = np.diff(strengths) if np.all(steps <= 0): currents, strengths = currents[::-1], strengths[::-1] @@ -434,8 +496,144 @@ def compute_hardware_from_mean_strength(self, strength: float) -> float: raise PyAMLException( f"Average strength of the series on {self._powerconverter} is not monotonic, it cannot be inverted" ) + tolerance = _SPREAD_RTOL * max(abs(strengths[0]), abs(strengths[-1]), 1.0) + if strength < strengths[0] - tolerance or strength > strengths[-1] + tolerance: + 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)) + 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 ``I``. It solves ``sum((f_i(I) - k_i) / span_i) == 0``, where ``span_i`` + is the signed strength range ``f_i(I_max) - f_i(I_min)`` of the excitation curve of magnet ``i``. Each term is + dimensionless and increasing, so every magnet weighs the same whatever its unit or its polarity, and the + equation always has a solution, exact for a request reachable with one hardware value. 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, if an excitation curve has no strength + range, or if an excitation curve is not monotonic. + """ + 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 for the requested strengths of 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._solve_normalized_residual(strengths, indices) + if not clamped: + reached = np.array([self.__sub_models[idx].compute_strengths([current])[0] for idx in indices]) + spans = np.array([self._strength_span(idx) for idx in indices]) + relative_gaps = (reached - strengths) / np.abs(spans) + if np.any(np.abs(relative_gaps) > _SPREAD_RTOL): + worst = int(np.argmax(np.abs(relative_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 {reached[worst] - strengths[worst]:.3g} " + f"(magnet {indices[worst]}, {100 * relative_gaps[worst]:.3g} % of its curve range)" + ) + return current + + def _solve_normalized_residual(self, strengths: np.ndarray, indices: list[int]) -> tuple[float, bool]: + """ + Solve ``g(I) = sum((f_i(I) - k_i) / span_i) == 0`` over some magnets of the series. + + ``g`` is piecewise linear on the union of the curve breakpoints and increasing, so its zero is found exactly by + interpolation. When ``g`` is zero on an interval (for example two saturated magnets), the middle of the interval + is returned. When the zero is outside the curves, the hardware value is clamped to the end of the curves and a + warning is logged. + + 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 + ------- + tuple[float, bool] + Hardware value, and whether it was clamped to the end of the excitation curves. + + Raises + ------ + PyAMLException + If an excitation curve has no strength range, or if an excitation curve is not monotonic. + """ + spans = np.array([self._strength_span(idx) for idx in indices]) + if np.any(spans == 0.0): + flat = [idx for idx, span in zip(indices, spans, strict=True) if span == 0.0] + raise PyAMLException( + f"Magnets {flat} of the series on {self._powerconverter} have an excitation curve with no strength " + f"range, the hardware value cannot be computed from their strength" + ) + currents = np.unique(np.concatenate([self.__curves[idx].get_curve()[:, 0] for idx in indices])) + residuals = np.array( + [ + sum( + (self.__sub_models[idx].compute_strengths([current])[0] - strength) / span + for idx, strength, span in zip(indices, strengths, spans, strict=True) + ) + for current in currents + ] + ) + if np.any(np.diff(residuals) < -_SPREAD_RTOL): + raise PyAMLException( + f"An excitation curve of the series on {self._powerconverter} is not monotonic, " + f"the hardware value cannot be computed" + ) + residuals[np.abs(residuals) < _ZERO_ATOL] = 0.0 + if residuals[0] > _SPREAD_RTOL or residuals[-1] < -_SPREAD_RTOL: + current = currents[0] if residuals[0] > 0.0 else currents[-1] + logger.warning( + f"Strengths {strengths} requested for the series on {self._powerconverter} are outside the range " + f"reachable by the excitation curves, the hardware value is clamped to the end of the curves ({current})" + ) + return float(current), True + return float((_first_zero(currents, residuals) + _last_zero(currents, residuals)) / 2.0), False + + def _strength_span(self, idx: int) -> float: + """Return the signed strength range ``f_i(I_max) - f_i(I_min)`` of the excitation curve of magnet ``idx``.""" + currents = self.__curves[idx].get_curve()[:, 0] + model = self.__sub_models[idx] + return float(model.compute_strengths([currents.max()])[0] - model.compute_strengths([currents.min()])[0]) + + def _polarities(self) -> np.ndarray: + """Return ``+1`` for each magnet with the polarity of the first magnet, ``-1`` otherwise.""" + signs = np.array([1.0 if self._strength_span(idx) >= 0.0 else -1.0 for idx in range(self.__nbMagnets)]) + return signs * signs[0] + def shares_power_supply(self) -> bool: """Return True, all the magnets of the series share one power supply.""" return True diff --git a/pyaml/magnet/serialized_magnet.py b/pyaml/magnet/serialized_magnet.py index 44cf55c1..665a7fd6 100644 --- a/pyaml/magnet/serialized_magnet.py +++ b/pyaml/magnet/serialized_magnet.py @@ -118,7 +118,12 @@ class ReadWriteSerializedStrength(ReadWriteSerializedHardware): Read/write access to the average strength of a serialized-magnet family. Magnet ``i`` of the series has the strength ``k_i = f_i(I)``. The family - strength is the average strength ``sum(f_i(I)) / N``. + strength is the average strength ``sum(p_i * f_i(I)) / N``, each strength being + taken with the polarity of the first magnet (``p_i = -1`` for a magnet wired in + reverse, ``+1`` otherwise). + + A single strength value for the whole family is only defined when every magnet has the same function: for a + series mixing functions (e.g. quadrupoles and sextupoles), use the strength of each magnet instead. Parameters ---------- @@ -126,6 +131,8 @@ class ReadWriteSerializedStrength(ReadWriteSerializedHardware): Per-magnet hardware accessors. Writing one of them moves the whole series. model : MagnetModel Serialized magnet model converting between average strength and hardware value. + homogeneous : bool, optional + Whether every magnet of the series has the same function. Default is ``True``. Methods ------- @@ -135,10 +142,23 @@ class ReadWriteSerializedStrength(ReadWriteSerializedHardware): Set the hardware value giving the requested average strength. unit() Return the physical strength unit. + + Raises + ------ + PyAMLException + On ``get()`` or ``set()`` when the series mixes magnet functions. """ + def __init__(self, hardwares: list[abstract.ReadWriteFloatScalar], model: MagnetModel, homogeneous: bool = True): + """ + Initialize a family strength accessor for serialized magnets. + """ + super().__init__(hardwares, model) + self._homogeneous = homogeneous + def get(self) -> float: """Return the average strength of the series.""" + self._check_homogeneous() return self._model.compute_mean_strength(self._hardwares[0].get()) def set(self, value: float): @@ -150,12 +170,178 @@ def set(self, value: float): value : float Average strength of the series. """ + self._check_homogeneous() self._hardwares[0].set(self._model.compute_hardware_from_mean_strength(value)) def unit(self) -> str: """Return the physical strength unit.""" return self._model.get_strength_units()[0] + def _check_homogeneous(self): + if not self._homogeneous: + raise PyAMLException( + "A single strength is not defined for a serialized family mixing magnet functions, " + "use the strength of each magnet (family.strengths or magnet.strength)" + ) + + +class ReadWriteSerializedHardwareArray(abstract.ReadWriteFloatArray): + """ + Read/write access to the hardware value of each magnet of a serialized-magnet family. + + Every magnet of the series reads the hardware value ``I`` of the shared power supply. + + Parameters + ---------- + hardwares : list[abstract.ReadWriteFloatScalar] + Per-magnet hardware accessors. Writing one of them moves the whole series. + model : MagnetModel + Serialized magnet model used for units. + + Methods + ------- + get() + Return the hardware value of each magnet. + set(value) + Set the shared hardware value to the average of the requested values. + set_and_wait(value) + Set the hardware values and wait for convergence. + unit() + Return the hardware-value units. + """ + + def __init__(self, hardwares: list[abstract.ReadWriteFloatScalar], model: MagnetModel): + """ + Initialize a per-magnet hardware accessor for serialized magnets. + """ + self._hardwares = hardwares + self._model = model + + def get(self) -> np.ndarray: + """Return the hardware value of each magnet.""" + return np.array([hardware.get() for hardware in self._hardwares]) + + def set(self, value: np.ndarray): + """ + Set the shared hardware value to the average of the requested values. + + Parameters + ---------- + value : np.ndarray + Requested hardware value of each magnet, in declaration order. + + Raises + ------ + PyAMLException + If ``value`` is not a 1-D sequence holding one value per magnet. A single value for the whole series is + written with ``family.hardware``. + """ + if np.ndim(value) != 1 or len(value) != len(self._hardwares): + raise PyAMLException( + f"One hardware value per magnet is expected ({len(self._hardwares)} values), got {np.shape(value)}" + ) + self._hardwares[0].set(float(np.mean(value))) + + def set_and_wait(self, value: np.ndarray): + """ + Set the hardware values and wait for convergence. + + Parameters + ---------- + value : np.ndarray + Requested hardware value of each magnet, in declaration order. + + Raises + ------ + NotImplementedError + Waiting for readback convergence is not implemented for this accessor. + """ + raise NotImplementedError("Not implemented yet.") + + def unit(self) -> list[str]: + """Return the hardware-value units.""" + return self._model.get_hardware_units() + + +class ReadWriteSerializedStrengthArray(abstract.ReadWriteFloatArray): + """ + Read/write access to the strength of each magnet of a serialized-magnet family. + + Magnet ``i`` of the series has the strength ``k_i = f_i(I)``. Writing the strengths ``k_i`` sets the hardware + value ``I`` solving ``sum((f_i(I) - k_i) / span_i) == 0``, ``span_i`` being the signed strength range of the + excitation curve of magnet ``i``. It is exact when the requested strengths are reachable with one hardware value, + and every magnet weighs the same whatever its unit or polarity, so this works for a series mixing magnet functions + or polarities. + + Parameters + ---------- + strengths : list[abstract.ReadWriteFloatScalar] + Per-magnet strength accessors. + hardwares : list[abstract.ReadWriteFloatScalar] + Per-magnet hardware accessors. Writing one of them moves the whole series. + model : MagnetModel + Serialized magnet model converting between strengths and hardware value. + + Methods + ------- + get() + Return the strength of each magnet. + set(value) + Set the hardware value giving the requested strengths. + set_and_wait(value) + Set the strengths and wait for convergence. + unit() + Return the physical strength units. + """ + + def __init__( + self, + strengths: list[abstract.ReadWriteFloatScalar], + hardwares: list[abstract.ReadWriteFloatScalar], + model: MagnetModel, + ): + """ + Initialize a per-magnet strength accessor for serialized magnets. + """ + self._strengths = strengths + self._hardwares = hardwares + self._model = model + + def get(self) -> np.ndarray: + """Return the strength of each magnet.""" + return np.array([strength.get() for strength in self._strengths]) + + def set(self, value: np.ndarray): + """ + Set the hardware value giving the requested strengths. + + Parameters + ---------- + value : np.ndarray + Requested strength of each magnet, in declaration order. + """ + self._hardwares[0].set(self._model.compute_hardware_from_strengths(value)) + + def set_and_wait(self, value: np.ndarray): + """ + Set the strengths and wait for convergence. + + Parameters + ---------- + value : np.ndarray + Requested strength of each magnet, in declaration order. + + Raises + ------ + NotImplementedError + Waiting for readback convergence is not implemented for this accessor. + """ + raise NotImplementedError("Not implemented yet.") + + def unit(self) -> list[str]: + """Return the physical strength units.""" + return self._model.get_strength_units() + @register_schema class SerializedMagnets(Element, DynamicValidation): @@ -171,8 +357,9 @@ class SerializedMagnets(Element, DynamicValidation): ---------- name : str Name of the serialized magnet group. - function : str - Magnet function identifier used to select the concrete virtual magnet type. + function : list[str] | str + Magnet function identifier selecting the concrete virtual magnet type. A single value applies to every + element, a list gives one function per element for a series mixing magnet functions. elements : list[str] | str Names of the individual magnets in the group. model : MagnetModel | None, optional @@ -185,9 +372,14 @@ class SerializedMagnets(Element, DynamicValidation): Attributes ---------- strength - Average strength of the magnets of the series, in physics unit + Average strength of the magnets of the series, in physics unit. Only defined when every magnet has the same + function. hardware Hardware value of the power supply shared by the series + strengths + Strength of each magnet of the series, in physics unit + hardwares + Hardware value of each magnet of the series Methods ------- @@ -195,6 +387,8 @@ class SerializedMagnets(Element, DynamicValidation): Return the number of magnets in the serialized group. get_magnets() Return the group's virtual single-function magnets. + is_homogeneous() + Return True when every magnet of the series has the same function. attach(peer, strengths, hardwares) Attach the group and its virtual magnets to a runtime peer. set_energy(energy) @@ -221,11 +415,20 @@ class SerializedMagnets(Element, DynamicValidation): Accessor ``get()`` ``set(x)`` ====================== =============================== ================================================== ``family.hardware`` ``I`` write ``I = x`` - ``family.strength`` average ``sum(f_i(I)) / N`` solve ``sum(f_i(I)) / N = x`` for ``I``, write it + ``family.strength`` ``sum(p_i * f_i(I)) / N`` solve ``sum(p_i * f_i(I)) / N = x`` for ``I``, write it + ``family.hardwares`` ``[I, ..., I]`` write ``I = mean(x)`` + ``family.strengths`` ``[f_1(I), ..., f_N(I)]`` solve ``sum((f_i(I) - x_i) / span_i) = 0``, write it ``magnet[i].hardware`` ``I`` write ``I = x``, all the magnets follow ``magnet[i].strength`` ``k_i = f_i(I)`` write ``I = f_i^-1(x)``, the others follow ====================== =============================== ================================================== + ``p_i`` is ``-1`` for a magnet wired in reverse compared to the first magnet, ``+1`` otherwise. ``span_i`` is the + signed strength range of the excitation curve of magnet ``i``, so every magnet weighs the same whatever its unit or + polarity (see :meth:`LinearSerializedMagnetModel.compute_hardware_from_strengths`). + + ``family.strength`` raises a ``PyAMLException`` for a series mixing magnet functions: a single strength value + cannot be both, say, a quadrupole and a sextupole strength. + Configuration ------------- The model can be defined inline or in a separate file. See @@ -245,7 +448,7 @@ class SerializedMagnets(Element, DynamicValidation): def __init__( self, name: str, - function: str, + function: list[str] | str, elements: list[str] | str, model: MagnetModel | None = None, description: str | None = None, @@ -256,25 +459,32 @@ def __init__( """ super().__init__(name, None, description) - self.function = function self.model = model - self.polynom = None + self.polynoms = None self.__strengths = None self.__hardwares = None + self.__strengths_array = None + self.__hardwares_array = None self.__virtuals: list[Magnet] = [] self.__elements = elements if isinstance(elements, list) else [elements] + self.functions: list[str] = function if isinstance(function, list) else [function] * len(self.__elements) + if len(self.functions) != len(self.__elements): + raise PyAMLException( + f"{name}: function must be a single value or one value per element " + f"({len(self.__elements)} elements, {len(self.functions)} functions)" + ) if peer is None: # Configuration part + for func in self.functions: + if func not in function_map: + raise PyAMLException(func + " not implemented for serialized magnet") # The submodels are built once here: the attached copies share them with the per-magnet accessors self.model.set_number_of_magnets(len(self.__elements)) - self.polynom = function_map[self.function].polynom - if self.function not in function_map: - raise PyAMLException(self.function + " not implemented for serialized magnet") - for element in self.__elements: - # Check mapping validity + self.polynoms = [function_map[func].polynom for func in self.functions] + for element, func in zip(self.__elements, self.functions, strict=True): # Create the virtual magnet for the corresponding magnet - vm = self.__create_virtual_magnet(element) + vm = self.__create_virtual_magnet(element, func) self.__virtuals.append(vm) # Register the virtual element in the factory to have a coherent factory and improve error reporting ELEMENT_REGISTRY.register(vm) @@ -285,17 +495,19 @@ def __init__( def _fill_device(self, holder: "ElementHolder") -> None: holder._fill_serialized_magnets(self) - def __create_virtual_magnet(self, name: str) -> Magnet: + def __create_virtual_magnet(self, name: str, function: str) -> Magnet: """ Create a virtual magnet for one serialized element. - The configured function mapping selects the concrete magnet class, and + The function of the element selects the concrete magnet class, and the serialized group's model is shared with the new virtual magnet. Parameters ---------- name : str Name assigned to the virtual magnet. + function : str + Function of the element, selecting the virtual magnet class. Returns ------- @@ -303,10 +515,14 @@ def __create_virtual_magnet(self, name: str) -> Magnet: Newly created virtual magnet linked to the serialized group. """ args = {"name": name, "model": self.model} - virtual: Magnet = function_map[self.function](**args) + virtual: Magnet = function_map[function](**args) virtual.set_model_name(self.get_name()) return virtual + def is_homogeneous(self) -> bool: + """Return True when every magnet of the series has the same function.""" + return len(set(self.functions)) == 1 + def get_nb_magnets(self) -> int: """Return the number of magnets in the serialized group.""" return len(self.__elements) @@ -339,11 +555,13 @@ def attach( Virtual magnets of the group, each bound to ``peer``. """ l = [] - n_ser_mag = SerializedMagnets(self._name, self.function, self.__elements, self.model, self.description, peer) + n_ser_mag = SerializedMagnets(self._name, self.functions, self.__elements, self.model, self.description, peer) if self.model.has_hardware(): n_ser_mag.__hardwares = ReadWriteSerializedHardware(hardwares, self.model) + n_ser_mag.__hardwares_array = ReadWriteSerializedHardwareArray(hardwares, self.model) if self.model.has_physics(): - n_ser_mag.__strengths = ReadWriteSerializedStrength(hardwares, self.model) + n_ser_mag.__strengths = ReadWriteSerializedStrength(hardwares, self.model, self.is_homogeneous()) + n_ser_mag.__strengths_array = ReadWriteSerializedStrengthArray(strengths, hardwares, self.model) l.append(n_ser_mag) # Construct a single magnet for each magnet. sub_magnets: list[Magnet] = [] @@ -375,6 +593,26 @@ def hardware(self) -> abstract.ReadWriteFloatScalar: raise PyAMLException(f"{str(self)} has no model that supports hardware units") return self.__hardwares + @property + def strengths(self) -> abstract.ReadWriteFloatArray: + """ + Strength of each magnet of the series, in physics unit + """ + self.check_peer() + if self.__strengths_array is None: + raise PyAMLException(f"{str(self)} has no model that supports physics units") + return self.__strengths_array + + @property + def hardwares(self) -> abstract.ReadWriteFloatArray: + """ + Hardware value of each magnet of the series, all equal to the shared power-supply value + """ + self.check_peer() + if self.__hardwares_array is None: + raise PyAMLException(f"{str(self)} has no model that supports hardware units") + return self.__hardwares_array + def set_energy(self, energy: float): """ Set beam energy for serialized-magnet strength conversion. diff --git a/tests/config/sr_serialized_magnets.yaml b/tests/config/sr_serialized_magnets.yaml index a8a3c486..bbb7fb88 100644 --- a/tests/config/sr_serialized_magnets.yaml +++ b/tests/config/sr_serialized_magnets.yaml @@ -37,6 +37,22 @@ devices: type: pyaml.magnet.csvcurve file: sr/magnet_models/serialized_strength.csv powerconverter: srmag/ps-corr-sh1/c01-a-ch1/current + - type: pyaml.magnet.serialized_magnet + name: mixedFunctionSeries + function: [B1, B2] + elements: + - QF6D-C01 + - SD1A-C01 + model: + type: pyaml.magnet.linear_serialized_model + calibration_factors: 1.0 + calibration_offsets: 0.0 + unit: m-1 + hardware_unit: A + curves: + type: pyaml.magnet.csvcurve + file: sr/magnet_models/serialized_strength.csv + powerconverter: srmag/ps-mixed-series/current - type: pyaml.magnet.serialized_magnet name: QF1A function: B1 diff --git a/tests/config/sr_serialized_magnets_range.yaml b/tests/config/sr_serialized_magnets_range.yaml new file mode 100644 index 00000000..0e18339f --- /dev/null +++ b/tests/config/sr_serialized_magnets_range.yaml @@ -0,0 +1,48 @@ +type: pyaml.accelerator +facility: ESRF +machine: sr +energy: 6e9 +simulators: + - type: pyaml.lattice.simulator + lattice: sr/lattices/ebs.mat + name: design +controls: + - type: tango.pyaml.controlsystem + tango_host: ebs-simu-3:10000 + name: live + catalog: + type: tango.pyaml.static_catalog + entries: + - type: tango.pyaml.static_catalog_entry + key: srmag/ps-qf1a/series/current + device: + type: tango.pyaml.attribute + attribute: srmag/ps-qf1a/series/current + unit: A + range: [0, 90] +data_folder: /data/store +arrays: + - type: pyaml.arrays.magnet + name: QF1A_PART + elements: + - QF1A-C06 + - QF1A-C07 +devices: + - type: pyaml.magnet.serialized_magnet + name: QF1A + function: B1 + elements: + - QF1A-C05 + - QF1A-C06 + - QF1A-C07 + - QF1A-C08 + model: + type: pyaml.magnet.linear_serialized_model + calibration_factors: [1.0, 1.01, 0.99, 1.02] + calibration_offsets: 0.0 + unit: m-1 + hardware_unit: A + curves: + type: pyaml.magnet.csvcurve + file: sr/magnet_models/QF1_strength.csv + powerconverter: srmag/ps-qf1a/series/current diff --git a/tests/magnet/test_serialized_magnets.py b/tests/magnet/test_serialized_magnets.py index 564ab6c1..2aae7685 100644 --- a/tests/magnet/test_serialized_magnets.py +++ b/tests/magnet/test_serialized_magnets.py @@ -1,8 +1,11 @@ import numpy as np import pytest +from pyaml import PyAMLException from pyaml.accelerator import Accelerator +from pyaml.magnet.quadrupole import Quadrupole from pyaml.magnet.serialized_magnet import SerializedMagnets +from pyaml.magnet.sextupole import Sextupole def check_no_diff(array: list[np.float64]) -> bool: @@ -202,3 +205,100 @@ def test_strength_computation(sr_file): # The family strength is the average strength of the magnets assert abs(np.mean(magnets_from_lattice_strengths) - 0.8) < 1e-6 assert magnets_strengths == magnets_from_lattice_strengths + + +SR_FILE = "tests/config/sr_serialized_magnets.yaml" +# A quadrupole (B1) and a sextupole (B2) on the same power supply +MIXED = "mixedFunctionSeries" + + +def _effective_length(element) -> float: + return element.Length if element.Length else 1.0 + + +def test_partial_strength_request_keeps_one_hardware_value_per_power_supply(): + sr: Accelerator = Accelerator.load(SR_FILE, include_locations=False, ignore_external=True) + family: SerializedMagnets = sr.design.serialized_magnet.get("QF1A") + model = family.model + current = np.array([family.hardware.get()]) + + assert model.compute_hardware_values_partial(current, {1: 0.8}).shape == current.shape + with pytest.raises(PyAMLException, match="current_hardware"): + model.compute_hardware_values_partial(np.array([80.0, 80.0]), {1: 0.8}) + + +def test_series_with_one_function_is_homogeneous(): + sr: Accelerator = Accelerator.load(SR_FILE, include_locations=False, ignore_external=True) + assert sr.design.serialized_magnet.get("QF1A").is_homogeneous() + assert not sr.design.serialized_magnet.get(MIXED).is_homogeneous() + + +def test_mixed_function_series_has_no_single_family_strength(): + sr: Accelerator = Accelerator.load(SR_FILE, include_locations=False, ignore_external=True) + family: SerializedMagnets = sr.design.serialized_magnet.get(MIXED) + with pytest.raises(PyAMLException, match="mixing magnet functions"): + family.strength.get() + with pytest.raises(PyAMLException, match="mixing magnet functions"): + family.strength.set(0.5) + + +def test_mixed_function_series_keeps_hardware_and_per_magnet_access(): + sr: Accelerator = Accelerator.load(SR_FILE, include_locations=False, ignore_external=True) + family: SerializedMagnets = sr.design.serialized_magnet.get(MIXED) + quad, sextu = family.get_magnets() + assert isinstance(quad, Quadrupole) + assert isinstance(sextu, Sextupole) + + family.hardware.set(50.0) + + assert family.hardware.get() == pytest.approx(50.0) + assert np.allclose(family.hardwares.get(), 50.0) + assert np.allclose(family.strengths.get(), [quad.strength.get(), sextu.strength.get()]) + + +def test_mixed_function_series_writes_each_magnet_in_its_own_polynom(): + sr: Accelerator = Accelerator.load(SR_FILE, include_locations=False, ignore_external=True) + family: SerializedMagnets = sr.design.serialized_magnet.get(MIXED) + ring = sr.design.get_lattice() + quad_element = ring.get_elements("QF6D-C01")[0] + sextu_element = ring.get_elements("SD1A-C01")[0] + sextu_b1_before = sextu_element.PolynomB[1] + + family.hardware.set(50.0) + + quad_strength = family.model.get_sub_model(0).compute_strengths([50.0])[0] + sextu_strength = family.model.get_sub_model(1).compute_strengths([50.0])[0] + assert quad_element.PolynomB[1] * _effective_length(quad_element) == pytest.approx(quad_strength) + assert sextu_element.PolynomB[2] * _effective_length(sextu_element) == pytest.approx(sextu_strength) + assert sextu_element.PolynomB[1] == sextu_b1_before + + +def test_mixed_function_series_strengths_converge_for_a_consistent_request(caplog): + sr: Accelerator = Accelerator.load(SR_FILE, include_locations=False, ignore_external=True) + family: SerializedMagnets = sr.design.serialized_magnet.get(MIXED) + family.hardware.set(50.0) + requested = family.strengths.get() + family.hardware.set(30.0) + + with caplog.at_level("WARNING"): + family.strengths.set(requested) + + assert family.hardware.get() == pytest.approx(50.0) + assert np.allclose(family.strengths.get(), requested) + assert "Inconsistent strengths" not in caplog.text + + +def test_mixed_function_series_strengths_warn_for_an_inconsistent_request(caplog): + sr: Accelerator = Accelerator.load(SR_FILE, include_locations=False, ignore_external=True) + family: SerializedMagnets = sr.design.serialized_magnet.get(MIXED) + family.hardware.set(50.0) + requested = family.strengths.get() * np.array([1.0, 1.1]) + + with caplog.at_level("WARNING"): + family.strengths.set(requested) + + # One power supply: the series gets the hardware value solving sum((f_i(I) - k_i) / span_i) == 0. Both magnets + # of this family share the same curve and calibration, so their spans are equal and the obtained strengths + # average to the requested ones + assert np.mean(family.strengths.get()) == pytest.approx(np.mean(requested)) + assert "Inconsistent strengths" in caplog.text diff --git a/tests/magnet/test_serialized_magnets_cs.py b/tests/magnet/test_serialized_magnets_cs.py index 1cc7a3e3..d33d7f51 100644 --- a/tests/magnet/test_serialized_magnets_cs.py +++ b/tests/magnet/test_serialized_magnets_cs.py @@ -1,6 +1,7 @@ import numpy as np import pytest +from pyaml import PyAMLException from pyaml.accelerator import Accelerator from pyaml.magnet.serialized_magnet import SerializedMagnets @@ -12,6 +13,20 @@ def _sub_strengths(family: SerializedMagnets, current: float) -> np.ndarray: return np.array([family.model.get_sub_model(i).compute_strengths([current])[0] for i in range(family.get_nb_magnets())]) +# Current range of the excitation curve sr/magnet_models/QF1_strength.csv +QF1_CURRENT_RANGE = (0.0, 110.0) + + +def _normalized_residual(family: SerializedMagnets, current: float, requested, indices) -> float: + """Return sum((f_i(I) - k_i) / span_i) over the requested magnets, span_i being the strength range of their curve.""" + total = 0.0 + for idx, strength in zip(indices, requested, strict=True): + model = family.model.get_sub_model(idx) + span = model.compute_strengths([QF1_CURRENT_RANGE[1]])[0] - model.compute_strengths([QF1_CURRENT_RANGE[0]])[0] + total += (model.compute_strengths([current])[0] - strength) / span + return total + + @pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) def test_serialized_magnets_attach_to_control_system(install_test_package): sr: Accelerator = Accelerator.load(CONFIG) @@ -139,17 +154,22 @@ def test_design_and_live_agree(install_test_package): @pytest.mark.parametrize("mode", MODES) @pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) -def test_array_with_one_magnet_of_a_series_moves_the_whole_series(install_test_package, mode): +def test_array_with_one_magnet_of_a_series_drives_the_series_with_its_own_submodel(install_test_package, mode, caplog): from pyaml.arrays.magnet_array import MagnetArray sr: Accelerator = Accelerator.load(CONFIG) family = _family(sr, mode) family.hardware.set(80.0) array = MagnetArray("ONE", [family.get_magnets()[1]]) - array.strengths.set(np.array([0.5])) + + with caplog.at_level("WARNING"): + array.strengths.set(np.array([0.5])) + + # Same as magnets[1].strength.set(0.5): I = f_1^-1(0.5), the other magnets follow as f_j(I) current = family.model.get_sub_model(1).compute_hardware_values([0.5])[0] assert family.hardware.get() == pytest.approx(current) - assert array.strengths.get() == pytest.approx([0.5]) + assert np.allclose(family.strengths.get(), _sub_strengths(family, current)) + assert "Inconsistent strengths" not in caplog.text @pytest.mark.parametrize("mode", MODES) @@ -178,8 +198,141 @@ def test_array_with_inconsistent_magnets_of_a_series(install_test_package, mode, strengths[:4] = requested with caplog.at_level("WARNING"): array.strengths.set(strengths) - # The power supply gets the mean of the currents asked by the magnets of the series - expected = np.mean([qf1a.model.get_sub_model(i).compute_hardware_values([s])[0] for i, s in enumerate(requested)]) - assert qf1a.hardware.get() == pytest.approx(expected) + # Every magnet of QF1A is in the array: the power supply gets the I solving sum((f_i(I) - k_i) / span_i) == 0 + assert _normalized_residual(qf1a, qf1a.hardware.get(), requested, range(4)) == pytest.approx(0.0, abs=1e-9) assert qd2a.hardware.get() == pytest.approx(20.0) assert "Inconsistent strengths" in caplog.text + + +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_array_with_part_of_a_series_drives_the_series_from_the_requested_magnets_only(install_test_package, mode): + from pyaml.arrays.magnet_array import MagnetArray + + sr: Accelerator = Accelerator.load(CONFIG) + family = _family(sr, mode) + family.hardware.set(80.0) + magnets = family.get_magnets() + array = MagnetArray("PART", [magnets[1], magnets[2]]) + requested = np.array([0.78, 0.82]) + + array.strengths.set(requested) + + # Only magnets 1 and 2 define I: sum((f_i(I) - k_i) / span_i) == 0 over them. Magnets 0 and 3 follow as f_j(I) + current = family.hardware.get() + assert _normalized_residual(family, current, requested, [1, 2]) == pytest.approx(0.0, abs=1e-9) + assert np.allclose(family.strengths.get(), _sub_strengths(family, current)) + + +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_strength_outside_the_excitation_curve_is_clamped_with_a_warning(install_test_package, mode, caplog): + sr: Accelerator = Accelerator.load(CONFIG) + family = _family(sr, mode) + family.hardware.set(80.0) + # Highest current of the QF1 excitation curve + curve_end = family.model.get_sub_model(0).compute_hardware_values([1e3])[0] + + with caplog.at_level("WARNING"): + family.strengths.set(np.full(family.get_nb_magnets(), 1e3)) + + assert family.hardware.get() == pytest.approx(curve_end) + assert "outside the range" in caplog.text + assert "Inconsistent strengths" not in caplog.text + + +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_family_strengths_and_hardwares_read_every_magnet(install_test_package, mode): + sr: Accelerator = Accelerator.load(CONFIG) + family = _family(sr, mode) + family.hardware.set(80.0) + + # One call reads every magnet of the series: k_i = f_i(I), and the shared current I + assert np.allclose(family.strengths.get(), [m.strength.get() for m in family.get_magnets()]) + assert np.allclose(family.strengths.get(), _sub_strengths(family, 80.0)) + assert np.allclose(family.hardwares.get(), 80.0) + + +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_family_strengths_set_balances_the_normalized_errors(install_test_package, mode, caplog): + sr: Accelerator = Accelerator.load(CONFIG) + family = _family(sr, mode) + family.hardware.set(80.0) + requested = np.array([0.78, 0.79, 0.80, 0.81]) + + with caplog.at_level("WARNING"): + family.strengths.set(requested) + + # One power supply: the series gets the I solving sum((f_i(I) - k_i) / span_i) == 0 + current = family.hardware.get() + assert _normalized_residual(family, current, requested, range(4)) == pytest.approx(0.0, abs=1e-9) + # QF1A has a different calibration factor per magnet: these strengths cannot all be reached with one current + assert "Inconsistent strengths" in caplog.text + + +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_family_hardwares_set_applies_the_average_current(install_test_package, mode): + sr: Accelerator = Accelerator.load(CONFIG) + family = _family(sr, mode) + + family.hardwares.set(np.array([79.0, 80.0, 81.0, 84.0])) + + assert family.hardware.get() == pytest.approx(81.0) + assert np.allclose(family.hardwares.get(), 81.0) + + +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_family_array_set_of_get_does_not_change_the_machine(install_test_package, mode, caplog): + sr: Accelerator = Accelerator.load(CONFIG) + family = _family(sr, mode) + family.hardware.set(80.0) + with caplog.at_level("WARNING"): + for accessor in [family.strengths, family.hardwares]: + accessor.set(accessor.get()) + assert family.hardware.get() == pytest.approx(80.0) + assert "Inconsistent strengths" not in caplog.text + + +# QF1A on a power supply limited to [0, 90] A, with an array holding magnets C06 and C07 only +RANGE_CONFIG = "tests/config/sr_serialized_magnets_range.yaml" + + +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_power_supply_range_protects_every_write_to_a_series(install_test_package): + sr: Accelerator = Accelerator.load(RANGE_CONFIG) + family: SerializedMagnets = sr.live.serialized_magnet.get("QF1A") + magnets = family.get_magnets() + family.hardware.set(80.0) + # Strengths that need 100 A, outside the range of the power supply + too_strong = _sub_strengths(family, 100.0) + + writes = [ + lambda: sr.live.magnets.get("QF1A_PART").strengths.set(too_strong[[1, 2]]), + lambda: family.strengths.set(too_strong), + lambda: family.strength.set(np.mean(too_strong)), + lambda: magnets[1].strength.set(too_strong[1]), + lambda: family.hardware.set(100.0), + ] + for write in writes: + with pytest.raises(PyAMLException): + write() + # Nothing is written: the series stays at 80 A + assert family.hardware.get() == pytest.approx(80.0) + + +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_array_with_part_of_a_series_within_the_power_supply_range(install_test_package): + sr: Accelerator = Accelerator.load(RANGE_CONFIG) + family: SerializedMagnets = sr.live.serialized_magnet.get("QF1A") + family.hardware.set(80.0) + within = _sub_strengths(family, 70.0) + + sr.live.magnets.get("QF1A_PART").strengths.set(within[[1, 2]]) + + # Magnets C06 and C07 ask for 70 A, magnets C05 and C08 follow + assert family.hardware.get() == pytest.approx(70.0) + assert np.allclose(family.strengths.get(), within) diff --git a/tests/magnet/test_serialized_magnets_expected_behaviour.py b/tests/magnet/test_serialized_magnets_expected_behaviour.py new file mode 100644 index 00000000..e1166cb4 --- /dev/null +++ b/tests/magnet/test_serialized_magnets_expected_behaviour.py @@ -0,0 +1,405 @@ +""" +Expected behaviour of a series of magnets, written as a basis for discussion. + +What is a series? + N magnets wired on ONE power supply. There is a single current I for all of them. + Magnet i has its own excitation curve f_i: its strength is k_i = f_i(I). + +What happens when we ask for strengths? + We ask for one strength k_i per magnet. The power supply can only take one current, so the + code looks for the current I that best matches all the requested strengths: + + sum over i of (f_i(I) - k_i) / span_i == 0 + + span_i is the strength range of the curve of magnet i, f_i(I_max) - f_i(I_min), with its sign. + Dividing by span_i makes every term dimensionless and increasing with I, so: + - every magnet weighs the same, whatever its unit (a sextupole does not outweigh a quadrupole), + - a magnet wired in reverse (negative calibration) is handled like the others, + - the equation always has a solution. + + A request is "consistent" when one single current gives all the requested strengths: the code + must then find this current exactly. Otherwise the request is "inconsistent": the code returns a + compromise and logs a warning. + +How to read and modify these tests? + Each test builds a small series with simple curves (current in A, strength = curve value x + calibration), asks for strengths, and checks the current. To discuss another expected behaviour, + change the curves, the calibrations, the request or the expected current in the test. +""" + +import numpy as np +import pytest + +from pyaml import PyAMLException +from pyaml.accelerator import Accelerator +from pyaml.magnet.inline_curve import InlineCurve +from pyaml.magnet.linear_serialized_model import LinearSerializedMagnetModel + +# Excitation curves: each row is [current (A), strength]. +# Straight line, 1.0 at 100 A. +LINEAR = np.array([[0.0, 0.0], [100.0, 1.0]]) +# Steep up to 50 A, then flatter: the magnet starts to saturate. +SATURATING = np.array([[0.0, 0.0], [50.0, 0.6], [100.0, 0.8]]) +# Straight line up to 110 A, used for the quadrupole and the sextupole example. +QUAD_SEXTU = np.array([[0.0, 0.0], [110.0, 1.1]]) +# Fully saturated above 80 A: the strength stays at 0.8 from 80 A to 110 A. +PLATEAU = np.array([[0.0, 0.0], [80.0, 0.8], [110.0, 0.8]]) +# Two slopes: steep up to 50 A, flatter after. +TWO_SLOPES = np.array([[0.0, 0.0], [50.0, 1.0], [100.0, 1.5]]) + + +def make_series(curves: list[np.ndarray], calibrations: list[float]) -> LinearSerializedMagnetModel: + """Build the model of a series: one curve and one calibration factor per magnet.""" + series = LinearSerializedMagnetModel( + curves=[InlineCurve(mat=curve) for curve in curves], + calibration_factors=calibrations, + calibration_offsets=0.0, + unit="m-1", + hardware_unit="A", + powerconverter="the/power/supply", + ) + series.set_number_of_magnets(len(curves)) + # Magnetic rigidity of 1: the strength is directly curve value x calibration + series.set_magnet_rigidity(1.0) + return series + + +def relative_errors(series: LinearSerializedMagnetModel, current: float, requested: list[float]) -> list[str]: + """Return, for each magnet, how far the strength obtained at `current` is from the requested one.""" + obtained = series.compute_strengths([current]) + return [f"{100 * (got - wanted) / wanted:+.1f} %" for got, wanted in zip(obtained, requested, strict=True)] + + +# --------------------------------------------------------------------------------------------------------------------- +# Magnets wired in reverse (opposite polarities) +# --------------------------------------------------------------------------------------------------------------------- + + +def test_opposite_polarities_with_different_saturations_reach_a_consistent_request(): + # Step 1: a series of two magnets on ONE power supply. + # - magnet 0 follows the LINEAR curve, calibration +1 + # - magnet 1 follows the SATURATING curve, calibration -1: it is wired in reverse, + # its strength goes down when the current goes up + series = make_series([LINEAR, SATURATING], calibrations=[1.0, -1.0]) + + # Step 2: ask for the strengths the two magnets really have at 40 A. + # This request is "consistent": one single current (40 A) gives both strengths. + requested = series.compute_strengths([40.0]) + + # Step 3: the power-supply current computed from these strengths. + current = series.compute_hardware_from_strengths(requested) + + # Expected: we find 40 A back, even though the two magnets have opposite polarities. + assert current == pytest.approx(40.0) + + +def test_opposite_polarities_with_identical_curves_reach_a_consistent_request(): + # Step 1: two magnets on ONE power supply, same LINEAR curve, the second one wired in reverse. + # At any current, their strengths are exactly opposite: k_1 = -k_0. + series = make_series([LINEAR, LINEAR], calibrations=[1.0, -1.0]) + + # Step 2: ask for the strengths they have at 40 A: +0.4 and -0.4. + requested = series.compute_strengths([40.0]) + + # Step 3: the power-supply current computed from these strengths. + current = series.compute_hardware_from_strengths(requested) + + # Expected: 40 A. The two strengths add up to 0, but this must not give 0 A. + assert current == pytest.approx(40.0) + + +def test_family_strength_of_opposite_polarities_uses_the_polarity_of_the_first_magnet(): + # Step 1: the same two magnets, the second one wired in reverse. + series = make_series([LINEAR, LINEAR], calibrations=[1.0, -1.0]) + + # Step 2: read the family strength at 40 A. + # The raw strengths are +0.4 and -0.4. The family strength takes each strength with the + # polarity of the first magnet, so the reversed magnet counts as +0.4, not -0.4. + family_strength = series.compute_mean_strength(40.0) + + # Expected: 0.4 (and not 0.0, the plain average of +0.4 and -0.4). + assert family_strength == pytest.approx(0.4) + + # Step 3: set this family strength back. + # Expected: 40 A. + assert series.compute_hardware_from_mean_strength(0.4) == pytest.approx(40.0) + + +def test_family_strength_of_opposite_polarities_with_different_saturations(): + # Step 1: LINEAR magnet, plus a SATURATING magnet wired in reverse. + series = make_series([LINEAR, SATURATING], calibrations=[1.0, -1.0]) + + # Step 2: at 40 A, the strengths are +0.40 and -0.48. + # With the polarity of the first magnet, the family strength is (0.40 + 0.48) / 2 = 0.44. + assert series.compute_mean_strength(40.0) == pytest.approx(0.44) + + # Step 3: set this family strength back. + # Expected: 40 A. + assert series.compute_hardware_from_mean_strength(0.44) == pytest.approx(40.0) + + +# --------------------------------------------------------------------------------------------------------------------- +# Magnets of different natures on one power supply (strengths with very different magnitudes) +# --------------------------------------------------------------------------------------------------------------------- + + +def test_quadrupole_and_sextupole_share_the_error_when_the_sextupole_asks_for_more(): + # Step 1: a "quadrupole" (calibration 1) and a "sextupole" (calibration 60) on the same curve. + # The sextupole strengths are 60 times bigger numbers than the quadrupole strengths. + series = make_series([QUAD_SEXTU, QUAD_SEXTU], calibrations=[1.0, 60.0]) + + # Step 2: an inconsistent request. + # - the quadrupole alone would need 50 A (0.5) + # - the sextupole alone would need 60 A (36 = 60 x 0.6) + # No single current gives both: the result is a compromise. + requested = [0.5, 36.0] + + # Step 3: the power-supply current computed from these strengths. + current = series.compute_hardware_from_strengths(np.array(requested)) + + # Expected: 55 A, halfway. Both magnets miss their target by about the same amount, the bigger + # numbers of the sextupole do not make it win. + assert current == pytest.approx(55.0), f"relative errors: {relative_errors(series, current, requested)}" + + +def test_quadrupole_and_sextupole_share_the_error_when_the_quadrupole_asks_for_more(): + # Step 1: the same quadrupole (calibration 1) and sextupole (calibration 60). + series = make_series([QUAD_SEXTU, QUAD_SEXTU], calibrations=[1.0, 60.0]) + + # Step 2: an inconsistent request, the other way round. + # - the quadrupole alone would need 60 A (0.6) + # - the sextupole alone would need 50 A (30 = 60 x 0.5) + requested = [0.6, 30.0] + + # Step 3: the power-supply current computed from these strengths. + current = series.compute_hardware_from_strengths(np.array(requested)) + + # Expected: 55 A, halfway again. + assert current == pytest.approx(55.0), f"relative errors: {relative_errors(series, current, requested)}" + + +# --------------------------------------------------------------------------------------------------------------------- +# Particular situations +# --------------------------------------------------------------------------------------------------------------------- + + +def test_two_saturated_magnets_give_the_middle_of_the_plateau(): + # Step 1: two magnets with the PLATEAU curve: from 80 A to 110 A the strength stays at 0.8. + series = make_series([PLATEAU, PLATEAU], calibrations=[1.0, 1.0]) + + # Step 2: ask for 0.8 on both magnets. Every current between 80 A and 110 A gives it. + current = series.compute_hardware_from_strengths(np.array([0.8, 0.8])) + + # Expected: the middle of the interval, 95 A, rather than one of its ends. + assert current == pytest.approx(95.0) + + +def test_strengths_outside_the_curves_are_clamped_with_a_warning(caplog): + # Step 1: two magnets with the LINEAR curve: the strongest they can be is 1.0, at 100 A. + series = make_series([LINEAR, LINEAR], calibrations=[1.0, 1.0]) + + # Step 2: ask for 5.0, far more than the curves can give. + with caplog.at_level("WARNING"): + current = series.compute_hardware_from_strengths(np.array([5.0, 5.0])) + + # Expected: the current stops at the end of the curves (100 A), and a warning says so. + assert current == pytest.approx(100.0) + assert "outside the range reachable by the excitation curves" in caplog.text + + +def test_a_curve_with_no_strength_range_cannot_give_a_current(): + # Step 1: magnet 1 has a flat curve: its strength is 0.3 whatever the current. + flat = np.array([[0.0, 0.3], [100.0, 0.3]]) + series = make_series([LINEAR, flat], calibrations=[1.0, 1.0]) + + # Step 2: ask for strengths. The flat magnet says nothing about the current. + # Expected: an error, rather than a current computed from a meaningless strength. + with pytest.raises(PyAMLException, match="no strength range"): + series.compute_hardware_from_strengths(np.array([0.4, 0.3])) + + +def test_the_model_has_a_single_rule_to_turn_strengths_into_a_current(): + # Step 1: two magnets with the TWO_SLOPES curve, calibrations 1.0 and 1.2. + series = make_series([TWO_SLOPES, TWO_SLOPES], calibrations=[1.0, 1.2]) + + # Step 2: the same request through the two methods of the model that turn strengths into a current. + requested = [0.8, 1.6] + from_hardware_values = series.compute_hardware_values(requested)[0] + from_strengths = series.compute_hardware_from_strengths(requested) + + # Expected: the same current. + assert from_hardware_values == pytest.approx(from_strengths) + + +# --------------------------------------------------------------------------------------------------------------------- +# One hardware value per magnet, on a real series +# --------------------------------------------------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "wrong_value", + [ + pytest.param([10.0, 20.0], id="2 values for 31 magnets"), + pytest.param(7.0, id="a single number"), + pytest.param([80.0] * 34, id="34 values for 31 magnets"), + ], +) +def test_family_hardwares_need_one_value_per_magnet(wrong_value): + # Step 1: the QF1A series of the test lattice: 31 magnets on one power supply. + sr = Accelerator.load("tests/config/sr_serialized_magnets.yaml", include_locations=False, ignore_external=True) + family = sr.design.serialized_magnet.get("QF1A") + family.hardware.set(80.0) + + # Step 2: write a list of hardware values that does not have one value per magnet. + # Expected: an error, and the power supply is not moved. + with pytest.raises(PyAMLException, match="One hardware value per magnet"): + family.hardwares.set(wrong_value) + assert family.hardware.get() == pytest.approx(80.0) + + # Step 3: with one value per magnet it works (a single value is written with family.hardware). + family.hardwares.set([75.0] * family.get_nb_magnets()) + assert family.hardware.get() == pytest.approx(75.0) + + +# --------------------------------------------------------------------------------------------------------------------- +# Review of the normalized residual (PR #485): expected behaviours that the current code does not have yet. +# These tests fail on purpose, to illustrate the cases to discuss. +# --------------------------------------------------------------------------------------------------------------------- + +# The LINEAR curve tabulated up to 200 A instead of 100 A: same magnet, same slope, longer table. +LINEAR_UP_TO_200 = np.array([[0.0, 0.0], [200.0, 2.0]]) +# The LINEAR curve with one more point at 300 A, on the same straight line: nothing changes below 100 A. +LINEAR_EXTENDED = np.array([[0.0, 0.0], [100.0, 1.0], [300.0, 3.0]]) + +TANGO_DUMMY = [{"name": "tango-pyaml", "path": "tests/dummy_cs/tango-pyaml"}] + + +def plateau_series_fragment(calibrations: list[float]) -> dict: + """Return a configuration fragment adding a series of two PLATEAU magnets, and an array holding them.""" + return { + "arrays": [{"type": "pyaml.arrays.magnet", "name": "PLATEAU_MAGNETS", "elements": ["QF4A-C10", "QF4B-C10"]}], + "devices": [ + { + "type": "pyaml.magnet.serialized_magnet", + "name": "plateauSeries", + "function": "B1", + "elements": ["QF4A-C10", "QF4B-C10"], + "model": { + "type": "pyaml.magnet.linear_serialized_model", + "calibration_factors": calibrations, + "calibration_offsets": 0.0, + "unit": "m-1", + "hardware_unit": "A", + "powerconverter": "srmag/ps-plateau/series/current", + "curves": {"type": "pyaml.magnet.inline_curve", "mat": PLATEAU.tolist()}, + }, + } + ], + } + + +def load_with_plateau_series(accelerator_from_fragments, calibrations: list[float]) -> Accelerator: + """Load the live test configuration plus a series of two PLATEAU magnets on one power supply.""" + return accelerator_from_fragments("tests/config/sr_serialized_magnets_cs.yaml", plateau_series_fragment(calibrations)) + + +# Point 1: the weight of a magnet must not depend on how far its curve is tabulated + + +@pytest.mark.parametrize( + "other_curve", + [ + pytest.param(LINEAR_UP_TO_200, id="same line tabulated up to 200 A"), + pytest.param(LINEAR_EXTENDED, id="same line with one more point at 300 A"), + ], +) +def test_the_compromise_does_not_depend_on_how_far_a_curve_is_tabulated(other_curve): + # Step 1: two magnets with the same slope (1.0 at 100 A). + # In the reference series, both curves stop at 100 A. + # In the other series, the second curve is the same straight line, only tabulated further. + reference = make_series([LINEAR, LINEAR], calibrations=[1.0, 1.0]) + other = make_series([LINEAR, other_curve], calibrations=[1.0, 1.0]) + + # Step 2: an inconsistent request: magnet 0 alone would need 40 A, magnet 1 alone 60 A. + requested = np.array([0.4, 0.6]) + + # Step 3: the power-supply current computed from these strengths, for both series. + reference_current = reference.compute_hardware_from_strengths(requested) + other_current = other.compute_hardware_from_strengths(requested) + + # Expected: the same current (50 A). The two magnets are physically identical, so the end of a table, which + # anybody can extend when a curve is measured again, must not change the compromise. + assert reference_current == pytest.approx(50.0) + assert other_current == pytest.approx(reference_current) + + +# Point 2: on a plateau, writing back the strengths that were read must not move the power supply + + +@pytest.mark.parametrize("current", [85.0, 100.0]) +def test_writing_back_the_strengths_read_on_a_plateau_keeps_the_current(current): + # Step 1: two PLATEAU magnets: from 80 A to 110 A the strength stays at 0.8. + series = make_series([PLATEAU, PLATEAU], calibrations=[1.0, 1.0]) + + # Step 2: the power supply is at `current`, on the plateau. Read the strengths and write them back, as an array + # holding both magnets does. + requested = series.compute_strengths([current]) + new_current = series.compute_hardware_values_partial(np.array([current]), {0: requested[0], 1: requested[1]})[0] + + # Expected: the current does not move. Every current of the plateau gives the requested strengths, and the + # present one is among them: there is no reason to move the power supply (today it goes to 95 A). + assert new_current == pytest.approx(current) + + +def test_writing_back_the_strength_of_one_magnet_read_on_a_plateau_keeps_the_current(): + # Step 1: the same two PLATEAU magnets, the power supply at 85 A. + series = make_series([PLATEAU, PLATEAU], calibrations=[1.0, 1.0]) + current = 85.0 + + # Step 2: an array holding only magnet 0 writes back the strength it read. + requested = series.compute_strengths([current])[0] + new_current = series.compute_hardware_values_partial(np.array([current]), {0: requested})[0] + + # Expected: the current does not move (today it goes to 95 A, the middle of the plateau). + assert new_current == pytest.approx(current) + + +@pytest.mark.parametrize("accessor", ["family.strengths", "family.strength", "magnet.strength", "array.strengths"]) +@pytest.mark.parametrize("install_test_package", TANGO_DUMMY, indirect=True) +def test_every_accessor_keeps_the_current_when_writing_back_on_a_plateau( + install_test_package, accelerator_from_fragments, accessor +): + # Step 1: in live, a series of two PLATEAU magnets on one power supply, set at 85 A (on the plateau). + sr = load_with_plateau_series(accelerator_from_fragments, calibrations=[1.0, 1.0]) + family = sr.live.serialized_magnet.get("plateauSeries") + family.hardware.set(85.0) + + # Step 2: read a strength and write it back, through one accessor. + if accessor == "family.strengths": + family.strengths.set(family.strengths.get()) + elif accessor == "family.strength": + family.strength.set(family.strength.get()) + elif accessor == "magnet.strength": + magnet = sr.live.magnet.get("QF4A-C10") + magnet.strength.set(magnet.strength.get()) + else: + array = sr.live.magnets.get("PLATEAU_MAGNETS") + array.strengths.set(array.strengths.get()) + + # Expected: the power supply stays at 85 A whatever the accessor. Today it goes to 95 A or 110 A depending on + # the accessor, so the same request does not give the same current. + assert family.hardware.get() == pytest.approx(85.0) + + +# Point 3: the family strength of a series mixing polarities + + +def test_family_strength_does_not_depend_on_the_order_of_the_magnets(): + # Step 1: the same two LINEAR magnets, one of them wired in reverse, declared in both orders. + declared_one_way = make_series([LINEAR, LINEAR], calibrations=[1.0, -1.0]) + declared_the_other_way = make_series([LINEAR, LINEAR], calibrations=[-1.0, 1.0]) + + # Step 2: read the family strength at 40 A. + # Expected: the same value. Swapping two lines of the configuration changes nothing in the machine, so it must + # not change the sign of the family strength (today +0.4 one way, -0.4 the other way). + assert declared_one_way.compute_mean_strength(40.0) == pytest.approx(declared_the_other_way.compute_mean_strength(40.0))