source: sasview/src/sas/fit/BumpsFitting.py @ 386ffe1

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 386ffe1 was 386ffe1, checked in by pkienzle, 9 years ago

remove scipy levenburg marquardt and park from ui

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