From c09443eba78a0a693913bb35d57f45cf49363a8f Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 31 Jul 2024 17:15:10 -0600 Subject: [PATCH 01/26] coil current scale attribute --- desc/coils.py | 71 +++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 60 insertions(+), 11 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index aacb4f19a3..f482e324f2 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -139,12 +139,17 @@ class _Coil(_MagneticField, Optimizable, ABC): ---------- current : float Current through the coil, in Amperes. + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). + """ - _io_attrs_ = _MagneticField._io_attrs_ + ["_current"] + _io_attrs_ = _MagneticField._io_attrs_ + ["_current", "_current_scale"] - def __init__(self, current, *args, **kwargs): + def __init__(self, current, current_scale=1, *args, **kwargs): self._current = float(np.squeeze(current)) + self._current_scale = float(np.squeeze(current)) super().__init__(*args, **kwargs) @optimizable_parameter @@ -158,6 +163,16 @@ def current(self, new): assert jnp.isscalar(new) or new.size == 1 self._current = float(np.squeeze(new)) + @property + def current_scale(self): + """float: Current passing through the coil, in Amperes.""" + return self._current_scale + + @current_scale.setter + def current_scale(self, new): + assert jnp.isscalar(new) or new.size == 1 + self._current_scale = float(np.squeeze(new)) + @property def num_coils(self): """int: Number of coils.""" @@ -257,7 +272,10 @@ def compute_magnetic_field( data["x"] = rpz2xyz(data["x"]) B = biot_savart_quad( - coords, data["x"], data["x_s"] * data["ds"][:, None], current + coords, + data["x"], + data["x_s"] * data["ds"][:, None], + current * self.current_scale, ) if basis.lower() == "rpz": @@ -899,20 +917,27 @@ def from_values( def _check_type(coil0, coil): + errorif( + coil0.current_scale != coil.current_scale, + ValueError, + "Coils in a CoilSet must all have the same current scale, got scales " + + f"{coil0.current_scale} and {coil.current_scale}. " + + "Consider using a MixedCoilSet.", + ) errorif( not isinstance(coil, coil0.__class__), TypeError, ( - "coils in a CoilSet must all be the same type, got types " - + f"{type(coil0)}, {type(coil)}. Consider using a MixedCoilSet" + "Coils in a CoilSet must all be the same type, got types " + + f"{type(coil0)} and {type(coil)}. Consider using a MixedCoilSet." ), ) errorif( isinstance(coil0, CoilSet), TypeError, ( - "coils in a CoilSet must all be base Coil types, not CoilSet. " - + "Consider using a MixedCoilSet" + "Coils in a CoilSet must all be base Coil types, not CoilSet. " + + "Consider using a MixedCoilSet." ), ) attrs = { @@ -921,7 +946,6 @@ def _check_type(coil0, coil): FourierPlanarCoil: ["r_basis"], SplineXYZCoil: ["method", "N", "knots"], } - for attr in attrs[coil0.__class__]: a0 = getattr(coil0, attr) a1 = getattr(coil, attr) @@ -929,9 +953,9 @@ def _check_type(coil0, coil): not equals(a0, a1), ValueError, ( - "coils in a CoilSet must have the same parameterization, got a " + "Coils in a CoilSet must have the same parameterization, got a " + f"mismatch between attr {attr}, with values {a0} and {a1}." - + " Consider using a MixedCoilSet" + + " Consider using a MixedCoilSet." ), ) @@ -961,6 +985,7 @@ class CoilSet(OptimizableCollection, _Coil, MutableSequence): _io_attrs_ = _Coil._io_attrs_ + ["_coils", "_NFP", "_sym"] _io_attrs_.remove("_current") + _io_attrs_.remove("_current_scale") def __init__( self, @@ -1022,6 +1047,16 @@ def current(self, new): for coil, cur in zip(self.coils, new): coil.current = cur + @property + def current_scale(self): + """list: current scale of each coil.""" + return self.coils[0].current_scale + + @current_scale.setter + def current_scale(self, new): + for coil in self.coils: + coil.current_scale = new + def _make_arraylike(self, x): if isinstance(x, dict): x = [x] * len(self) @@ -1552,7 +1587,9 @@ def flatten_coils(coilset): contour_Y = np.asarray(coords[0:, 1]) contour_Z = np.asarray(coords[0:, 2]) - currents = np.ones_like(contour_X) * float(coil.current) + currents = np.ones_like(contour_X) * float( + coil.current * coil.current_scale + ) if endpoint: currents[-1] = 0 # this last point must have 0 current else: # close the curves if needed @@ -1875,6 +1912,18 @@ def __init__(self, *coils, name="", check_intersection=True): if check_intersection: self.is_self_intersecting() + @property + def current_scale(self): + """list: current scale of each coil.""" + return [coil.current_scale for coil in self.coils] + + @current_scale.setter + def current_scale(self, new): + if jnp.isscalar(new): + new = [new] * len(self) + for coil, scale in zip(self.coils, new): + coil.current_scale = scale + @property def num_coils(self): """int: Number of coils.""" From e0d8b9b1da19c74b50586cb5ef9aeff4bb44a8a2 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 31 Jul 2024 17:39:07 -0600 Subject: [PATCH 02/26] fixing init methods --- desc/coils.py | 82 +++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 69 insertions(+), 13 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index f482e324f2..9b01c2fc80 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -149,7 +149,7 @@ class _Coil(_MagneticField, Optimizable, ABC): def __init__(self, current, current_scale=1, *args, **kwargs): self._current = float(np.squeeze(current)) - self._current_scale = float(np.squeeze(current)) + self._current_scale = float(np.squeeze(current_scale)) super().__init__(*args, **kwargs) @optimizable_parameter @@ -449,6 +449,9 @@ class FourierRZCoil(_Coil, FourierRZCurve): number of field periods sym : bool whether to enforce stellarator symmetry + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). name : str name for this coil @@ -495,12 +498,25 @@ def __init__( modes_Z=None, NFP=1, sym="auto", + current_scale=1, name="", ): - super().__init__(current, R_n, Z_n, modes_R, modes_Z, NFP, sym, name) + super().__init__( + current, current_scale, R_n, Z_n, modes_R, modes_Z, NFP, sym, name + ) @classmethod - def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name=""): + def from_values( + cls, + current, + coords, + N=10, + NFP=1, + current_scale=1, + basis="rpz", + sym=False, + name="", + ): """Fit coordinates to FourierRZCoil representation. Parameters @@ -515,6 +531,9 @@ def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name= NFP : int Number of field periods, the curve will have a discrete toroidal symmetry according to NFP. + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). basis : {"rpz", "xyz"} basis for input coordinates. Defaults to "rpz" sym : bool @@ -522,7 +541,6 @@ def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name= name : str name for this coil - Returns ------- coil : FourierRZCoil @@ -534,6 +552,7 @@ def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name= ) return FourierRZCoil( current=current, + current_scale=current_scale, R_n=curve.R_n, Z_n=curve.Z_n, modes_R=curve.R_basis.modes[:, 2], @@ -555,6 +574,9 @@ class FourierXYZCoil(_Coil, FourierXYZCurve): fourier coefficients for X, Y, Z modes : array-like mode numbers associated with X_n etc. + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). name : str name for this coil @@ -602,12 +624,15 @@ def __init__( Y_n=[0, 0, 0], Z_n=[-2, 0, 0], modes=None, + current_scale=1, name="", ): - super().__init__(current, X_n, Y_n, Z_n, modes, name) + super().__init__(current, current_scale, X_n, Y_n, Z_n, modes, name) @classmethod - def from_values(cls, current, coords, N=10, s=None, basis="xyz", name=""): + def from_values( + cls, current, coords, N=10, s=None, current_scale=1, basis="xyz", name="" + ): """Fit coordinates to FourierXYZCoil representation. Parameters @@ -624,6 +649,9 @@ def from_values(cls, current, coords, N=10, s=None, basis="xyz", name=""): Should be monotonic, 1D array of same length as coords if None, defaults to normalized arclength + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). basis : {"rpz", "xyz"} basis for input coordinates. Defaults to "xyz" Returns @@ -635,6 +663,7 @@ def from_values(cls, current, coords, N=10, s=None, basis="xyz", name=""): curve = super().from_values(coords=coords, N=N, s=s, basis=basis, name=name) return FourierXYZCoil( current=current, + current_scale=current_scale, X_n=curve.X_n, Y_n=curve.Y_n, Z_n=curve.Z_n, @@ -664,6 +693,9 @@ class FourierPlanarCoil(_Coil, FourierPlanarCurve): mode numbers associated with r_n basis : {'xyz', 'rpz'} Coordinate system for center and normal vectors. Default = 'xyz'. + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). name : str Name for this coil. @@ -712,12 +744,15 @@ def __init__( r_n=2, modes=None, basis="xyz", + current_scale=1, name="", ): - super().__init__(current, center, normal, r_n, modes, basis, name) + super().__init__( + current, current_scale, center, normal, r_n, modes, basis, name + ) @classmethod - def from_values(cls, current, coords, N=10, basis="xyz", name=""): + def from_values(cls, current, coords, N=10, current_scale=1, basis="xyz", name=""): """Fit coordinates to FourierPlanarCoil representation. Parameters @@ -729,6 +764,9 @@ def from_values(cls, current, coords, N=10, basis="xyz", name=""): corresponding to xyz or rpz depending on the basis argument. N : int Fourier resolution of the new r representation. + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). basis : {"rpz", "xyz"} Basis for input coordinates. Defaults to "xyz". name : str @@ -743,6 +781,7 @@ def from_values(cls, current, coords, N=10, basis="xyz", name=""): curve = super().from_values(coords=coords, N=N, basis=basis, name=name) return FourierPlanarCoil( current=current, + current_scale=current_scale, center=curve.center, normal=curve.normal, r_n=curve.r_n, @@ -782,6 +821,9 @@ class SplineXYZCoil(_Coil, SplineXYZCurve): - ``'monotonic-0'``: same as `'monotonic'` but with 0 first derivatives at both endpoints + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). name : str name for this curve @@ -797,9 +839,10 @@ def __init__( Z, knots=None, method="cubic", + current_scale=1, name="", ): - super().__init__(current, X, Y, Z, knots, method, name) + super().__init__(current, current_scale, X, Y, Z, knots, method, name) def compute_magnetic_field( self, coords, params=None, basis="rpz", source_grid=None, transforms=None @@ -855,7 +898,9 @@ def compute_magnetic_field( # coils curvature which is a 2nd derivative of the position, and doing that # with only possibly c1 cubic splines is inaccurate, so we don't do it # (for now, maybe in the future?) - B = biot_savart_hh(coords, coil_pts_start, coil_pts_end, current) + B = biot_savart_hh( + coords, coil_pts_start, coil_pts_end, current * self.current_scale + ) if basis == "rpz": B = xyz2rpz_vec(B, x=coords[:, 0], y=coords[:, 1]) @@ -863,7 +908,14 @@ def compute_magnetic_field( @classmethod def from_values( - cls, current, coords, knots=None, method="cubic", name="", basis="xyz" + cls, + current, + coords, + knots=None, + method="cubic", + current_scale=1, + basis="xyz", + name="", ): """Create SplineXYZCoil from coordinate values. @@ -891,10 +943,13 @@ def from_values( - `'cubic2'`: C2 cubic splines (aka natural splines) - `'catmull-rom'`: C1 cubic centripetal "tension" splines - name : str - name for this curve + current_scale : float + Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. + Default = 1 (current in Amperes with no scaling). basis : {"rpz", "xyz"} basis for input coordinates. Defaults to "xyz" + name : str + name for this curve Returns ------- @@ -907,6 +962,7 @@ def from_values( ) return SplineXYZCoil( current=current, + current_scale=current_scale, X=curve.X, Y=curve.Y, Z=curve.Z, From c674795c92fa4cdc945d173da1d12407773466fc Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 31 Jul 2024 17:39:34 -0600 Subject: [PATCH 03/26] very simple test --- tests/test_coils.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/test_coils.py b/tests/test_coils.py index ffd6558e8b..6d72199d02 100644 --- a/tests/test_coils.py +++ b/tests/test_coils.py @@ -348,6 +348,26 @@ def test_from_symmetry(self): )[0] np.testing.assert_allclose(B_true, B_approx, rtol=1e-3, atol=1e-10) + @pytest.mark.unit + def test_current_scale(self): + """Test different current scales give the same magnetic field.""" + coil = FourierPlanarCoil() + coil.current = 1e7 # Amps + coils = CoilSet.linspaced_angular(coil, n=5) + grid = LinearGrid(N=32, endpoint=False) + transforms = get_transforms(["x", "x_s", "ds"], coil, grid=grid) + B_A = coils.compute_magnetic_field( + [10, 0, 0], basis="rpz", source_grid=grid, transforms=transforms + )[0] + + coils.current = 1 # Mega-Amps + coils.current_scale = 1e7 + B_MA = coils.compute_magnetic_field( + [10, 0, 0], basis="rpz", source_grid=grid, transforms=transforms + )[0] + + np.testing.assert_allclose(B_A, B_MA) + @pytest.mark.unit def test_is_self_intersecting_warnings(self): """Test warning in from_symmetry for self-intersection.""" From 73c3007234c88aa56adac7a50d6bb96227e054cd Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Fri, 9 Aug 2024 15:10:26 -0600 Subject: [PATCH 04/26] revert all commits LOL --- desc/coils.py | 151 +++++++------------------------------------- tests/test_coils.py | 20 ------ 2 files changed, 23 insertions(+), 148 deletions(-) diff --git a/desc/coils.py b/desc/coils.py index 9b01c2fc80..aacb4f19a3 100644 --- a/desc/coils.py +++ b/desc/coils.py @@ -139,17 +139,12 @@ class _Coil(_MagneticField, Optimizable, ABC): ---------- current : float Current through the coil, in Amperes. - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). - """ - _io_attrs_ = _MagneticField._io_attrs_ + ["_current", "_current_scale"] + _io_attrs_ = _MagneticField._io_attrs_ + ["_current"] - def __init__(self, current, current_scale=1, *args, **kwargs): + def __init__(self, current, *args, **kwargs): self._current = float(np.squeeze(current)) - self._current_scale = float(np.squeeze(current_scale)) super().__init__(*args, **kwargs) @optimizable_parameter @@ -163,16 +158,6 @@ def current(self, new): assert jnp.isscalar(new) or new.size == 1 self._current = float(np.squeeze(new)) - @property - def current_scale(self): - """float: Current passing through the coil, in Amperes.""" - return self._current_scale - - @current_scale.setter - def current_scale(self, new): - assert jnp.isscalar(new) or new.size == 1 - self._current_scale = float(np.squeeze(new)) - @property def num_coils(self): """int: Number of coils.""" @@ -272,10 +257,7 @@ def compute_magnetic_field( data["x"] = rpz2xyz(data["x"]) B = biot_savart_quad( - coords, - data["x"], - data["x_s"] * data["ds"][:, None], - current * self.current_scale, + coords, data["x"], data["x_s"] * data["ds"][:, None], current ) if basis.lower() == "rpz": @@ -449,9 +431,6 @@ class FourierRZCoil(_Coil, FourierRZCurve): number of field periods sym : bool whether to enforce stellarator symmetry - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). name : str name for this coil @@ -498,25 +477,12 @@ def __init__( modes_Z=None, NFP=1, sym="auto", - current_scale=1, name="", ): - super().__init__( - current, current_scale, R_n, Z_n, modes_R, modes_Z, NFP, sym, name - ) + super().__init__(current, R_n, Z_n, modes_R, modes_Z, NFP, sym, name) @classmethod - def from_values( - cls, - current, - coords, - N=10, - NFP=1, - current_scale=1, - basis="rpz", - sym=False, - name="", - ): + def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name=""): """Fit coordinates to FourierRZCoil representation. Parameters @@ -531,9 +497,6 @@ def from_values( NFP : int Number of field periods, the curve will have a discrete toroidal symmetry according to NFP. - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). basis : {"rpz", "xyz"} basis for input coordinates. Defaults to "rpz" sym : bool @@ -541,6 +504,7 @@ def from_values( name : str name for this coil + Returns ------- coil : FourierRZCoil @@ -552,7 +516,6 @@ def from_values( ) return FourierRZCoil( current=current, - current_scale=current_scale, R_n=curve.R_n, Z_n=curve.Z_n, modes_R=curve.R_basis.modes[:, 2], @@ -574,9 +537,6 @@ class FourierXYZCoil(_Coil, FourierXYZCurve): fourier coefficients for X, Y, Z modes : array-like mode numbers associated with X_n etc. - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). name : str name for this coil @@ -624,15 +584,12 @@ def __init__( Y_n=[0, 0, 0], Z_n=[-2, 0, 0], modes=None, - current_scale=1, name="", ): - super().__init__(current, current_scale, X_n, Y_n, Z_n, modes, name) + super().__init__(current, X_n, Y_n, Z_n, modes, name) @classmethod - def from_values( - cls, current, coords, N=10, s=None, current_scale=1, basis="xyz", name="" - ): + def from_values(cls, current, coords, N=10, s=None, basis="xyz", name=""): """Fit coordinates to FourierXYZCoil representation. Parameters @@ -649,9 +606,6 @@ def from_values( Should be monotonic, 1D array of same length as coords if None, defaults to normalized arclength - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). basis : {"rpz", "xyz"} basis for input coordinates. Defaults to "xyz" Returns @@ -663,7 +617,6 @@ def from_values( curve = super().from_values(coords=coords, N=N, s=s, basis=basis, name=name) return FourierXYZCoil( current=current, - current_scale=current_scale, X_n=curve.X_n, Y_n=curve.Y_n, Z_n=curve.Z_n, @@ -693,9 +646,6 @@ class FourierPlanarCoil(_Coil, FourierPlanarCurve): mode numbers associated with r_n basis : {'xyz', 'rpz'} Coordinate system for center and normal vectors. Default = 'xyz'. - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). name : str Name for this coil. @@ -744,15 +694,12 @@ def __init__( r_n=2, modes=None, basis="xyz", - current_scale=1, name="", ): - super().__init__( - current, current_scale, center, normal, r_n, modes, basis, name - ) + super().__init__(current, center, normal, r_n, modes, basis, name) @classmethod - def from_values(cls, current, coords, N=10, current_scale=1, basis="xyz", name=""): + def from_values(cls, current, coords, N=10, basis="xyz", name=""): """Fit coordinates to FourierPlanarCoil representation. Parameters @@ -764,9 +711,6 @@ def from_values(cls, current, coords, N=10, current_scale=1, basis="xyz", name=" corresponding to xyz or rpz depending on the basis argument. N : int Fourier resolution of the new r representation. - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). basis : {"rpz", "xyz"} Basis for input coordinates. Defaults to "xyz". name : str @@ -781,7 +725,6 @@ def from_values(cls, current, coords, N=10, current_scale=1, basis="xyz", name=" curve = super().from_values(coords=coords, N=N, basis=basis, name=name) return FourierPlanarCoil( current=current, - current_scale=current_scale, center=curve.center, normal=curve.normal, r_n=curve.r_n, @@ -821,9 +764,6 @@ class SplineXYZCoil(_Coil, SplineXYZCurve): - ``'monotonic-0'``: same as `'monotonic'` but with 0 first derivatives at both endpoints - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). name : str name for this curve @@ -839,10 +779,9 @@ def __init__( Z, knots=None, method="cubic", - current_scale=1, name="", ): - super().__init__(current, current_scale, X, Y, Z, knots, method, name) + super().__init__(current, X, Y, Z, knots, method, name) def compute_magnetic_field( self, coords, params=None, basis="rpz", source_grid=None, transforms=None @@ -898,9 +837,7 @@ def compute_magnetic_field( # coils curvature which is a 2nd derivative of the position, and doing that # with only possibly c1 cubic splines is inaccurate, so we don't do it # (for now, maybe in the future?) - B = biot_savart_hh( - coords, coil_pts_start, coil_pts_end, current * self.current_scale - ) + B = biot_savart_hh(coords, coil_pts_start, coil_pts_end, current) if basis == "rpz": B = xyz2rpz_vec(B, x=coords[:, 0], y=coords[:, 1]) @@ -908,14 +845,7 @@ def compute_magnetic_field( @classmethod def from_values( - cls, - current, - coords, - knots=None, - method="cubic", - current_scale=1, - basis="xyz", - name="", + cls, current, coords, knots=None, method="cubic", name="", basis="xyz" ): """Create SplineXYZCoil from coordinate values. @@ -943,13 +873,10 @@ def from_values( - `'cubic2'`: C2 cubic splines (aka natural splines) - `'catmull-rom'`: C1 cubic centripetal "tension" splines - current_scale : float - Scale of current, e.g. 1e3 for kilo-Amperes, 1e6 for Mega-Amperes, etc. - Default = 1 (current in Amperes with no scaling). - basis : {"rpz", "xyz"} - basis for input coordinates. Defaults to "xyz" name : str name for this curve + basis : {"rpz", "xyz"} + basis for input coordinates. Defaults to "xyz" Returns ------- @@ -962,7 +889,6 @@ def from_values( ) return SplineXYZCoil( current=current, - current_scale=current_scale, X=curve.X, Y=curve.Y, Z=curve.Z, @@ -973,27 +899,20 @@ def from_values( def _check_type(coil0, coil): - errorif( - coil0.current_scale != coil.current_scale, - ValueError, - "Coils in a CoilSet must all have the same current scale, got scales " - + f"{coil0.current_scale} and {coil.current_scale}. " - + "Consider using a MixedCoilSet.", - ) errorif( not isinstance(coil, coil0.__class__), TypeError, ( - "Coils in a CoilSet must all be the same type, got types " - + f"{type(coil0)} and {type(coil)}. Consider using a MixedCoilSet." + "coils in a CoilSet must all be the same type, got types " + + f"{type(coil0)}, {type(coil)}. Consider using a MixedCoilSet" ), ) errorif( isinstance(coil0, CoilSet), TypeError, ( - "Coils in a CoilSet must all be base Coil types, not CoilSet. " - + "Consider using a MixedCoilSet." + "coils in a CoilSet must all be base Coil types, not CoilSet. " + + "Consider using a MixedCoilSet" ), ) attrs = { @@ -1002,6 +921,7 @@ def _check_type(coil0, coil): FourierPlanarCoil: ["r_basis"], SplineXYZCoil: ["method", "N", "knots"], } + for attr in attrs[coil0.__class__]: a0 = getattr(coil0, attr) a1 = getattr(coil, attr) @@ -1009,9 +929,9 @@ def _check_type(coil0, coil): not equals(a0, a1), ValueError, ( - "Coils in a CoilSet must have the same parameterization, got a " + "coils in a CoilSet must have the same parameterization, got a " + f"mismatch between attr {attr}, with values {a0} and {a1}." - + " Consider using a MixedCoilSet." + + " Consider using a MixedCoilSet" ), ) @@ -1041,7 +961,6 @@ class CoilSet(OptimizableCollection, _Coil, MutableSequence): _io_attrs_ = _Coil._io_attrs_ + ["_coils", "_NFP", "_sym"] _io_attrs_.remove("_current") - _io_attrs_.remove("_current_scale") def __init__( self, @@ -1103,16 +1022,6 @@ def current(self, new): for coil, cur in zip(self.coils, new): coil.current = cur - @property - def current_scale(self): - """list: current scale of each coil.""" - return self.coils[0].current_scale - - @current_scale.setter - def current_scale(self, new): - for coil in self.coils: - coil.current_scale = new - def _make_arraylike(self, x): if isinstance(x, dict): x = [x] * len(self) @@ -1643,9 +1552,7 @@ def flatten_coils(coilset): contour_Y = np.asarray(coords[0:, 1]) contour_Z = np.asarray(coords[0:, 2]) - currents = np.ones_like(contour_X) * float( - coil.current * coil.current_scale - ) + currents = np.ones_like(contour_X) * float(coil.current) if endpoint: currents[-1] = 0 # this last point must have 0 current else: # close the curves if needed @@ -1968,18 +1875,6 @@ def __init__(self, *coils, name="", check_intersection=True): if check_intersection: self.is_self_intersecting() - @property - def current_scale(self): - """list: current scale of each coil.""" - return [coil.current_scale for coil in self.coils] - - @current_scale.setter - def current_scale(self, new): - if jnp.isscalar(new): - new = [new] * len(self) - for coil, scale in zip(self.coils, new): - coil.current_scale = scale - @property def num_coils(self): """int: Number of coils.""" diff --git a/tests/test_coils.py b/tests/test_coils.py index 6d72199d02..ffd6558e8b 100644 --- a/tests/test_coils.py +++ b/tests/test_coils.py @@ -348,26 +348,6 @@ def test_from_symmetry(self): )[0] np.testing.assert_allclose(B_true, B_approx, rtol=1e-3, atol=1e-10) - @pytest.mark.unit - def test_current_scale(self): - """Test different current scales give the same magnetic field.""" - coil = FourierPlanarCoil() - coil.current = 1e7 # Amps - coils = CoilSet.linspaced_angular(coil, n=5) - grid = LinearGrid(N=32, endpoint=False) - transforms = get_transforms(["x", "x_s", "ds"], coil, grid=grid) - B_A = coils.compute_magnetic_field( - [10, 0, 0], basis="rpz", source_grid=grid, transforms=transforms - )[0] - - coils.current = 1 # Mega-Amps - coils.current_scale = 1e7 - B_MA = coils.compute_magnetic_field( - [10, 0, 0], basis="rpz", source_grid=grid, transforms=transforms - )[0] - - np.testing.assert_allclose(B_A, B_MA) - @pytest.mark.unit def test_is_self_intersecting_warnings(self): """Test warning in from_symmetry for self-intersection.""" From 6a76eb0d594d865228586439e18812255ff7967f Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Fri, 9 Aug 2024 15:11:34 -0600 Subject: [PATCH 05/26] add scaling matrix to factorize_linear_constraints --- desc/objectives/utils.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 924e8cfb43..323dbd9071 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -141,17 +141,20 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 b = jnp.asarray(b) xp = put(xp, unfixed_idx, Ainv_full @ b) xp = jnp.asarray(xp) + xp_scale = np.where(xp == 0, 1, xp) + D = jnp.diag(xp_scale) + Dinv = jnp.diag(1 / xp_scale) @jit def project(x): """Project a full state vector into the reduced optimization vector.""" - x_reduced = Z.T @ ((x - xp)[unfixed_idx]) + x_reduced = Z.T @ Dinv @ ((x - xp)[unfixed_idx]) return jnp.atleast_1d(jnp.squeeze(x_reduced)) @jit def recover(x_reduced): """Recover the full state vector from the reduced optimization vector.""" - dx = put(jnp.zeros(objective.dim_x), unfixed_idx, Z @ x_reduced) + dx = put(jnp.zeros(objective.dim_x), unfixed_idx, D @ Z @ x_reduced) return jnp.atleast_1d(jnp.squeeze(xp + dx)) # check that all constraints are actually satisfiable From 82ac97f0ef63c6d4d3881d421759ac0a5f5a619e Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Fri, 9 Aug 2024 15:18:38 -0600 Subject: [PATCH 06/26] return diag scaling matrix D --- desc/objectives/utils.py | 4 +++- desc/optimize/_constraint_wrappers.py | 3 ++- desc/perturbations.py | 8 ++++---- desc/vmec.py | 2 +- tests/test_linear_objectives.py | 8 ++++---- 5 files changed, 14 insertions(+), 11 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 323dbd9071..2c1aa75d92 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -33,6 +33,8 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 Combined RHS vector. Z : ndarray Null space operator for full combined A such that A @ Z == 0. + D : ndarray + Diagonal scaling matrix, based on the magnitude of the particular solution xp. unfixed_idx : ndarray Indices of x that correspond to non-fixed values. project, recover : function @@ -200,7 +202,7 @@ def recover(x_reduced): "or be due to floating point error.", ) - return xp, A, b, Z, unfixed_idx, project, recover + return xp, A, b, Z, D, unfixed_idx, project, recover def softmax(arr, alpha): diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index 93bf5385ae..7f46aef4ff 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -103,6 +103,7 @@ def build(self, use_jit=None, verbose=1): self._A, self._b, self._Z, + self._D, self._unfixed_idx, self._project, self._recover, @@ -533,7 +534,7 @@ def _set_eq_state_vector(self): self._args.remove(arg) linear_constraint = ObjectiveFunction(self._linear_constraints) linear_constraint.build() - _, A, _, self._Z, self._unfixed_idx, _, _ = factorize_linear_constraints( + _, A, _, self._Z, _, self._unfixed_idx, _, _ = factorize_linear_constraints( self._constraint, linear_constraint ) diff --git a/desc/perturbations.py b/desc/perturbations.py index 3fa622b248..a9950100e9 100644 --- a/desc/perturbations.py +++ b/desc/perturbations.py @@ -185,7 +185,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces if verbose > 0: print("Factorizing linear constraints") timer.start("linear constraint factorize") - xp, _, _, Z, unfixed_idx, project, recover = factorize_linear_constraints( + xp, _, _, Z, _, unfixed_idx, project, recover = factorize_linear_constraints( objective, constraint ) timer.stop("linear constraint factorize") @@ -388,7 +388,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces con.update_target(eq_new) constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - xp, _, _, Z, unfixed_idx, project, recover = factorize_linear_constraints( + xp, _, _, Z, _, unfixed_idx, project, recover = factorize_linear_constraints( objective, constraint ) @@ -547,7 +547,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - _, _, _, Z, unfixed_idx, project, recover = factorize_linear_constraints( + _, _, _, Z, _, unfixed_idx, project, recover = factorize_linear_constraints( objective_f, constraint ) @@ -752,7 +752,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces con.update_target(eq_new) constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - _, _, _, Z, unfixed_idx, project, recover = factorize_linear_constraints( + _, _, _, Z, _, unfixed_idx, project, recover = factorize_linear_constraints( objective_f, constraint ) diff --git a/desc/vmec.py b/desc/vmec.py index 7a68f52de9..b136027a96 100644 --- a/desc/vmec.py +++ b/desc/vmec.py @@ -192,7 +192,7 @@ def load( constraints = maybe_add_self_consistency(eq, constraints) objective = ObjectiveFunction(constraints) objective.build(verbose=0) - _, _, _, _, _, project, recover = factorize_linear_constraints( + _, _, _, _, _, _, project, recover = factorize_linear_constraints( objective, objective ) args = objective.unpack_state(recover(project(objective.x(eq))), False)[0] diff --git a/tests/test_linear_objectives.py b/tests/test_linear_objectives.py index 167765cd83..7f9b57b217 100644 --- a/tests/test_linear_objectives.py +++ b/tests/test_linear_objectives.py @@ -451,7 +451,7 @@ def test_correct_indexing_passed_modes(): constraint = ObjectiveFunction(constraints, use_jit=False) constraint.build() - xp, A, b, Z, unfixed_idx, project, recover = factorize_linear_constraints( + xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( objective, constraint ) @@ -514,7 +514,7 @@ def test_correct_indexing_passed_modes_and_passed_target(): constraint = ObjectiveFunction(constraints, use_jit=False) constraint.build() - xp, A, b, Z, unfixed_idx, project, recover = factorize_linear_constraints( + xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( objective, constraint ) @@ -574,7 +574,7 @@ def test_correct_indexing_passed_modes_axis(): constraint = ObjectiveFunction(constraints, use_jit=False) constraint.build() - xp, A, b, Z, unfixed_idx, project, recover = factorize_linear_constraints( + xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( objective, constraint ) @@ -703,7 +703,7 @@ def test_correct_indexing_passed_modes_and_passed_target_axis(): constraint = ObjectiveFunction(constraints, use_jit=False) constraint.build() - xp, A, b, Z, unfixed_idx, project, recover = factorize_linear_constraints( + xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( objective, constraint ) From 964c72cd961805d39b77a983ece95a80dc021baa Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Fri, 9 Aug 2024 15:24:39 -0600 Subject: [PATCH 07/26] improved scaling --- desc/objectives/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 2c1aa75d92..72ba584c2b 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -143,7 +143,7 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 b = jnp.asarray(b) xp = put(xp, unfixed_idx, Ainv_full @ b) xp = jnp.asarray(xp) - xp_scale = np.where(xp == 0, 1, xp) + xp_scale = np.where(jnp.abs(xp) < 1, 1, jnp.abs(xp)) D = jnp.diag(xp_scale) Dinv = jnp.diag(1 / xp_scale) From 2a612021e67640d488a85e961e4b2b4ca9f514be Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Mon, 12 Aug 2024 10:48:04 -0600 Subject: [PATCH 08/26] change D to rescale x_reduced instead of x_full --- desc/objectives/utils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 72ba584c2b..1789c05246 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -143,20 +143,20 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 b = jnp.asarray(b) xp = put(xp, unfixed_idx, Ainv_full @ b) xp = jnp.asarray(xp) - xp_scale = np.where(jnp.abs(xp) < 1, 1, jnp.abs(xp)) - D = jnp.diag(xp_scale) - Dinv = jnp.diag(1 / xp_scale) + x_scale = jnp.abs(Z.T @ xp[unfixed_idx]) + x_scale = np.where(x_scale < 1, 1, x_scale) + D = jnp.diag(x_scale) @jit def project(x): """Project a full state vector into the reduced optimization vector.""" - x_reduced = Z.T @ Dinv @ ((x - xp)[unfixed_idx]) + x_reduced = jnp.diag(1 / jnp.diagonal(D)) @ Z.T @ (x - xp)[unfixed_idx] return jnp.atleast_1d(jnp.squeeze(x_reduced)) @jit def recover(x_reduced): """Recover the full state vector from the reduced optimization vector.""" - dx = put(jnp.zeros(objective.dim_x), unfixed_idx, D @ Z @ x_reduced) + dx = put(jnp.zeros(objective.dim_x), unfixed_idx, Z @ D @ x_reduced) return jnp.atleast_1d(jnp.squeeze(xp + dx)) # check that all constraints are actually satisfiable From 5bf24c71c63acbd35d8caf0bd2cd021524491eb0 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Mon, 12 Aug 2024 10:48:28 -0600 Subject: [PATCH 09/26] add D matrix to perturbation logic --- desc/optimize/_constraint_wrappers.py | 30 +++++++++++++++++---------- desc/perturbations.py | 22 ++++++++------------ 2 files changed, 28 insertions(+), 24 deletions(-) diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index 7f46aef4ff..18133acdd9 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -114,9 +114,9 @@ def build(self, use_jit=None, verbose=1): self._dim_x = self._objective.dim_x self._dim_x_reduced = self._Z.shape[1] - # equivalent matrix for A[unfixed_idx]@Z == A@unfixed_idx_mat + # equivalent matrix for A[unfixed_idx] @ Z @ D == A @ unfixed_idx_mat self._unfixed_idx_mat = ( - jnp.eye(self._objective.dim_x)[:, self._unfixed_idx] @ self._Z + jnp.eye(self._objective.dim_x)[:, self._unfixed_idx] @ self._Z @ self._D ) self._built = True @@ -262,7 +262,7 @@ def grad(self, x_reduced, constants=None): """ x = self.recover(x_reduced) df = self._objective.grad(x, constants) - return df[self._unfixed_idx] @ self._Z + return df[self._unfixed_idx] @ self._Z @ self._D def hess(self, x_reduced, constants=None): """Compute Hessian of self.compute_scalar. @@ -282,13 +282,19 @@ def hess(self, x_reduced, constants=None): """ x = self.recover(x_reduced) df = self._objective.hess(x, constants) - return self._Z.T @ df[self._unfixed_idx, :][:, self._unfixed_idx] @ self._Z + return ( + jnp.diag(1 / jnp.diagonal(self._D)) + @ self._Z.T + @ df[self._unfixed_idx, :][:, self._unfixed_idx] + @ self._Z + @ self._D + ) def _jac(self, x_reduced, constants=None, op="scaled"): x = self.recover(x_reduced) if self._objective._deriv_mode == "blocked": fun = getattr(self._objective, "jac_" + op) - return fun(x, constants)[:, self._unfixed_idx] @ self._Z + return fun(x, constants)[:, self._unfixed_idx] @ self._Z @ self._D v = self._unfixed_idx_mat df = getattr(self._objective, "jvp_" + op)(v.T, x, constants) @@ -402,7 +408,7 @@ def jvp_unscaled(self, v, x_reduced, constants=None): def _vjp(self, v, x_reduced, constants=None, op="vjp_scaled"): x = self.recover(x_reduced) df = getattr(self._objective, op)(v, x, constants) - return df[self._unfixed_idx] @ self._Z + return df[self._unfixed_idx] @ self._Z @ self._D def vjp_scaled(self, v, x_reduced, constants=None): """Compute vector-Jacobian product of self.compute_scaled. @@ -534,8 +540,8 @@ def _set_eq_state_vector(self): self._args.remove(arg) linear_constraint = ObjectiveFunction(self._linear_constraints) linear_constraint.build() - _, A, _, self._Z, _, self._unfixed_idx, _, _ = factorize_linear_constraints( - self._constraint, linear_constraint + _, _, _, self._Z, self._D, self._unfixed_idx, _, _ = ( + factorize_linear_constraints(self._constraint, linear_constraint) ) # dx/dc - goes from the full state to optimization variables for eq @@ -619,13 +625,15 @@ def build(self, use_jit=None, verbose=1): # noqa: C901 ) self._dimx_per_thing = [t.dim_x for t in self.things] - # equivalent matrix for A[unfixed_idx]@Z == A@unfixed_idx_mat + # equivalent matrix for A[unfixed_idx] @ Z @ D == A @ unfixed_idx_mat self._unfixed_idx_mat = jnp.eye(self._objective.dim_x) self._unfixed_idx_mat = jnp.split( self._unfixed_idx_mat, np.cumsum([t.dim_x for t in self.things]), axis=-1 ) self._unfixed_idx_mat[self._eq_idx] = ( - self._unfixed_idx_mat[self._eq_idx][:, self._unfixed_idx] @ self._Z + self._unfixed_idx_mat[self._eq_idx][:, self._unfixed_idx] + @ self._Z + @ self._D ) self._unfixed_idx_mat = np.concatenate( [np.atleast_2d(foo) for foo in self._unfixed_idx_mat], axis=-1 @@ -1019,7 +1027,7 @@ def jvp_unscaled(self, v, x, constants=None): @functools.partial(jit, static_argnames=("self", "op")) def _jvp_f(self, xf, dc, constants, op): Fx = getattr(self._constraint, "jac_" + op)(xf, constants) - Fx_reduced = Fx[:, self._unfixed_idx] @ self._Z + Fx_reduced = Fx[:, self._unfixed_idx] @ self._Z @ self._D Fc = Fx @ (self._dxdc @ dc) Fxh = Fx_reduced cutoff = jnp.finfo(Fxh.dtype).eps * max(Fxh.shape) diff --git a/desc/perturbations.py b/desc/perturbations.py index a9950100e9..8f86cd826c 100644 --- a/desc/perturbations.py +++ b/desc/perturbations.py @@ -185,7 +185,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces if verbose > 0: print("Factorizing linear constraints") timer.start("linear constraint factorize") - xp, _, _, Z, _, unfixed_idx, project, recover = factorize_linear_constraints( + xp, _, _, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( objective, constraint ) timer.stop("linear constraint factorize") @@ -282,7 +282,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces weight = jnp.diag(weight) else: weight = weight[unfixed_idx, unfixed_idx] - W = Z.T @ weight @ Z + W = jnp.diag(1 / jnp.diagonal(D)) @ Z.T @ weight @ Z @ D scale_inv = W scale = jnp.linalg.inv(scale_inv) @@ -291,7 +291,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces print("Computing df") timer.start("df computation") Jx = objective.jac_scaled_error(x) - Jx_reduced = Jx[:, unfixed_idx] @ Z @ scale + Jx_reduced = Jx[:, unfixed_idx] @ Z @ D @ scale RHS1 = objective.jvp_scaled(tangents, x) if include_f: f = objective.compute_scaled_error(x) @@ -388,9 +388,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces con.update_target(eq_new) constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - xp, _, _, Z, _, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint - ) + _, _, _, _, _, _, _, recover = factorize_linear_constraints(objective, constraint) # update other attributes dx_reduced = dx1_reduced + dx2_reduced + dx3_reduced @@ -547,7 +545,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - _, _, _, Z, _, unfixed_idx, project, recover = factorize_linear_constraints( + _, _, _, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( objective_f, constraint ) @@ -564,7 +562,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces dx2_reduced = 0 # dx/dx_reduced - dxdx_reduced = jnp.eye(eq.dim_x)[:, unfixed_idx] @ Z + dxdx_reduced = jnp.eye(eq.dim_x)[:, unfixed_idx] @ Z @ D # dx/dc dxdc = [] @@ -612,8 +610,8 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces timer.disp("dg computation") # projections onto optimization space - Fx_reduced = Fx[:, unfixed_idx] @ Z - Gx_reduced = Gx[:, unfixed_idx] @ Z + Fx_reduced = Fx[:, unfixed_idx] @ Z @ D + Gx_reduced = Gx[:, unfixed_idx] @ Z @ D Fc = Fx @ dxdc Gc = Gx @ dxdc @@ -752,9 +750,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces con.update_target(eq_new) constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - _, _, _, Z, _, unfixed_idx, project, recover = factorize_linear_constraints( - objective_f, constraint - ) + _, _, _, _, _, _, _, recover = factorize_linear_constraints(objective_f, constraint) # update other attributes dx_reduced = dx1_reduced + dx2_reduced From 869f408de20f77d34b3f22c77aeffbda4cf36e3c Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Mon, 12 Aug 2024 22:14:27 -0600 Subject: [PATCH 10/26] improved scaling --- desc/objectives/utils.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 05a9a438af..356b3115b7 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -143,9 +143,8 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 b = jnp.asarray(b) xp = put(xp, unfixed_idx, Ainv_full @ b) xp = jnp.asarray(xp) - x_scale = jnp.abs(Z.T @ xp[unfixed_idx]) - x_scale = np.where(x_scale < 1, 1, x_scale) - D = jnp.diag(x_scale) + x_scale = jnp.abs(Z.T) @ jnp.abs(xp[unfixed_idx]) / jnp.sum(jnp.abs(Z.T), axis=1) + D = jnp.diag(np.where(x_scale < 1, 1, x_scale)) @jit def project(x): From 240752c4c4ffbf28ef601935a008e9a59f639245 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Mon, 12 Aug 2024 22:15:01 -0600 Subject: [PATCH 11/26] added test --- tests/test_optimizer.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/tests/test_optimizer.py b/tests/test_optimizer.py index 7c5f810d98..e4c06766bb 100644 --- a/tests/test_optimizer.py +++ b/tests/test_optimizer.py @@ -5,6 +5,7 @@ import numpy as np import pytest from numpy.random import default_rng +from scipy.constants import mu_0 from scipy.optimize import ( BFGS, NonlinearConstraint, @@ -32,6 +33,7 @@ FixParameters, FixPressure, FixPsi, + FixSumCoilCurrent, ForceBalance, GenericObjective, MagneticWell, @@ -1350,3 +1352,28 @@ def test_quad_flux_with_surface_current_field(): (field_modular_opt,), result = opt.optimize( field, objective=obj, constraints=constraints, maxiter=1, copy=True ) + + +@pytest.mark.unit +def test_optimize_coil_currents(DummyCoilSet): + """Tests optimization takes step sizes proportional to variable scales.""" + eq = desc.examples.get("precise_QH") + coils = load(load_from=str(DummyCoilSet["output_path_sym"]), file_format="hdf5") + grid = LinearGrid(rho=1.0, M=eq.M_grid, N=eq.N_grid, NFP=eq.NFP, sym=eq.sym) + current = 2 * np.pi * eq.compute("G", grid=grid)["G"][0] / mu_0 + for coil in coils: + coil.current = current / coils.num_coils + + objective = ObjectiveFunction(QuadraticFlux(eq=eq, field=coils, vacuum=True)) + constraints = FixSumCoilCurrent(coils) + optimizer = Optimizer("lsq-exact") + [coils_opt], _ = optimizer.optimize( + things=coils, + objective=objective, + constraints=constraints, + verbose=2, + ) + # check that optimized coil currents changed sufficiently from initial values + np.testing.assert_array_less( + 1e5, np.abs(np.asarray(coils_opt.current) - np.asarray(coils.current)) + ) From 127d453c84fa6026516d26094d12efab3b5c96de Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Tue, 13 Aug 2024 17:20:24 -0600 Subject: [PATCH 12/26] looks like it's working! --- desc/objectives/utils.py | 41 ++++++++++++++++++--------- desc/optimize/_constraint_wrappers.py | 36 +++++++++++++++-------- desc/perturbations.py | 10 +++---- tests/test_optimizer.py | 6 ++-- 4 files changed, 61 insertions(+), 32 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 356b3115b7..a50c6ff915 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -132,34 +132,49 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 ) A = A[unfixed_rows][:, unfixed_idx] b = b[unfixed_rows] + + # unscaled particular solution to get scale of x unfixed_idx = indices_idx if A.size: - Ainv_full, Z = svd_inv_null(A) + A_inv, Z = svd_inv_null(A) else: - Ainv_full = A.T + A_inv = A.T Z = np.eye(A.shape[1]) - Ainv_full = jnp.asarray(Ainv_full) - Z = jnp.asarray(Z) - b = jnp.asarray(b) - xp = put(xp, unfixed_idx, Ainv_full @ b) + xp = put(xp, unfixed_idx, A_inv @ b) + x_scale = np.where(np.abs(xp) < 1, 1, np.abs(xp)) # TODO: adjust threshold? + D = np.diag(x_scale) + + # scaled system + A = A @ D[unfixed_idx][:, unfixed_idx] + if A.size: + A_inv, Z = svd_inv_null(A) + else: + A_inv = A.T + Z = np.eye(A.shape[1]) + xp = put(xp, unfixed_idx, A_inv @ b) + + # cast to jnp arrays xp = jnp.asarray(xp) - x_scale = jnp.abs(Z.T) @ jnp.abs(xp[unfixed_idx]) / jnp.sum(jnp.abs(Z.T), axis=1) - D = jnp.diag(np.where(x_scale < 1, 1, x_scale)) + A = jnp.asarray(A) + b = jnp.asarray(b) + Z = jnp.asarray(Z) + D = jnp.asarray(D) @jit - def project(x): + def project(x_full): """Project a full state vector into the reduced optimization vector.""" - x_reduced = jnp.diag(1 / jnp.diagonal(D)) @ Z.T @ (x - xp)[unfixed_idx] + x_reduced = Z.T @ (jnp.diag(1 / jnp.diagonal(D)) @ x_full - xp)[unfixed_idx] return jnp.atleast_1d(jnp.squeeze(x_reduced)) @jit def recover(x_reduced): """Recover the full state vector from the reduced optimization vector.""" - dx = put(jnp.zeros(objective.dim_x), unfixed_idx, Z @ D @ x_reduced) - return jnp.atleast_1d(jnp.squeeze(xp + dx)) + dx = put(jnp.zeros(objective.dim_x), unfixed_idx, Z @ x_reduced) + x_full = D @ (xp + dx) + return jnp.atleast_1d(jnp.squeeze(x_full)) # check that all constraints are actually satisfiable - params = objective.unpack_state(xp, False) + params = objective.unpack_state(D @ xp, False) for con in constraint.objectives: xpi = [params[i] for i, t in enumerate(objective.things) if t in con.things] y1 = con.compute_unscaled(*xpi) diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index 18133acdd9..9f220a9ea5 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -114,9 +114,11 @@ def build(self, use_jit=None, verbose=1): self._dim_x = self._objective.dim_x self._dim_x_reduced = self._Z.shape[1] - # equivalent matrix for A[unfixed_idx] @ Z @ D == A @ unfixed_idx_mat + # equivalent matrix for A[unfixed_idx] @ D @ Z == A @ unfixed_idx_mat self._unfixed_idx_mat = ( - jnp.eye(self._objective.dim_x)[:, self._unfixed_idx] @ self._Z @ self._D + jnp.eye(self._objective.dim_x)[:, self._unfixed_idx] + @ self._D[self._unfixed_idx][:, self._unfixed_idx] + @ self._Z ) self._built = True @@ -262,7 +264,7 @@ def grad(self, x_reduced, constants=None): """ x = self.recover(x_reduced) df = self._objective.grad(x, constants) - return df[self._unfixed_idx] @ self._Z @ self._D + return df @ self._unfixed_idx_mat def hess(self, x_reduced, constants=None): """Compute Hessian of self.compute_scalar. @@ -283,18 +285,20 @@ def hess(self, x_reduced, constants=None): x = self.recover(x_reduced) df = self._objective.hess(x, constants) return ( - jnp.diag(1 / jnp.diagonal(self._D)) - @ self._Z.T + self._Z.T + @ jnp.diag( + 1 / jnp.diagonal(self._D[self._unfixed_idx][:, self._unfixed_idx]) + ) @ df[self._unfixed_idx, :][:, self._unfixed_idx] - @ self._Z - @ self._D + @ self._D[self._unfixed_idx][:, self._unfixed_idx] + @ self._Z # TODO: replace with self._unfixed_idx_mat ) def _jac(self, x_reduced, constants=None, op="scaled"): x = self.recover(x_reduced) if self._objective._deriv_mode == "blocked": fun = getattr(self._objective, "jac_" + op) - return fun(x, constants)[:, self._unfixed_idx] @ self._Z @ self._D + return fun(x, constants) @ self._unfixed_idx_mat v = self._unfixed_idx_mat df = getattr(self._objective, "jvp_" + op)(v.T, x, constants) @@ -408,7 +412,11 @@ def jvp_unscaled(self, v, x_reduced, constants=None): def _vjp(self, v, x_reduced, constants=None, op="vjp_scaled"): x = self.recover(x_reduced) df = getattr(self._objective, op)(v, x, constants) - return df[self._unfixed_idx] @ self._Z @ self._D + return ( + df[self._unfixed_idx] + @ self._D[self._unfixed_idx][:, self._unfixed_idx] + @ self._Z + ) def vjp_scaled(self, v, x_reduced, constants=None): """Compute vector-Jacobian product of self.compute_scaled. @@ -625,15 +633,15 @@ def build(self, use_jit=None, verbose=1): # noqa: C901 ) self._dimx_per_thing = [t.dim_x for t in self.things] - # equivalent matrix for A[unfixed_idx] @ Z @ D == A @ unfixed_idx_mat + # equivalent matrix for A[unfixed_idx] @ D @ Z == A @ unfixed_idx_mat self._unfixed_idx_mat = jnp.eye(self._objective.dim_x) self._unfixed_idx_mat = jnp.split( self._unfixed_idx_mat, np.cumsum([t.dim_x for t in self.things]), axis=-1 ) self._unfixed_idx_mat[self._eq_idx] = ( self._unfixed_idx_mat[self._eq_idx][:, self._unfixed_idx] + @ self._D[self._unfixed_idx][:, self._unfixed_idx] @ self._Z - @ self._D ) self._unfixed_idx_mat = np.concatenate( [np.atleast_2d(foo) for foo in self._unfixed_idx_mat], axis=-1 @@ -1027,7 +1035,11 @@ def jvp_unscaled(self, v, x, constants=None): @functools.partial(jit, static_argnames=("self", "op")) def _jvp_f(self, xf, dc, constants, op): Fx = getattr(self._constraint, "jac_" + op)(xf, constants) - Fx_reduced = Fx[:, self._unfixed_idx] @ self._Z @ self._D + Fx_reduced = ( + Fx[:, self._unfixed_idx] + @ self._D[self._unfixed_idx][:, self._unfixed_idx] + @ self._Z # TODO: replace with self._unfixed_idx_mat + ) Fc = Fx @ (self._dxdc @ dc) Fxh = Fx_reduced cutoff = jnp.finfo(Fxh.dtype).eps * max(Fxh.shape) diff --git a/desc/perturbations.py b/desc/perturbations.py index 8f86cd826c..bfde718233 100644 --- a/desc/perturbations.py +++ b/desc/perturbations.py @@ -282,7 +282,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces weight = jnp.diag(weight) else: weight = weight[unfixed_idx, unfixed_idx] - W = jnp.diag(1 / jnp.diagonal(D)) @ Z.T @ weight @ Z @ D + W = Z.T @ weight @ Z scale_inv = W scale = jnp.linalg.inv(scale_inv) @@ -291,7 +291,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces print("Computing df") timer.start("df computation") Jx = objective.jac_scaled_error(x) - Jx_reduced = Jx[:, unfixed_idx] @ Z @ D @ scale + Jx_reduced = Jx[:, unfixed_idx] @ D @ Z @ scale RHS1 = objective.jvp_scaled(tangents, x) if include_f: f = objective.compute_scaled_error(x) @@ -562,7 +562,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces dx2_reduced = 0 # dx/dx_reduced - dxdx_reduced = jnp.eye(eq.dim_x)[:, unfixed_idx] @ Z @ D + dxdx_reduced = jnp.eye(eq.dim_x)[:, unfixed_idx] @ D @ Z # dx/dc dxdc = [] @@ -610,8 +610,8 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces timer.disp("dg computation") # projections onto optimization space - Fx_reduced = Fx[:, unfixed_idx] @ Z @ D - Gx_reduced = Gx[:, unfixed_idx] @ Z @ D + Fx_reduced = Fx[:, unfixed_idx] @ D @ Z + Gx_reduced = Gx[:, unfixed_idx] @ D @ Z Fc = Fx @ dxdc Gc = Gx @ dxdc diff --git a/tests/test_optimizer.py b/tests/test_optimizer.py index e4c06766bb..71cc515a58 100644 --- a/tests/test_optimizer.py +++ b/tests/test_optimizer.py @@ -1372,8 +1372,10 @@ def test_optimize_coil_currents(DummyCoilSet): objective=objective, constraints=constraints, verbose=2, + copy=True, ) - # check that optimized coil currents changed sufficiently from initial values + # check that optimized coil currents changed by more than 15% from initial values np.testing.assert_array_less( - 1e5, np.abs(np.asarray(coils_opt.current) - np.asarray(coils.current)) + np.asarray(coils.current) * 0.15, + np.abs(np.asarray(coils_opt.current) - np.asarray(coils.current)), ) From a30fc3f923f0364845514340d58cc6425ae837fe Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Tue, 13 Aug 2024 17:42:05 -0600 Subject: [PATCH 13/26] fixed indexing --- desc/optimize/_constraint_wrappers.py | 24 ++++++------------------ desc/perturbations.py | 8 ++++---- 2 files changed, 10 insertions(+), 22 deletions(-) diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index 9f220a9ea5..cf4c036698 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -116,9 +116,7 @@ def build(self, use_jit=None, verbose=1): # equivalent matrix for A[unfixed_idx] @ D @ Z == A @ unfixed_idx_mat self._unfixed_idx_mat = ( - jnp.eye(self._objective.dim_x)[:, self._unfixed_idx] - @ self._D[self._unfixed_idx][:, self._unfixed_idx] - @ self._Z + jnp.eye(self._objective.dim_x) @ self._D[:, self._unfixed_idx] @ self._Z ) self._built = True @@ -286,11 +284,9 @@ def hess(self, x_reduced, constants=None): df = self._objective.hess(x, constants) return ( self._Z.T - @ jnp.diag( - 1 / jnp.diagonal(self._D[self._unfixed_idx][:, self._unfixed_idx]) - ) - @ df[self._unfixed_idx, :][:, self._unfixed_idx] - @ self._D[self._unfixed_idx][:, self._unfixed_idx] + @ jnp.diag(1 / jnp.diagonal(self._D[self._unfixed_idx, :])) + @ df + @ self._D[:, self._unfixed_idx] @ self._Z # TODO: replace with self._unfixed_idx_mat ) @@ -412,11 +408,7 @@ def jvp_unscaled(self, v, x_reduced, constants=None): def _vjp(self, v, x_reduced, constants=None, op="vjp_scaled"): x = self.recover(x_reduced) df = getattr(self._objective, op)(v, x, constants) - return ( - df[self._unfixed_idx] - @ self._D[self._unfixed_idx][:, self._unfixed_idx] - @ self._Z - ) + return df @ self._unfixed_idx_mat def vjp_scaled(self, v, x_reduced, constants=None): """Compute vector-Jacobian product of self.compute_scaled. @@ -1035,11 +1027,7 @@ def jvp_unscaled(self, v, x, constants=None): @functools.partial(jit, static_argnames=("self", "op")) def _jvp_f(self, xf, dc, constants, op): Fx = getattr(self._constraint, "jac_" + op)(xf, constants) - Fx_reduced = ( - Fx[:, self._unfixed_idx] - @ self._D[self._unfixed_idx][:, self._unfixed_idx] - @ self._Z # TODO: replace with self._unfixed_idx_mat - ) + Fx_reduced = Fx @ self._D[:, self._unfixed_idx] @ self._Z Fc = Fx @ (self._dxdc @ dc) Fxh = Fx_reduced cutoff = jnp.finfo(Fxh.dtype).eps * max(Fxh.shape) diff --git a/desc/perturbations.py b/desc/perturbations.py index bfde718233..5ce3e28fc4 100644 --- a/desc/perturbations.py +++ b/desc/perturbations.py @@ -291,7 +291,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces print("Computing df") timer.start("df computation") Jx = objective.jac_scaled_error(x) - Jx_reduced = Jx[:, unfixed_idx] @ D @ Z @ scale + Jx_reduced = Jx @ D[:, unfixed_idx] @ Z @ scale RHS1 = objective.jvp_scaled(tangents, x) if include_f: f = objective.compute_scaled_error(x) @@ -562,7 +562,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces dx2_reduced = 0 # dx/dx_reduced - dxdx_reduced = jnp.eye(eq.dim_x)[:, unfixed_idx] @ D @ Z + dxdx_reduced = jnp.eye(eq.dim_x) @ D[:, unfixed_idx] @ Z # dx/dc dxdc = [] @@ -610,8 +610,8 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces timer.disp("dg computation") # projections onto optimization space - Fx_reduced = Fx[:, unfixed_idx] @ D @ Z - Gx_reduced = Gx[:, unfixed_idx] @ D @ Z + Fx_reduced = Fx @ D[:, unfixed_idx] @ Z + Gx_reduced = Gx @ D[:, unfixed_idx] @ Z Fc = Fx @ dxdc Gc = Gx @ dxdc From 0ac5456797b2bfff9fb0e82737d2c87125af1f06 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 14 Aug 2024 11:08:39 -0600 Subject: [PATCH 14/26] also rescale fixed values of xp --- desc/objectives/utils.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index a50c6ff915..a5a0dafbf3 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -133,8 +133,10 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 A = A[unfixed_rows][:, unfixed_idx] b = b[unfixed_rows] - # unscaled particular solution to get scale of x unfixed_idx = indices_idx + fixed_idx = np.delete(np.arange(xp.size), unfixed_idx) + + # unscaled particular solution to get scale of x if A.size: A_inv, Z = svd_inv_null(A) else: @@ -152,6 +154,7 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 A_inv = A.T Z = np.eye(A.shape[1]) xp = put(xp, unfixed_idx, A_inv @ b) + xp = put(xp, fixed_idx, (jnp.diag(1 / jnp.diagonal(D)) @ xp)[fixed_idx]) # cast to jnp arrays xp = jnp.asarray(xp) From bc6fe03861916f006a9a590711e52d6ba6902716 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 14 Aug 2024 15:42:15 -0600 Subject: [PATCH 15/26] making improvements --- desc/examples/HELIOTRON_output.h5 | Bin 338380 -> 324656 bytes desc/objectives/utils.py | 15 +++++++-------- desc/optimize/_constraint_wrappers.py | 12 +++++------- desc/perturbations.py | 8 ++++---- tests/test_linear_objectives.py | 16 ++++++++-------- 5 files changed, 24 insertions(+), 27 deletions(-) diff --git a/desc/examples/HELIOTRON_output.h5 b/desc/examples/HELIOTRON_output.h5 index aee9e5bab248ddacb4e54638555ecef8a73b1c34..a1cb51c45686a70cc32a8c68d7c9169f261c1c11 100644 GIT binary patch delta 14232 zcmbVS2UrwWw4OV6z+LK6gjEDqMJX1*La^;(jV%#H6h#e26R~3hVgq9d8Vk6}n4+%6 z80-*~!M1||ji@M=C`OZDjD}buu_8s^xp!tZD|z4h@Xhz_?0?%i=bm%!xoyr}cF(SG zOAo*yPq~~nJyd>PjP!)G5AD!Cv`uKIF6}3UcMR#$`QwknyR@G;dG_2-7~?Y4n1)rX zRZRMLOM6yS&JIn(Oy(F z<(dXcGuUNXgGmphI;?bQAn7r#3maznr>=t=R&k75!-XA1gkA3*)8O8XHRd>Y8q?iF z8Yat4_lXUwIHqyKg^o8ick*2C(LfOfdo-N#Orwjlr?Ju_05SZD*~k-l<+^BtRW2k} zX;e}Kw1GcljL`+C-J}f0!8oo}lN968z2CGe+4g(n7w5K%m3M3E2K~BVy4%2jvaR2n zPR{#t^~xHr#O2c*@0Z*Z_f*aJE%9HCf!ir;T^sd-hC@m$>KbMx1vdpGshp4O&JTE6J{5b=7)-){6M-EwQms2O_$ z4)jRD4({4h$ESMSoLM^YNnV>TPrf>OzfrpW==;ZytiN~o!h@Iz_Pg&S8nPz@RO~yV zdpgq>7gj%dn7ieAyPL;;B5=iOigm$sZ(-J9;Q@nOIY@pvK#=_!Ce6_YS5+UVo8Nu_9{2qNVRH|Kr@Li8Z;w z4v$x-)aBhy`}5eDE1kxtd~{`RP1mXCuZJGgPXBs-x7h|73)>ePvb)8l)La^uUcJbnY1Hes2@`KzC#gUFn{XRU5gwXe4Wlch&}4eBvuK0Bescj2EJz8So7kQm*~@aMMPe>2*hV)(>W_FG>i%PwdE!C0 zc9ks-2>Mr>%O-^H>?}IX?WkUyJV7jFbB~07J5_&q`M`kwGhUbrF8FY5?tRcroOu6U z=@&g#>U*cPN{OClYv}t$xv!~>uej&$-(Ho?&N7=GWi|bI>NEYVh|SG^{p+FmT~nLi zH;j9zKRu{u!v6D@%-fwOAKu{glXvi=`FWh3zpkw_liU}tjb3p&zl!sV z&|1|dExn8${!@VR693313T2&DogiB0AOAM)9#*;S03&#EeFGM-Z37r0iQ%#!=s_;vlfi7B#1DyX{j)9 zw4lZo7a?4}O5NefscZ35{u=4BsHzE zpp6!^8Cua^=+>1X_FKQ)(Md6O7`0M(xU`aa1Z-}}p%Hw78Rr>S*?JhY+XL+$Nha1l zM!dKKPr}WivWuE4%uW6VLSaL#p%t}3 z+3#THVYxYf0E>B0EU{Rwuqdw_jrdpJ;t;&i(snZx0_CH|>_}a1%+pIO{uXk!dD|T2 z3zo5}+{zaq$RP`@lTi3i9CwFr8qSV4Uas-v&(Gx+8xv}^>`$%v6R)vHuBo=EiWc*3 zHb^Q*ZZK$!X|)=m%m$5C;U^N)#>VnPDgj3D$E=XALBP8P3P}HjNOuH$?fn%6cma|W z(9(k1SWriQNdbK6ckBXTju!$=cM~o)L!gI#P)pLo2&hFyR9+yOZG%6C6lE12OjspC zHOZtaEM(VOP^JZKfo32`3EzG&2|Npq_xLsUN zR$_XTk4Qi|A)p6p%HazND^)Kd|8Eu<9zg+^=(3R77wHiD+F#h-4^`S>`vAGaF~d+e zGV1VX)-;?i-Hs0;FMeATS`L5y!x79A!qnl&8#!99W_t3Qqfr|h8r=E)qmfq%%$E7+ z-B(xLLsw+_`bE4tal^5E)+2FZ{4AgNqbJ2puY0U|>UG3Wm3Vi+%8Q%D$RD!4d-v4_ zF~84^X2XUS8K&)c-u4CmgIIpUC2i7K7?H`&*brd#;+9)a>Ajb&yUq z-V;CkUC*S39~9RVLisXw@vRc`)jQ`}gB9Z;CWLjLZq`hxG&Bh}t~7+^i!E@|GUq!#@Rrb~gvd6_%J?i&uo%@v>O2le$y+IL?szO!qCVaDNghEwAP z{cZopHgV?Sz~%E={B35!pT>OR;A&_ydw-J~>Pe!98O9BWXeb7(8-AuY%st-)-)Z?o^q72bp4o7Z^ z7-c9taC&6NVaLothrVijGqQ;ov+m^DFWR&;Fc-bcvwX*kv5T$;bsDEK9QyhGBcESf z#gWhUZJik(Azr<=+GpGrA4A-WkM2CdUWOFL)V0~*o4x2*La&WtolTId`j1k442YuHnMKh#ZANN7JZ(n7Dt6HbxS$X zeSwx4~82wk&9E;PT6S)>Gs5X=H$wNFl|)#8%OPn@px86$o;a>(^L zp6kU2>m6-Px1O167Y_9K;LCxAGu{7elD9KKOwpbFsQSZ4=A|b(OwL(lvVyA0CSr@i+ z7GJ(PwB84tNZFnB0qn(ZAC2ztRgWM{BHlCxMWNwB$WH8p*n9;PDB$#PK5iGjnKvF? zM|?&y@)RmRM7`~x1*~YGg>5vSD&ituH4_D*D1K-F_7NgyqT3ir1wz#vWJD-hs*b~= zkrbp8l2CJ!y*?TJV+VC19v2DymqA@N8c8z}J3m?y3x5h(|FpINb^uTS0BAwP5-L}L zTtMc`YMApBE^Mj7dy^G!wo(y-8G*=5UnnGOCgmgKuILmHGJ;MIuwE9*w;?yoj!-0u z8zCzyUI7Uscu?IjKD!6j@l{!^Y8%5z|$fdyoqohBIQl`H(&6D_8|WFOLW- z_M*+el?Mf0eXFUGH#=f)Re%%|tX6hp8^CT-q|c(1Z_-UbK1fbC0c?RJr@RMXPb>15 zO8I05kSU-hYRl48Xd$2kgv&rzWN|1Y9451B`3PbC5u{hM-^2))3JLHSEx4ATa0hl} ztYk#Z&tzWyB39Aau2@-TISM$efD^Hj74yqc8B!O=O3Pw6TPR1Kxglq}+Oado$i_(- zL;XxB_?t}PuxMjVPlTpbD2%mFn2f1Pvf8ttYWwofr6jW!^Xs;7@o0=BS9T4^#a2C| zwtZugA3$dIR(!8YR03-O%|GZB{5d~HzWn;fNXx(f7dr1*QtFEvR<%c{&-BOfX*l;qo@Ke&HE@zb{*g9=f zFHPE7mp7`HCcAR$ay6u(O;Za)yLuJu`4`WT?*K&}2Ph_bQPPfpiHu6xsb1O{($Z6D zkd|v(U(n^+!z(%J0v*Fed{%7h6c?7H@|?x>oIdz~=I~S>^I|=h+U?~29W$M4x<5Mo zCEs-Mj90tHV&=zA1OIB&EP33+E1`~z&Ve~GzZY}KWr>Sl`rHQtwO`&`6I|YsshGIq z_>|Ag=l9wyt{lQ7C!xe1w)xQJ#1bM%q3LCmZ|hqSvLI3Lj%$Agkw~n_uK`TX6Uc*% z} zPboHrOi>Bh8a4vkr|y${!*l=}@B8}LU67Iy7cnp1fid%38Z(iNSs2S;7bzm6MndpT z?2M4J{O>9M`|=;`!oLVrUPL_q7yl6tof2l6NyDpriil6<1+kYfn_EkVq8uqTIOHM6 z5Ao-`+4M!mo}6z+!4R@o$9G>Kx;YjiLqHhxT%yFaAlW@vH(U*&4GR zRkHehWPn>&$+1rr&O+>KUnK-}VOQAUK2=CKYlG?<3E6|#7}kE%BT0Mt|6t=+CLcRB z@q31evFG>$&28dQm88u0Aw-!ULSwq&VrMH1=PL|9g)%zNK@&Q?X-0rT+4WVzy^(B; z9a~T(8@aSfiiX;WY&up~`AH15yvzHnuYLYdDOnc9(c=ZSToynznFvs65rBCq_gJL> zjRLeUB>}Q$5&?8Cg;leN0Pnq&1sMBM5+HvOEAOq@=}OyqFL^R03YsNEIWbcuX_BXN zBs#~XNy*ezl5lv@OF9`J6=2)_C1Tge9Pnnyx9qt2J!$N}TW-AXCCdH0Hdx+n-_MZC z=yn_Dzogcgq+-ufK|GY*_7={*>I+#@UsGW=b0*fb7Gh z-}-8av*sA#{C>5}8IL3A9$9Gt$hZ#xDV!6liJs3}?a0~13fHwV*R)e$pSN0M&@s~N zzl#hyMz&Cv$At1?_9BZq!!Y&XOV2@o=|rI}0;w`f*r6zxCKF1aB$UDeC}FiB6lcF+ zW05Rh`FZFrnKB$q6GqXlAPWclyh+|XaRB06vYMopOZ>^Dq>R(HQI#44_ zsbhmxsGYF?1h!>SKEJJmg>V%?&lZYms3B22Kg$RE@)2sRM$JgSjg=i#Y??P!K|0Cx z0SL{u(q>m^vx9s>Al5ecwuvP-dk7nI>^%`@JmKt#Z0OsTz8lfEioVt4&1VN;oiRqE zK@P@ks&`;!lO8h7lIQ=rye8#N+l(+j)^yYL2!=_+0J+Ubnp>VVM3gcLjt_;oydJA1RLbx9~$ zwtC9lK)FZIaoCf-d(n4q`tC#Ded)U&bA#CvxJ`jr9(g2k@Ma7!h zNKq+?Mx`jF-2@rnIknQ>nN95|MY7FEB_y=KQxLmUscclhH_{4E-blM*Tv)3_?#s1u z}x7My(bR9n#Z(AvqQy~_r%-Tnp*y?x9W5j z_5SZHO1J9vb+R?Q>Ih|NtKOzg*{XM}ql!?r>fP#OP4$-Y$yS|@8-&$?kL87|tOKi! zFU?7i^D9S?{=r%m_da^)?*n zQcw7V$-=2^AQAgp9i2bfh-AsX<{~xvguu)ru$q`p5YS)~764LaPaKhe$)!dv-PS8R zEEj&raqK0y=iy5>J%~_Bfh%C^tuNUwT5rwVA#J_&Y|O(}(!Yz+oZRJFH}|a9ob2kX zn>V4&EgRdW^_r8a?ct}8Z&t6lGFK_(KK$8ZctR5;NpG#Kn(Iih7a#W;tBl$HoJH#9 zP{&$oa(`{z>6Z1_`cXQ@#Gv8qBiRFSWCE6Nqhc&54iH?zjX-R~V&U*^9E;fa1m0sW zn0)*kLBAJYRnu);*F-E|77a|q^gfC$n?|f@jN@Po&ZH1FsBG#&3Sl_1X)7rtaVa?s zAlN&Joc*|_t$WkS!uLvRF*byksu_5b7mj6t6S?bR%QGk;QMsB;KUi zC^g&4r!GgaC{aia;%pI~4>C))Z=V8^P9OnDCypMofxNjl zHy(N`Yk0S!;pDZnvqZzC=@b$T>oO?>4Wj`5M0?iPNHnb7E$M)0xVcNtn+h%~kriB; zFV!+oqaINdoSI8X5(TFpq%bvFjSdk+^qZY0wc5&?T7j+<{r+sx@9!3LPeJTsS;MJ^ zWew*SN*dmLifGvA96cef0FHsQ3`!n1q5fGL%|eYEg)E1YD?1hQ5y_CNlM$-uj(($I zm2Q)f)@wQ?gH4FGg&lQJ6t#8ESvGT+Tw6_N%D}^W)qY1rClIue^31&WQrOI6p;~q`W%I+ zka-r6JfNGqUXk1_2KbOicNYtfIBpV)pA>O4vS$=z(G^M4BwnNDida(_Aw+%aIweW? z{7g|~lnR%%u*~6e?Kn+St0Nl9Adok{p1c_kDjRXaffLwse3~d_iFW$*q#wN!#DNBq`*wrvNx|a8j5gq&rFfCTI93gvqxa{Dxp|fqGMo ztO`N5kL!f^!J(Xvz5JapBb0-$aPC?x^jz`FkjIpQ>{aVvL_mYET+4FS)Ccudi`bfpMICt%dS%4bbo-;ss5iSN?>hv^ObSjUHTOct8K@Je#kYp ze5ixx6ytwS4cVA5e%4gESUNadT#U)F?(0VU3>R?SK6AM_>?Q+$crIr|F0jyLGm81O zMi>$c=Wz>QaWa6?7n~Mz?GgDBNWHyS@=7MR`PmnZr~}1Czv)1lDT)KB#zF@qrR^)0 z9VmY!?H|md_;d_~uuQ>aqX;ApR52RRaib<778T3>QvyOuhWp8aezBlCK$(u3gBjpI z`D2C9^~8T377Nf3@E`7s;y>+#3W@(PlPM(rQ#lFXM&9&3 zYVMlwW=E@*DIKj$f=Xk_Q6~>aozRIu(w@X3Topt~kxusSOdwxc17oKx&4rkZBUS+1 zNUiJ#0+PsrO|xk+QBw9C+Q>8v8qFn;D5-QFAZe6!IwNbT5A;*g(qIc3t{`^I8KGbg z2WvS~ve0B7QP9LQLgId+poJC%tyUBiK9lk%3QC_v`NMs;noXp+L_rlxX|9RnehQGt z-{S4hV139c51#sS#MdAQ8KCr|C^mA9;A7@ISk|^Ra zRakSMVsb#sH#-4C%yMVoFWD_YUUybLq7@UMB2K#_;SwN_^(V#T92al8OE9&{4n->w zt7#t5%BHonYWgyy=HDyioV~10-D!{Hjxky5Xd7b5{B(elp$}SAUIb{{Nv$CFsdU_qkSJk(|te?SsaC@i#Vve-K69Gw| z$jybSEvFH-9jLlig5T2np;Y8b7(FZU?JH6H%gX{5M)c29U_*_$IBFo zA$hbdncUcX+7>#-cYg?Gwz_Y4L-m^u(R?%Lue#M}!%jWJ+|sj?$i`S%L3Rghy@KrV-_@w)DjJ%#z-=|1(+eIX62-+!;7`_-%Pu8{JYE3>5$JpIC(J{IwZ;$z|0Kyi|f z)fLiJY9X$;NcB%#y6O^z#K%-u2qeKt`x77tPFG{uB2Kzo!rp+S;MB^3+F4L%SOC)y zPF<|xn&e_^k0Jrg0TF5r)Ra7IDb%DMmP#&{q$NOEDTTzt0B2>t23W@l;cUcrBgtZl9Y}Cln^8b7?3W5kWvODhHeBzx($%- zu5X68-|yc4|9jtio`>0M@4fo0z4l)FJe2&U=pLYq;Q@E^VwLjC;WAKSTs$H|qM|~g z{MRLfo>@pdc_MK`>V^gE4X*J&BM{=*|1(1Bzaa?V{|kp2eBXZ-VN7uUpAotV{zZf_ z;S1J(@&npKtT0h38W<^&$iHH8z{v1uU;#vG|E$$A(Z6WrOZ+c(EffD62a?ER_E) z)Zkob38nr&3l*UHXaB-dsh<3+R!`1lDNz5rR@SgUVj5TxKR$yJU+&dd0t- zj8+(z7A=Mp_hJkS)c}{%um>zWm&4);us{|%XtlWF9n(0amWHz1`+v5K*7n4&>m&*$ z@{FAIB_K}vQ(qsEQ2LgvJZYck6-M5j&C%FD#a`0HH?Vqs1owe{E#vSAl2Zye_@!BZ zT(Pr$VDiD7rVeSj`IY#OlKC$8iB(?RV#Z|KW70Z{2l3PX<;~J`#l$$H;@aEPua;*@ zo3kMU!^^gHhh@)!rHis#QF7%|vB zTJ?IGrf8MBM6@mK?ZO?4X);9ovi+oGppnqv=aqg+B{iK*XC1qQ#J7-KJO1TM-RjFG zv6AKxu+lT9{AL50UfE=f_f6(U*K+kukMkueo&^L_8MNuqIWO<6#^U!~CjsR<;(VM& zShte(uRlsXx>lEK@2cK6(i`LC*%`P5Sy^SAC9j759;lmKV*1M)0sUerm04JEyuSza zD?5=FmWa<7X{*ke71_tlTqj0x_)MTca=;zWT7z@!~LAGVEKdh?pWfVFd;`i+fP{&_G?>d)-{7{|GvC zkWs>R*_r@K*jM&y)UP`S`d5sT2PV!%0~gu{NMK%^9H^mm&R3{kb}pCy8<=Taw5Xv9 zEzrJOL27ywzbzzFjuV10sOc%Enu zs0X+FNB{@G3t&0#Vc^2~vw$JkhVUpTjHO(lg%|*0%Ev|ZK;R@647i{qXuO>C1}*^n zjV+}CGRs9Jxbc8b@T^-P+Hz+I1P_oXFI8(G057m%z|}{ArgABrH~=^(@6yEqKzQuu z3_)a=+0leN2>fC_s??3NSIopSyb5&$r{Xspf~ z@N>>fpNkZ*0DxQhgThdJAiR9EVgnzXgqk~daZDph>kpS~Dlh2kFX-DjQ7-ucpgh5e zRyxrf@F_3vFady(a2a3*8r|_zoLSLew`Tv>lBm%Wpb5c3ZY)m2NYX54$Sd#{q64{ zpro9CPY(pDzrY9gp=E}G{=g1@>chs4Xkq5I%(3;FaFN#_5qLWmKw6%3%!Z4eh0|dH zSJQD4v;a9e0+$dSSbz;K_=F1;1Sz1GAVei8@Eruk;>LhSQi7Th4=4eYEmYhWvBJrX zflwbnhX@KV`GF1*?u#IoL(tO!D05U0{R{xg9xd(*dmbZXcmPRkuoH?Q+S>>zVF2as zfD4A8m5CJLx3jiYJL1e1}>-eCU zNAUZjKr=%BAs_{iQKNWXB?mziObO4@27&=fIJXX<1W+REbO1(x>OXf=MHJR~sKOI; z0e-v-ox~N$;lKev1E0|aaG-kv)r8pW9102oNl(q5uOYjczMRRlw;-c_T`P{?uYfSP zj^?tS%YkiVV<$Ae&h@pvV}d3E>ozm&m(? zMYig(FOg1s3qI_I!*v(|g1+DRF(4Q-MbtX|DIhfaIk_W`fsuq}D33ukW~iPmi!{hff&u_aP^C#J3&4O9 z81h6YehrDC-!y*cC4`z+7?%8jHh}iQC=w4+_8oykr9p-? z=)e0YPlXbRk-6){l*Ww2pcNPxfQpepFj&CGCY88IjCwH&>j#VSPg<&24TQl+3=ApN zZj)-1LLYABe6Os-n91PY`Lc=XW5x%~ch<4_mA^-Az19HybSnR0NMFD@q~>{lwW>)4 za%9hR^92SR!uD=;{fh9pcMX&RxmmylEs!FZS`u~f5J_Zzx2J|-QRf^4K} zEpmOF3Bh3Ss8`o!gSaC;;#`wdsBXggIyh)(4( zT&;X0Z0I0R+0~QZ5gzVH7Hht#z^XTp`j6yjMkw#m2=qy>)fG`Tv}A55LMC7IJ>P$& z1vzoVx413$PQE_As@AiK1sU+knM9wN55hujXk1xfhr|frf&?~vuY2!6ojUx8AF1E9 zi`PPb9~l6l$q?nDfnYc#zbMPxkjIF2*yP)ifI=`b5yAr0*pN>t@hh-Rb_m7~KTd8N zd{q2Ia`cfGkjg(ed+b!;n4DSzH!JGN|NFf`7Fq<6l*#+ zTPVbnl#N8gzY%$y$V4>PD1cNoulno3Xp2lJ*W-81=Y(*1CDUcry+DREh(c$-1lK`o zA_&=uIx!$jv)iQt_ z;#Xuz!AZnqfQS6-@4Al_6V%87PnIy>OdTYA6<^zsA{bfo?x}UfqdNK2UyGeD8h7Ne zxyI5w0W~Dgo#YedJOyN=KlDxTPJ13iEPv$u<5p#Z^Sey7vK!F9vQy zCg(DSPexslEo-Ir8=Hbip2VrUJN3E{>)|veVM{xRP>*#x*iQm^P4of+^u|jetJ@Im`*m z0-X5>DDTn50>C9GK(K`f09RfC^;baqBLH?XJ6>T9i~_(T#FzzoOF&H$<)P5AwE>vl zP!k;25#R-k;18t0Ob9JU;1?KlF@e~31z-T6jd){w zsKN7et@i2yb0Qc(osYCH^$EHH!mofh6lat^i1l!^LfVMVH)w@Wv*!x+UCB+ll>7Ev zE+RDr)q->YnBZPDTtN(si*QOqVW37X=vb}*5gh>Gw~W~X25X^)&YklZ1)!Y+verQm zzDM)ZK`?wkZ5tpqaA`p;08@C5n^*vlfw;i2MvQQg0w4v2(Om$b_Vz0tNodh6ljsvD zIP*$Y32KliZ~XiWTCa0)Ji0uCwqJGs>kJC+yHaCB2QFIyOyd35feOKf0Gt5Oj1F8t z5;TEuszUFFe{~RXUjcbMP=OxedlMR9GC+{`08)gYSVIKL<+>nnmf?kopi)Cb%_Q0~ zLqzu!+Eb{BbCWv^5fjtEXAoY$;e~z}@P-+Hrd*`!CO+81=+b-pYv*Brs9!~Gxp1g4 zz(Vw&JL(HrYD{y!>+!!6?E`I-M-L71Yk|AX2rwv-;uvJ2qDvQWsV=f4*hPf;V3aKs^S2ZaqB0a2(X<_J*!Z^k;qXO93~=#`!> zC>m_k-&Nlg1!Mk4|Hc34yZw*8`~T>B{Exop|LA+6`Y3sH7*4$C66oS6gBEzz2f_P} z0cPk`7mbAW%vB$adgTN3+|@Yxj^RK1Xqzscb=2mKe#FmfHZ(+t{y{w#3SeNMbI=tZ z92AU_hyPy#xNr}5FgNejJo>za=6%)2N5TI|PXJ$A1NeX=1l>9y0U{vGxDfI|2!vXD z1ilRy-P9)Vd69v;j#%DAM>_5LWu%81A(nrEP!WzAIL~+@mjFm=9pMXrP@%4{e(3eK50B25W>>RbolX{DIWz~Oneg9ng;@POlU7$#ziizp!lDS!zWeqHDx`v;v8 z2x>6VZpFZuz{WgR>ueqIh7yEQi)7<6$1tOSb9K0{0O3mjB)NfLphxQk+4%qJHHE@L zCqB^X211q>l!yr)+eNKKaQS6YA_aG|VZrQh7y$uzLm?PdPiVndz=q({4K!92%pnIK zp(O%BovVpAFBdZt3OLVuLMY%|Cy6US{t^JGZX!^r<`@8JY$EuiK_U3yCz}AG3kx{G zL6qDB>4WgDb6%LN4KIBLVkY<@f6?7mxbT3>=0(=G)9yyB6L9}0H|yev8V^q#|O=9BHo~rJZSSG2hlzO#efOo4KL(3Us)LFs`fiC4Z%GqtcYXbi^^~__P=JYDnszd*7>`6 zh6`HULJ+2aYA^}ndoTUdf3C&#g^RBC(uIFEFa5)lg;oQ!eHjttDB#>LbXNe|6~MO* zAQx?K1M2;N#kK+X+7Up7sQ(DMg%6fP;hj55efz>uBVW+3>D9r7+}tC~*wP#D z_2-zRh{`IIOpwVo;(j%Hi9AP9o`<8~m8jq=Ao3CbC2S)m8bCuJa2jglTujlWn6Os# z`=DzRMcDvDZ@}n6vt|2YA>pD6O~gzGs2_y76i~!ug7>%u?&~|a}T-qgx0?wTZxdLuq0wB#Dbj*JSf%JARq8%K9nmRWyaR-1gt~VoQw?Mp@ zpm#e6!hJN?ipzOW>&}IEwPUpSp&f+n2}l?h#0W;%@c~#MP$enK7*sldpltvx?ZAz3 z!7Tv6xrGXNU_5{j4Z_yNa6w^2=c6&`QQ`}#jC0G6P{Ph%4Z-fE2^3cV{Vsr5p+?`f z>>`$#!1@4)42+1<1>FSTG4mi65dZF_L53)bbAxQrc?)%v>=Oh(K;0ao--dz^N+QF# z9s#zvKqj0O51qk0&xLzm5r$s@Z&4gBk9T;e>n=2hm^1)&(h~#Q0^s~Bw69k{CrU@u zQO^|&w;y(nK!^|CLe<{n-i- zKbLm>e@S!WT$;qWG)eTQDHscB!1=0?KVLQQK1whyXdkD%|Mx1+`4U}3t+ey}aEJoV zQyggT0s!If!BHjD1iViUOe!vrb#Yn}xzt4e-bMdv8O7IMEXDQ9d2J@F^STS2xybk6 z0t(<``u{j_onr-Do|RI$uwqLYVVOL*pg%9qQ%T7%E<8pN%mXguKA$4Z#f8r(g0sPF zDCN%8727-CE8uU#z|v@q&-2q=)ZkUIW_Ve8T^4H=SKa?xtl3@~3U)z}O^SB$JwCip z=u$0=|2V0h>+5uB+le(=Nz@?aK9H^(M z3f$}Z;P_D+WKJzPm|tTqc_O!YDl~@B>h8zCY9UW-wbVG1`DcH+JZf!T;d=6>dQX)9 zhvJjFq4#g|JxUeZmYF^>nLZjnTmCZaJ+vijx3=dkw`xfYk^LY&JeAItAr zOh^s4`RP#p^0>^RA}nd4`It(2#IzN95dKC=WFy_aM@_G#xiC`umMC6|>4t@weSD+x z-B(uje{>o4-TwA)Ly0DqTPKwT(^VQ5CU!+vCX(7IU8nR|FkZ)6`*?Ibk1?v)W}D{X z!r!BDnYd>JQnt3LRq2+^p6*6Eshn{qRkGFuJ3He%c`sI5m`iK@w|_t=%R)j8nj4Te z65

9rU!sYK4v^@J&Io+=n~H15Yk0|MN`$mrszvcpscy=?)i;L09OM7bry@nYR% zZmSXD9`g{d%n+3P_=oA>Sos-!7{uD>or6^9+(d%>WOBkw<^5{avI#K%K;%VZ=%$ zewMK(SDxS4+B^A+qa0DB@d1yBl=rFql3WIX*!s7-X5LImr``|N8WFbIlLNN~ zWN!hvcZm7h#G@W37pXVTG?>;*Keko0E#~j`&xGEbh#(Ql^ZivvyAu4`?j%I;d-ZHb zRJDHs-FwlYfrm_HX5cC3@7^24*d+@;fME_%E~msT3OQ#>8ZIA4-ckqX^`@pQ#IZ)+ zVq?2ZnZu@GV~M3}sCrnMyu5+Lk$mPod~W8lhc_s80)BkR)2)^43cGF6c)Q9@GWYJu z#*t?s)AFir(#hCYBG~4*A3mFm!qWY-(X!!i{jnGKlW0wz*7{^c#8A;lc*rzZ*7lN- zm!>q;+*t|B_irgMcu8UC0L51k9#RoyJ^ktb?i0@o`BdI^%KpGW-`0F!Wha2vlb!_` z|HhpkE6+&#uEN^SGBF%EIZqj{&?Rsj@oASD*RK;!xBJJAdDA;Dd^B-h&8lZ%73nCy zc$2TUV^zNGQqr7Tz)V)X)06+O%K63`({u8pQaO^Qpa>se@do@PC^u+!JAIq3()}|zaSHs2H$)qVsie$;m9)nCV67W?11G9 zf8aIgVI17=A4yVo97zY|$By^!GvmiEb6uOH1C)~zdzAGgaUu8Kp7j)xl{jI^-oFVX z`&mEZ%!_P^q?jjp_MTO2d~G#sDzdShGU5(vB+o?G^1-#PTc0=6E1+KS$}dMDFZGDZ zzge!>ev!wYm*ZKVHH}qJ6L_Oz7aKiI7kKAu{6s*V+J~fu$B?%R_?F+b$71a`xK2mt z<8Lm5(;i)$&w$)$MdH%97!`&nd7JDisb^UW9O?e@Cb1FE)!K%;D+=2SFWO~(Sc#mE zzu74T9$@9wckWV_;XZ}>4)T_d+jg7! zhlr2c+x}POjJ`j_be*Ofi5=CQ>?_;dX~=)L>OQB>v0gYfPfJs=5V4+an5$WlNHS*n zZhV_GY4a9Tgs+(>%Kg`Fs>0^Ji5HpokLI+kpzLy~t?e=n&;wkq zV*iEi+`}{E)UK2ae{(eWCHtS@qPk7X-7kUVFa!LO$IWp@b9qefCV1ac8f8BoO^0sg zRE%G!N5j#=q|z|MU>OPRe_`Q^%>S_5`1~ZQ1FL521*UH|a@- zcXboYKw)T*v7ol3eX+cug8v@qDgQ@AH}eygwD16sBfboW{V=r0G`DOnnve!2vqYjW zne>H{#|jj55CGXAi||24`##3?D3P-9FIj^7QhuBEo_8jTbfeDSPu?Xb)y6rRm8RBQ zcsuaO&*bkgoM?LdGiCd3da>sq^d#u)&@16rXaCb0|H3g3;XFbY2XWW1WS)A}V7pg+ zZQ25gImLHZn){&E8oXvTS;TOoCu(lOBTdrNH+lURscMb|1e0zHF%sua;Agye6m;Ex zl%fPds)yw$NfzA|8cI{y8HB^>G+rR3IIJRW{5S3WJwJDX@sRmDn%dTib@P9OEOoi| zc?Tc7Rn-k?){UH(FlN${w3rC(7H%yH|^l5|E`ZM4N+}0lP{~WAzBbA96x%BH=!MEqh@9dLAv}!-n+eA}{uR%9wu7~8r zWM|F(!Px2!Z+^27GO`ualY>3@yyU~fth?`CV8rw)M@ZR8nVX1xEr2`;!jXw<&Rx}? z+47%ezjZ&Jl@qG5q;bs4CB-$-AaI=l>y`3vDjI@&BV3zzU0Y9hrYx9#yCzLjPCuug zNFB_Ks`=Rw@~gka?iw4>jG0{Z;%$HD_@>!uW+a)KPOG#V;FzQPhn4lwqn%nG zvOIlAbynV2T#aJ{F>oQ+>Kt2eouV(Wb<^`eO%nOEj3iY{rQc|lW|&naKYiwGYfTq#xn#HB&d`0(Z6j*)@d$^;6xT$5pllR_%e;Tk9%A6U-C;mv zpyDZ3xoSCVC2IV5b!L&W5>K;>TCPD-S0T`=e)lP{uJPE$pD&gi_hqf=8moQ<{8jg_ zp_7f|u{-u{h3?(cpU~85KT~1g_KeyfChTT!m{=5glzoj1&!I8iJ8>;J0{d?4R@i=- zG_gAr;|aT8f>Gql!lOPqCebvE&pxHDI8|KzUB9~~0-5#q?lN@7BxaY*1xP;{^GR|e zrkCKn8K`20l`z!1NupeNi@>7`y1O(mT&^|tZo2hY-39ts@mniasn2zWgJ-8}*R{gR z?RjcfEB0pj^-eaP-&&Abcz3Frv`FUr&iAN0d4*R?T<@EPFfZf$J|o=x7hXdF8Rd-| zf25Y7r!milc$s}w`zv?vy=JignNMG5ME?qgq`0H1*E>_s)Zmw~|I&O>&?Y@yu(dD6 zEB1ki8?j^t)ag;DG(ASS%%k_|y|FW12?0evsc@VkEmR_nk)rr9T*mLzf zON3mZ&r-71Bm0Pj(wnueIwtMbHZ5o0H`1ic`VccZnbK)lr_ix-HC0}8^e6$XO4_4U-}_EUur1y91diuKARCb z*kP`{J4?S^x95HPQN9Z6^NGgPCQkxlc9N;87j__0@zG9*80NZDl4mNj5@?Zl*KDim zThyL40y!1O@rJ&e)6{HIoP_V|kc%laA&|zG`<4lm$3;I&T%H&kzy9soIMT#6LnXC$ z*CbNnsk7v|sz&6Zk1`y(=adBl@I8rNaPT@?gg>co=s2(;9g)b^N!g~8c0SoTY#U{^ z9vPKjv8`d$lA+gOn$HE;Qn<4mG8u!-`UECJAz7B(y95+&<@G~ezvG`+(BRTP`B8L3 zFgq0bQjR`KPynewvZB10L008>uyUsWfpP8cwFO;yEKCFGF(({Jw)UDg@lh`IFx$1s zg*}-OlNUnGT#q4jPJ0c^Fq`>xryOaY;ixoEw!RzLQyvh)#hpk!wHk~zs2yEd-D6yW zGTx#shfgB;ch%I(I*9z)80AN-975IK6a1WoT6xg!CR)yGKEIa=lcSY&S0cifEgGim z|J+QHHuhyjPQyIq5U12xw?bsr-obah5HFtJtXRUAGK_@V^-k8yY^{4?>!xE_geDO? z49(LYe>G0CS-9E^?Qb$(cUBmphT|9`xems4Kd>!IiLyKkNeO?;#(7WvrO6^T#aODC z6x8!>l|)|1tLbo>G2*F1EXlCJ@SW62@JIVG6FT0tzWWYt`{4m(goDU_qkE46NQkv* zW!@KKu%AJCka-t8?Q=5tc_r;*vqJkd1`Nz zvh4pYrV00~oCtN!%Er<<5l(BpP0{dP@&}Y#{F+aeQWSgfvus@B&xu1fTZALey6q54 zW9$J3anJ8q346|%1@Gx_4AP$t?!5DXI_~Wq1NdvApKIpt^xa!+j_mLCO)_~X{Nbi8 zadX0)0e!DKn{7SBO)K7KIYvmgL55=Z$1^aN^Fez0{%}!{)=!1!iQ_C}d|3_)>i@=zflb{>_iysN3lkYGzlKdB>DkBtNpv z-=3${XV?`NJ~F~gA+-MZ9PepRf%B2pQzL^pa_^(+#$I357g& zmGo7kuo?e0)4p5g&0=jzZQ*uFTbdCvdjjo-tC)@q(f3(e?9lzi{(gw#UZa?br9CJf z2f?Tea`f^Wz?V-$oc=!>a{;4Lf$KMPLly7B6TS<1i6gwYC-MTiZb1w#XJ zYd1nCw-)=4N|*bOXbxtQSFspZK5%lZZFu!oYf-mN5wtuQpR2#upX9&0gQYygP)hs( zi8rNMV34pnLKsT-E43n{@OD%}i~mi(ni1W%Wz%6*gk=!su(){rTh}{GH2ifY&kkdJ zGAbaCkEd{>?iNZHR~2W@ymz-Ro1TLPDrb|UM(4x312Fm^fqI^*r6ZC zCqzz(>9^2~m3LZRLVI6Ubl3xO;DQMPr}*k^bKEZ1c3BRbTL+ey-;MpgJ{vk) zHQ&8IScRt`&uzHG(-)q&vD23j04U!z%#$x~tFKrxbXVM_SJ;>`R`P{LhAA;jPaQo8 zhz+~>-CnXw5|h~B&4w~Yw(%6gMV^_8W0={LaUX+cjai3-*g|>IDs9Phm!W5J>Oma8 z5ec@kY2ESDwoMw*iTB+7uG^g*qU;YxKH|s7r-@0VDF;Y#VJj`bffkixD-BO_w&6AM z*F1QN&Cx~Ez8W}zXW;D)Vklo$u zZ1ZSWdqi@rSd%5E4FAEwN&57j&#BCO^wQz4nyKL(u6vSn%UOgk<>uSoeO*SFk2~HI zr7nH5o;+&D{kt=u1Fz7aD0Zp+Q+QMDLrdY*9q3T(Oo}aNFr&qw^UG;~j~&1bR?RJC z)r>-Pc1nCm`pzYz+XCc>VS4EF?KT06GiKh6JMFv9Q^d(~3?drkyZjqk)EYc4f$d${ zTAwZ2rMq5JY3zJY^sU5NFRe%ol`tkxI%r6XCD9ylp6@q*z>KF+bJo5^vU)o1ECFpi zfbwd(`!;R%?F&W{X|K9UhOs{*-iIM#Ox{bcC(TGOdZ$S_WbjjoKbi@4l+hlfU3Caf z@KD<^-w?RGtl_=W-W5s>`2-YS+dlg7jb_lNHP&gylB&A8HZ1=bDsR0kZx^^YCUEoZ*-I$g zN{@yo-fm!eVHC65CCy{7MVKY~dDB;cH3~j$s=1NxsbJ*^NkbB|G*$TR2hAUyG1yX)L)@L^{InBV*}oF}-sSs)BjA72(Kp?EzKQ?aB7ifWqE^8D^6=@4)2U z`BdJ3l~%6mTKeN@js8B+8}u-`EM(e>wWC$djH{ARpV`j1@Q-t z1m!ZJr{FTB!DYpLP~UVZo$bjt>_}JPBfh^k0zMpxeAHVY!u~43um?FNl(%3tYJzsZ zsONV}>)3wEpnu(4BsD8PVlv!eMWmzZO)o7xMb`CuYq;<3v@-UK*&HJO9sEZgxYIaJ z_Gd1H%+y@{1zv@-7AG$0Utbcv5y}@4#MP~*`GaZoC|$qtRetcAP$7-eS&Q-?`d4D_ zr0vje;>3}NeIJa9!TP{gsi5=nkx6&Fe=`p!&&cm<)dkxcwg&v zoM0on3vGqzlwS6l69cvW5qC8_?6FD6jk8ZNe(Sn~N^d7-I1G2M^+^)pl$egyuR#}@ zEIkVcdke;m$#m0S`<1L8vb^DmfA)3e+v}zlxhh%%+xicNepYV;uFa*^)rgGv5t`Z* zhHP`F>ocaQ9r*+Ru!s9&-yFLGf36wXFMoT={-s@O?zcPnXD*-g$f|7Il=smF&&tH= zkJ1iPd#e%V3;c;0JjMh8dhVGSGuFpY;o4sRDe=^rwLW>d#^u7DaOnn(-~HePOclvAtIe$VBj$&Wv8 zqWv>sK^e2!12X=}rCF$N6>2cgpUqPz+V0S!sBBe4V#Y@FXXVN7n+&O#Cf<#{5iq*- z)8l57*a_as!|xjG8S9?bHZ%{rrM&m==XRi>w{GQFl^=Ye2^uR zkE|W+#TBM-bgEkBr*-m>muu`Ks_qnZIHGW)!nudF|51lW(XC4cMs8adlYd=Cam)uy zyk)3TdvK?rO!FO~tWV)%%`98!+dHpCggy;qrr&%TZaMX9CzsBM!7lbwz2Cg-#L~j= z>_k7x;;syDQDfD5$2qcy<6$|c^}dA3hQg;#Vs~FHlNJ$@R2axf*A8I9r{~BvMB`os zN~bd>c0cD)6g8Y1$Y+1ou1_7FpFGn=Tt?2UtugzKBXSP>m9J0XQEyPG9rRl8bcvIV z@f2{g46O@s9ebd z4$Xw%M?`IqO=|&ZUd;aLm|&C%fyw?IX+%~ zHj{1G!ZLeY&q-2mX;XE%nUv7m**%qW?WYP+-ca%xpYsTK=yF4}{4?DOH0CJnPYq8< zF(b_vLr;-t-y-VhWS{;(yTu)ck0&pbyWbejuBZq6CH8^;rNawL{LGm82=;hXP3o@B zsm#oZsv^fJ*^^ROQ2EYRKr!vx-4m1O50uV~6N5m#ucBLPIfBLb=Q_xy@)~KyTd}&R z1eU&%Kjq60@vYOo{}d)l3W2gug5}d`ou>$b>-cY^cbm>oHBVdOqC-0 zWaIapo>0dmKDva%S2HO!y~lB_{rma(PI|etME)xe>T-X1F-S;qwn@JoZ%P^;W-$@Q z+cT|U4Jq#W#oWR;j?fiidFI6%wDwvku~wZ~4MG!nryQN?7?QZ{P zi=gC{DxI9IBDUURF_X8wb0_3lYE-x%*z6jKDhA9crfzyptHJid)9`-$)wC%Xey6ym z_Z|u<;8i`_W)zs-Qv^XhG5rmT3vLDKZD&^yaF^y4u4#U{`)CBdKCGB28n(|suy4_b z9pj*bnf&oCHGhKYMA)JMXRHu!%lZ=>2Ui{Nbfo98-?J zu$@2d$GV&~72dz8o^ii=u@X9kCmj9j8+4~(Y{d6YkJ*x1fcGqudK9y6x>~9GqluSZ zImx3>J+yuucFFuvD z_iY?_A=%7Gz`UGF8Qs@LLqctOli!2j-Waj)p`mPG*V9#;lt|5b){iCH>Y5bxvpMfp zAI(S6yaly=l6QJ3uWI#&E96NbN4jrN?J3- zuKg;030d=G=xZ334QpOug3m~Ye;edEAbf3_jAx*Y{!63jWO5FUng*QE8U1 zDl0ll+EpOUDdi_ka@I5+7WbXK4I5<2zb6t8;e5?L${X%xOlb{G6X;oPe;4A2%k-_n zbsVL5K{T}>R-J!NiDGDE=De-%kw|qO=di0DUOV+k%Ss&^7I(Yh#7Y&EVc;v&+4E*k zPSV2!+w#okC35knaFrYm-OfE%VU_Prt*?tH5raXh=3ico%Gv#8Fw;^!dI0Q~J`-^f zOC{Pmr4kykW5a{u1+wu{_>jxiWs-I>mIN&Z$v@ccMSux6;dJM*gcnvMAWk1{&ZT^)BH54ZzZE9S!e3zaYYa(njI^1s_rTtmn#HpC~iMY>1 z$;slO@in^F)e<8G28wIP$fF&J@*nc43=)>k#N`C3cYh8@MMkGAj(sLzsCRiVj5Ljs zdb4N8ur%#;8wX0!Wq;>a#Bx&QSlL$78Wu}xN>}8t`pgeu(c;^8^qfZAzQL!ale4-z ztLG)xoa!^}x^)(g`Sybz{3r}V;AtDU^XL~F-+IO5EqSB;RRptRt>12`$hPSm-j4%c ze0R-HYZziUEvBrQ)9>1^82 zk^1;Z^@qpl%P&1-AadAvSz#nX!SreM*b@6)Qrkf%{kZ5`?~cMK9~W;<`x4NkXqui} zk4h|dU-D<7$NbUujn%zhE{sl_{nHCS$sYxW*-b!Zq(|>2G zc8)Nq)8TRcax1n~v|6R~UqA4xR(HDUV+OVfQ_AHui&tryjLJ&8TS^d7ih&gm<=L(} zsPa-XHo)0xy>>jMsC<(A^CoO8xE+6^@aZZNFn3l^l!Mxr3tPQ>+sN@&F(-3aT0y7%=|_Ku3_O+Udb&5eZE9I*ch4fu&ea)FFLkGlJ!=EBt8L=6Z&O*}lpx&%B z2-hnPDqsm@OPG1A7BMUS*zfZC2sex=a%=pttxCWae$Bt@kdrlA`)gW13%EV>I%Z4a z<}Y#E&p1#^hrd5<#}EGgZb8SilK zGB0qeC8chvZ{Keap<7z7Tzz{d$YzPCAw>8WhO>&*aCA1S)1K*mkIEG5lm64VI$Hkl zX&g<#Cs5PFNa9MOlL7IjvuR)U)nGmN#~Uo22_YL~8_WS>w(e4_PE@`O-%Yhhjc!+| z&n>u1*ww!&J%#R&v36UQXHG|GyqtsVeTri5TpsJ(Js|#Ue0JO3@*QKlh*jusOi>Fy z$LWc3!KxI$92sbsAVsWF3nfDYoNP3~N5*>McL+4@?IHJ9FxiWaQwo~DEPp+F5C=ck zGW#{+;T^{`i<0TJS(#tMh!HpcI5+jk(?1eb>zX(=i+Yh2|65k1sarnBzliF=^IxU} zUeVd>Be87I;odOPuIt6>rj!x{Piy8R&<-jks* z1_wB^8<_b-FDER^xkP1$m9g=sw41(WBYz!M$U&B{9BYzl`IKsuUG@oOoIPF58pP>r zzv`3B_02}s*&fmli|~}EzaFm_-tZ$*=^+=Z=JU`0*(OxN z6vX}Ba=hWrGa}+=gVXi57NCf?n^U-%TnG6Ev2A!{GHX~pWU3#$*WIrNyVBq#&~jFj z>NbQ9OiZ!r8f6;aD3(?VBZ62d5=8`;Nv%+%AU$?)Kj@g%J@tA?QJy7r2FO^XW{wHF zrYuQq3vfotQGI#IN1snDODk0ZPcKd3B_ca6Dkm%AtPUYe$CTVlV}Uk?2d~ju2@uCj z-G0DtBx6|kHOFq)C6S?Su7IyU$Gkfs(GX5)XeA$c>ztw!*wWXcRLJmus*G9E&DZdUCA`Hmh zi)t9osSex;mb{N&(>xM99^%{Q#LP%rG-``9P?Rw{di~cd3W(VfWm9qh(@O6{Su-YY z8J{^DzNhh5QJTkrvJtCtHH5pQxdn9LYJEK5t!drF^VFo8wuDF$%ZKdpmI(V|y7k8h>6drPkA<%W%x}*svCO9wf95Bn)6J%hIeryMDuP@wqBDrqYtd(sp1GCP-MCdc=IdH{j3tiZ!EQ20!J@UK2Et;k zlZ{S=58QvP7Ns4-;np`lD8a5jZu*$7k+0~Snm49MKD!o1vzmS^%fS+Td(78|@wSGd z24k75DZA*0>(8O%yFJm&W~t_P8RNkF-S64FQ=VHU{VrHbw>Ti6BB&qZBIdWgLw_h* z`+!(Rx14m-`w0s?yzPxzEa8;U7PoH`mp$brRBz=o#sEKJdpB5kLMK7_V2({8)LZv>Ewy;7~fow1tB7 z#V(akmeZ*`O}ByHq8+H9BxpExVu4;}`2w9O~rY~hl9ZQXXQ zrSD7oUASD{yxhm;Pz@h8VA?u)x_g>;0+ya~wG`>fd(|Kl^4Y}!5;$nOwIHsfZk@9T zs`BLlQo%NZW6P)i3BeP73HxdX~UaH84It*N?P;QX48?5DPyhe73HVP-|^{$Rx7uSgNViR zAE>kvvOqIm=wu)M>ZA@D7S=f~H38vab@1()+-|>tEvWx)C2+rP3%kF^EH*rIAlKMb z+$loW(d_Rv#p!SAg9IwivvD_vnZD1DYfaItJa2i`R$|T_oNbXq88SfGO%RW54qxQ< zYyETJx%7a*kA6W$c02u9R=EsOwC-Q*BFUkdIUg+O=IyR3sYmUqb>5?l&+Z)kr73*% zEC0E9LmcVteI>S-E{O(Q=Hw%b*Jh4c%ZU`R=??{#8R49m77)lqWH&qyPukoeMU}e<-5iY#;o+! zR~xP)tql#N{-_eIZXiA5*N|NDp;_Qq`~HdRoe<$U0Wxn<1)J~`2N+dw!d64atVhP+ z*^QtIPn^-TQxS&I{PyeJ+oU)_-Khn45=SarvZg}=Y3`*Z(b-RAHc8)$!b)}7QVV~0 zmTP!KrN5h@lx$suui-7M8vjY*mEOm(_b=?*ol)zL-e{`e`StxoM?Z%xUfM7F)iidN zfoTcd*`Eu+vsguVW{z@#mCm1SM%x9$qSFuMozJ|&n(xkWOv^1{Q3YlE`K@gBQsh%k zDW#bB@uQ@r(7m@@zMHqSKeFDw=KZ=@o0nRE{cR-|wriRqNB!N8tE%l_?w;u*#6rGO z$SV_^Q7E;Y@#omuo6GzGCnlrQ9jVfod=$SJe_q>C`x$TdwzbMVo@tG`F2=bbLR%$s7;jetmW7-7OD|!K9em z94c6p;0V-ldu=}?K2YD5(sZv%hk_(_ca$(@vxyiw_9B^%-(R2fc`Dqd;x*->0N*>o zrww*WT!G&$Z0>(g<08I|nVAJmXWmDit*A}d@CA@?K`O0_nR;A_3G+smK=@COj#BrU zlEwmMDBFErd;20;sn^C&i(5-l(*v_8U+7MIo7GFvQyWL+*NSy*>#p4w~((jrxX5s!nyWZ=@!)6;zsa8il1q7N zztw$~w0T@Z{6u2zIkD#501pLE#L5nSUrE-!_bdliO=-Egj9DkqpyQiQuEr5TaYAuZ z-}v+Tn#xYJmIgI{oBlqmcTUC@F74@#qUoxi;(zFDSd*WGNzb%#`07^V9Q~lUnSH%t zw$souGx>yedu-{hquXQ>-|*K0;0s)tvEMG4O>X=)Sg3fkBBbp&nLW_ zM~uGwF02(4V=CdfN>2udbDdorKTwl!NA-)j(M_jh7^pk$KC?zkhl<^N0#$#7j-8*@+9 z!*AZ_y3Rv0lTDSmS z6Wat_jq)|oX>Fy>LX&D`#ik^qw+`>N?B)6e1P@etX66(dtdHRj(g`Jy*O9wf*3eB#VSw4u>zNt=LS?SmzLYF_UaCo1>TXqD&>t3L$|gW?}|+-$e8Qe9$Gd6&5J?z);j zbYH2k??XXry@tf%4Q;AM>Q)!89viWKO#jZT8LvwxjJWT-mN&1zJ9(k-$+&+S8@Bj9 z9wWSw(ld<6Zock%xuVg&$VHI1sk36#48k;2UTxCGCgNP)w~$R_3fFkFG=(Ll($j)&VQHVxWm4C zz_pe;(keak*M9KhV}#0W2d~Yo+vdx5WbBB!(Ohp`e6H#O*XihJ_;n)b>JFC``9FDf zZf2cYg55D&x@zvv>S$Yddx`Fgxhcb3eI=$C9ydI0{^cBHSrw~P7OmWLr&sZsm6YQ; zV)OXLV_MJmJUMN9b$iZP)~iOh?)|DoED4tlHg;#)*}pcOmS3VW@qxE->EHUZ)Y|Rq z{oceZ4f#4b)3f2~8S6gg_Q(`R-WVYlD@B4{n(>tmd?j0Q)~4g)$2$6R~kN# zm3`VDw_>r^t6AFxIo;og>7z8i_~j;dXp>91rG8mg4_t6H(Ri0Q_k-{DQ}d~Eg%quv zmUXAo4z`r;I(&Df+zxIg`eN)KEMBXxt! z4>{lM?r0Tw$L+}XwO8+v^`7WukRx2Q#Mj$5jGUWrwd%C#(t8C(ZcRcv*=)IimeUym zxkK9$mS^W^3ga)I@0RZ?8M!H3Wvg9 zeXiQ{LbCo&ih4+w&Etg5?^>7I?EHE49#Pw#DfLL(`VBW6ZY}Rx*Yl%Bp`)|yXu;fK z!;Z3<2P7VuNEiraWF={o90n(XHe7V^dXrPseR~^!)5XHY4k^ifDjjNttG@qQWV9k{ z;gch#9tTEe^H!;=b<1=J5+o|@CT%t-v>REFra_o4@*hbebCVC1suE%U{xW1H2s{o;=3t8f%oBx9Fx@ z;$hH0S%990uOCLd3CIfE@b2sJo(Ww&wY#Uc^lJ*-Om58c=Dj(Tbu4nGLie&FxtquK zw2!R=I_(#JPKtaVeDRiDqmWQ;`nsqlxo26)!QN@vw{~$4PV7*)}2nu~SB0X!dy9Po@KX0~s~!x4!B?)yrs_+{7PYdKwp!oVM*ODy>&= zx)la*-Ia0a&TjIdJnf#y^z*&JYv0GdpW)aYXds*Qg)zZa2N8kTY z3eb8`)s?%5RneUmkZUHGJuYat%)Jzqc`gMDuQ(rO>*TG`o%UCrOpw~rg3GTx4r%@A zP;@8qe#eK#lMeA7A;x+OLnE4vrWtm{Hgws?-Bmjo)xtNIwan{j;N=`TrE9V`?W_5Z zxeC>(R^zrtrk}ggymjC~N9WtyaV9-uUat3#`W1Xh|Fx8|i{SIkz%CUDvv*fo|5_08 zV{2l`)=ie)UaA`m5-SzU_APzXl(I5&-}Tr@4#pDR&${InXnPiA>Q?CmrId5d$KMm4 zt_a?mweBjv+x$VJpF&h$Nigq-gl~Y%-QLD^I-OBLPshGC80DwG@q*m-S3gImwwbY0 zV@mFxsrWT3Wa=UEMYh`6gm=2FRz-@t^G<5D#F)A5NInqi@2A758WW$Z`9d6&>fk(mG)FU$+c4jEWAM)IgsBvxdsmifDZrEA2 z@w{Jaw&`z^I<@~+#lL&^Xl{P#bE$C5s@X}=RL`ma&EO5eP2nb`rd^hE%q}hIXTQ5` z@nxF8CbW2xFxXpFJE*o;eY_0^Y~7REK2zOq$wa@;`O+DCo!iZ;DmObESZtj8Msm^a z;I7^;r=O`j4X!n{@x0&<@=quSyj08K9V)qCwYl@L_n2;foe5#RdxcgpcAWZEwz~XD zvRP*~vRwP0B_FzTq`dR3J%%U`FEY%`qxf=)_EK6 z6Mx-!*uG+;Cb9j`N=1o?KZB2(<+%Llo_@CGpLg6(dlhY(ie(Jd^jZ4V2S>P^j{h)a ztZe4M(m7*)a0~x}b*rO3wyo5bRaN5MRk;$JSk0e zdI1CFA$4zZXZRDmQ8img4=%wc>e6sk-dQIpqU}6}sOGru0)jEcQiR+kh+x$(U2+`5 z)aa!!1$LJiaDxH4u;dQ@`I0&t2=)SBsPvX89%NA#zt|d-ryrQ3&Ag9<&cw&#qFid9FBmhE>4G#B#h0J~+u|OK z9AW}j=VB9sLC_jvLL0W}*Y7ed%#=Dv?*bITt-(cO{a;ibZjEYBi8ms81>VPo z5WONEV2J3IREi*?S9%$3{S0M1nOIHa0lEaV8vjOtjK(2RoEf9-lVAtVxq#fET*kd6P^ zSoao{e&6MH8()uJ*3ey@_h}$ra@hUFhqE%C=hX!5DONORHHRz4js$K%9WD=tcpO|F z3bASl_7;WM`yPfU#Da$yq7aw7#t?-#s{z8n5SP6FJ;XkM6p03Sf2R9k(Ewm709Xf* z&}zY!P>c;fVu)gF(TX67@q#vb7@fj%;7W0s?J_2!YXyLRwSa))DiGus(%`=sYqg7F z?A%Mo*sKr5I9Aa*&wa*B|Bs*BeVgM7(vN5b0R%M zaa#tkJ($7}Gx|bI|MBpN_!?2|II4kW-k(T4DFeEL8E-hGJ1QIJP?^qbUCQ_&nC`(0 zL#K@RGA6_JDEPy?kIJsar$t5N;zXSS+lZ*F1B%^l;f<1U=Pg;RMoVH7Cw8ZVJMd+g z7e+IWokcCFhT0g2d2$~x9T4-@e8Lbi4|oz#?L#{jG7mVjA&f_v%_)(ltGrjh>Cj7{ zNw;JFFyLkiOWBn56A(an4Dba{KsFmBpxI~z8(^V;3__~07f_vYBo`_1fHg{Z7Dp#3 z;Ijfj>Twm|XbcL|=+ovM>DgG`)GT?DOO2`q98EE`8vHJ^r1X*4j_?sYDIitrry)Tg z5uJ&90)E>!05!A({Jw7bx7>I58+O2Au~f%Pup0ypC#j8f;4<-Ez#m=@qNSusZF(Wg zW81PqL6#I=H4GkGj!M=gxe1flQb0=%KYe{miXBSbae*GfWCMos@M%1hd7*-v0f&&u z=bvZ>I#Mci2DC+Mk!uk%MF|W-jo}DaT%124pUz%jG7;NG25k*7gxL$6Cm{$!3OtP< z%BO;Bz-l62l)CdFz{lfK2CNbR!Xuxr)eDxhfyo~v-}eVt#g^V)FZ#%(*gt?sP-zS7 zkI$=se0*Mka(z5EbR6hSzz}Vs&twb{YXU7WM6B^PhY+(yCD7(~Uk#TtzUs|~ooHXf z8DFbxu?^&FlO2Z0*ROLBMD!u4KlI6U(Ly&Iq7!h9_@uhbczJ{IGB=-}u1vzs&gTn; zlgrse)qK91D!H01z4jI3Y_bX|a<+_d_Jzn{_81~(wU%LsoRxEd zFdi<*013r`kI)v<6(YLwDGNQaj$^TwjsT;U5OI!f$7UX+^Tzlf!2v!&&Ylb1CpaDCXj?l!%#=h2lp8@9`+GaF2a50iRk6IAN7g7{#^=+Om#Y;_$lG| z**A#xcx3|^RJP?1SswHPQJk1XLntnH3pT@{jHVHX`07QEK;mHkcLbwzMEErD!%fL( zDZ*qArENvVQkhn44iVPCFS9~V^>}>OdE{(1TNldZ5{d)-Tze8`OJ6u`+!`83uk1(z^H*z-Z@_q<>@vK7XGhr;n2sR|O=u2= zC^SB|A&iebyK*TQ0f#dZQC$292-!BMMJ&};L`RPPJgQ|=Zo446EQHjFa54c1bn*8d zC?F*T(nn#86r`x&Bczh1m>L@A&kgXuESBwEJR0qy^F8c5vOSGO(m|Zd zB2lhBrJhHTy9v+x)WT@;3Run}`9GscZ7CItr@s@Zt(hTE(+g;gK4Nx4BS<0|#8oqs zu~Vt!Cb%n*2=IoW`!C5wL@0R<-=~Ri0ImR^w$Tva6iH-2vIr2DSbXnllq+#B`PHc? zS29H##iHC(nOtdx-7z70L;*9b#H&EL5=wPeh;pT>7C(iAfI>975%4J_1pFHD?Vf5cVICoMr^WZ*AV>A0;`0Mw5sMC<3Z|R^)4l;(`2cLe?K7>kQ z7NN=Fq%cGYnJSGS+J`I|i1M}2mq_MX_$S2T*TUEXPO4PjPk1ev@(`E~3P?_c-NbK~ zrm&4x#+aElq&!4mr)%bNF|y7)cGix7E(u3 zOfwMr&SI?)rF`Eov>eK*Gnb7B5K9O6fo Date: Wed, 14 Aug 2024 19:02:04 -0600 Subject: [PATCH 16/26] making fixes --- desc/optimize/_constraint_wrappers.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index 7fd0d991b1..d6426059c9 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -281,11 +281,9 @@ def hess(self, x_reduced, constants=None): x = self.recover(x_reduced) df = self._objective.hess(x, constants) return ( - self._Z.T - * (1 / self._D[None, self._unfixed_idx]) - @ df - * self._D[None, self._unfixed_idx] - @ self._Z # TODO: replace with self._unfixed_idx_mat + (self._Z.T * (1 / self._D)[None, self._unfixed_idx]) + @ df[self._unfixed_idx, :][:, self._unfixed_idx] + @ (self._Z * self._D[self._unfixed_idx, None]) ) def _jac(self, x_reduced, constants=None, op="scaled"): @@ -628,11 +626,9 @@ def build(self, use_jit=None, verbose=1): # noqa: C901 self._unfixed_idx_mat = jnp.split( self._unfixed_idx_mat, np.cumsum([t.dim_x for t in self.things]), axis=-1 ) - self._unfixed_idx_mat[self._eq_idx] = ( - self._unfixed_idx_mat[self._eq_idx][:, self._unfixed_idx] - * self._D[None, self._unfixed_idx] - @ self._Z - ) + self._unfixed_idx_mat[self._eq_idx] = self._unfixed_idx_mat[self._eq_idx][ + :, self._unfixed_idx + ] @ (self._Z * self._D[self._unfixed_idx, None]) self._unfixed_idx_mat = np.concatenate( [np.atleast_2d(foo) for foo in self._unfixed_idx_mat], axis=-1 ) @@ -1025,7 +1021,8 @@ def jvp_unscaled(self, v, x, constants=None): @functools.partial(jit, static_argnames=("self", "op")) def _jvp_f(self, xf, dc, constants, op): Fx = getattr(self._constraint, "jac_" + op)(xf, constants) - Fx_reduced = Fx * self._D[None, self._unfixed_idx] @ self._Z + # TODO: replace with self._unfixed_idx_mat? + Fx_reduced = Fx @ jnp.diag(self._D)[:, self._unfixed_idx] @ self._Z Fc = Fx @ (self._dxdc @ dc) Fxh = Fx_reduced cutoff = jnp.finfo(Fxh.dtype).eps * max(Fxh.shape) From f846496b13271d6bb1a5175e4a2e19e805833dff Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 14 Aug 2024 19:25:12 -0600 Subject: [PATCH 17/26] remove changes to HELIOTRON example --- desc/examples/HELIOTRON_output.h5 | Bin 324656 -> 338380 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/desc/examples/HELIOTRON_output.h5 b/desc/examples/HELIOTRON_output.h5 index a1cb51c45686a70cc32a8c68d7c9169f261c1c11..aee9e5bab248ddacb4e54638555ecef8a73b1c34 100644 GIT binary patch delta 26871 zcma&N1z1$ww=m9`p}Rr4yHP@7Ktc>crBgtZl9Y}Cln^8b7?3W5kWvODhHeBzx($%- zu5X68-|yc4|9jtio`>0M@4fo0z4l)FJe2&U=pLYq;Q@E^VwLjC;WAKSTs$H|qM|~g z{MRLfo>@pdc_MK`>V^gE4X*J&BM{=*|1(1Bzaa?V{|kp2eBXZ-VN7uUpAotV{zZf_ z;S1J(@&npKtT0h38W<^&$iHH8z{v1uU;#vG|E$$A(Z6WrOZ+c(EffD62a?ER_E) z)Zkob38nr&3l*UHXaB-dsh<3+R!`1lDNz5rR@SgUVj5TxKR$yJU+&dd0t- zj8+(z7A=Mp_hJkS)c}{%um>zWm&4);us{|%XtlWF9n(0amWHz1`+v5K*7n4&>m&*$ z@{FAIB_K}vQ(qsEQ2LgvJZYck6-M5j&C%FD#a`0HH?Vqs1owe{E#vSAl2Zye_@!BZ zT(Pr$VDiD7rVeSj`IY#OlKC$8iB(?RV#Z|KW70Z{2l3PX<;~J`#l$$H;@aEPua;*@ zo3kMU!^^gHhh@)!rHis#QF7%|vB zTJ?IGrf8MBM6@mK?ZO?4X);9ovi+oGppnqv=aqg+B{iK*XC1qQ#J7-KJO1TM-RjFG zv6AKxu+lT9{AL50UfE=f_f6(U*K+kukMkueo&^L_8MNuqIWO<6#^U!~CjsR<;(VM& zShte(uRlsXx>lEK@2cK6(i`LC*%`P5Sy^SAC9j759;lmKV*1M)0sUerm04JEyuSza zD?5=FmWa<7X{*ke71_tlTqj0x_)MTca=;zWT7z@!~LAGVEKdh?pWfVFd;`i+fP{&_G?>d)-{7{|GvC zkWs>R*_r@K*jM&y)UP`S`d5sT2PV!%0~gu{NMK%^9H^mm&R3{kb}pCy8<=Taw5Xv9 zEzrJOL27ywzbzzFjuV10sOc%Enu zs0X+FNB{@G3t&0#Vc^2~vw$JkhVUpTjHO(lg%|*0%Ev|ZK;R@647i{qXuO>C1}*^n zjV+}CGRs9Jxbc8b@T^-P+Hz+I1P_oXFI8(G057m%z|}{ArgABrH~=^(@6yEqKzQuu z3_)a=+0leN2>fC_s??3NSIopSyb5&$r{Xspf~ z@N>>fpNkZ*0DxQhgThdJAiR9EVgnzXgqk~daZDph>kpS~Dlh2kFX-DjQ7-ucpgh5e zRyxrf@F_3vFady(a2a3*8r|_zoLSLew`Tv>lBm%Wpb5c3ZY)m2NYX54$Sd#{q64{ zpro9CPY(pDzrY9gp=E}G{=g1@>chs4Xkq5I%(3;FaFN#_5qLWmKw6%3%!Z4eh0|dH zSJQD4v;a9e0+$dSSbz;K_=F1;1Sz1GAVei8@Eruk;>LhSQi7Th4=4eYEmYhWvBJrX zflwbnhX@KV`GF1*?u#IoL(tO!D05U0{R{xg9xd(*dmbZXcmPRkuoH?Q+S>>zVF2as zfD4A8m5CJLx3jiYJL1e1}>-eCU zNAUZjKr=%BAs_{iQKNWXB?mziObO4@27&=fIJXX<1W+REbO1(x>OXf=MHJR~sKOI; z0e-v-ox~N$;lKev1E0|aaG-kv)r8pW9102oNl(q5uOYjczMRRlw;-c_T`P{?uYfSP zj^?tS%YkiVV<$Ae&h@pvV}d3E>ozm&m(? zMYig(FOg1s3qI_I!*v(|g1+DRF(4Q-MbtX|DIhfaIk_W`fsuq}D33ukW~iPmi!{hff&u_aP^C#J3&4O9 z81h6YehrDC-!y*cC4`z+7?%8jHh}iQC=w4+_8oykr9p-? z=)e0YPlXbRk-6){l*Ww2pcNPxfQpepFj&CGCY88IjCwH&>j#VSPg<&24TQl+3=ApN zZj)-1LLYABe6Os-n91PY`Lc=XW5x%~ch<4_mA^-Az19HybSnR0NMFD@q~>{lwW>)4 za%9hR^92SR!uD=;{fh9pcMX&RxmmylEs!FZS`u~f5J_Zzx2J|-QRf^4K} zEpmOF3Bh3Ss8`o!gSaC;;#`wdsBXggIyh)(4( zT&;X0Z0I0R+0~QZ5gzVH7Hht#z^XTp`j6yjMkw#m2=qy>)fG`Tv}A55LMC7IJ>P$& z1vzoVx413$PQE_As@AiK1sU+knM9wN55hujXk1xfhr|frf&?~vuY2!6ojUx8AF1E9 zi`PPb9~l6l$q?nDfnYc#zbMPxkjIF2*yP)ifI=`b5yAr0*pN>t@hh-Rb_m7~KTd8N zd{q2Ia`cfGkjg(ed+b!;n4DSzH!JGN|NFf`7Fq<6l*#+ zTPVbnl#N8gzY%$y$V4>PD1cNoulno3Xp2lJ*W-81=Y(*1CDUcry+DREh(c$-1lK`o zA_&=uIx!$jv)iQt_ z;#Xuz!AZnqfQS6-@4Al_6V%87PnIy>OdTYA6<^zsA{bfo?x}UfqdNK2UyGeD8h7Ne zxyI5w0W~Dgo#YedJOyN=KlDxTPJ13iEPv$u<5p#Z^Sey7vK!F9vQy zCg(DSPexslEo-Ir8=Hbip2VrUJN3E{>)|veVM{xRP>*#x*iQm^P4of+^u|jetJ@Im`*m z0-X5>DDTn50>C9GK(K`f09RfC^;baqBLH?XJ6>T9i~_(T#FzzoOF&H$<)P5AwE>vl zP!k;25#R-k;18t0Ob9JU;1?KlF@e~31z-T6jd){w zsKN7et@i2yb0Qc(osYCH^$EHH!mofh6lat^i1l!^LfVMVH)w@Wv*!x+UCB+ll>7Ev zE+RDr)q->YnBZPDTtN(si*QOqVW37X=vb}*5gh>Gw~W~X25X^)&YklZ1)!Y+verQm zzDM)ZK`?wkZ5tpqaA`p;08@C5n^*vlfw;i2MvQQg0w4v2(Om$b_Vz0tNodh6ljsvD zIP*$Y32KliZ~XiWTCa0)Ji0uCwqJGs>kJC+yHaCB2QFIyOyd35feOKf0Gt5Oj1F8t z5;TEuszUFFe{~RXUjcbMP=OxedlMR9GC+{`08)gYSVIKL<+>nnmf?kopi)Cb%_Q0~ zLqzu!+Eb{BbCWv^5fjtEXAoY$;e~z}@P-+Hrd*`!CO+81=+b-pYv*Brs9!~Gxp1g4 zz(Vw&JL(HrYD{y!>+!!6?E`I-M-L71Yk|AX2rwv-;uvJ2qDvQWsV=f4*hPf;V3aKs^S2ZaqB0a2(X<_J*!Z^k;qXO93~=#`!> zC>m_k-&Nlg1!Mk4|Hc34yZw*8`~T>B{Exop|LA+6`Y3sH7*4$C66oS6gBEzz2f_P} z0cPk`7mbAW%vB$adgTN3+|@Yxj^RK1Xqzscb=2mKe#FmfHZ(+t{y{w#3SeNMbI=tZ z92AU_hyPy#xNr}5FgNejJo>za=6%)2N5TI|PXJ$A1NeX=1l>9y0U{vGxDfI|2!vXD z1ilRy-P9)Vd69v;j#%DAM>_5LWu%81A(nrEP!WzAIL~+@mjFm=9pMXrP@%4{e(3eK50B25W>>RbolX{DIWz~Oneg9ng;@POlU7$#ziizp!lDS!zWeqHDx`v;v8 z2x>6VZpFZuz{WgR>ueqIh7yEQi)7<6$1tOSb9K0{0O3mjB)NfLphxQk+4%qJHHE@L zCqB^X211q>l!yr)+eNKKaQS6YA_aG|VZrQh7y$uzLm?PdPiVndz=q({4K!92%pnIK zp(O%BovVpAFBdZt3OLVuLMY%|Cy6US{t^JGZX!^r<`@8JY$EuiK_U3yCz}AG3kx{G zL6qDB>4WgDb6%LN4KIBLVkY<@f6?7mxbT3>=0(=G)9yyB6L9}0H|yev8V^q#|O=9BHo~rJZSSG2hlzO#efOo4KL(3Us)LFs`fiC4Z%GqtcYXbi^^~__P=JYDnszd*7>`6 zh6`HULJ+2aYA^}ndoTUdf3C&#g^RBC(uIFEFa5)lg;oQ!eHjttDB#>LbXNe|6~MO* zAQx?K1M2;N#kK+X+7Up7sQ(DMg%6fP;hj55efz>uBVW+3>D9r7+}tC~*wP#D z_2-zRh{`IIOpwVo;(j%Hi9AP9o`<8~m8jq=Ao3CbC2S)m8bCuJa2jglTujlWn6Os# z`=DzRMcDvDZ@}n6vt|2YA>pD6O~gzGs2_y76i~!ug7>%u?&~|a}T-qgx0?wTZxdLuq0wB#Dbj*JSf%JARq8%K9nmRWyaR-1gt~VoQw?Mp@ zpm#e6!hJN?ipzOW>&}IEwPUpSp&f+n2}l?h#0W;%@c~#MP$enK7*sldpltvx?ZAz3 z!7Tv6xrGXNU_5{j4Z_yNa6w^2=c6&`QQ`}#jC0G6P{Ph%4Z-fE2^3cV{Vsr5p+?`f z>>`$#!1@4)42+1<1>FSTG4mi65dZF_L53)bbAxQrc?)%v>=Oh(K;0ao--dz^N+QF# z9s#zvKqj0O51qk0&xLzm5r$s@Z&4gBk9T;e>n=2hm^1)&(h~#Q0^s~Bw69k{CrU@u zQO^|&w;y(nK!^|CLe<{n-i- zKbLm>e@S!WT$;qWG)eTQDHscB!1=0?KVLQQK1whyXdkD%|Mx1+`4U}3t+ey}aEJoV zQyggT0s!If!BHjD1iViUOe!vrb#Yn}xzt4e-bMdv8O7IMEXDQ9d2J@F^STS2xybk6 z0t(<``u{j_onr-Do|RI$uwqLYVVOL*pg%9qQ%T7%E<8pN%mXguKA$4Z#f8r(g0sPF zDCN%8727-CE8uU#z|v@q&-2q=)ZkUIW_Ve8T^4H=SKa?xtl3@~3U)z}O^SB$JwCip z=u$0=|2V0h>+5uB+le(=Nz@?aK9H^(M z3f$}Z;P_D+WKJzPm|tTqc_O!YDl~@B>h8zCY9UW-wbVG1`DcH+JZf!T;d=6>dQX)9 zhvJjFq4#g|JxUeZmYF^>nLZjnTmCZaJ+vijx3=dkw`xfYk^LY&JeAItAr zOh^s4`RP#p^0>^RA}nd4`It(2#IzN95dKC=WFy_aM@_G#xiC`umMC6|>4t@weSD+x z-B(uje{>o4-TwA)Ly0DqTPKwT(^VQ5CU!+vCX(7IU8nR|FkZ)6`*?Ibk1?v)W}D{X z!r!BDnYd>JQnt3LRq2+^p6*6Eshn{qRkGFuJ3He%c`sI5m`iK@w|_t=%R)j8nj4Te z65

9rU!sYK4v^@J&Io+=n~H15Yk0|MN`$mrszvcpscy=?)i;L09OM7bry@nYR% zZmSXD9`g{d%n+3P_=oA>Sos-!7{uD>or6^9+(d%>WOBkw<^5{avI#K%K;%VZ=%$ zewMK(SDxS4+B^A+qa0DB@d1yBl=rFql3WIX*!s7-X5LImr``|N8WFbIlLNN~ zWN!hvcZm7h#G@W37pXVTG?>;*Keko0E#~j`&xGEbh#(Ql^ZivvyAu4`?j%I;d-ZHb zRJDHs-FwlYfrm_HX5cC3@7^24*d+@;fME_%E~msT3OQ#>8ZIA4-ckqX^`@pQ#IZ)+ zVq?2ZnZu@GV~M3}sCrnMyu5+Lk$mPod~W8lhc_s80)BkR)2)^43cGF6c)Q9@GWYJu z#*t?s)AFir(#hCYBG~4*A3mFm!qWY-(X!!i{jnGKlW0wz*7{^c#8A;lc*rzZ*7lN- zm!>q;+*t|B_irgMcu8UC0L51k9#RoyJ^ktb?i0@o`BdI^%KpGW-`0F!Wha2vlb!_` z|HhpkE6+&#uEN^SGBF%EIZqj{&?Rsj@oASD*RK;!xBJJAdDA;Dd^B-h&8lZ%73nCy zc$2TUV^zNGQqr7Tz)V)X)06+O%K63`({u8pQaO^Qpa>se@do@PC^u+!JAIq3()}|zaSHs2H$)qVsie$;m9)nCV67W?11G9 zf8aIgVI17=A4yVo97zY|$By^!GvmiEb6uOH1C)~zdzAGgaUu8Kp7j)xl{jI^-oFVX z`&mEZ%!_P^q?jjp_MTO2d~G#sDzdShGU5(vB+o?G^1-#PTc0=6E1+KS$}dMDFZGDZ zzge!>ev!wYm*ZKVHH}qJ6L_Oz7aKiI7kKAu{6s*V+J~fu$B?%R_?F+b$71a`xK2mt z<8Lm5(;i)$&w$)$MdH%97!`&nd7JDisb^UW9O?e@Cb1FE)!K%;D+=2SFWO~(Sc#mE zzu74T9$@9wckWV_;XZ}>4)T_d+jg7! zhlr2c+x}POjJ`j_be*Ofi5=CQ>?_;dX~=)L>OQB>v0gYfPfJs=5V4+an5$WlNHS*n zZhV_GY4a9Tgs+(>%Kg`Fs>0^Ji5HpokLI+kpzLy~t?e=n&;wkq zV*iEi+`}{E)UK2ae{(eWCHtS@qPk7X-7kUVFa!LO$IWp@b9qefCV1ac8f8BoO^0sg zRE%G!N5j#=q|z|MU>OPRe_`Q^%>S_5`1~ZQ1FL521*UH|a@- zcXboYKw)T*v7ol3eX+cug8v@qDgQ@AH}eygwD16sBfboW{V=r0G`DOnnve!2vqYjW zne>H{#|jj55CGXAi||24`##3?D3P-9FIj^7QhuBEo_8jTbfeDSPu?Xb)y6rRm8RBQ zcsuaO&*bkgoM?LdGiCd3da>sq^d#u)&@16rXaCb0|H3g3;XFbY2XWW1WS)A}V7pg+ zZQ25gImLHZn){&E8oXvTS;TOoCu(lOBTdrNH+lURscMb|1e0zHF%sua;Agye6m;Ex zl%fPds)yw$NfzA|8cI{y8HB^>G+rR3IIJRW{5S3WJwJDX@sRmDn%dTib@P9OEOoi| zc?Tc7Rn-k?){UH(FlN${w3rC(7H%yH|^l5|E`ZM4N+}0lP{~WAzBbA96x%BH=!MEqh@9dLAv}!-n+eA}{uR%9wu7~8r zWM|F(!Px2!Z+^27GO`ualY>3@yyU~fth?`CV8rw)M@ZR8nVX1xEr2`;!jXw<&Rx}? z+47%ezjZ&Jl@qG5q;bs4CB-$-AaI=l>y`3vDjI@&BV3zzU0Y9hrYx9#yCzLjPCuug zNFB_Ks`=Rw@~gka?iw4>jG0{Z;%$HD_@>!uW+a)KPOG#V;FzQPhn4lwqn%nG zvOIlAbynV2T#aJ{F>oQ+>Kt2eouV(Wb<^`eO%nOEj3iY{rQc|lW|&naKYiwGYfTq#xn#HB&d`0(Z6j*)@d$^;6xT$5pllR_%e;Tk9%A6U-C;mv zpyDZ3xoSCVC2IV5b!L&W5>K;>TCPD-S0T`=e)lP{uJPE$pD&gi_hqf=8moQ<{8jg_ zp_7f|u{-u{h3?(cpU~85KT~1g_KeyfChTT!m{=5glzoj1&!I8iJ8>;J0{d?4R@i=- zG_gAr;|aT8f>Gql!lOPqCebvE&pxHDI8|KzUB9~~0-5#q?lN@7BxaY*1xP;{^GR|e zrkCKn8K`20l`z!1NupeNi@>7`y1O(mT&^|tZo2hY-39ts@mniasn2zWgJ-8}*R{gR z?RjcfEB0pj^-eaP-&&Abcz3Frv`FUr&iAN0d4*R?T<@EPFfZf$J|o=x7hXdF8Rd-| zf25Y7r!milc$s}w`zv?vy=JignNMG5ME?qgq`0H1*E>_s)Zmw~|I&O>&?Y@yu(dD6 zEB1ki8?j^t)ag;DG(ASS%%k_|y|FW12?0evsc@VkEmR_nk)rr9T*mLzf zON3mZ&r-71Bm0Pj(wnueIwtMbHZ5o0H`1ic`VccZnbK)lr_ix-HC0}8^e6$XO4_4U-}_EUur1y91diuKARCb z*kP`{J4?S^x95HPQN9Z6^NGgPCQkxlc9N;87j__0@zG9*80NZDl4mNj5@?Zl*KDim zThyL40y!1O@rJ&e)6{HIoP_V|kc%laA&|zG`<4lm$3;I&T%H&kzy9soIMT#6LnXC$ z*CbNnsk7v|sz&6Zk1`y(=adBl@I8rNaPT@?gg>co=s2(;9g)b^N!g~8c0SoTY#U{^ z9vPKjv8`d$lA+gOn$HE;Qn<4mG8u!-`UECJAz7B(y95+&<@G~ezvG`+(BRTP`B8L3 zFgq0bQjR`KPynewvZB10L008>uyUsWfpP8cwFO;yEKCFGF(({Jw)UDg@lh`IFx$1s zg*}-OlNUnGT#q4jPJ0c^Fq`>xryOaY;ixoEw!RzLQyvh)#hpk!wHk~zs2yEd-D6yW zGTx#shfgB;ch%I(I*9z)80AN-975IK6a1WoT6xg!CR)yGKEIa=lcSY&S0cifEgGim z|J+QHHuhyjPQyIq5U12xw?bsr-obah5HFtJtXRUAGK_@V^-k8yY^{4?>!xE_geDO? z49(LYe>G0CS-9E^?Qb$(cUBmphT|9`xems4Kd>!IiLyKkNeO?;#(7WvrO6^T#aODC z6x8!>l|)|1tLbo>G2*F1EXlCJ@SW62@JIVG6FT0tzWWYt`{4m(goDU_qkE46NQkv* zW!@KKu%AJCka-t8?Q=5tc_r;*vqJkd1`Nz zvh4pYrV00~oCtN!%Er<<5l(BpP0{dP@&}Y#{F+aeQWSgfvus@B&xu1fTZALey6q54 zW9$J3anJ8q346|%1@Gx_4AP$t?!5DXI_~Wq1NdvApKIpt^xa!+j_mLCO)_~X{Nbi8 zadX0)0e!DKn{7SBO)K7KIYvmgL55=Z$1^aN^Fez0{%}!{)=!1!iQ_C}d|3_)>i@=zflb{>_iysN3lkYGzlKdB>DkBtNpv z-=3${XV?`NJ~F~gA+-MZ9PepRf%B2pQzL^pa_^(+#$I357g& zmGo7kuo?e0)4p5g&0=jzZQ*uFTbdCvdjjo-tC)@q(f3(e?9lzi{(gw#UZa?br9CJf z2f?Tea`f^Wz?V-$oc=!>a{;4Lf$KMPLly7B6TS<1i6gwYC-MTiZb1w#XJ zYd1nCw-)=4N|*bOXbxtQSFspZK5%lZZFu!oYf-mN5wtuQpR2#upX9&0gQYygP)hs( zi8rNMV34pnLKsT-E43n{@OD%}i~mi(ni1W%Wz%6*gk=!su(){rTh}{GH2ifY&kkdJ zGAbaCkEd{>?iNZHR~2W@ymz-Ro1TLPDrb|UM(4x312Fm^fqI^*r6ZC zCqzz(>9^2~m3LZRLVI6Ubl3xO;DQMPr}*k^bKEZ1c3BRbTL+ey-;MpgJ{vk) zHQ&8IScRt`&uzHG(-)q&vD23j04U!z%#$x~tFKrxbXVM_SJ;>`R`P{LhAA;jPaQo8 zhz+~>-CnXw5|h~B&4w~Yw(%6gMV^_8W0={LaUX+cjai3-*g|>IDs9Phm!W5J>Oma8 z5ec@kY2ESDwoMw*iTB+7uG^g*qU;YxKH|s7r-@0VDF;Y#VJj`bffkixD-BO_w&6AM z*F1QN&Cx~Ez8W}zXW;D)Vklo$u zZ1ZSWdqi@rSd%5E4FAEwN&57j&#BCO^wQz4nyKL(u6vSn%UOgk<>uSoeO*SFk2~HI zr7nH5o;+&D{kt=u1Fz7aD0Zp+Q+QMDLrdY*9q3T(Oo}aNFr&qw^UG;~j~&1bR?RJC z)r>-Pc1nCm`pzYz+XCc>VS4EF?KT06GiKh6JMFv9Q^d(~3?drkyZjqk)EYc4f$d${ zTAwZ2rMq5JY3zJY^sU5NFRe%ol`tkxI%r6XCD9ylp6@q*z>KF+bJo5^vU)o1ECFpi zfbwd(`!;R%?F&W{X|K9UhOs{*-iIM#Ox{bcC(TGOdZ$S_WbjjoKbi@4l+hlfU3Caf z@KD<^-w?RGtl_=W-W5s>`2-YS+dlg7jb_lNHP&gylB&A8HZ1=bDsR0kZx^^YCUEoZ*-I$g zN{@yo-fm!eVHC65CCy{7MVKY~dDB;cH3~j$s=1NxsbJ*^NkbB|G*$TR2hAUyG1yX)L)@L^{InBV*}oF}-sSs)BjA72(Kp?EzKQ?aB7ifWqE^8D^6=@4)2U z`BdJ3l~%6mTKeN@js8B+8}u-`EM(e>wWC$djH{ARpV`j1@Q-t z1m!ZJr{FTB!DYpLP~UVZo$bjt>_}JPBfh^k0zMpxeAHVY!u~43um?FNl(%3tYJzsZ zsONV}>)3wEpnu(4BsD8PVlv!eMWmzZO)o7xMb`CuYq;<3v@-UK*&HJO9sEZgxYIaJ z_Gd1H%+y@{1zv@-7AG$0Utbcv5y}@4#MP~*`GaZoC|$qtRetcAP$7-eS&Q-?`d4D_ zr0vje;>3}NeIJa9!TP{gsi5=nkx6&Fe=`p!&&cm<)dkxcwg&v zoM0on3vGqzlwS6l69cvW5qC8_?6FD6jk8ZNe(Sn~N^d7-I1G2M^+^)pl$egyuR#}@ zEIkVcdke;m$#m0S`<1L8vb^DmfA)3e+v}zlxhh%%+xicNepYV;uFa*^)rgGv5t`Z* zhHP`F>ocaQ9r*+Ru!s9&-yFLGf36wXFMoT={-s@O?zcPnXD*-g$f|7Il=smF&&tH= zkJ1iPd#e%V3;c;0JjMh8dhVGSGuFpY;o4sRDe=^rwLW>d#^u7DaOnn(-~HePOclvAtIe$VBj$&Wv8 zqWv>sK^e2!12X=}rCF$N6>2cgpUqPz+V0S!sBBe4V#Y@FXXVN7n+&O#Cf<#{5iq*- z)8l57*a_as!|xjG8S9?bHZ%{rrM&m==XRi>w{GQFl^=Ye2^uR zkE|W+#TBM-bgEkBr*-m>muu`Ks_qnZIHGW)!nudF|51lW(XC4cMs8adlYd=Cam)uy zyk)3TdvK?rO!FO~tWV)%%`98!+dHpCggy;qrr&%TZaMX9CzsBM!7lbwz2Cg-#L~j= z>_k7x;;syDQDfD5$2qcy<6$|c^}dA3hQg;#Vs~FHlNJ$@R2axf*A8I9r{~BvMB`os zN~bd>c0cD)6g8Y1$Y+1ou1_7FpFGn=Tt?2UtugzKBXSP>m9J0XQEyPG9rRl8bcvIV z@f2{g46O@s9ebd z4$Xw%M?`IqO=|&ZUd;aLm|&C%fyw?IX+%~ zHj{1G!ZLeY&q-2mX;XE%nUv7m**%qW?WYP+-ca%xpYsTK=yF4}{4?DOH0CJnPYq8< zF(b_vLr;-t-y-VhWS{;(yTu)ck0&pbyWbejuBZq6CH8^;rNawL{LGm82=;hXP3o@B zsm#oZsv^fJ*^^ROQ2EYRKr!vx-4m1O50uV~6N5m#ucBLPIfBLb=Q_xy@)~KyTd}&R z1eU&%Kjq60@vYOo{}d)l3W2gug5}d`ou>$b>-cY^cbm>oHBVdOqC-0 zWaIapo>0dmKDva%S2HO!y~lB_{rma(PI|etME)xe>T-X1F-S;qwn@JoZ%P^;W-$@Q z+cT|U4Jq#W#oWR;j?fiidFI6%wDwvku~wZ~4MG!nryQN?7?QZ{P zi=gC{DxI9IBDUURF_X8wb0_3lYE-x%*z6jKDhA9crfzyptHJid)9`-$)wC%Xey6ym z_Z|u<;8i`_W)zs-Qv^XhG5rmT3vLDKZD&^yaF^y4u4#U{`)CBdKCGB28n(|suy4_b z9pj*bnf&oCHGhKYMA)JMXRHu!%lZ=>2Ui{Nbfo98-?J zu$@2d$GV&~72dz8o^ii=u@X9kCmj9j8+4~(Y{d6YkJ*x1fcGqudK9y6x>~9GqluSZ zImx3>J+yuucFFuvD z_iY?_A=%7Gz`UGF8Qs@LLqctOli!2j-Waj)p`mPG*V9#;lt|5b){iCH>Y5bxvpMfp zAI(S6yaly=l6QJ3uWI#&E96NbN4jrN?J3- zuKg;030d=G=xZ334QpOug3m~Ye;edEAbf3_jAx*Y{!63jWO5FUng*QE8U1 zDl0ll+EpOUDdi_ka@I5+7WbXK4I5<2zb6t8;e5?L${X%xOlb{G6X;oPe;4A2%k-_n zbsVL5K{T}>R-J!NiDGDE=De-%kw|qO=di0DUOV+k%Ss&^7I(Yh#7Y&EVc;v&+4E*k zPSV2!+w#okC35knaFrYm-OfE%VU_Prt*?tH5raXh=3ico%Gv#8Fw;^!dI0Q~J`-^f zOC{Pmr4kykW5a{u1+wu{_>jxiWs-I>mIN&Z$v@ccMSux6;dJM*gcnvMAWk1{&ZT^)BH54ZzZE9S!e3zaYYa(njI^1s_rTtmn#HpC~iMY>1 z$;slO@in^F)e<8G28wIP$fF&J@*nc43=)>k#N`C3cYh8@MMkGAj(sLzsCRiVj5Ljs zdb4N8ur%#;8wX0!Wq;>a#Bx&QSlL$78Wu}xN>}8t`pgeu(c;^8^qfZAzQL!ale4-z ztLG)xoa!^}x^)(g`Sybz{3r}V;AtDU^XL~F-+IO5EqSB;RRptRt>12`$hPSm-j4%c ze0R-HYZziUEvBrQ)9>1^82 zk^1;Z^@qpl%P&1-AadAvSz#nX!SreM*b@6)Qrkf%{kZ5`?~cMK9~W;<`x4NkXqui} zk4h|dU-D<7$NbUujn%zhE{sl_{nHCS$sYxW*-b!Zq(|>2G zc8)Nq)8TRcax1n~v|6R~UqA4xR(HDUV+OVfQ_AHui&tryjLJ&8TS^d7ih&gm<=L(} zsPa-XHo)0xy>>jMsC<(A^CoO8xE+6^@aZZNFn3l^l!Mxr3tPQ>+sN@&F(-3aT0y7%=|_Ku3_O+Udb&5eZE9I*ch4fu&ea)FFLkGlJ!=EBt8L=6Z&O*}lpx&%B z2-hnPDqsm@OPG1A7BMUS*zfZC2sex=a%=pttxCWae$Bt@kdrlA`)gW13%EV>I%Z4a z<}Y#E&p1#^hrd5<#}EGgZb8SilK zGB0qeC8chvZ{Keap<7z7Tzz{d$YzPCAw>8WhO>&*aCA1S)1K*mkIEG5lm64VI$Hkl zX&g<#Cs5PFNa9MOlL7IjvuR)U)nGmN#~Uo22_YL~8_WS>w(e4_PE@`O-%Yhhjc!+| z&n>u1*ww!&J%#R&v36UQXHG|GyqtsVeTri5TpsJ(Js|#Ue0JO3@*QKlh*jusOi>Fy z$LWc3!KxI$92sbsAVsWF3nfDYoNP3~N5*>McL+4@?IHJ9FxiWaQwo~DEPp+F5C=ck zGW#{+;T^{`i<0TJS(#tMh!HpcI5+jk(?1eb>zX(=i+Yh2|65k1sarnBzliF=^IxU} zUeVd>Be87I;odOPuIt6>rj!x{Piy8R&<-jks* z1_wB^8<_b-FDER^xkP1$m9g=sw41(WBYz!M$U&B{9BYzl`IKsuUG@oOoIPF58pP>r zzv`3B_02}s*&fmli|~}EzaFm_-tZ$*=^+=Z=JU`0*(OxN z6vX}Ba=hWrGa}+=gVXi57NCf?n^U-%TnG6Ev2A!{GHX~pWU3#$*WIrNyVBq#&~jFj z>NbQ9OiZ!r8f6;aD3(?VBZ62d5=8`;Nv%+%AU$?)Kj@g%J@tA?QJy7r2FO^XW{wHF zrYuQq3vfotQGI#IN1snDODk0ZPcKd3B_ca6Dkm%AtPUYe$CTVlV}Uk?2d~ju2@uCj z-G0DtBx6|kHOFq)C6S?Su7IyU$Gkfs(GX5)XeA$c>ztw!*wWXcRLJmus*G9E&DZdUCA`Hmh zi)t9osSex;mb{N&(>xM99^%{Q#LP%rG-``9P?Rw{di~cd3W(VfWm9qh(@O6{Su-YY z8J{^DzNhh5QJTkrvJtCtHH5pQxdn9LYJEK5t!drF^VFo8wuDF$%ZKdpmI(V|y7k8h>6drPkA<%W%x}*svCO9wf95Bn)6J%hIeryMDuP@wqBDrqYtd(sp1GCP-MCdc=IdH{j3tiZ!EQ20!J@UK2Et;k zlZ{S=58QvP7Ns4-;np`lD8a5jZu*$7k+0~Snm49MKD!o1vzmS^%fS+Td(78|@wSGd z24k75DZA*0>(8O%yFJm&W~t_P8RNkF-S64FQ=VHU{VrHbw>Ti6BB&qZBIdWgLw_h* z`+!(Rx14m-`w0s?yzPxzEa8;U7PoH`mp$brRBz=o#sEKJdpB5kLMK7_V2({8)LZv>Ewy;7~fow1tB7 z#V(akmeZ*`O}ByHq8+H9BxpExVu4;}`2w9O~rY~hl9ZQXXQ zrSD7oUASD{yxhm;Pz@h8VA?u)x_g>;0+ya~wG`>fd(|Kl^4Y}!5;$nOwIHsfZk@9T zs`BLlQo%NZW6P)i3BeP73HxdX~UaH84It*N?P;QX48?5DPyhe73HVP-|^{$Rx7uSgNViR zAE>kvvOqIm=wu)M>ZA@D7S=f~H38vab@1()+-|>tEvWx)C2+rP3%kF^EH*rIAlKMb z+$loW(d_Rv#p!SAg9IwivvD_vnZD1DYfaItJa2i`R$|T_oNbXq88SfGO%RW54qxQ< zYyETJx%7a*kA6W$c02u9R=EsOwC-Q*BFUkdIUg+O=IyR3sYmUqb>5?l&+Z)kr73*% zEC0E9LmcVteI>S-E{O(Q=Hw%b*Jh4c%ZU`R=??{#8R49m77)lqWH&qyPukoeMU}e<-5iY#;o+! zR~xP)tql#N{-_eIZXiA5*N|NDp;_Qq`~HdRoe<$U0Wxn<1)J~`2N+dw!d64atVhP+ z*^QtIPn^-TQxS&I{PyeJ+oU)_-Khn45=SarvZg}=Y3`*Z(b-RAHc8)$!b)}7QVV~0 zmTP!KrN5h@lx$suui-7M8vjY*mEOm(_b=?*ol)zL-e{`e`StxoM?Z%xUfM7F)iidN zfoTcd*`Eu+vsguVW{z@#mCm1SM%x9$qSFuMozJ|&n(xkWOv^1{Q3YlE`K@gBQsh%k zDW#bB@uQ@r(7m@@zMHqSKeFDw=KZ=@o0nRE{cR-|wriRqNB!N8tE%l_?w;u*#6rGO z$SV_^Q7E;Y@#omuo6GzGCnlrQ9jVfod=$SJe_q>C`x$TdwzbMVo@tG`F2=bbLR%$s7;jetmW7-7OD|!K9em z94c6p;0V-ldu=}?K2YD5(sZv%hk_(_ca$(@vxyiw_9B^%-(R2fc`Dqd;x*->0N*>o zrww*WT!G&$Z0>(g<08I|nVAJmXWmDit*A}d@CA@?K`O0_nR;A_3G+smK=@COj#BrU zlEwmMDBFErd;20;sn^C&i(5-l(*v_8U+7MIo7GFvQyWL+*NSy*>#p4w~((jrxX5s!nyWZ=@!)6;zsa8il1q7N zztw$~w0T@Z{6u2zIkD#501pLE#L5nSUrE-!_bdliO=-Egj9DkqpyQiQuEr5TaYAuZ z-}v+Tn#xYJmIgI{oBlqmcTUC@F74@#qUoxi;(zFDSd*WGNzb%#`07^V9Q~lUnSH%t zw$souGx>yedu-{hquXQ>-|*K0;0s)tvEMG4O>X=)Sg3fkBBbp&nLW_ zM~uGwF02(4V=CdfN>2udbDdorKTwl!NA-)j(M_jh7^pk$KC?zkhl<^N0#$#7j-8*@+9 z!*AZ_y3Rv0lTDSmS z6Wat_jq)|oX>Fy>LX&D`#ik^qw+`>N?B)6e1P@etX66(dtdHRj(g`Jy*O9wf*3eB#VSw4u>zNt=LS?SmzLYF_UaCo1>TXqD&>t3L$|gW?}|+-$e8Qe9$Gd6&5J?z);j zbYH2k??XXry@tf%4Q;AM>Q)!89viWKO#jZT8LvwxjJWT-mN&1zJ9(k-$+&+S8@Bj9 z9wWSw(ld<6Zock%xuVg&$VHI1sk36#48k;2UTxCGCgNP)w~$R_3fFkFG=(Ll($j)&VQHVxWm4C zz_pe;(keak*M9KhV}#0W2d~Yo+vdx5WbBB!(Ohp`e6H#O*XihJ_;n)b>JFC``9FDf zZf2cYg55D&x@zvv>S$Yddx`Fgxhcb3eI=$C9ydI0{^cBHSrw~P7OmWLr&sZsm6YQ; zV)OXLV_MJmJUMN9b$iZP)~iOh?)|DoED4tlHg;#)*}pcOmS3VW@qxE->EHUZ)Y|Rq z{oceZ4f#4b)3f2~8S6gg_Q(`R-WVYlD@B4{n(>tmd?j0Q)~4g)$2$6R~kN# zm3`VDw_>r^t6AFxIo;og>7z8i_~j;dXp>91rG8mg4_t6H(Ri0Q_k-{DQ}d~Eg%quv zmUXAo4z`r;I(&Df+zxIg`eN)KEMBXxt! z4>{lM?r0Tw$L+}XwO8+v^`7WukRx2Q#Mj$5jGUWrwd%C#(t8C(ZcRcv*=)IimeUym zxkK9$mS^W^3ga)I@0RZ?8M!H3Wvg9 zeXiQ{LbCo&ih4+w&Etg5?^>7I?EHE49#Pw#DfLL(`VBW6ZY}Rx*Yl%Bp`)|yXu;fK z!;Z3<2P7VuNEiraWF={o90n(XHe7V^dXrPseR~^!)5XHY4k^ifDjjNttG@qQWV9k{ z;gch#9tTEe^H!;=b<1=J5+o|@CT%t-v>REFra_o4@*hbebCVC1suE%U{xW1H2s{o;=3t8f%oBx9Fx@ z;$hH0S%990uOCLd3CIfE@b2sJo(Ww&wY#Uc^lJ*-Om58c=Dj(Tbu4nGLie&FxtquK zw2!R=I_(#JPKtaVeDRiDqmWQ;`nsqlxo26)!QN@vw{~$4PV7*)}2nu~SB0X!dy9Po@KX0~s~!x4!B?)yrs_+{7PYdKwp!oVM*ODy>&= zx)la*-Ia0a&TjIdJnf#y^z*&JYv0GdpW)aYXds*Qg)zZa2N8kTY z3eb8`)s?%5RneUmkZUHGJuYat%)Jzqc`gMDuQ(rO>*TG`o%UCrOpw~rg3GTx4r%@A zP;@8qe#eK#lMeA7A;x+OLnE4vrWtm{Hgws?-Bmjo)xtNIwan{j;N=`TrE9V`?W_5Z zxeC>(R^zrtrk}ggymjC~N9WtyaV9-uUat3#`W1Xh|Fx8|i{SIkz%CUDvv*fo|5_08 zV{2l`)=ie)UaA`m5-SzU_APzXl(I5&-}Tr@4#pDR&${InXnPiA>Q?CmrId5d$KMm4 zt_a?mweBjv+x$VJpF&h$Nigq-gl~Y%-QLD^I-OBLPshGC80DwG@q*m-S3gImwwbY0 zV@mFxsrWT3Wa=UEMYh`6gm=2FRz-@t^G<5D#F)A5NInqi@2A758WW$Z`9d6&>fk(mG)FU$+c4jEWAM)IgsBvxdsmifDZrEA2 z@w{Jaw&`z^I<@~+#lL&^Xl{P#bE$C5s@X}=RL`ma&EO5eP2nb`rd^hE%q}hIXTQ5` z@nxF8CbW2xFxXpFJE*o;eY_0^Y~7REK2zOq$wa@;`O+DCo!iZ;DmObESZtj8Msm^a z;I7^;r=O`j4X!n{@x0&<@=quSyj08K9V)qCwYl@L_n2;foe5#RdxcgpcAWZEwz~XD zvRP*~vRwP0B_FzTq`dR3J%%U`FEY%`qxf=)_EK6 z6Mx-!*uG+;Cb9j`N=1o?KZB2(<+%Llo_@CGpLg6(dlhY(ie(Jd^jZ4V2S>P^j{h)a ztZe4M(m7*)a0~x}b*rO3wyo5bRaN5MRk;$JSk0e zdI1CFA$4zZXZRDmQ8img4=%wc>e6sk-dQIpqU}6}sOGru0)jEcQiR+kh+x$(U2+`5 z)aa!!1$LJiaDxH4u;dQ@`I0&t2=)SBsPvX89%NA#zt|d-ryrQ3&Ag9<&cw&#qFid9FBmhE>4G#B#h0J~+u|OK z9AW}j=VB9sLC_jvLL0W}*Y7ed%#=Dv?*bITt-(cO{a;ibZjEYBi8ms81>VPo z5WONEV2J3IREi*?S9%$3{S0M1nOIHa0lEaV8vjOtjK(2RoEf9-lVAtVxq#fET*kd6P^ zSoao{e&6MH8()uJ*3ey@_h}$ra@hUFhqE%C=hX!5DONORHHRz4js$K%9WD=tcpO|F z3bASl_7;WM`yPfU#Da$yq7aw7#t?-#s{z8n5SP6FJ;XkM6p03Sf2R9k(Ewm709Xf* z&}zY!P>c;fVu)gF(TX67@q#vb7@fj%;7W0s?J_2!YXyLRwSa))DiGus(%`=sYqg7F z?A%Mo*sKr5I9Aa*&wa*B|Bs*BeVgM7(vN5b0R%M zaa#tkJ($7}Gx|bI|MBpN_!?2|II4kW-k(T4DFeEL8E-hGJ1QIJP?^qbUCQ_&nC`(0 zL#K@RGA6_JDEPy?kIJsar$t5N;zXSS+lZ*F1B%^l;f<1U=Pg;RMoVH7Cw8ZVJMd+g z7e+IWokcCFhT0g2d2$~x9T4-@e8Lbi4|oz#?L#{jG7mVjA&f_v%_)(ltGrjh>Cj7{ zNw;JFFyLkiOWBn56A(an4Dba{KsFmBpxI~z8(^V;3__~07f_vYBo`_1fHg{Z7Dp#3 z;Ijfj>Twm|XbcL|=+ovM>DgG`)GT?DOO2`q98EE`8vHJ^r1X*4j_?sYDIitrry)Tg z5uJ&90)E>!05!A({Jw7bx7>I58+O2Au~f%Pup0ypC#j8f;4<-Ez#m=@qNSusZF(Wg zW81PqL6#I=H4GkGj!M=gxe1flQb0=%KYe{miXBSbae*GfWCMos@M%1hd7*-v0f&&u z=bvZ>I#Mci2DC+Mk!uk%MF|W-jo}DaT%124pUz%jG7;NG25k*7gxL$6Cm{$!3OtP< z%BO;Bz-l62l)CdFz{lfK2CNbR!Xuxr)eDxhfyo~v-}eVt#g^V)FZ#%(*gt?sP-zS7 zkI$=se0*Mka(z5EbR6hSzz}Vs&twb{YXU7WM6B^PhY+(yCD7(~Uk#TtzUs|~ooHXf z8DFbxu?^&FlO2Z0*ROLBMD!u4KlI6U(Ly&Iq7!h9_@uhbczJ{IGB=-}u1vzs&gTn; zlgrse)qK91D!H01z4jI3Y_bX|a<+_d_Jzn{_81~(wU%LsoRxEd zFdi<*013r`kI)v<6(YLwDGNQaj$^TwjsT;U5OI!f$7UX+^Tzlf!2v!&&Ylb1CpaDCXj?l!%#=h2lp8@9`+GaF2a50iRk6IAN7g7{#^=+Om#Y;_$lG| z**A#xcx3|^RJP?1SswHPQJk1XLntnH3pT@{jHVHX`07QEK;mHkcLbwzMEErD!%fL( zDZ*qArENvVQkhn44iVPCFS9~V^>}>OdE{(1TNldZ5{d)-Tze8`OJ6u`+!`83uk1(z^H*z-Z@_q<>@vK7XGhr;n2sR|O=u2= zC^SB|A&iebyK*TQ0f#dZQC$292-!BMMJ&};L`RPPJgQ|=Zo446EQHjFa54c1bn*8d zC?F*T(nn#86r`x&Bczh1m>L@A&kgXuESBwEJR0qy^F8c5vOSGO(m|Zd zB2lhBrJhHTy9v+x)WT@;3Run}`9GscZ7CItr@s@Zt(hTE(+g;gK4Nx4BS<0|#8oqs zu~Vt!Cb%n*2=IoW`!C5wL@0R<-=~Ri0ImR^w$Tva6iH-2vIr2DSbXnllq+#B`PHc? zS29H##iHC(nOtdx-7z70L;*9b#H&EL5=wPeh;pT>7C(iAfI>975%4J_1pFHD?Vf5cVICoMr^WZ*AV>A0;`0Mw5sMC<3Z|R^)4l;(`2cLe?K7>kQ z7NN=Fq%cGYnJSGS+J`I|i1M}2mq_MX_$S2T*TUEXPO4PjPk1ev@(`E~3P?_c-NbK~ zrm&4x#+aElq&!4mr)%bNF|y7)cGix7E(u3 zOfwMr&SI?)rF`Eov>eK*Gnb7B5K9O6foPjP!)G5AD!Cv`uKIF6}3UcMR#$`QwknyR@G;dG_2-7~?Y4n1)rX zRZRMLOM6yS&JIn(Oy(F z<(dXcGuUNXgGmphI;?bQAn7r#3maznr>=t=R&k75!-XA1gkA3*)8O8XHRd>Y8q?iF z8Yat4_lXUwIHqyKg^o8ick*2C(LfOfdo-N#Orwjlr?Ju_05SZD*~k-l<+^BtRW2k} zX;e}Kw1GcljL`+C-J}f0!8oo}lN968z2CGe+4g(n7w5K%m3M3E2K~BVy4%2jvaR2n zPR{#t^~xHr#O2c*@0Z*Z_f*aJE%9HCf!ir;T^sd-hC@m$>KbMx1vdpGshp4O&JTE6J{5b=7)-){6M-EwQms2O_$ z4)jRD4({4h$ESMSoLM^YNnV>TPrf>OzfrpW==;ZytiN~o!h@Iz_Pg&S8nPz@RO~yV zdpgq>7gj%dn7ieAyPL;;B5=iOigm$sZ(-J9;Q@nOIY@pvK#=_!Ce6_YS5+UVo8Nu_9{2qNVRH|Kr@Li8Z;w z4v$x-)aBhy`}5eDE1kxtd~{`RP1mXCuZJGgPXBs-x7h|73)>ePvb)8l)La^uUcJbnY1Hes2@`KzC#gUFn{XRU5gwXe4Wlch&}4eBvuK0Bescj2EJz8So7kQm*~@aMMPe>2*hV)(>W_FG>i%PwdE!C0 zc9ks-2>Mr>%O-^H>?}IX?WkUyJV7jFbB~07J5_&q`M`kwGhUbrF8FY5?tRcroOu6U z=@&g#>U*cPN{OClYv}t$xv!~>uej&$-(Ho?&N7=GWi|bI>NEYVh|SG^{p+FmT~nLi zH;j9zKRu{u!v6D@%-fwOAKu{glXvi=`FWh3zpkw_liU}tjb3p&zl!sV z&|1|dExn8${!@VR693313T2&DogiB0AOAM)9#*;S03&#EeFGM-Z37r0iQ%#!=s_;vlfi7B#1DyX{j)9 zw4lZo7a?4}O5NefscZ35{u=4BsHzE zpp6!^8Cua^=+>1X_FKQ)(Md6O7`0M(xU`aa1Z-}}p%Hw78Rr>S*?JhY+XL+$Nha1l zM!dKKPr}WivWuE4%uW6VLSaL#p%t}3 z+3#THVYxYf0E>B0EU{Rwuqdw_jrdpJ;t;&i(snZx0_CH|>_}a1%+pIO{uXk!dD|T2 z3zo5}+{zaq$RP`@lTi3i9CwFr8qSV4Uas-v&(Gx+8xv}^>`$%v6R)vHuBo=EiWc*3 zHb^Q*ZZK$!X|)=m%m$5C;U^N)#>VnPDgj3D$E=XALBP8P3P}HjNOuH$?fn%6cma|W z(9(k1SWriQNdbK6ckBXTju!$=cM~o)L!gI#P)pLo2&hFyR9+yOZG%6C6lE12OjspC zHOZtaEM(VOP^JZKfo32`3EzG&2|Npq_xLsUN zR$_XTk4Qi|A)p6p%HazND^)Kd|8Eu<9zg+^=(3R77wHiD+F#h-4^`S>`vAGaF~d+e zGV1VX)-;?i-Hs0;FMeATS`L5y!x79A!qnl&8#!99W_t3Qqfr|h8r=E)qmfq%%$E7+ z-B(xLLsw+_`bE4tal^5E)+2FZ{4AgNqbJ2puY0U|>UG3Wm3Vi+%8Q%D$RD!4d-v4_ zF~84^X2XUS8K&)c-u4CmgIIpUC2i7K7?H`&*brd#;+9)a>Ajb&yUq z-V;CkUC*S39~9RVLisXw@vRc`)jQ`}gB9Z;CWLjLZq`hxG&Bh}t~7+^i!E@|GUq!#@Rrb~gvd6_%J?i&uo%@v>O2le$y+IL?szO!qCVaDNghEwAP z{cZopHgV?Sz~%E={B35!pT>OR;A&_ydw-J~>Pe!98O9BWXeb7(8-AuY%st-)-)Z?o^q72bp4o7Z^ z7-c9taC&6NVaLothrVijGqQ;ov+m^DFWR&;Fc-bcvwX*kv5T$;bsDEK9QyhGBcESf z#gWhUZJik(Azr<=+GpGrA4A-WkM2CdUWOFL)V0~*o4x2*La&WtolTId`j1k442YuHnMKh#ZANN7JZ(n7Dt6HbxS$X zeSwx4~82wk&9E;PT6S)>Gs5X=H$wNFl|)#8%OPn@px86$o;a>(^L zp6kU2>m6-Px1O167Y_9K;LCxAGu{7elD9KKOwpbFsQSZ4=A|b(OwL(lvVyA0CSr@i+ z7GJ(PwB84tNZFnB0qn(ZAC2ztRgWM{BHlCxMWNwB$WH8p*n9;PDB$#PK5iGjnKvF? zM|?&y@)RmRM7`~x1*~YGg>5vSD&ituH4_D*D1K-F_7NgyqT3ir1wz#vWJD-hs*b~= zkrbp8l2CJ!y*?TJV+VC19v2DymqA@N8c8z}J3m?y3x5h(|FpINb^uTS0BAwP5-L}L zTtMc`YMApBE^Mj7dy^G!wo(y-8G*=5UnnGOCgmgKuILmHGJ;MIuwE9*w;?yoj!-0u z8zCzyUI7Uscu?IjKD!6j@l{!^Y8%5z|$fdyoqohBIQl`H(&6D_8|WFOLW- z_M*+el?Mf0eXFUGH#=f)Re%%|tX6hp8^CT-q|c(1Z_-UbK1fbC0c?RJr@RMXPb>15 zO8I05kSU-hYRl48Xd$2kgv&rzWN|1Y9451B`3PbC5u{hM-^2))3JLHSEx4ATa0hl} ztYk#Z&tzWyB39Aau2@-TISM$efD^Hj74yqc8B!O=O3Pw6TPR1Kxglq}+Oado$i_(- zL;XxB_?t}PuxMjVPlTpbD2%mFn2f1Pvf8ttYWwofr6jW!^Xs;7@o0=BS9T4^#a2C| zwtZugA3$dIR(!8YR03-O%|GZB{5d~HzWn;fNXx(f7dr1*QtFEvR<%c{&-BOfX*l;qo@Ke&HE@zb{*g9=f zFHPE7mp7`HCcAR$ay6u(O;Za)yLuJu`4`WT?*K&}2Ph_bQPPfpiHu6xsb1O{($Z6D zkd|v(U(n^+!z(%J0v*Fed{%7h6c?7H@|?x>oIdz~=I~S>^I|=h+U?~29W$M4x<5Mo zCEs-Mj90tHV&=zA1OIB&EP33+E1`~z&Ve~GzZY}KWr>Sl`rHQtwO`&`6I|YsshGIq z_>|Ag=l9wyt{lQ7C!xe1w)xQJ#1bM%q3LCmZ|hqSvLI3Lj%$Agkw~n_uK`TX6Uc*% z} zPboHrOi>Bh8a4vkr|y${!*l=}@B8}LU67Iy7cnp1fid%38Z(iNSs2S;7bzm6MndpT z?2M4J{O>9M`|=;`!oLVrUPL_q7yl6tof2l6NyDpriil6<1+kYfn_EkVq8uqTIOHM6 z5Ao-`+4M!mo}6z+!4R@o$9G>Kx;YjiLqHhxT%yFaAlW@vH(U*&4GR zRkHehWPn>&$+1rr&O+>KUnK-}VOQAUK2=CKYlG?<3E6|#7}kE%BT0Mt|6t=+CLcRB z@q31evFG>$&28dQm88u0Aw-!ULSwq&VrMH1=PL|9g)%zNK@&Q?X-0rT+4WVzy^(B; z9a~T(8@aSfiiX;WY&up~`AH15yvzHnuYLYdDOnc9(c=ZSToynznFvs65rBCq_gJL> zjRLeUB>}Q$5&?8Cg;leN0Pnq&1sMBM5+HvOEAOq@=}OyqFL^R03YsNEIWbcuX_BXN zBs#~XNy*ezl5lv@OF9`J6=2)_C1Tge9Pnnyx9qt2J!$N}TW-AXCCdH0Hdx+n-_MZC z=yn_Dzogcgq+-ufK|GY*_7={*>I+#@UsGW=b0*fb7Gh z-}-8av*sA#{C>5}8IL3A9$9Gt$hZ#xDV!6liJs3}?a0~13fHwV*R)e$pSN0M&@s~N zzl#hyMz&Cv$At1?_9BZq!!Y&XOV2@o=|rI}0;w`f*r6zxCKF1aB$UDeC}FiB6lcF+ zW05Rh`FZFrnKB$q6GqXlAPWclyh+|XaRB06vYMopOZ>^Dq>R(HQI#44_ zsbhmxsGYF?1h!>SKEJJmg>V%?&lZYms3B22Kg$RE@)2sRM$JgSjg=i#Y??P!K|0Cx z0SL{u(q>m^vx9s>Al5ecwuvP-dk7nI>^%`@JmKt#Z0OsTz8lfEioVt4&1VN;oiRqE zK@P@ks&`;!lO8h7lIQ=rye8#N+l(+j)^yYL2!=_+0J+Ubnp>VVM3gcLjt_;oydJA1RLbx9~$ zwtC9lK)FZIaoCf-d(n4q`tC#Ded)U&bA#CvxJ`jr9(g2k@Ma7!h zNKq+?Mx`jF-2@rnIknQ>nN95|MY7FEB_y=KQxLmUscclhH_{4E-blM*Tv)3_?#s1u z}x7My(bR9n#Z(AvqQy~_r%-Tnp*y?x9W5j z_5SZHO1J9vb+R?Q>Ih|NtKOzg*{XM}ql!?r>fP#OP4$-Y$yS|@8-&$?kL87|tOKi! zFU?7i^D9S?{=r%m_da^)?*n zQcw7V$-=2^AQAgp9i2bfh-AsX<{~xvguu)ru$q`p5YS)~764LaPaKhe$)!dv-PS8R zEEj&raqK0y=iy5>J%~_Bfh%C^tuNUwT5rwVA#J_&Y|O(}(!Yz+oZRJFH}|a9ob2kX zn>V4&EgRdW^_r8a?ct}8Z&t6lGFK_(KK$8ZctR5;NpG#Kn(Iih7a#W;tBl$HoJH#9 zP{&$oa(`{z>6Z1_`cXQ@#Gv8qBiRFSWCE6Nqhc&54iH?zjX-R~V&U*^9E;fa1m0sW zn0)*kLBAJYRnu);*F-E|77a|q^gfC$n?|f@jN@Po&ZH1FsBG#&3Sl_1X)7rtaVa?s zAlN&Joc*|_t$WkS!uLvRF*byksu_5b7mj6t6S?bR%QGk;QMsB;KUi zC^g&4r!GgaC{aia;%pI~4>C))Z=V8^P9OnDCypMofxNjl zHy(N`Yk0S!;pDZnvqZzC=@b$T>oO?>4Wj`5M0?iPNHnb7E$M)0xVcNtn+h%~kriB; zFV!+oqaINdoSI8X5(TFpq%bvFjSdk+^qZY0wc5&?T7j+<{r+sx@9!3LPeJTsS;MJ^ zWew*SN*dmLifGvA96cef0FHsQ3`!n1q5fGL%|eYEg)E1YD?1hQ5y_CNlM$-uj(($I zm2Q)f)@wQ?gH4FGg&lQJ6t#8ESvGT+Tw6_N%D}^W)qY1rClIue^31&WQrOI6p;~q`W%I+ zka-r6JfNGqUXk1_2KbOicNYtfIBpV)pA>O4vS$=z(G^M4BwnNDida(_Aw+%aIweW? z{7g|~lnR%%u*~6e?Kn+St0Nl9Adok{p1c_kDjRXaffLwse3~d_iFW$*q#wN!#DNBq`*wrvNx|a8j5gq&rFfCTI93gvqxa{Dxp|fqGMo ztO`N5kL!f^!J(Xvz5JapBb0-$aPC?x^jz`FkjIpQ>{aVvL_mYET+4FS)Ccudi`bfpMICt%dS%4bbo-;ss5iSN?>hv^ObSjUHTOct8K@Je#kYp ze5ixx6ytwS4cVA5e%4gESUNadT#U)F?(0VU3>R?SK6AM_>?Q+$crIr|F0jyLGm81O zMi>$c=Wz>QaWa6?7n~Mz?GgDBNWHyS@=7MR`PmnZr~}1Czv)1lDT)KB#zF@qrR^)0 z9VmY!?H|md_;d_~uuQ>aqX;ApR52RRaib<778T3>QvyOuhWp8aezBlCK$(u3gBjpI z`D2C9^~8T377Nf3@E`7s;y>+#3W@(PlPM(rQ#lFXM&9&3 zYVMlwW=E@*DIKj$f=Xk_Q6~>aozRIu(w@X3Topt~kxusSOdwxc17oKx&4rkZBUS+1 zNUiJ#0+PsrO|xk+QBw9C+Q>8v8qFn;D5-QFAZe6!IwNbT5A;*g(qIc3t{`^I8KGbg z2WvS~ve0B7QP9LQLgId+poJC%tyUBiK9lk%3QC_v`NMs;noXp+L_rlxX|9RnehQGt z-{S4hV139c51#sS#MdAQ8KCr|C^mA9;A7@ISk|^Ra zRakSMVsb#sH#-4C%yMVoFWD_YUUybLq7@UMB2K#_;SwN_^(V#T92al8OE9&{4n->w zt7#t5%BHonYWgyy=HDyioV~10-D!{Hjxky5Xd7b5{B(elp$}SAUIb{{Nv$CFsdU_qkSJk(|te?SsaC@i#Vve-K69Gw| z$jybSEvFH-9jLlig5T2np;Y8b7(FZU?JH6H%gX{5M)c29U_*_$IBFo zA$hbdncUcX+7>#-cYg?Gwz_Y4L-m^u(R?%Lue#M}!%jWJ+|sj?$i`S%L3Rghy@KrV-_@w)DjJ%#z-=|1(+eIX62-+!;7`_-%Pu8{JYE3>5$JpIC(J{IwZ;$z|0Kyi|f z)fLiJY9X$;NcB%#y6O^z#K%-u2qeKt`x77tPFG{uB2Kzo!rp+S;MB^3+F4L%SOC)y zPF<|xn&e_^k0Jrg0TF5r)Ra7IDb%DMmP#&{q$NOEDTTzt0B2>t23W@l;cU Date: Wed, 14 Aug 2024 19:25:57 -0600 Subject: [PATCH 18/26] add custom/auto x_scale --- desc/objectives/utils.py | 29 ++++++++++++++++++--------- desc/optimize/_constraint_wrappers.py | 4 +++- desc/perturbations.py | 12 +++++++---- tests/test_linear_objectives.py | 8 ++++---- 4 files changed, 34 insertions(+), 19 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index bc2c23fbee..ba80ddc75e 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -9,7 +9,9 @@ from desc.utils import Index, errorif, flatten_list, svd_inv_null, unique_list, warnif -def factorize_linear_constraints(objective, constraint): # noqa: C901 +def factorize_linear_constraints( # noqa: C901 + objective, constraint, things=None, x_scale="auto" +): """Compute and factorize A to get pseudoinverse and nullspace. Given constraints of the form Ax=b, factorize A to find a particular solution xp @@ -22,6 +24,13 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 Objective function to optimize. constraint : ObjectiveFunction Objective function of linear constraints to enforce. + things : Optimizable or tuple/list of Optimizable + Things to optimize. Defaults to ``objective.things``. + Only used if ``x_scale='auto'``. + x_scale : array_like or ``'auto'``, optional + Characteristic scale of each variable. Setting ``x_scale`` is equivalent + to reformulating the problem in scaled variables ``xs = x / x_scale``. + If set to ``'auto'``, the scale is determined from the initial state vector. Returns ------- @@ -136,16 +145,16 @@ def factorize_linear_constraints(objective, constraint): # noqa: C901 unfixed_idx = indices_idx fixed_idx = np.delete(np.arange(xp.size), unfixed_idx) - # unscaled particular solution to get scale of x - if A.size: - A_inv, Z = svd_inv_null(A) - else: - A_inv = A.T - Z = np.eye(A.shape[1]) - xp = put(xp, unfixed_idx, A_inv @ b) - D = np.where(np.abs(xp) < 1, 1, np.abs(xp)) # TODO: adjust threshold? + # get x_scale if not provided + if x_scale == "auto": + if things is None: + things = objective.things + elif not isinstance(things, list): + things = [things] + x_scale = objective.x(*[things[things.index(t)] for t in objective.things]) + D = np.where(np.abs(x_scale) < 1e1, 1, np.abs(x_scale)) # TODO: adjust threshold? - # scaled system + # null space & particular solution A = A * D[None, unfixed_idx] if A.size: A_inv, Z = svd_inv_null(A) diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index d6426059c9..b3b4e71261 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -537,7 +537,9 @@ def _set_eq_state_vector(self): linear_constraint = ObjectiveFunction(self._linear_constraints) linear_constraint.build() _, _, _, self._Z, self._D, self._unfixed_idx, _, _ = ( - factorize_linear_constraints(self._constraint, linear_constraint) + factorize_linear_constraints( + self._constraint, linear_constraint, things=[self._eq] + ) ) # dx/dc - goes from the full state to optimization variables for eq diff --git a/desc/perturbations.py b/desc/perturbations.py index 03529a1ed1..397270586f 100644 --- a/desc/perturbations.py +++ b/desc/perturbations.py @@ -186,7 +186,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces print("Factorizing linear constraints") timer.start("linear constraint factorize") xp, _, _, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint + objective, constraint, things=eq ) timer.stop("linear constraint factorize") if verbose > 1: @@ -388,7 +388,9 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces con.update_target(eq_new) constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - _, _, _, _, _, _, _, recover = factorize_linear_constraints(objective, constraint) + _, _, _, _, _, _, _, recover = factorize_linear_constraints( + objective, constraint, things=[eq_new] + ) # update other attributes dx_reduced = dx1_reduced + dx2_reduced + dx3_reduced @@ -546,7 +548,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces constraint.build(verbose=verbose) _, _, _, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective_f, constraint + objective_f, constraint, things=eq ) # state vector @@ -750,7 +752,9 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces con.update_target(eq_new) constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - _, _, _, _, _, _, _, recover = factorize_linear_constraints(objective_f, constraint) + _, _, _, _, _, _, _, recover = factorize_linear_constraints( + objective_f, constraint, things=[eq_new] + ) # update other attributes dx_reduced = dx1_reduced + dx2_reduced diff --git a/tests/test_linear_objectives.py b/tests/test_linear_objectives.py index f67f934a55..139488af7a 100644 --- a/tests/test_linear_objectives.py +++ b/tests/test_linear_objectives.py @@ -452,7 +452,7 @@ def test_correct_indexing_passed_modes(): constraint.build() xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint + objective, constraint, things=eq ) x1 = objective.x(eq) @@ -515,7 +515,7 @@ def test_correct_indexing_passed_modes_and_passed_target(): constraint.build() xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint + objective, constraint, things=eq ) x1 = objective.x(eq) @@ -575,7 +575,7 @@ def test_correct_indexing_passed_modes_axis(): constraint.build() xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint + objective, constraint, things=eq ) x1 = objective.x(eq) @@ -704,7 +704,7 @@ def test_correct_indexing_passed_modes_and_passed_target_axis(): constraint.build() xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint + objective, constraint, things=eq ) x1 = objective.x(eq) From 7624f34d7f486dbfb5ed58df90667d26b5aeafc1 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 14 Aug 2024 22:17:25 -0600 Subject: [PATCH 19/26] fix remaining tests --- desc/objectives/utils.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index ba80ddc75e..4467903708 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -6,7 +6,15 @@ import numpy as np from desc.backend import cond, jit, jnp, logsumexp, put -from desc.utils import Index, errorif, flatten_list, svd_inv_null, unique_list, warnif +from desc.utils import ( + Index, + errorif, + flatten_list, + get_instance, + svd_inv_null, + unique_list, + warnif, +) def factorize_linear_constraints( # noqa: C901 @@ -145,14 +153,20 @@ def factorize_linear_constraints( # noqa: C901 unfixed_idx = indices_idx fixed_idx = np.delete(np.arange(xp.size), unfixed_idx) - # get x_scale if not provided + # compute x_scale if not provided if x_scale == "auto": if things is None: things = objective.things - elif not isinstance(things, list): - things = [things] - x_scale = objective.x(*[things[things.index(t)] for t in objective.things]) - D = np.where(np.abs(x_scale) < 1e1, 1, np.abs(x_scale)) # TODO: adjust threshold? + else: + things = [things] if not isinstance(things, list) else things + things = [get_instance(things, type(t)) for t in objective.things] + x_scale = objective.x(*things) + errorif( + x_scale.shape != xp.shape, + ValueError, + "x_scale must be the same size as the full state vector.", + ) + D = np.where(np.abs(x_scale) < 1e1, 1, np.abs(x_scale)) # null space & particular solution A = A * D[None, unfixed_idx] From 8ed20d7cff36da778a2bccb9e38ab2e2601bf9ca Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Wed, 14 Aug 2024 22:25:13 -0600 Subject: [PATCH 20/26] update docstring --- desc/objectives/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 4467903708..10f3cf9bdd 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -51,7 +51,7 @@ def factorize_linear_constraints( # noqa: C901 Z : ndarray Null space operator for full combined A such that A @ Z == 0. D : ndarray - Scale of the full state vector x, based on the particular solution xp. + Scale of the full state vector x, as set by the parameter ``x_scale``. unfixed_idx : ndarray Indices of x that correspond to non-fixed values. project, recover : function From 461c0ef85fd46d47cce6e910b53ebdad67c3ab8a Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Thu, 15 Aug 2024 10:36:31 -0600 Subject: [PATCH 21/26] increase threshold to 20 --- desc/objectives/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 10f3cf9bdd..a7c2e7e8b1 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -166,7 +166,7 @@ def factorize_linear_constraints( # noqa: C901 ValueError, "x_scale must be the same size as the full state vector.", ) - D = np.where(np.abs(x_scale) < 1e1, 1, np.abs(x_scale)) + D = np.where(np.abs(x_scale) < 2e1, 1, np.abs(x_scale)) # null space & particular solution A = A * D[None, unfixed_idx] From 56103e6dbce9851eda7d74fae5da89a72ff9eb48 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Thu, 15 Aug 2024 13:27:56 -0600 Subject: [PATCH 22/26] all tests hopefully passing now? --- tests/test_bootstrap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_bootstrap.py b/tests/test_bootstrap.py index 34ffe08e22..811480bbaf 100644 --- a/tests/test_bootstrap.py +++ b/tests/test_bootstrap.py @@ -1589,7 +1589,7 @@ def test_bootstrap_optimization_comparison_qa(): objective=objective, constraints=constraints, optimizer="proximal-lsq-exact", - maxiter=4, + maxiter=5, gtol=1e-16, verbose=3, ) @@ -1622,5 +1622,5 @@ def test_bootstrap_optimization_comparison_qa(): grid.compress(data2[""]), grid.compress(data2[" Redl"]), rtol=1.8e-2 ) np.testing.assert_allclose( - grid.compress(data1[""]), grid.compress(data2[""]), rtol=1.8e-2 + grid.compress(data1[""]), grid.compress(data2[""]), rtol=1.9e-2 ) From 6187f566f3249cb428c6b9ce50d9485ea745f66f Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Fri, 16 Aug 2024 10:47:52 -0600 Subject: [PATCH 23/26] making requested changes --- desc/objectives/utils.py | 5 +++-- desc/optimize/_constraint_wrappers.py | 8 +++++--- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index a7c2e7e8b1..58882c3ba9 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -164,9 +164,10 @@ def factorize_linear_constraints( # noqa: C901 errorif( x_scale.shape != xp.shape, ValueError, - "x_scale must be the same size as the full state vector.", + "x_scale must be the same size as the full state vector. " + + f"Got size {x_scale.size} for state vector of size {xp.size}.", ) - D = np.where(np.abs(x_scale) < 2e1, 1, np.abs(x_scale)) + D = np.where(np.abs(x_scale) < 1e2, 1, np.abs(x_scale)) # null space & particular solution A = A * D[None, unfixed_idx] diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index b3b4e71261..b138ef2a3d 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -260,7 +260,7 @@ def grad(self, x_reduced, constants=None): """ x = self.recover(x_reduced) df = self._objective.grad(x, constants) - return df @ self._unfixed_idx_mat + return df[self._unfixed_idx] @ (self._Z * self._D[self._unfixed_idx, None]) def hess(self, x_reduced, constants=None): """Compute Hessian of self.compute_scalar. @@ -290,7 +290,9 @@ def _jac(self, x_reduced, constants=None, op="scaled"): x = self.recover(x_reduced) if self._objective._deriv_mode == "blocked": fun = getattr(self._objective, "jac_" + op) - return fun(x, constants) @ self._unfixed_idx_mat + return fun(x, constants)[:, self._unfixed_idx] @ ( + self._Z * self._D[self._unfixed_idx, None] + ) v = self._unfixed_idx_mat df = getattr(self._objective, "jvp_" + op)(v.T, x, constants) @@ -404,7 +406,7 @@ def jvp_unscaled(self, v, x_reduced, constants=None): def _vjp(self, v, x_reduced, constants=None, op="vjp_scaled"): x = self.recover(x_reduced) df = getattr(self._objective, op)(v, x, constants) - return df @ self._unfixed_idx_mat + return df[:, self._unfixed_idx] @ (self._Z * self._D[self._unfixed_idx, None]) def vjp_scaled(self, v, x_reduced, constants=None): """Compute vector-Jacobian product of self.compute_scaled. From 6b9233aa8383ab7a1e826844a70930065b6c9515 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Fri, 16 Aug 2024 11:17:56 -0600 Subject: [PATCH 24/26] fix bug from last commit --- desc/optimize/_constraint_wrappers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index b138ef2a3d..4021dd42f6 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -406,7 +406,7 @@ def jvp_unscaled(self, v, x_reduced, constants=None): def _vjp(self, v, x_reduced, constants=None, op="vjp_scaled"): x = self.recover(x_reduced) df = getattr(self._objective, op)(v, x, constants) - return df[:, self._unfixed_idx] @ (self._Z * self._D[self._unfixed_idx, None]) + return df[self._unfixed_idx] @ (self._Z * self._D[self._unfixed_idx, None]) def vjp_scaled(self, v, x_reduced, constants=None): """Compute vector-Jacobian product of self.compute_scaled. From cc1ddee55a8dc1ae74413504bacf01bc28bfa9db Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Fri, 16 Aug 2024 16:29:10 -0600 Subject: [PATCH 25/26] always use default things --- desc/optimize/_constraint_wrappers.py | 4 +--- desc/perturbations.py | 12 ++++-------- tests/test_linear_objectives.py | 8 ++++---- 3 files changed, 9 insertions(+), 15 deletions(-) diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py index 4021dd42f6..8d4e480161 100644 --- a/desc/optimize/_constraint_wrappers.py +++ b/desc/optimize/_constraint_wrappers.py @@ -539,9 +539,7 @@ def _set_eq_state_vector(self): linear_constraint = ObjectiveFunction(self._linear_constraints) linear_constraint.build() _, _, _, self._Z, self._D, self._unfixed_idx, _, _ = ( - factorize_linear_constraints( - self._constraint, linear_constraint, things=[self._eq] - ) + factorize_linear_constraints(self._constraint, linear_constraint) ) # dx/dc - goes from the full state to optimization variables for eq diff --git a/desc/perturbations.py b/desc/perturbations.py index 397270586f..03529a1ed1 100644 --- a/desc/perturbations.py +++ b/desc/perturbations.py @@ -186,7 +186,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces print("Factorizing linear constraints") timer.start("linear constraint factorize") xp, _, _, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint, things=eq + objective, constraint ) timer.stop("linear constraint factorize") if verbose > 1: @@ -388,9 +388,7 @@ def perturb( # noqa: C901 - FIXME: break this up into simpler pieces con.update_target(eq_new) constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - _, _, _, _, _, _, _, recover = factorize_linear_constraints( - objective, constraint, things=[eq_new] - ) + _, _, _, _, _, _, _, recover = factorize_linear_constraints(objective, constraint) # update other attributes dx_reduced = dx1_reduced + dx2_reduced + dx3_reduced @@ -548,7 +546,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces constraint.build(verbose=verbose) _, _, _, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective_f, constraint, things=eq + objective_f, constraint ) # state vector @@ -752,9 +750,7 @@ def optimal_perturb( # noqa: C901 - FIXME: break this up into simpler pieces con.update_target(eq_new) constraint = ObjectiveFunction(constraints) constraint.build(verbose=verbose) - _, _, _, _, _, _, _, recover = factorize_linear_constraints( - objective_f, constraint, things=[eq_new] - ) + _, _, _, _, _, _, _, recover = factorize_linear_constraints(objective_f, constraint) # update other attributes dx_reduced = dx1_reduced + dx2_reduced diff --git a/tests/test_linear_objectives.py b/tests/test_linear_objectives.py index 139488af7a..f67f934a55 100644 --- a/tests/test_linear_objectives.py +++ b/tests/test_linear_objectives.py @@ -452,7 +452,7 @@ def test_correct_indexing_passed_modes(): constraint.build() xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint, things=eq + objective, constraint ) x1 = objective.x(eq) @@ -515,7 +515,7 @@ def test_correct_indexing_passed_modes_and_passed_target(): constraint.build() xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint, things=eq + objective, constraint ) x1 = objective.x(eq) @@ -575,7 +575,7 @@ def test_correct_indexing_passed_modes_axis(): constraint.build() xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint, things=eq + objective, constraint ) x1 = objective.x(eq) @@ -704,7 +704,7 @@ def test_correct_indexing_passed_modes_and_passed_target_axis(): constraint.build() xp, A, b, Z, D, unfixed_idx, project, recover = factorize_linear_constraints( - objective, constraint, things=eq + objective, constraint ) x1 = objective.x(eq) From ddb63d299fb08fde9ebeb9fd86a3bf4c71511822 Mon Sep 17 00:00:00 2001 From: daniel-dudt Date: Mon, 19 Aug 2024 10:04:37 -0600 Subject: [PATCH 26/26] always use objective.things for auto scale --- desc/objectives/utils.py | 28 +++++----------------------- 1 file changed, 5 insertions(+), 23 deletions(-) diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py index 58882c3ba9..5167573f51 100644 --- a/desc/objectives/utils.py +++ b/desc/objectives/utils.py @@ -6,20 +6,10 @@ import numpy as np from desc.backend import cond, jit, jnp, logsumexp, put -from desc.utils import ( - Index, - errorif, - flatten_list, - get_instance, - svd_inv_null, - unique_list, - warnif, -) - - -def factorize_linear_constraints( # noqa: C901 - objective, constraint, things=None, x_scale="auto" -): +from desc.utils import Index, errorif, flatten_list, svd_inv_null, unique_list, warnif + + +def factorize_linear_constraints(objective, constraint, x_scale="auto"): # noqa: C901 """Compute and factorize A to get pseudoinverse and nullspace. Given constraints of the form Ax=b, factorize A to find a particular solution xp @@ -32,9 +22,6 @@ def factorize_linear_constraints( # noqa: C901 Objective function to optimize. constraint : ObjectiveFunction Objective function of linear constraints to enforce. - things : Optimizable or tuple/list of Optimizable - Things to optimize. Defaults to ``objective.things``. - Only used if ``x_scale='auto'``. x_scale : array_like or ``'auto'``, optional Characteristic scale of each variable. Setting ``x_scale`` is equivalent to reformulating the problem in scaled variables ``xs = x / x_scale``. @@ -155,12 +142,7 @@ def factorize_linear_constraints( # noqa: C901 # compute x_scale if not provided if x_scale == "auto": - if things is None: - things = objective.things - else: - things = [things] if not isinstance(things, list) else things - things = [get_instance(things, type(t)) for t in objective.things] - x_scale = objective.x(*things) + x_scale = objective.x(*objective.things) errorif( x_scale.shape != xp.shape, ValueError,