Changeset 2c4a190 in sasmodels


Ignore:
Timestamp:
Dec 13, 2018 11:08:26 AM (5 weeks ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master
Children:
8803a38
Parents:
c6084f1
Message:

simplify use of multiple scattering in bumps fits

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • example/multiscatfit.py

    r49d1f8b8 r2c4a190  
    1515 
    1616    # Show the model without fitting 
    17     PYTHONPATH=..:../explore:../../bumps:../../sasview/src python multiscatfit.py 
     17    PYTHONPATH=..:../../bumps:../../sasview/src python multiscatfit.py 
    1818 
    1919    # Run the fit 
    20     PYTHONPATH=..:../explore:../../bumps:../../sasview/src ../../bumps/run.py \ 
     20    PYTHONPATH=..:../../bumps:../../sasview/src ../../bumps/run.py \ 
    2121    multiscatfit.py --store=/tmp/t1 
    2222 
     
    5555    ) 
    5656 
     57# Tie the model to the data 
     58M = Experiment(data=data, model=model) 
     59 
     60# Stack mulitple scattering on top of the existing resolution function. 
     61M.resolution = MultipleScattering(resolution=M.resolution, probability=0.) 
     62 
    5763# SET THE FITTING PARAMETERS 
    5864model.radius_polar.range(15, 3000) 
     
    6571model.scale.range(0, 0.1) 
    6672 
    67 # Mulitple scattering probability parameter 
    68 # HACK: the probability is stuffed in as an extra parameter to the experiment. 
    69 probability = Parameter(name="probability", value=0.0) 
    70 probability.range(0.0, 0.9) 
     73# The multiple scattering probability parameter is in the resolution function 
     74# instead of the scattering function, so access it through M.resolution 
     75M.scattering_probability.range(0.0, 0.9) 
    7176 
    72 M = Experiment(data=data, model=model, extra_pars={'probability': probability}) 
    73  
    74 # Stack mulitple scattering on top of the existing resolution function. 
    75 # Because resolution functions in sasview don't have fitting parameters, 
    76 # we instead allow the multiple scattering calculator to take a function 
    77 # instead of a probability.  This function returns the current value of 
    78 # the parameter. ** THIS IS TEMPORARY ** when multiple scattering is 
    79 # properly integrated into sasmodels and sasview, its fittable parameter 
    80 # will be treated like the model parameters. 
    81 M.resolution = MultipleScattering(resolution=M.resolution, 
    82                                   probability=lambda: probability.value, 
    83                                   ) 
    84 M._kernel_inputs = M.resolution.q_calc 
     77# Let bumps know that we are fitting this experiment 
    8578problem = FitProblem(M) 
    8679 
  • sasmodels/bumps_model.py

    r49d1f8b8 r2c4a190  
    3535    # when bumps is not on the path. 
    3636    from bumps.names import Parameter # type: ignore 
     37    from bumps.parameter import Reference # type: ignore 
    3738except ImportError: 
    3839    pass 
     
    139140    def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): 
    140141        # type: (Data, Model, float) -> None 
     142        # Allow resolution function to define fittable parameters.  We do this 
     143        # by creating reference parameters within the resolution object rather 
     144        # than modifying the object itself to use bumps parameters.  We need 
     145        # to reset the parameters each time the object has changed.  These 
     146        # additional parameters need to be returned from the fitting engine. 
     147        # To make them available to the user, they are added as top-level 
     148        # attributes to the experiment object.  The only change to the 
     149        # resolution function is that it needs an optional 'fittable' attribute 
     150        # which maps the internal name to the user visible name for the 
     151        # for the parameter. 
     152        self._resolution = None 
     153        self._resolution_pars = {} 
    141154        # remember inputs so we can inspect from outside 
    142155        self.name = data.filename if name is None else name 
     
    145158        self._interpret_data(data, model.sasmodel) 
    146159        self._cache = {} 
     160        # CRUFT: no longer need extra parameters 
     161        # Multiple scattering probability is now retrieved directly from the 
     162        # multiple scattering resolution function. 
    147163        self.extra_pars = extra_pars 
    148164 
     
    162178        return len(self.Iq) 
    163179 
     180    @property 
     181    def resolution(self): 
     182        return self._resolution 
     183 
     184    @resolution.setter 
     185    def resolution(self, value): 
     186        self._resolution = value 
     187 
     188        # Remove old resolution fitting parameters from experiment 
     189        for name in self._resolution_pars: 
     190            delattr(self, name) 
     191 
     192        # Create new resolution fitting parameters 
     193        res_pars = getattr(self._resolution, 'fittable', {}) 
     194        self._resolution_pars = { 
     195            name: Reference(self._resolution, refname, name=name) 
     196            for refname, name in res_pars.items() 
     197        } 
     198 
     199        # Add new resolution fitting parameters as experiment attributes 
     200        for name, ref in self._resolution_pars.items(): 
     201            setattr(self, name, ref) 
     202 
    164203    def parameters(self): 
    165204        # type: () -> Dict[str, Parameter] 
     
    168207        """ 
    169208        pars = self.model.parameters() 
    170         if self.extra_pars: 
     209        if self.extra_pars is not None: 
    171210            pars.update(self.extra_pars) 
     211        pars.update(self._resolution_pars) 
    172212        return pars 
    173213 
  • sasmodels/direct_model.py

    r7b9e4dd r2c4a190  
    242242            else: 
    243243                Iq, dIq = None, None 
    244             #self._theory = np.zeros_like(q) 
    245             q_vectors = [res.q_calc] 
    246244        elif self.data_type == 'Iqxy': 
    247245            #if not model.info.parameters.has_2d: 
     
    260258            res = resolution2d.Pinhole2D(data=data, index=index, 
    261259                                         nsigma=3.0, accuracy=accuracy) 
    262             #self._theory = np.zeros_like(self.Iq) 
    263             q_vectors = res.q_calc 
    264260        elif self.data_type == 'Iq': 
    265261            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     
    286282            else: 
    287283                res = resolution.Perfect1D(data.x[index]) 
    288  
    289             #self._theory = np.zeros_like(self.Iq) 
    290             q_vectors = [res.q_calc] 
    291284        elif self.data_type == 'Iq-oriented': 
    292285            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     
    304297                                      qx_width=data.dxw[index], 
    305298                                      qy_width=data.dxl[index]) 
    306             q_vectors = res.q_calc 
    307299        else: 
    308300            raise ValueError("Unknown data type") # never gets here 
     
    310302        # Remember function inputs so we can delay loading the function and 
    311303        # so we can save/restore state 
    312         self._kernel_inputs = q_vectors 
    313304        self._kernel = None 
    314305        self.Iq, self.dIq, self.index = Iq, dIq, index 
     
    347338        # type: (ParameterSet, float) -> np.ndarray 
    348339        if self._kernel is None: 
    349             self._kernel = self._model.make_kernel(self._kernel_inputs) 
     340            # TODO: change interfaces so that resolution returns kernel inputs 
     341            # Maybe have resolution always return a tuple, or maybe have 
     342            # make_kernel accept either an ndarray or a pair of ndarrays. 
     343            kernel_inputs = self.resolution.q_calc 
     344            if isinstance(kernel_inputs, np.ndarray): 
     345                kernel_inputs = (kernel_inputs,) 
     346            self._kernel = self._model.make_kernel(kernel_inputs) 
    350347 
    351348        # Need to pull background out of resolution for multiple scattering 
  • sasmodels/multiscat.py

    rb3703f5 r2c4a190  
    342342 
    343343    *probability* is related to the expected number of scattering 
    344     events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a 
    345     hack to allow probability to be a fitted parameter, the "value" 
    346     can be a function that takes no parameters and returns the current 
    347     value of the probability.  *coverage* determines how many scattering 
    348     steps to consider.  The default is 0.99, which sets $n$ such that 
    349     $1 \ldots n$ covers 99% of the Poisson probability mass function. 
     344    events in the sample $\lambda$ as $p = 1 - e^{-\lambda}$. 
     345    *coverage* determines how many scattering steps to consider.  The 
     346    default is 0.99, which sets $n$ such that $1 \ldots n$ covers 99% 
     347    of the Poisson probability mass function. 
    350348 
    351349    *is2d* is True then 2D scattering is used, otherwise it accepts 
     
    399397        self.qmin = qmin 
    400398        self.nq = nq 
    401         self.probability = probability 
     399        self.probability = 0. if probability is None else probability 
    402400        self.coverage = coverage 
    403401        self.is2d = is2d 
     
    456454        self.Iqxy = None # type: np.ndarray 
    457455 
     456        # Label probability as a fittable parameter, and give its external name 
     457        # Note that the external name must be a valid python identifier, since 
     458        # is will be set as an experiment attribute. 
     459        self.fittable = {'probability': 'scattering_probability'} 
     460 
    458461    def apply(self, theory): 
    459462        if self.is2d: 
     
    463466        Iq_calc = Iq_calc.reshape(self.nq, self.nq) 
    464467 
     468        # CRUFT: don't need probability as a function anymore 
    465469        probability = self.probability() if callable(self.probability) else self.probability 
    466470        coverage = self.coverage 
Note: See TracChangeset for help on using the changeset viewer.