source: sasview/src/sas/sascalc/fit/BumpsFitting.py @ 2b3eb3d

ticket-1094-headlessticket-1243
Last change on this file since 2b3eb3d was 1386b2f, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

lint

  • Property mode set to 100644
File size: 14.1 KB
RevLine 
[6fe5100]1"""
2BumpsFitting module runs the bumps optimizer.
3"""
[249a7c6]4import os
[35086c3]5from datetime import timedelta, datetime
[1a5d5f2]6import traceback
[35086c3]7
[9a5097c]8import numpy as np
[6fe5100]9
10from bumps import fitters
[7945367]11try:
12    from bumps.options import FIT_CONFIG
13    # Default bumps to use the Levenberg-Marquardt optimizer
14    FIT_CONFIG.selected_id = fitters.LevenbergMarquardtFit.id
15    def get_fitter():
16        return FIT_CONFIG.selected_fitter, FIT_CONFIG.selected_values
[1386b2f]17except ImportError:
[7945367]18    # CRUFT: Bumps changed its handling of fit options around 0.7.5.6
19    # Default bumps to use the Levenberg-Marquardt optimizer
20    fitters.FIT_DEFAULT = 'lm'
21    def get_fitter():
22        fitopts = fitters.FIT_OPTIONS[fitters.FIT_DEFAULT]
23        return fitopts.fitclass, fitopts.options.copy()
24
25
[249a7c6]26from bumps.mapper import SerialMapper, MPMapper
[e3efa6b3]27from bumps import parameter
28from bumps.fitproblem import FitProblem
[7945367]29
[345e7e4]30
[b699768]31from sas.sascalc.fit.AbstractFitEngine import FitEngine
32from sas.sascalc.fit.AbstractFitEngine import FResult
33from sas.sascalc.fit.expression import compile_constraints
[6fe5100]34
[35086c3]35class Progress(object):
36    def __init__(self, history, max_step, pars, dof):
37        remaining_time = int(history.time[0]*(float(max_step)/history.step[0]-1))
38        # Depending on the time remaining, either display the expected
39        # time of completion, or the amount of time remaining.  Use precision
40        # appropriate for the duration.
41        if remaining_time >= 1800:
42            completion_time = datetime.now() + timedelta(seconds=remaining_time)
43            if remaining_time >= 36000:
44                time = completion_time.strftime('%Y-%m-%d %H:%M')
45            else:
46                time = completion_time.strftime('%H:%M')
47        else:
48            if remaining_time >= 3600:
49                time = '%dh %dm'%(remaining_time//3600, (remaining_time%3600)//60)
50            elif remaining_time >= 60:
51                time = '%dm %ds'%(remaining_time//60, remaining_time%60)
52            else:
53                time = '%ds'%remaining_time
54        chisq = "%.3g"%(2*history.value[0]/dof)
55        step = "%d of %d"%(history.step[0], max_step)
56        header = "=== Steps: %s  chisq: %s  ETA: %s\n"%(step, chisq, time)
57        parameters = ["%15s: %-10.3g%s"%(k,v,("\n" if i%3==2 else " | "))
[1386b2f]58                      for i, (k, v) in enumerate(zip(pars, history.point[0]))]
[35086c3]59        self.msg = "".join([header]+parameters)
60
61    def __str__(self):
62        return self.msg
63
64
[85f17f6]65class BumpsMonitor(object):
[35086c3]66    def __init__(self, handler, max_step, pars, dof):
[85f17f6]67        self.handler = handler
68        self.max_step = max_step
[35086c3]69        self.pars = pars
70        self.dof = dof
[ed4aef2]71
[85f17f6]72    def config_history(self, history):
73        history.requires(time=1, value=2, point=1, step=1)
[ed4aef2]74
[85f17f6]75    def __call__(self, history):
[e3efa6b3]76        if self.handler is None: return
[35086c3]77        self.handler.set_result(Progress(history, self.max_step, self.pars, self.dof))
[85f17f6]78        self.handler.progress(history.step[0], self.max_step)
[1386b2f]79        if len(history.step) > 1 and history.step[1] > history.step[0]:
[85f17f6]80            self.handler.improvement()
81        self.handler.update_fit()
82
[ed4aef2]83class ConvergenceMonitor(object):
84    """
85    ConvergenceMonitor contains population summary statistics to show progress
86    of the fit.  This is a list [ (best, 0%, 25%, 50%, 75%, 100%) ] or
87    just a list [ (best, ) ] if population size is 1.
88    """
89    def __init__(self):
90        self.convergence = []
91
92    def config_history(self, history):
93        history.requires(value=1, population_values=1)
94
95    def __call__(self, history):
96        best = history.value[0]
97        try:
98            p = history.population_values[0]
[1386b2f]99            n, p = len(p), np.sort(p)
100            QI, Qmid = int(0.2*n), int(0.5*n)
101            self.convergence.append((best, p[0], p[QI], p[Qmid], p[-1-QI], p[-1]))
102        except Exception:
103            self.convergence.append((best, best, best, best, best, best))
[ed4aef2]104
[e3efa6b3]105
[4e9f227]106# Note: currently using bumps parameters for each parameter object so that
107# a SasFitness can be used directly in bumps with the usual semantics.
108# The disadvantage of this technique is that we need to copy every parameter
109# back into the model each time the function is evaluated.  We could instead
[fd5ac0d]110# define reference parameters for each sas parameter, but then we would not
[4e9f227]111# be able to express constraints using python expressions in the usual way
112# from bumps, and would instead need to use string expressions.
[e3efa6b3]113class SasFitness(object):
[6fe5100]114    """
[e3efa6b3]115    Wrap SAS model as a bumps fitness object
[6fe5100]116    """
[5044543]117    def __init__(self, model, data, fitted=[], constraints={},
118                 initial_values=None, **kw):
[4e9f227]119        self.name = model.name
120        self.model = model.model
[6fe5100]121        self.data = data
[9f7fbd9]122        if self.data.smearer is not None:
123            self.data.smearer.model = self.model
[e3efa6b3]124        self._define_pars()
125        self._init_pars(kw)
[5044543]126        if initial_values is not None:
127            self._reset_pars(fitted, initial_values)
[4e9f227]128        self.constraints = dict(constraints)
[e3efa6b3]129        self.set_fitted(fitted)
[5044543]130        self.update()
131
132    def _reset_pars(self, names, values):
[1386b2f]133        for k, v in zip(names, values):
[5044543]134            self._pars[k].value = v
[e3efa6b3]135
136    def _define_pars(self):
137        self._pars = {}
138        for k in self.model.getParamList():
[1386b2f]139            name = ".".join((self.name, k))
[e3efa6b3]140            value = self.model.getParam(k)
[1386b2f]141            bounds = self.model.details.get(k, ["", None, None])[1:3]
[e3efa6b3]142            self._pars[k] = parameter.Parameter(value=value, bounds=bounds,
143                                                fixed=True, name=name)
[4e9f227]144        #print parameter.summarize(self._pars.values())
[e3efa6b3]145
146    def _init_pars(self, kw):
[1386b2f]147        for k, v in kw.items():
[e3efa6b3]148            # dispersion parameters initialized with _field instead of .field
[1386b2f]149            if k.endswith('_width'):
150                k = k[:-6]+'.width'
151            elif k.endswith('_npts'):
152                k = k[:-5]+'.npts'
153            elif k.endswith('_nsigmas'):
154                k = k[:-7]+'.nsigmas'
155            elif k.endswith('_type'):
156                k = k[:-5]+'.type'
[e3efa6b3]157            if k not in self._pars:
158                formatted_pars = ", ".join(sorted(self._pars.keys()))
159                raise KeyError("invalid parameter %r for %s--use one of: %s"
160                               %(k, self.model, formatted_pars))
161            if '.' in k and not k.endswith('.width'):
162                self.model.setParam(k, v)
163            elif isinstance(v, parameter.BaseParameter):
164                self._pars[k] = v
[1386b2f]165            elif isinstance(v, (tuple, list)):
[e3efa6b3]166                low, high = v
167                self._pars[k].value = (low+high)/2
[1386b2f]168                self._pars[k].range(low, high)
[95d58d3]169            else:
[e3efa6b3]170                self._pars[k].value = v
171
172    def set_fitted(self, param_list):
[6fe5100]173        """
[e3efa6b3]174        Flag a set of parameters as fitted parameters.
[6fe5100]175        """
[1386b2f]176        for k, p in self._pars.items():
[4e9f227]177            p.fixed = (k not in param_list or k in self.constraints)
[4a0dc427]178        self.fitted_par_names = [k for k in param_list if k not in self.constraints]
[bf5e985]179        self.computed_par_names = [k for k in param_list if k in self.constraints]
180        self.fitted_pars = [self._pars[k] for k in self.fitted_par_names]
181        self.computed_pars = [self._pars[k] for k in self.computed_par_names]
[6fe5100]182
[e3efa6b3]183    # ===== Fitness interface ====
184    def parameters(self):
185        return self._pars
[6fe5100]186
[e3efa6b3]187    def update(self):
[1386b2f]188        for k, v in self._pars.items():
[4e9f227]189            #print "updating",k,v,v.value
[1386b2f]190            self.model.setParam(k, v.value)
[e3efa6b3]191        self._dirty = True
[6fe5100]192
[e3efa6b3]193    def _recalculate(self):
194        if self._dirty:
[9f7fbd9]195            self._residuals, self._theory \
196                = self.data.residuals(self.model.evalDistribution)
[e3efa6b3]197            self._dirty = False
[6fe5100]198
[e3efa6b3]199    def numpoints(self):
[9a5097c]200        return np.sum(self.data.idx) # number of fitted points
[6fe5100]201
[e3efa6b3]202    def nllf(self):
[9a5097c]203        return 0.5*np.sum(self.residuals()**2)
[e3efa6b3]204
205    def theory(self):
206        self._recalculate()
207        return self._theory
208
209    def residuals(self):
210        self._recalculate()
211        return self._residuals
212
213    # Not implementing the data methods for now:
214    #
215    #     resynth_data/restore_data/save/plot
[6fe5100]216
[191c648]217class ParameterExpressions(object):
218    def __init__(self, models):
219        self.models = models
220        self._setup()
221
222    def _setup(self):
223        exprs = {}
224        for M in self.models:
225            exprs.update((".".join((M.name, k)), v) for k, v in M.constraints.items())
226        if exprs:
227            symtab = dict((".".join((M.name, k)), p)
228                          for M in self.models
[1386b2f]229                          for k, p in M.parameters().items())
[191c648]230            self.update = compile_constraints(symtab, exprs)
231        else:
232            self.update = lambda: 0
233
234    def __call__(self):
235        self.update()
236
237    def __getstate__(self):
238        return self.models
239
240    def __setstate__(self, state):
241        self.models = state
242        self._setup()
243
[6fe5100]244class BumpsFit(FitEngine):
245    """
246    Fit a model using bumps.
247    """
248    def __init__(self):
249        """
250        Creates a dictionary (self.fit_arrange_dict={})of FitArrange elements
251        with Uid as keys
252        """
253        FitEngine.__init__(self)
254        self.curr_thread = None
255
256    def fit(self, msg_q=None,
257            q=None, handler=None, curr_thread=None,
258            ftol=1.49012e-8, reset_flag=False):
[e3efa6b3]259        # Build collection of bumps fitness calculators
[bf5e985]260        models = [SasFitness(model=M.get_model(),
261                             data=M.get_data(),
262                             constraints=M.constraints,
[5044543]263                             fitted=M.pars,
264                             initial_values=M.vals if reset_flag else None)
[bf5e985]265                  for M in self.fit_arrange_dict.values()
266                  if M.get_to_fit()]
[233c121]267        if len(models) == 0:
268            raise RuntimeError("Nothing to fit")
[e3efa6b3]269        problem = FitProblem(models)
270
[191c648]271        # TODO: need better handling of parameter expressions and bounds constraints
272        # so that they are applied during polydispersity calculations.  This
273        # will remove the immediate need for the setp_hook in bumps, though
274        # bumps may still need something similar, such as a sane class structure
275        # which allows a subclass to override setp.
276        problem.setp_hook = ParameterExpressions(models)
[4e9f227]277
[e3efa6b3]278        # Run the fit
279        result = run_bumps(problem, handler, curr_thread)
[6fe5100]280        if handler is not None:
281            handler.update_fit(last=True)
[e3efa6b3]282
[eff93b8]283        # TODO: shouldn't reference internal parameters of fit problem
[e3efa6b3]284        varying = problem._parameters
285        # collect the results
286        all_results = []
287        for M in problem.models:
288            fitness = M.fitness
289            fitted_index = [varying.index(p) for p in fitness.fitted_pars]
[e1442d4]290            param_list = fitness.fitted_par_names + fitness.computed_par_names
[e3efa6b3]291            R = FResult(model=fitness.model, data=fitness.data,
[e1442d4]292                        param_list=param_list)
[e3efa6b3]293            R.theory = fitness.theory()
294            R.residuals = fitness.residuals()
[5044543]295            R.index = fitness.data.idx
[e3efa6b3]296            R.fitter_id = self.fitter_id
[eff93b8]297            # TODO: should scale stderr by sqrt(chisq/DOF) if dy is unknown
[e3efa6b3]298            R.success = result['success']
[e1442d4]299            if R.success:
[1a5d5f2]300                if result['stderr'] is None:
[9a5097c]301                    R.stderr = np.NaN*np.ones(len(param_list))
[1a5d5f2]302                else:
[9a5097c]303                    R.stderr = np.hstack((result['stderr'][fitted_index],
304                                          np.NaN*np.ones(len(fitness.computed_pars))))
305                R.pvec = np.hstack((result['value'][fitted_index],
[1386b2f]306                                    [p.value for p in fitness.computed_pars]))
[9a5097c]307                R.fitness = np.sum(R.residuals**2)/(fitness.numpoints() - len(fitted_index))
[e1442d4]308            else:
[9a5097c]309                R.stderr = np.NaN*np.ones(len(param_list))
[1386b2f]310                R.pvec = np.asarray([p.value for p in fitness.fitted_pars+fitness.computed_pars])
[9a5097c]311                R.fitness = np.NaN
[e3efa6b3]312            R.convergence = result['convergence']
313            if result['uncertainty'] is not None:
314                R.uncertainty_state = result['uncertainty']
315            all_results.append(R)
[1a5d5f2]316        all_results[0].mesg = result['errors']
[e3efa6b3]317
[6fe5100]318        if q is not None:
[e3efa6b3]319            q.put(all_results)
[6fe5100]320            return q
[e3efa6b3]321        else:
322            return all_results
[6fe5100]323
[e3efa6b3]324def run_bumps(problem, handler, curr_thread):
[85f17f6]325    def abort_test():
326        if curr_thread is None: return False
327        try: curr_thread.isquit()
328        except KeyboardInterrupt:
329            if handler is not None:
330                handler.stop("Fitting: Terminated!!!")
331            return True
332        return False
333
[7945367]334    fitclass, options = get_fitter()
335    steps = options.get('steps', 0)
336    if steps == 0:
[1386b2f]337        pop = options.get('pop', 0)*len(problem._parameters)
[7945367]338        samples = options.get('samples', 0)
339        steps = (samples+pop-1)/pop if pop != 0 else samples
340    max_step = steps + options.get('burn', 0)
[35086c3]341    pars = [p.name for p in problem._parameters]
[9a5097c]342    #x0 = np.asarray([p.value for p in problem._parameters])
[e3efa6b3]343    options['monitors'] = [
[35086c3]344        BumpsMonitor(handler, max_step, pars, problem.dof),
[e3efa6b3]345        ConvergenceMonitor(),
346        ]
[95d58d3]347    fitdriver = fitters.FitDriver(fitclass, problem=problem,
[042f065]348                                  abort_test=abort_test, **options)
[1386b2f]349    omp_threads = int(os.environ.get('OMP_NUM_THREADS', '0'))
[e1442d4]350    mapper = MPMapper if omp_threads == 1 else SerialMapper
[6fe5100]351    fitdriver.mapper = mapper.start_mapper(problem, None)
[233c121]352    #import time; T0 = time.time()
[6fe5100]353    try:
354        best, fbest = fitdriver.fit()
[1a5d5f2]355        errors = []
356    except Exception as exc:
[9a5097c]357        best, fbest = None, np.NaN
[1a30720]358        errors = [str(exc), traceback.format_exc()]
[95d58d3]359    finally:
360        mapper.stop_mapper(fitdriver.mapper)
[e3efa6b3]361
362
363    convergence_list = options['monitors'][-1].convergence
[9a5097c]364    convergence = (2*np.asarray(convergence_list)/problem.dof
[1386b2f]365                   if convergence_list else np.empty((0, 1), 'd'))
[e1442d4]366
367    success = best is not None
[1a5d5f2]368    try:
369        stderr = fitdriver.stderr() if success else None
370    except Exception as exc:
371        errors.append(str(exc))
372        errors.append(traceback.format_exc())
373        stderr = None
[e3efa6b3]374    return {
[e1442d4]375        'value': best if success else None,
[1a5d5f2]376        'stderr': stderr,
[e1442d4]377        'success': success,
[e3efa6b3]378        'convergence': convergence,
379        'uncertainty': getattr(fitdriver.fitter, 'state', None),
[1a5d5f2]380        'errors': '\n'.join(errors),
[e3efa6b3]381        }
Note: See TracBrowser for help on using the repository browser.