FunctionalLayer can now accept lists of parameters

[hoogerheide]

For certain types of fits (splines, mixtures) with similar parameters, there is some advantage to being able to pass in lists of parameters to FunctionalLayer. This allows much cleaner scripting; the size of the list can be set dynamically without each parameter needing a separate keyword argument in the FunctionalLayer call.

This functionality requires the ability of _set_vars and FunctionalLayer._fpars to be able to unpack and repack parameter lists, respectively.

Note that the magnetism _fpars is also changed but untested.

[pkienzle]

An example would be useful, both for documentation and for testing.

I think it would be useful to keep the parameters in a list. Something like:

def poly(z, p):
    return np.polyval(p, z)
substrate = SLD('Si', rho=2.07)
profile = FunctionalProfile(100, 0, poly, p=[1, 0, 0])
sample = silicon | profile | air
# access in a loop
for coeff in profile.p:
    coeff.range(-1, 1)
# or address individually
profile.p[0].range(0, 10)

[hoogerheide]

Here's an example with the current code. The initialization value of a sinusoid is only so there's something to see when previewing the profile. This doesn't have the ability to call the parameter from the list as @pkienzle suggests. Would this require a subclass of Parameter that can handle lists, or is there an easier way?

from refl1d.probe import NeutronProbe
import numpy as np
import scipy.interpolate as si
from refl1d.names import SLD, FitProblem, Parameter, Slab, Experiment
from refl1d.flayer import FunctionalProfile
import matplotlib.pyplot as plt

def spline(z, yctrl):
    # create equally spaced control points in z
    zctrl=np.linspace(0, max(z), len(yctrl), endpoint=True)

    return si.pchip_interpolate(zctrl, yctrl, z)

n_controlpoints = 10
dimension = 100
spline_yvals = [Parameter(name='y%i' % i, value=np.sin(i/(2*np.pi*0.05))).range(-2, 2) for i in range(n_controlpoints)]

silicon = SLD('Si', rho=2.07)
air = SLD('air', rho=0)
silicon_layer = Slab(material=silicon)
layer = FunctionalProfile(100, 0, spline, yctrl=spline_yvals)
air_layer = Slab(material=air)

sample = silicon | layer | air

spline_yvals[0]= silicon.rho
spline_yvals[-1] = air.rho

probe=NeutronProbe(T=np.linspace(0, 5, 101), L=5)

model = Experiment(sample=sample, probe=probe, dz=0.5, step_interfaces=True)

problem=FitProblem(model)

