""" BumpsFitting module runs the bumps optimizer. """ import os from datetime import timedelta, datetime import numpy as np from bumps import fitters try: from bumps.options import FIT_CONFIG # Default bumps to use the Levenberg-Marquardt optimizer FIT_CONFIG.selected_id = fitters.LevenbergMarquardtFit.id def get_fitter(): return FIT_CONFIG.selected_fitter, FIT_CONFIG.selected_values except: # CRUFT: Bumps changed its handling of fit options around 0.7.5.6 # Default bumps to use the Levenberg-Marquardt optimizer fitters.FIT_DEFAULT = 'lm' def get_fitter(): fitopts = fitters.FIT_OPTIONS[fitters.FIT_DEFAULT] return fitopts.fitclass, fitopts.options.copy() from bumps.mapper import SerialMapper, MPMapper from bumps import parameter from bumps.fitproblem import FitProblem from sas.sascalc.fit.AbstractFitEngine import FitEngine from sas.sascalc.fit.AbstractFitEngine import FResult from sas.sascalc.fit.expression import compile_constraints class Progress(object): def __init__(self, history, max_step, pars, dof): remaining_time = int(history.time[0]*(float(max_step)/history.step[0]-1)) # Depending on the time remaining, either display the expected # time of completion, or the amount of time remaining. Use precision # appropriate for the duration. if remaining_time >= 1800: completion_time = datetime.now() + timedelta(seconds=remaining_time) if remaining_time >= 36000: time = completion_time.strftime('%Y-%m-%d %H:%M') else: time = completion_time.strftime('%H:%M') else: if remaining_time >= 3600: time = '%dh %dm'%(remaining_time//3600, (remaining_time%3600)//60) elif remaining_time >= 60: time = '%dm %ds'%(remaining_time//60, remaining_time%60) else: time = '%ds'%remaining_time chisq = "%.3g"%(2*history.value[0]/dof) step = "%d of %d"%(history.step[0], max_step) header = "=== Steps: %s chisq: %s ETA: %s\n"%(step, chisq, time) parameters = ["%15s: %-10.3g%s"%(k,v,("\n" if i%3==2 else " | ")) for i,(k,v) in enumerate(zip(pars,history.point[0]))] self.msg = "".join([header]+parameters) def __str__(self): return self.msg class BumpsMonitor(object): def __init__(self, handler, max_step, pars, dof): self.handler = handler self.max_step = max_step self.pars = pars self.dof = dof def config_history(self, history): history.requires(time=1, value=2, point=1, step=1) def __call__(self, history): if self.handler is None: return self.handler.set_result(Progress(history, self.max_step, self.pars, self.dof)) self.handler.progress(history.step[0], self.max_step) if len(history.step)>1 and history.step[1] > history.step[0]: self.handler.improvement() self.handler.update_fit() class ConvergenceMonitor(object): """ ConvergenceMonitor contains population summary statistics to show progress of the fit. This is a list [ (best, 0%, 25%, 50%, 75%, 100%) ] or just a list [ (best, ) ] if population size is 1. """ def __init__(self): self.convergence = [] def config_history(self, history): history.requires(value=1, population_values=1) def __call__(self, history): best = history.value[0] try: p = history.population_values[0] n,p = len(p), np.sort(p) QI,Qmid, = int(0.2*n),int(0.5*n) self.convergence.append((best, p[0],p[QI],p[Qmid],p[-1-QI],p[-1])) except: self.convergence.append((best, best,best,best,best,best)) # Note: currently using bumps parameters for each parameter object so that # a SasFitness can be used directly in bumps with the usual semantics. # The disadvantage of this technique is that we need to copy every parameter # back into the model each time the function is evaluated. We could instead # define reference parameters for each sas parameter, but then we would not # be able to express constraints using python expressions in the usual way # from bumps, and would instead need to use string expressions. class SasFitness(object): """ Wrap SAS model as a bumps fitness object """ def __init__(self, model, data, fitted=[], constraints={}, initial_values=None, **kw): self.name = model.name self.model = model.model self.data = data if self.data.smearer is not None: self.data.smearer.model = self.model self._define_pars() self._init_pars(kw) if initial_values is not None: self._reset_pars(fitted, initial_values) self.constraints = dict(constraints) self.set_fitted(fitted) self.update() def _reset_pars(self, names, values): for k,v in zip(names, values): self._pars[k].value = v def _define_pars(self): self._pars = {} for k in self.model.getParamList(): name = ".".join((self.name,k)) value = self.model.getParam(k) bounds = self.model.details.get(k,["",None,None])[1:3] self._pars[k] = parameter.Parameter(value=value, bounds=bounds, fixed=True, name=name) #print parameter.summarize(self._pars.values()) def _init_pars(self, kw): for k,v in kw.items(): # dispersion parameters initialized with _field instead of .field if k.endswith('_width'): k = k[:-6]+'.width' elif k.endswith('_npts'): k = k[:-5]+'.npts' elif k.endswith('_nsigmas'): k = k[:-7]+'.nsigmas' elif k.endswith('_type'): k = k[:-5]+'.type' if k not in self._pars: formatted_pars = ", ".join(sorted(self._pars.keys())) raise KeyError("invalid parameter %r for %s--use one of: %s" %(k, self.model, formatted_pars)) if '.' in k and not k.endswith('.width'): self.model.setParam(k, v) elif isinstance(v, parameter.BaseParameter): self._pars[k] = v elif isinstance(v, (tuple,list)): low, high = v self._pars[k].value = (low+high)/2 self._pars[k].range(low,high) else: self._pars[k].value = v def set_fitted(self, param_list): """ Flag a set of parameters as fitted parameters. """ for k,p in self._pars.items(): p.fixed = (k not in param_list or k in self.constraints) self.fitted_par_names = [k for k in param_list if k not in self.constraints] self.computed_par_names = [k for k in param_list if k in self.constraints] self.fitted_pars = [self._pars[k] for k in self.fitted_par_names] self.computed_pars = [self._pars[k] for k in self.computed_par_names] # ===== Fitness interface ==== def parameters(self): return self._pars def update(self): for k,v in self._pars.items(): #print "updating",k,v,v.value self.model.setParam(k,v.value) self._dirty = True def _recalculate(self): if self._dirty: self._residuals, self._theory \ = self.data.residuals(self.model.evalDistribution) self._dirty = False def numpoints(self): return np.sum(self.data.idx) # number of fitted points def nllf(self): return 0.5*np.sum(self.residuals()**2) def theory(self): self._recalculate() return self._theory def residuals(self): self._recalculate() return self._residuals # Not implementing the data methods for now: # # resynth_data/restore_data/save/plot class ParameterExpressions(object): def __init__(self, models): self.models = models self._setup() def _setup(self): exprs = {} for M in self.models: exprs.update((".".join((M.name, k)), v) for k, v in M.constraints.items()) if exprs: symtab = dict((".".join((M.name, k)), p) for M in self.models for k,p in M.parameters().items()) self.update = compile_constraints(symtab, exprs) else: self.update = lambda: 0 def __call__(self): self.update() def __getstate__(self): return self.models def __setstate__(self, state): self.models = state self._setup() class BumpsFit(FitEngine): """ Fit a model using bumps. """ def __init__(self): """ Creates a dictionary (self.fit_arrange_dict={})of FitArrange elements with Uid as keys """ FitEngine.__init__(self) self.curr_thread = None def fit(self, msg_q=None, q=None, handler=None, curr_thread=None, ftol=1.49012e-8, reset_flag=False): # Build collection of bumps fitness calculators models = [SasFitness(model=M.get_model(), data=M.get_data(), constraints=M.constraints, fitted=M.pars, initial_values=M.vals if reset_flag else None) for M in self.fit_arrange_dict.values() if M.get_to_fit()] if len(models) == 0: raise RuntimeError("Nothing to fit") problem = FitProblem(models) # TODO: need better handling of parameter expressions and bounds constraints # so that they are applied during polydispersity calculations. This # will remove the immediate need for the setp_hook in bumps, though # bumps may still need something similar, such as a sane class structure # which allows a subclass to override setp. problem.setp_hook = ParameterExpressions(models) # Run the fit result = run_bumps(problem, handler, curr_thread) if handler is not None: handler.update_fit(last=True) # TODO: shouldn't reference internal parameters of fit problem varying = problem._parameters # collect the results all_results = [] for M in problem.models: fitness = M.fitness fitted_index = [varying.index(p) for p in fitness.fitted_pars] param_list = fitness.fitted_par_names + fitness.computed_par_names R = FResult(model=fitness.model, data=fitness.data, param_list=param_list) R.theory = fitness.theory() R.residuals = fitness.residuals() R.index = fitness.data.idx R.fitter_id = self.fitter_id # TODO: should scale stderr by sqrt(chisq/DOF) if dy is unknown R.success = result['success'] if R.success: R.stderr = np.hstack((result['stderr'][fitted_index], np.NaN*np.ones(len(fitness.computed_pars)))) R.pvec = np.hstack((result['value'][fitted_index], [p.value for p in fitness.computed_pars])) R.fitness = np.sum(R.residuals**2)/(fitness.numpoints() - len(fitted_index)) else: R.stderr = np.NaN*np.ones(len(param_list)) R.pvec = np.asarray( [p.value for p in fitness.fitted_pars+fitness.computed_pars]) R.fitness = np.NaN R.convergence = result['convergence'] if result['uncertainty'] is not None: R.uncertainty_state = result['uncertainty'] all_results.append(R) if q is not None: q.put(all_results) return q else: return all_results def run_bumps(problem, handler, curr_thread): def abort_test(): if curr_thread is None: return False try: curr_thread.isquit() except KeyboardInterrupt: if handler is not None: handler.stop("Fitting: Terminated!!!") return True return False fitclass, options = get_fitter() steps = options.get('steps', 0) if steps == 0: pop = options.get('pop',0)*len(problem._parameters) samples = options.get('samples', 0) steps = (samples+pop-1)/pop if pop != 0 else samples max_step = steps + options.get('burn', 0) pars = [p.name for p in problem._parameters] #x0 = np.asarray([p.value for p in problem._parameters]) options['monitors'] = [ BumpsMonitor(handler, max_step, pars, problem.dof), ConvergenceMonitor(), ] fitdriver = fitters.FitDriver(fitclass, problem=problem, abort_test=abort_test, **options) omp_threads = int(os.environ.get('OMP_NUM_THREADS','0')) mapper = MPMapper if omp_threads == 1 else SerialMapper fitdriver.mapper = mapper.start_mapper(problem, None) #import time; T0 = time.time() try: best, fbest = fitdriver.fit() except: import traceback; traceback.print_exc() raise finally: mapper.stop_mapper(fitdriver.mapper) convergence_list = options['monitors'][-1].convergence convergence = (2*np.asarray(convergence_list)/problem.dof if convergence_list else np.empty((0,1),'d')) success = best is not None return { 'value': best if success else None, 'stderr': fitdriver.stderr() if success else None, 'success': success, 'convergence': convergence, 'uncertainty': getattr(fitdriver.fitter, 'state', None), }