Merge FreeMagnetismInterface feature enhancement (#232)

* ENH: add free interface in addition to free layer for magnetism. * FIX: length of vectors not equal when dthetaM is not defined Also fixed issue with magnetic profile calculation when drhoM[-1] = 0. * FIX: Make thetaM behave the same as rhoM Removed conditions at lines 514 and 521. Also added same statement at line 520 as is used for rhoM. This should now behave in the same way as FreeInterface. * Lint * FIX: nans in mono.py as result of divide by zero Removed divide by zero when z[-1]=0. * FIX: changed FreeMagnetismInterface thetaM above or below default value when not specified If only tabove or tbelow are specified then they are made equal. If neither are specified they are made the same as DEFAULT_THETA_M. This is to get around the issue when interfacing with a non-magnetic layer where the thetaM value should not change. * Additional Edits for FreeMagnetismInterface Additional edits as suggested by Paul K. * create self-extracting version of unstable windows .exe.zip * numba required for readthedocs build * set max initial ylim on spin asymmetry plot * add link to unstable release page * ENH: add free interface in addition to free layer for magnetism. * FIX: length of vectors not equal when dthetaM is not defined Also fixed issue with magnetic profile calculation when drhoM[-1] = 0. * FIX: Make thetaM behave the same as rhoM Removed conditions at lines 514 and 521. Also added same statement at line 520 as is used for rhoM. This should now behave in the same way as FreeInterface. * Lint * FIX: nans in mono.py as result of divide by zero Removed divide by zero when z[-1]=0. * FIX: changed FreeMagnetismInterface thetaM above or below default value when not specified If only tabove or tbelow are specified then they are made equal. If neither are specified they are made the same as DEFAULT_THETA_M. This is to get around the issue when interfacing with a non-magnetic layer where the thetaM value should not change. * Additional Edits for FreeMagnetismInterface Additional edits as suggested by Paul K. * ENH: add free interface in addition to free layer for magnetism. * FIX: length of vectors not equal when dthetaM is not defined Also fixed issue with magnetic profile calculation when drhoM[-1] = 0. * FIX: Make thetaM behave the same as rhoM Removed conditions at lines 514 and 521. Also added same statement at line 520 as is used for rhoM. This should now behave in the same way as FreeInterface. * Lint * FIX: nans in mono.py as result of divide by zero Removed divide by zero when z[-1]=0. * FIX: changed FreeMagnetismInterface thetaM above or below default value when not specified If only tabove or tbelow are specified then they are made equal. If neither are specified they are made the same as DEFAULT_THETA_M. This is to get around the issue when interfacing with a non-magnetic layer where the thetaM value should not change. * Additional Edits for FreeMagnetismInterface Additional edits as suggested by Paul K. * ENH: add free interface in addition to free layer for magnetism. * FIX: Make thetaM behave the same as rhoM Removed conditions at lines 514 and 521. Also added same statement at line 520 as is used for rhoM. This should now behave in the same way as FreeInterface. * Additional Edits for FreeMagnetismInterface Additional edits as suggested by Paul K. * ENH: add free interface in addition to free layer for magnetism. * FIX: length of vectors not equal when dthetaM is not defined Also fixed issue with magnetic profile calculation when drhoM[-1] = 0. * FIX: Make thetaM behave the same as rhoM Removed conditions at lines 514 and 521. Also added same statement at line 520 as is used for rhoM. This should now behave in the same way as FreeInterface. * Lint * FIX: nans in mono.py as result of divide by zero Removed divide by zero when z[-1]=0. * FIX: changed FreeMagnetismInterface thetaM above or below default value when not specified If only tabove or tbelow are specified then they are made equal. If neither are specified they are made the same as DEFAULT_THETA_M. This is to get around the issue when interfacing with a non-magnetic layer where the thetaM value should not change. * Additional Edits for FreeMagnetismInterface Additional edits as suggested by Paul K. * ENH: add free interface in addition to free layer for magnetism. * FIX: Make thetaM behave the same as rhoM Removed conditions at lines 514 and 521. Also added same statement at line 520 as is used for rhoM. This should now behave in the same way as FreeInterface. * Additional Edits for FreeMagnetismInterface Additional edits as suggested by Paul K. * ENH: add free interface in addition to free layer for magnetism. * FIX: length of vectors not equal when dthetaM is not defined Also fixed issue with magnetic profile calculation when drhoM[-1] = 0. * FIX: changed FreeMagnetismInterface thetaM above or below default value when not specified If only tabove or tbelow are specified then they are made equal. If neither are specified they are made the same as DEFAULT_THETA_M. This is to get around the issue when interfacing with a non-magnetic layer where the thetaM value should not change. * Additional Edits for FreeMagnetismInterface Additional edits as suggested by Paul K. * Linting and bringing in line with master * Maint: Serialize `FreeMagnetismInterface` * Tidying up to minimise changes When updating the branch to the latest version from the master branch a lot of merge conflicts had to be gone through. Here I am correcting where I accepted minor changes that do not matter to the code but deviate from the master branch unnecessarily. * Fix: Serializing/De-serializing of FreeInterface * Tidy up imports * Sorting Merge conflicts * Further import tidy up * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Updated dataclass types for `FreeMagnetismInterface` * minor fixup on stored attribute types (they are always stored as Parameter after init, even ones that are supplied as float to the init function) * fixes to FreeInterface class, adding magnetism and setting last element of z to 1 if it is zero * add FreeMagnetismInterface class * add FreeMagnetismInterface to names * FIX: Limits of rhoM should be +- inf in `FreeMagnetism` * ENH: add magnetism attr to `FreeLayer` * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: Paul Kienzle <[email protected]> Co-authored-by: Brian Benjamin Maranville <[email protected]> Co-authored-by: Paul Kienzle <[email protected]> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
reflectometry · Feb 21, 2025 · beb331e · beb331e
1 parent b09d43b
commit beb331e
Show file tree

Hide file tree

Showing 3 changed files with 173 additions and 11 deletions.
diff --git a/refl1d/names.py b/refl1d/names.py
@@ -39,7 +39,7 @@
 from .sample.cheby import ChebyVF, FreeformCheby, cheby_approx, cheby_points
 from .sample.flayer import FunctionalMagnetism, FunctionalProfile
 from .sample.layers import Slab, Stack
-from .sample.magnetism import FreeMagnetism, Magnetism, MagnetismStack, MagnetismTwist
+from .sample.magnetism import FreeMagnetism, FreeMagnetismInterface, Magnetism, MagnetismStack, MagnetismTwist
 from .sample.material import SLD, Compound, Material, Mixture
 
 # Pull in common materials for reflectometry experiments.

diff --git a/refl1d/sample/magnetism.py b/refl1d/sample/magnetism.py
@@ -386,7 +386,7 @@ def parvec(vector, name, limits):
 
         self.rhoM, self.thetaM, self.z = [
             parvec(v, name + " " + part, limits)
-            for v, part, limits in zip((rhoM, thetaM, z), ("rhoM", "thetaM", "z"), ((0, None), (0, 360), (0, 1)))
+            for v, part, limits in zip((rhoM, thetaM, z), ("rhoM", "thetaM", "z"), ((None, None), (0, 360), (0, 1)))
         ]
         if len(self.z) != len(self.rhoM):
             raise ValueError("must have one position z for each rhoM")
@@ -437,3 +437,139 @@ def __str__(self):
 
     def __repr__(self):
         return "FreeMagnetism"
+
+
+@dataclass(init=False)
+class FreeMagnetismInterface(BaseMagnetism):
+    r"""
+    Spline change in magnetism throughout layer.
+
+    Defines monotonic splines for rhoM and thetaM with shared knot positions.
+
+    *dz* is the relative $z$ step between the knots, with position $z_k$
+    defined by $z_k = w \sum_{i=0}^k \delta z_i / \sum_{i=0}^n \delta z_i$,
+    where $n$ is the number of intervals. The resulting $z$ must be monotonic,
+    with $\delta z_i \ge 0$ for all intervals.
+
+    *drhoM* gives the relative $\rho_M$ step between knots. Unlike
+    $\rho_{Mk} = \sum_{i=0}^k \delta \rho_{Mi} / \sum_{i=0}^n \delta \delta \rho_{Mi}$.
+
+    *dthetaM* gives the magnetic angle for each knot.
+
+    *name* is the base name for the various layer parameters.
+
+    *dead_above* and *dead_below* define magnetically dead layers at the
+    nuclear boundaries.  These can be negative if magnetism extends beyond
+    the nuclear boundary.
+
+    *interface_above* and *interface_below* define the magnetic interface
+    at the boundaries, if it is different from the nuclear interface.
+
+    *mbelow* and *mabove* are the rhoM parameter values of
+    the above and below layers respectively. Do not specify if
+    layers either side are not magnetic.
+
+    *tbelow* and *tabove* are the thetaM parameter values of
+    the above and below layers respectively. Do not specify if
+    layers either side are not magnetic.
+    """
+
+    name: str
+    mbelow: Parameter
+    mabove: Parameter
+    tbelow: Parameter
+    tabove: Parameter
+    dz: List[Parameter]
+    drhoM: List[Parameter]
+    dthetaM: List[Parameter]
+
+    magnetic = True
+
+    def __init__(
+        self, dz=(), drhoM=(), dthetaM=(), mbelow=0, mabove=0, tbelow=None, tabove=None, name="MagInterface", **kw
+    ):
+        BaseMagnetism.__init__(self, name=name, **kw)
+
+        def parvec(vector, name, limits):
+            return [Parameter.default(p, name=name + "[%d]" % i, limits=limits) for i, p in enumerate(vector)]
+
+        self.drhoM, self.dthetaM, self.dz = [
+            parvec(v, name + " " + part, limits)
+            for v, part, limits in zip((drhoM, dthetaM, dz), ("drhoM", "dthetaM", "dz"), ((-1, 1), (-1, 1), (0, 1)))
+        ]
+        self.mbelow = Parameter.default(mbelow, name=name + " mbelow", limits=(-np.inf, np.inf))
+        self.mabove = Parameter.default(mabove, name=name + " mabove", limits=(-np.inf, np.inf))
+
+        # if only tabove or tbelow is defined then they are made equal
+        # this is to deal with the situation of a non-magnetic
+        # layer next to a magnetic one
+        if tabove is None and tbelow is None:
+            tbelow = tabove = DEFAULT_THETA_M
+        elif tbelow is None:
+            tbelow = tabove
+        elif tabove is None:
+            tabove = tbelow
+
+        self.tbelow = Parameter.default(tbelow, name=name + " tbelow", limits=(0, 360))
+        self.tabove = Parameter.default(tabove, name=name + " tabove", limits=(0, 360))
+        if len(self.dz) != len(self.drhoM):
+            raise ValueError("Need one dz for each drhoM")
+        if len(self.dthetaM) > 0 and len(self.drhoM) != len(self.dthetaM):
+            raise ValueError("Need one dthetaM for each drhoM")
+
+    def parameters(self):
+        parameters = BaseMagnetism.parameters(self)
+        parameters.update(
+            drhoM=self.drhoM,
+            dthetaM=self.dthetaM,
+            dz=self.dz,
+            mbelow=self.mbelow,
+            mabove=self.mabove,
+            tbelow=self.tbelow,
+            tabove=self.tabove,
+        )
+        return parameters
+
+    def profile(self, Pz, thickness):
+        z = np.hstack((0, np.cumsum(np.asarray([v.value for v in self.dz], "d"))))
+        if z[-1] == 0:
+            z[-1] = 1
+        z *= thickness / z[-1]
+
+        rhoM_fraction = np.hstack((0, np.cumsum(np.asarray([v.value for v in self.drhoM], "d"))))
+        # AJC added since without the line below FreeMagnetismInterface
+        # does not initialise properly - fixes strange behaviour at drho=0 on end point
+        if rhoM_fraction[-1] == 0:
+            rhoM_fraction[-1] = 1
+
+        rhoM_fraction *= 1 / rhoM_fraction[-1]
+        PrhoM = np.clip(monospline(z, rhoM_fraction, Pz), 0, 1)
+
+        if self.dthetaM:
+            thetaM_fraction = np.hstack((0, np.cumsum(np.asarray([v.value for v in self.dthetaM], "d"))))
+            if thetaM_fraction[-1] == 0:
+                thetaM_fraction[-1] = 1
+
+            thetaM_fraction *= 1 / thetaM_fraction[-1]
+            PthetaM = np.clip(monospline(z, thetaM_fraction, Pz), 0, 1)
+        else:
+            # AJC changed from len(z) to PrhoM - since PrhoM is the length of the vector
+            # we want PthetaM to match - otherwise slabs.add_magnetism throws an error
+            PthetaM = np.linspace(0.0, 1.0, len(PrhoM))
+
+        return PrhoM, PthetaM
+
+    def render(self, probe, slabs, thickness, anchor, sigma):
+        Pw, Pz = slabs.microslabs(thickness)
+        rhoM_profile, thetaM_profile = self.profile(Pz, thickness)
+        mbelow, mabove = self.mbelow.value, self.mabove.value
+        tbelow, tabove = self.tbelow.value, self.tabove.value
+        rhoM = (1 - rhoM_profile) * mbelow + rhoM_profile * mabove
+        thetaM = (1 - thetaM_profile) * tbelow + thetaM_profile * tabove
+        slabs.add_magnetism(anchor=anchor, w=Pw, rhoM=rhoM, thetaM=thetaM, sigma=sigma)
+
+    def __str__(self):
+        return "freemagint(%d)" % (len(self.drhoM))
+
+    def __repr__(self):
+        return "FreeMagnetismInterface"
diff --git a/refl1d/sample/mono.py b/refl1d/sample/mono.py
@@ -14,12 +14,15 @@
 
 from .. import utils
 from .layers import Layer
+from .magnetism import BaseMagnetism
+from .material import Scatterer
 
 
 # TODO: add left_sld, right_sld to all layers so that fresnel works
 # TODO: access left_sld, right_sld so freeform doesn't need left, right
 # TODO: restructure to use vector parameters
 # TODO: allow the number of layers to be adjusted by the fit
+@dataclass(init=False)
 class FreeLayer(Layer):
     """
     A freeform section of the sample modeled with splines.
@@ -35,7 +38,16 @@ class FreeLayer(Layer):
     with slabs.
     """
 
-    def __init__(self, below=None, above=None, thickness=0, z=(), rho=(), irho=(), name="Freeform"):
+    name: str
+    below: Scatterer
+    above: Scatterer
+    thickness: Par
+    z: List[Par]
+    rho: List[Par]
+    irho: List[Par]
+    magnetism: Optional[BaseMagnetism] = None
+
+    def __init__(self, below=None, above=None, thickness=0, z=(), rho=(), irho=(), name="Freeform", magnetism=None):
         self.name = name
         self.below, self.above = below, above
         self.thickness = Par.default(thickness, name=name + " thickness", limits=(0, inf))
@@ -53,6 +65,8 @@ def parvec(vector, name, limits):
         if len(self.irho) > 0 and len(self.z) != len(self.irho):
             raise ValueError("must have one z for each irho value")
 
+        self.magnetism = magnetism
+
     def parameters(self):
         return {
             "thickness": self.thickness,
@@ -61,6 +75,7 @@ def parameters(self):
             "z": self.z,
             "below": self.below.parameters(),
             "above": self.above.parameters(),
+            "magnetism": self.magnetism.parameters() if self.magnetism is not None else None,
         }
 
     def to_dict(self):
@@ -115,16 +130,19 @@ class FreeInterface(Layer):
     with slabs.
     """
 
-    name: Optional[str]
-    below: Optional[Any]
-    above: Optional[Any]
+    name: str
+    below: Scatterer
+    above: Scatterer
     thickness: Par
     interface: Par
-    dz: List[Union[float, Par]]
-    dp: List[Union[float, Par]]
+    dz: List[Par]
+    dp: List[Par]
+    magnetism: Optional[BaseMagnetism] = None
     # inflections: List[Any]
 
-    def __init__(self, thickness=0, interface=0, below=None, above=None, dz=None, dp=None, name="Interface"):
+    def __init__(
+        self, thickness=0, interface=0, below=None, above=None, dz=None, dp=None, name="Interface", magnetism=None
+    ):
         self.name = name
         self.below, self.above = below, above
         self.thickness = Par.default(thickness, limits=(0, inf), name=name + " thickness")
@@ -140,12 +158,11 @@ def __init__(self, thickness=0, interface=0, below=None, above=None, dz=None, dp
         if len(dz) != len(dp):
             raise ValueError("Need one dz for every dp")
 
-        # if len(z) != len(vf)+2:
-        #    raise ValueError("Only need vf for interior z, so len(z)=len(vf)+2")
         self.dz = [Par.default(p, name=name + " dz[%d]" % i, limits=(0, inf)) for i, p in enumerate(dz)]
         self.dp = [Par.default(p, name=name + " dp[%d]" % i, limits=(0, inf)) for i, p in enumerate(dp)]
         self.inflections = Par(name=name + " inflections")
         self.inflections.equals(ParFunction(inflections, dx=self.dz, dy=self.dp))
+        self.magnetism = magnetism
 
     def parameters(self):
         return {
@@ -156,6 +173,7 @@ def parameters(self):
             "below": self.below.parameters(),
             "above": self.above.parameters(),
             "inflections": self.inflections,
+            "magnetism": self.magnetism.parameters() if self.magnetism is not None else None,
         }
 
     def to_dict(self):
@@ -175,14 +193,22 @@ def profile(self, Pz):
         if p[-1] == 0:
             p[-1] = 1
         p *= 1 / p[-1]
+        # AJC included condition as if z[-1] == 0 then z *= thickness/z[-1] == [nan]*len(z)
+        # This then ends with bumps.mono.Monospline adding an extra element as a result of
+        # line 42 in bumps.mono.Monospline
+        if z[-1] == 0:
+            z[-1] = 1
         z *= thickness / z[-1]
         profile = clip(monospline(z, p, Pz), 0, 1)
         return profile
 
     def render(self, probe, slabs):
         thickness = self.thickness.value
+
+        # TODO: why is provided if it is ignored?
         # interface ignored for FreeInterface
         # interface = self.interface.value
+
         below_rho, below_irho = self.below.sld(probe)
         above_rho, above_irho = self.above.sld(probe)
         # Pz is the center, Pw is the width