"""

[pkienzle]

My inclination would be to extend the parameter class to support numpy arrays but there are some semantics to work out first. Meanwhile putting the support in the flayer class is a reasonable "temporary" solution.

I don't understand why your sample model works. Your implementation extracts the parameters from the list when you define FLayer and afterward references them as attribute name+k in the model. Setting spline_yvals[0] = silicon.rho after the fact should have no effect.

This is simple enough to work around by setting up the equality constraints prior to building FLayer object, but it goes against the spirit of bumps parameter handling, which allows you to set them afterward. You could write layer.yctrl0 = silicon.rho, but as a list layer.yctrl[0] is a little more natural.

Also DRY (don't repeat yourself): I restructured the code to isolate parameter access to functions.

-- model.py --

from refl1d.probe import NeutronProbe
import numpy as np
import scipy.interpolate as si
from refl1d.names import SLD, FitProblem, Parameter, Slab, Experiment
from refl1d.flayer import FunctionalProfile
import matplotlib.pyplot as plt

def spline(z, yctrl):
    # create equally spaced control points in z
    zctrl=np.linspace(0, max(z), len(yctrl), endpoint=True)

    return si.pchip_interpolate(zctrl, yctrl, z)

n_controlpoints = 10
dimension = 100
spline_yvals = [np.sin(i/(2*np.pi*0.05)) for i in range(n_controlpoints)]

silicon = SLD('Si', rho=2.07)
air = SLD('air', rho=0)
silicon_layer = Slab(material=silicon)
layer = FunctionalProfile(100, 0, spline, yctrl=spline_yvals)
air_layer = Slab(material=air)

sample = silicon | layer | air

layer.yctrl[0]= silicon.rho
for p in layer.yctrl[1:-1]:
    p.range(-1, 1)
layer.yctrl[-1] = air.rho

probe=NeutronProbe(T=np.linspace(0, 5, 101), L=5)

model = Experiment(sample=sample, probe=probe, dz=0.5, step_interfaces=True)

problem=FitProblem(model)

-- flayer.patch --

diff --git a/refl1d/flayer.py b/refl1d/flayer.py
index 4563e6e..b5712d7 100644
--- a/refl1d/flayer.py
+++ b/refl1d/flayer.py
@@ -78,10 +78,10 @@ class FunctionalProfile(Layer):
         self.start = SLD(name+" start", rho=rho_start, irho=irho_start)
         self.end = SLD(name+" end", rho=rho_end, irho=irho_end)
 
-        self._parameters = _set_vars(self, name, profile, kw, self.RESERVED)
+        _set_parameters(self, name, profile, kw, self.RESERVED)
 
     def parameters(self):
-        P = {k: getattr(self, k) for k in self._parameters}
+        P = _get_parameters(self)
         P['thickness'] = self.thickness
         #P['interface'] = self.interface
         return P
@@ -93,7 +93,7 @@ class FunctionalProfile(Layer):
             'thickness': self.thickness,
             'interface': self.interface,
             'profile': self.profile,
-            'parameters': {k: getattr(self, k) for k in self._parameters},
+            'parameters': _get_parameters(self),
             'tol': self.tol,
             'magnetism': self.magnetism,
         })
@@ -105,7 +105,7 @@ class FunctionalProfile(Layer):
         #print kw
         # TODO: always return rho, irho from profile function
         # return value may be a constant for rho or irho
-        phi = asarray(self.profile(Pz, **self._fpars()))
+        phi = asarray(self.profile(Pz, **_get_values(self)))
         if phi.shape != Pz.shape:
             raise TypeError("profile function '%s' did not return array phi(z)"
                             %self.profile.__name__)
@@ -114,10 +114,6 @@ class FunctionalProfile(Layer):
         slabs.extend(rho=[real(phi)], irho=[imag(phi)], w=Pw)
         #slabs.interface(self.interface.value)
 
-    def _fpars(self):
-        kw = dict((k, getattr(self, k).value) for k in self._parameters)
-        return  kw
-
 
 class FunctionalMagnetism(BaseMagnetism):
     """
@@ -153,7 +149,7 @@ class FunctionalMagnetism(BaseMagnetism):
         self.profile = profile
         self.tol = tol
 
-        self._parameters = _set_vars(self, name, profile, kw, self.RESERVED)
+        _set_parameters(self, name, profile, kw, self.RESERVED)
         rhoM_start = _MagnetismLimit(self, isend=False, isrhoM=True)
         rhoM_end = _MagnetismLimit(self, isend=True, isrhoM=True)
         thetaM_start = _MagnetismLimit(self, isend=False, isrhoM=False)
@@ -181,14 +177,14 @@ class FunctionalMagnetism(BaseMagnetism):
 
     def parameters(self):
         parameters = BaseMagnetism.parameters(self)
-        parameters.update((k, getattr(self, k)) for k in self._parameters)
+        parameters.update(_get_parameters(self))
         return parameters
 
     def to_dict(self):
         ret = BaseMagnetism.to_dict(self)
         ret.update(to_dict({
             'profile': self.profile,
-            'parameters': {k: getattr(self, k) for k in self._parameters},
+            'parameters': _get_parameters(self),
             'tol': self.tol,
         }))
 
@@ -196,8 +192,7 @@ class FunctionalMagnetism(BaseMagnetism):
         Pw, Pz = slabs.microslabs(thickness)
         if len(Pw) == 0:
             return
