""" Build a bumps model from a function and data. Example ------- Given a function *sin_model* which computes a sine wave at times *t*:: from numpy import sin def sin_model(t, freq, phase): return sin(2*pi*(freq*t + phase)) and given data *(y,dy)* measured at times *t*, we can define the fit problem as follows:: from bumps.names import * M = Curve(sin_model, t, y, dy, freq=20) The *freq* and *phase* keywords are optional initial values for the model parameters which otherwise default to zero. The model parameters can be accessed as attributes on the model to set fit range:: M.freq.range(2, 100) M.phase.range(0, 1) As usual, you can initialize or assign parameter expressions to the the parameters if you want to tie parameters together within or between models. Note: there is sometimes difficulty getting bumps to recognize the function during fits, which can be addressed by putting the definition in a separate file on the python path. With the windows binary distribution of bumps, this can be done in the problem definition file with the following code:: import os from bumps.names import * sys.path.insert(0, os.getcwd()) The model function can then be imported from the external module as usual:: from sin_model import sin_model """ __all__ = ["Curve", "PoissonCurve", "plot_err"] import inspect import numpy as np from numpy import log, pi, sqrt from .parameter import Parameter class Curve(object): r""" Model a measurement with a user defined function. The function *fn(x,p1,p2,...)* should return the expected value *y* for each point *x* given the parameters *p1*, *p2*, etc. *dy* is the uncertainty for each measured value *y*. If not specified, it defaults to 1. Initial values for the parameters can be set as *p=value* arguments to *Curve*. If no value is set, then the initial value will be taken from the default value given in the definition of *fn*, or set to 0 if the parameter is not defined with an initial value. Arbitrary non-fittable data can be passed to the function as parameters, but only if the parameter is given a default value of *None* in the function definition, and has the initial value set as an argument to *Curve*. Defining *state=dict(key=value, ...)* before *Curve*, and calling *Curve* as *Curve(..., \*\*state)* works pretty well. *Curve* takes two special keyword arguments: *name* and *plot*. *name* is added to each parameter name when the parameter is defined. The filename for the data is a good choice, since this allows you to keep the parameters straight when fitting multiple datasets simultaneously. Plotting defaults to a 1-D plot with error bars for the data, and a line for the function value. You can assign your own plot function with the *plot* keyword. The function should be defined as *plot(x,y,dy,fy,\*\*kw)*. The keyword arguments will be filled with the values of the parameters used to compute *fy*. It will be easiest to list the parameters you need to make your plot as positional arguments after *x,y,dy,fy* in the plot function declaration. For example, *plot(x,y,dy,fy,p3,\*\*kw)* will make the value of parameter *p3* available as a variable in your function. The special keyword *view* will be a string containing *linear*, *log*, *logx* or *loglog*. The data uncertainty is assumed to follow a gaussian distribution. If measurements draw from some other uncertainty distribution, then subclass Curve and replace nllf with the correct probability given the residuals. See the implementation of :class:`PoissonCurve` for an example. """ def __init__(self, fn, x, y, dy=None, name="", plot=None, **fnkw): self.x, self.y = np.asarray(x), np.asarray(y) if dy is None: self.dy = 1 else: self.dy = np.asarray(dy) if (self.dy <= 0).any(): raise ValueError("measurement uncertainty must be positive") self.fn = fn self.name = name # if name else fn.__name__ + " " # Make every name a parameter; initialize the parameters # with the default value if function is defined with keyword # initializers; override the initializers with any keyword # arguments specified in the fit function constructor. pnames, vararg, varkw, pvalues = inspect.getargspec(fn) if vararg or varkw: raise TypeError( "Function cannot have *args or **kwargs in declaration") # TODO: need "self" handling for passed methods # assume the first argument is x pnames = pnames[1:] # Parameters default to zero init = dict((p, 0) for p in pnames) # If the function provides default values, use those if pvalues: # ignore default value for "x" parameter if len(pvalues) > len(pnames): pvalues = pvalues[1:] init.update(zip(pnames[-len(pvalues):], pvalues)) # Non-fittable parameters need to be sent in as None state_vars = set(p for p, v in init.items() if v is None) # Regardless, use any values specified in the constructor, but first # check that they exist as function parameters. invalid = set(fnkw.keys()) - set(pnames) if invalid: raise TypeError("Invalid initializers: %s" % ", ".join(sorted(invalid))) init.update(fnkw) # Build parameters out of ranges and initial values # maybe: name=(p+name if name.startswith('_') else name+p) pars = dict((p, Parameter.default(init[p], name=name + p)) for p in pnames if p not in state_vars) # Make parameters accessible as model attributes for k, v in pars.items(): if hasattr(self, k): raise TypeError("Parameter cannot be named %s" % k) setattr(self, k, v) # Remember the function, parameters, and number of parameters self._function = fn self._pnames = [p for p in pnames if p not in state_vars] self._cached_theory = None self._plot = plot if plot is not None else plot_err self._state = dict((p, v) for p, v in init.items() if p in state_vars) def update(self): self._cached_theory = None def parameters(self): return dict((p, getattr(self, p)) for p in self._pnames) def numpoints(self): return np.prod(self.y.shape) def theory(self, x=None): if self._cached_theory is None: if x is None: x = self.x kw = dict((p, getattr(self, p).value) for p in self._pnames) kw.update(self._state) self._cached_theory = self._function(x, **kw) return self._cached_theory def simulate_data(self, noise=None): theory = self.theory() if noise is not None: if noise == 'data': pass elif noise < 0: self.dy = -theory*noise*0.01 else: self.dy = noise self.y = theory + np.random.randn(*theory.shape)*self.dy def residuals(self): return (self.theory() - self.y) / self.dy def nllf(self): r = self.residuals() return 0.5 * np.sum(r ** 2) def save(self, basename): # TODO: need header line with state vars as json # TODO: need to support nD x,y,dy data = np.vstack((self.x, self.y, self.dy, self.theory())) np.savetxt(basename + '.dat', data.T) def plot(self, view=None): import pylab kw = dict((p, getattr(self, p).value) for p in self._pnames) kw.update(self._state) #print "kw_plot",kw if view == 'residual': plot_resid(self.x, self.residuals()) else: plot_ratio = 4 h = pylab.subplot2grid((plot_ratio, 1), (0, 0), rowspan=plot_ratio-1) self._plot(self.x, self.y, self.dy, self.theory(), view=view, **kw) for tick_label in pylab.gca().get_xticklabels(): tick_label.set_visible(False) #pylab.gca().xaxis.set_visible(False) #pylab.gca().spines['bottom'].set_visible(False) #pylab.gca().set_xticks([]) pylab.subplot2grid((plot_ratio, 1), (plot_ratio-1, 0), sharex=h) plot_resid(self.x, self.residuals()) def plot_resid(x, resid): import pylab pylab.plot(x, resid, '.') pylab.gca().locator_params(axis='y', tight=True, nbins=4) pylab.axhline(y=1, ls='dotted') pylab.axhline(y=-1, ls='dotted') pylab.ylabel("Residuals") def plot_err(x, y, dy, fy, view=None, **kw): """ Plot data *y* and error *dy* against *x*. *view* is one of linear, log, logx or loglog. """ import pylab pylab.errorbar(x, y, yerr=dy, fmt='.') pylab.plot(x, fy, '-') if view == 'log': pylab.xscale('linear') pylab.yscale('log') elif view == 'logx': pylab.xscale('log') pylab.yscale('linear') elif view == 'loglog': pylab.xscale('log') pylab.yscale('log') else: # view == 'linear' pylab.xscale('linear') pylab.yscale('linear') _LOGFACTORIAL = np.array([log(np.prod(np.arange(1., k + 1))) for k in range(21)]) def logfactorial(n): """Compute the log factorial for each element of an array""" result = np.empty(n.shape, dtype='double') idx = (n <= 20) result[idx] = _LOGFACTORIAL[np.asarray(n[idx], 'int32')] n = n[~idx] result[~idx] = n * \ log(n) - n + log(n * (1 + 4 * n * (1 + 2 * n))) / 6 + log(pi) / 2 return result class PoissonCurve(Curve): r""" Model a measurement with Poisson uncertainty. The nllf is calculated using Poisson probabilities, but the curve itself is displayed using the approximation that $\sigma_y \approx \sqrt(y)$. See :class:`Curve` for details. """ def __init__(self, fn, x, y, name="", **fnkw): dy = sqrt(y) + (y == 0) if y is not None else None Curve.__init__(self, fn, x, y, dy, name=name, **fnkw) self._logfacty = logfactorial(y) if y is not None else None self._logfactysum = np.sum(self._logfacty) ## Assume gaussian residuals for now #def residuals(self): # # TODO: provide individual probabilities as residuals # # or perhaps the square roots --- whatever gives a better feel for # # which points are driving the fit # theory = self.theory() # return np.sqrt(self.y * log(theory) - theory - self._logfacty) def nllf(self): theory = self.theory() if (theory <= 0).any(): return 1e308 return -sum(self.y * log(theory) - theory) + self._logfactysum def simulate_data(self, noise=None): theory = self.theory() self.y = np.random.poisson(theory) self.dy = sqrt(self.y) + (self.y == 0) self._logfacty = logfactorial(self.y) self._logfactysum = np.sum(self._logfacty) def save(self, basename): # TODO: need header line with state vars as json # TODO: need to support nD x,y,dy data = np.vstack((self.x, self.y, self.theory())) np.savetxt(basename + '.dat', data.T)