Changeset 2c4a190 in sasmodels

Timestamp:

Dec 13, 2018 9:08:26 AM (6 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

8803a38

Parents:

c6084f1

Message:

simplify use of multiple scattering in bumps fits

Files:

• example/multiscatfit.py (modified) (3 diffs)
• sasmodels/bumps_model.py (modified) (5 diffs)
• sasmodels/direct_model.py (modified) (6 diffs)
• sasmodels/multiscat.py (modified) (4 diffs)

example/multiscatfit.py

@@ -15 +15 @@
 
     # Show the model without fitting
-    PYTHONPATH=..:../explore:../../bumps:../../sasview/src python multiscatfit.py
+    PYTHONPATH=..:../../bumps:../../sasview/src python multiscatfit.py
 
     # Run the fit
-    PYTHONPATH=..:../explore:../../bumps:../../sasview/src ../../bumps/run.py \
+    PYTHONPATH=..:../../bumps:../../sasview/src ../../bumps/run.py \
     multiscatfit.py --store=/tmp/t1
 
@@ -55 +55 @@
     )
 
+# Tie the model to the data
+M = Experiment(data=data, model=model)
+
+# Stack mulitple scattering on top of the existing resolution function.
+M.resolution = MultipleScattering(resolution=M.resolution, probability=0.)
+
 # SET THE FITTING PARAMETERS
 model.radius_polar.range(15, 3000)
@@ -65 +71 @@
 model.scale.range(0, 0.1)
 
-# Mulitple scattering probability parameter
-# HACK: the probability is stuffed in as an extra parameter to the experiment.
-probability = Parameter(name="probability", value=0.0)
-probability.range(0.0, 0.9)
+# The multiple scattering probability parameter is in the resolution function
+# instead of the scattering function, so access it through M.resolution
+M.scattering_probability.range(0.0, 0.9)
 
-M = Experiment(data=data, model=model, extra_pars={'probability': probability})
-
-# Stack mulitple scattering on top of the existing resolution function.
-# Because resolution functions in sasview don't have fitting parameters,
-# we instead allow the multiple scattering calculator to take a function
-# instead of a probability.  This function returns the current value of
-# the parameter. ** THIS IS TEMPORARY ** when multiple scattering is
-# properly integrated into sasmodels and sasview, its fittable parameter
-# will be treated like the model parameters.
-M.resolution = MultipleScattering(resolution=M.resolution,
-                                  probability=lambda: probability.value,
-                                  )
-M._kernel_inputs = M.resolution.q_calc
+# Let bumps know that we are fitting this experiment
 problem = FitProblem(M)
 

sasmodels/bumps_model.py

@@ -35 +35 @@
     # when bumps is not on the path.
     from bumps.names import Parameter # type: ignore
+    from bumps.parameter import Reference # type: ignore
 except ImportError:
     pass
@@ -139 +140 @@
     def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None):
         # type: (Data, Model, float) -> None
+        # Allow resolution function to define fittable parameters.  We do this
+        # by creating reference parameters within the resolution object rather
+        # than modifying the object itself to use bumps parameters.  We need
+        # to reset the parameters each time the object has changed.  These
+        # additional parameters need to be returned from the fitting engine.
+        # To make them available to the user, they are added as top-level
+        # attributes to the experiment object.  The only change to the
+        # resolution function is that it needs an optional 'fittable' attribute
+        # which maps the internal name to the user visible name for the
+        # for the parameter.
+        self._resolution = None
+        self._resolution_pars = {}
         # remember inputs so we can inspect from outside
         self.name = data.filename if name is None else name
@@ -145 +158 @@
         self._interpret_data(data, model.sasmodel)
         self._cache = {}
+        # CRUFT: no longer need extra parameters
+        # Multiple scattering probability is now retrieved directly from the
+        # multiple scattering resolution function.
         self.extra_pars = extra_pars
 
@@ -162 +178 @@
         return len(self.Iq)
 