-        kw = dict((k, getattr(self, k).value) for k in self._parameters)
-        P = self.profile(Pz, **kw)
+        P = self.profile(Pz, **_get_values(self))
 
         rhoM, thetaM = P if isinstance(P, tuple) else (P, DEFAULT_THETA_M)
         try:
@@ -212,15 +207,11 @@ class FunctionalMagnetism(BaseMagnetism):
         slabs.add_magnetism(
             anchor=anchor, w=Pw, rhoM=rhoM, thetaM=thetaM, sigma=sigma)
 
-    def _fpars(self):
-        kw = dict((k, getattr(self, k).value) for k in self._parameters)
-        return  kw
-
     def __repr__(self):
         return "FunctionalMagnetism(%s)"%self.name
 
 
-def _set_vars(self, name, profile, kw, reserved):
+def _set_parameters(self, name, profile, kw, reserved):
     # Query profile function for the list of arguments
     vars = inspect.getargspec(profile)[0]
     #print "vars", vars
@@ -237,9 +228,27 @@ def _set_vars(self, name, profile, kw, reserved):
     for k in vars:
         kw.setdefault(k, 0)
     for k, v in kw.items():
-        setattr(self, k, Parameter.default(v, name=name+" "+k))
+        try:
+            pv = [Parameter.default(vi, name=f"{name} {k}[{i}]")
+               for i, vi in enumerate(v)]
+        except TypeError:
+            pv = Parameter.default(v, name=f"{name} {k}")
+        setattr(self, k, pv)
+    self._parameters = vars
+
+def _get_parameters(self):
+    return {k: getattr(self, k) for k in self._parameters}
+
+def _get_values(self):
+    vals = {}
+    for k in self._parameters:
+        v = getattr(self, k)
+        if isinstance(v, list):
+            vals[k] = asarray([vk.value for vk in v], 'd')
+        else:
+            vals[k] = v.value
+    return vals
 
-    return vars
 
 class _LayerLimit(BaseParameter):
     def __init__(self, flayer, isend=True, isrho=True):
@@ -256,7 +265,7 @@ class _LayerLimit(BaseParameter):
     @property
     def value(self):
         z = asarray([0., self.flayer.thickness.value])
-        P = self.flayer.profile(asarray(z), **self.flayer._fpars())
+        P = self.flayer.profile(asarray(z), **_get_values(self.flayer))
         index = 1 if self.isend else 0
         return real(P[index]) if self.isrho else imag(P[index])
 
@@ -279,7 +288,7 @@ class _MagnetismLimit(BaseParameter):
     def value(self):
         zmax = self.flayer._calc_thickness()
         z = asarray([0., zmax])
-        P = self.flayer.profile(z, **self.flayer._fpars())
+        P = self.flayer.profile(z, **_get_values(self.flayer))
         rhoM, thetaM = P if isinstance(P, tuple) else (P, DEFAULT_THETA_M)
         rhoM, thetaM = [broadcast_to(v, z.shape) for v in (rhoM, thetaM)]
         index = -1 if self.isend else 0

[hoogerheide]

@pkienzle Yes you're right there was a logic problem with my model. It sort of worked with the preview and I didn't test it any further than that.

I'm happy with your b5712d7 commit but wouldn't mind testing it out. Where do I find it?

[pkienzle]

I put it in the branch list_aware_flayer on this repo. I don't think I can push to your branch or modify the PR to pull from mine.

[pkienzle]

I think @bmaranville fixed the errors with Test/publish. Can you please merge upstream master into this branch and verify?

[hoogerheide]

I think @bmaranville fixed the errors with Test/publish. Can you please merge upstream master into this branch and verify?

I removed the unnecessary section flagged by @pkienzle. Merge should be complete.

[pkienzle]

Feel free to merge if this is working for you.

I posted an issue to bumps requesting vector parameters: bumps/bumps#73

[hoogerheide]

I think I'd like to do some testing before we merge this in.

[hoogerheide]

Testing looks good. I think it can be merged, but I don't think I have permissions.

[bmaranville]

I think at some point we need to talk about making more concrete classes that support your use case, instead of using the generic FunctionalLayer, since there is application-specific organization of parameters (the functions are not changing radically from model to model, right?)

[hoogerheide]

the functions are not changing radically from model to model, right?

The content of the function and its required parameters can change quite a bit.

[bmaranville]

Is it a countable number of functions? Having different parameters is expected of course.

[hoogerheide]

I'm not sure what you mean by a countable number of functions. Typically the function will be built up from its various components (bilayer, proteins, other stuff). The objects all interact with each other inside the function, sometimes in nontrivial ways. Here are a couple typical functions (requires the molgroups "blm" and "prot" objects of course):

def bilayer(z, sigma, bulknsld, global_rough, rho_substrate,nf_tether, mult_tether, l_tether, l_lipid1, l_lipid2, vf_bilayer):
    """ Generic tBLM bilayer."""
    
    bulknsld = bulknsld * 1e-6
    rho_substrate = rho_substrate * 1e-6

    blm.fnSet(sigma, bulknsld, global_rough, rho_substrate, nf_tether, mult_tether, l_tether, l_lipid1, l_lipid2, vf_bilayer)
    normarea, area, nsl = blm.fnWriteProfile(z)

    # this replaces fnWriteCanvas2Model
    nsld = nsl / (normarea * np.gradient(z)) + (1.0 - area / normarea) * bulknsld

    return nsld * 1e6

def bilayer_protein(z, sigma, bulknsld, global_rough, rho_substrate,nf_tether, mult_tether, l_tether, l_lipid1, l_lipid2, vf_bilayer, dp, vf, nsld_prot_h2o, nsld_prot_d2o, protexchratio, nf_protein):
    """ Generic tBLM bilayer + protein"""
    
    bulknsld = bulknsld * 1e-6
    rho_substrate = rho_substrate * 1e-6
    #nsld_prot_h2o *= 1e-6
    #nsld_prot_d2o *= 1e-6

    blm.fnSet(sigma, bulknsld, global_rough, rho_substrate, nf_tether, mult_tether, l_tether, l_lipid1, l_lipid2, vf_bilayer)

    fracD = protexchratio * (bulknsld + 0.56e-6) / (6.36e-6 + 0.56e-6)
    nsld_prot = fracD * nsld_prot_d2o + (1 - fracD) * nsld_prot_h2o
    prot.fnSetRelative(spline_spacing, blm.lipid1.z - 0.5 * blm.lipid1.l, dp, vf, nsld_prot, nf_protein*vf_bilayer)
    prot.fnSetNormarea(blm.normarea)
    v1 = prot.fnGetVolume(blm.lipid1.z - 0.5*blm.lipid1.l, blm.methyl1.z + 0.5 * blm.methyl1.l)
    v2 = prot.fnGetVolume(blm.methyl2.z - 0.5 * blm.methyl2.l, blm.lipid2.z + 0.5*blm.lipid2.l)
#    print(v1, v2)
    blm.fnSet(sigma, bulknsld, global_rough, rho_substrate, nf_tether, mult_tether, l_tether, l_lipid1, l_lipid2, vf_bilayer, _hc_substitution_1=v1,
              _hc_substitution_2=v2)
    normarea, area, nsl = blm.fnWriteProfile(z)
    area, nsl = prot.fnOverlayProfile(z, area, nsl, normarea)

    nsld = nsl / (normarea * np.gradient(z)) + (1.0 - area / normarea) * bulknsld

    return nsld * 1e6

[bmaranville]

In the future I can imagine specifying models referring to these functions but from a library for membrane protein work... then you can compose layers with algebraic expressions using these well-defined functions.

[pkienzle]

Wasn't that a ship/surf project from last summer? I don't know its current status.