+    @property
+    def resolution(self):
+        return self._resolution
+
+    @resolution.setter
+    def resolution(self, value):
+        self._resolution = value
+
+        # Remove old resolution fitting parameters from experiment
+        for name in self._resolution_pars:
+            delattr(self, name)
+
+        # Create new resolution fitting parameters
+        res_pars = getattr(self._resolution, 'fittable', {})
+        self._resolution_pars = {
+            name: Reference(self._resolution, refname, name=name)
+            for refname, name in res_pars.items()
+        }
+
+        # Add new resolution fitting parameters as experiment attributes
+        for name, ref in self._resolution_pars.items():
+            setattr(self, name, ref)
+
     def parameters(self):
         # type: () -> Dict[str, Parameter]
@@ -168 +207 @@
         """
         pars = self.model.parameters()
-        if self.extra_pars:
+        if self.extra_pars is not None:
             pars.update(self.extra_pars)
+        pars.update(self._resolution_pars)
         return pars
 

sasmodels/direct_model.py

@@ -242 +242 @@
             else:
                 Iq, dIq = None, None
-            #self._theory = np.zeros_like(q)
-            q_vectors = [res.q_calc]
         elif self.data_type == 'Iqxy':
             #if not model.info.parameters.has_2d:
@@ -260 +258 @@
             res = resolution2d.Pinhole2D(data=data, index=index,
                                          nsigma=3.0, accuracy=accuracy)
-            #self._theory = np.zeros_like(self.Iq)
-            q_vectors = res.q_calc
         elif self.data_type == 'Iq':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
@@ -286 +282 @@
             else:
                 res = resolution.Perfect1D(data.x[index])
-
-            #self._theory = np.zeros_like(self.Iq)
-            q_vectors = [res.q_calc]
         elif self.data_type == 'Iq-oriented':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
@@ -304 +297 @@
                                       qx_width=data.dxw[index],
                                       qy_width=data.dxl[index])
-            q_vectors = res.q_calc
         else:
             raise ValueError("Unknown data type") # never gets here
@@ -310 +302 @@
         # Remember function inputs so we can delay loading the function and
         # so we can save/restore state
-        self._kernel_inputs = q_vectors
         self._kernel = None
         self.Iq, self.dIq, self.index = Iq, dIq, index
@@ -347 +338 @@
         # type: (ParameterSet, float) -> np.ndarray
         if self._kernel is None:
-            self._kernel = self._model.make_kernel(self._kernel_inputs)
+            # TODO: change interfaces so that resolution returns kernel inputs
+            # Maybe have resolution always return a tuple, or maybe have
+            # make_kernel accept either an ndarray or a pair of ndarrays.
+            kernel_inputs = self.resolution.q_calc
+            if isinstance(kernel_inputs, np.ndarray):
+                kernel_inputs = (kernel_inputs,)
+            self._kernel = self._model.make_kernel(kernel_inputs)
 
         # Need to pull background out of resolution for multiple scattering

sasmodels/multiscat.py

@@ -342 +342 @@
 
     *probability* is related to the expected number of scattering
-    events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a
-    hack to allow probability to be a fitted parameter, the "value"
-    can be a function that takes no parameters and returns the current
-    value of the probability.  *coverage* determines how many scattering
-    steps to consider.  The default is 0.99, which sets $n$ such that
-    $1 \ldots n$ covers 99% of the Poisson probability mass function.
+    events in the sample $\lambda$ as $p = 1 - e^{-\lambda}$.
+    *coverage* determines how many scattering steps to consider.  The
+    default is 0.99, which sets $n$ such that $1 \ldots n$ covers 99%
+    of the Poisson probability mass function.
 
     *is2d* is True then 2D scattering is used, otherwise it accepts
@@ -399 +397 @@
         self.qmin = qmin
         self.nq = nq
-        self.probability = probability
+        self.probability = 0. if probability is None else probability
         self.coverage = coverage
         self.is2d = is2d
@@ -456 +454 @@
         self.Iqxy = None # type: np.ndarray
 
+        # Label probability as a fittable parameter, and give its external name
+        # Note that the external name must be a valid python identifier, since
+        # is will be set as an experiment attribute.
+        self.fittable = {'probability': 'scattering_probability'}
+
     def apply(self, theory):
         if self.is2d:
@@ -463 +466 @@
         Iq_calc = Iq_calc.reshape(self.nq, self.nq)
 
+        # CRUFT: don't need probability as a function anymore
         probability = self.probability() if callable(self.probability) else self.probability
         coverage = self.coverage