Index: park-1.2.1/park/__init__.py
===================================================================
--- park-1.2.1/park/__init__.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/__init__.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,70 @@
+# This program is public domain
+"""
+PARK fitting service.
+
+The PARK fitting service is a set of python packages to support
+fitting in datasets.  Using common python infrastructure such
+as scipy optimizers, numpy arrays, and matplotlib plotting, park
+allows you to define models, associate them with datasets and
+simultaneously fit them.  Park provides a simple job queue to
+manage multiple fits using separate processes or running across
+the network on separate machines.
+
+Installation
+============
+
+The latest version of Park is available from 
+http://www.reflectometry.org/danse/park.
+
+
+Currently this is supplied as source from a zip file.  You can
+also retrieve the latest version from svn::
+
+    svn co svn://danse.us/park/branches/park-1.2
+
+If you are installing from source, you will need a python
+environment with numpy and scipy.  For the GUI version, you
+will also need matplotlib and wx.
+
+You will need a C compiler to build the resolution convolution 
+function.  On Windows this may require installing MinGW and
+adding distutils.cfg to your distutils directory::
+
+    [build]
+    compiler = mingw
+
+Once you have the required supporting packages, use the following
+to build and install::
+
+    python setup.py install
+    
+Usage
+=====
+
+To get started with park, you will need to first define the
+models that you are using.  These can be very basic models, 
+listing all possible fitting parameters.  When setting up the 
+fit later, you will be able to combine models, using
+parameter expressions to relate the values in one model with
+those in another.  See `park.model` for details.
+
+Once your models are constructed you can use them in a fit.
+See `park.fit` for details.
+
+Important classes and functions:
+`park.model.Model`, `park.parameter.ParameterSet`, `park.fit.Fit`
+
+:group models: model, assembly, parameter, data
+:group examples: peaks
+:group optimizer: fit, fitresult, fitmc, simplex
+:group support: expression, deps, version, setup, serial
+:group server: fitservice
+"""
+
+#from fitservice import *
+from parameter import *
+from data import *
+from model import *
+from assembly import *
+import fit
+
Index: park-1.2.1/park/assembly.py
===================================================================
--- park-1.2.1/park/assembly.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/assembly.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,639 @@
+# This program is public domain
+"""
+An assembly is a collection of fitting functions.  This provides
+the model representation that is the basis of the park fitting engine.
+
+Models can range from very simple one dimensional theory functions
+to complex assemblies of multidimensional datasets from different
+experimental techniques, each with their own theory function and
+a common underlying physical model.
+
+Usage
+=====
+
+First define the models you want to work with.  In the example
+below we will use an example of a simple multilayer system measured by
+specular reflection of xrays and neutrons.  The gold depth is the only
+fitting parameter, ranging from 10-30 A.  The interface depths are
+tied together using expressions.  In this case the expression is
+a simple copy, but any standard math functions can be used.  Some
+model developers may provide additional functions for use with the
+expression.
+
+Example models::
+
+    import reflectometry.model1d as refl
+    xray = refl.model('xray')
+    xray.incident('Air',rho=0)
+    xray.interface('iAu',sigma=5)
+    xray.layer('Au',rho=124.68,depth=[10,30])
+    xray.interface('iSi',sigma=5)
+    xray.substrate('Si',rho=20.07)
+    datax = refl.data('xray.dat')
+
+    neutron = refl.model('neutron')
+    neutron.incident('Air',rho=0)
+    neutron.interface('iAu',sigma='xray.iAu')
+    neutron.layer('Au',rho=4.66,depth='xray.Au.depth')
+    neutron.interface('iSi',sigma='xray.iSi')
+    neutron.substrate('Si',rho=2.07)
+    datan = refl.data('neutron.dat')
+
+As you can see from the above, parameters can be set to a value if
+the parameter is fixed, to a range if the parametemr is fitted, or
+to a string expression if the parameter is calculated from other
+parameters.  See park.Parameter.set for further details.
+
+Having constructed the models, we can now create an assembly::
+
+    import park
+    assembly = park.Assembly([(xray,datax), (neutron,datan)])
+
+Note: this would normally be done in the context of a fit
+using fit = park.Fit([(xray,datax), (neutron,datan)]), and later referenced
+using fit.assembly.
+
+Individual parts of the assembly are accessable using the
+model number 0, 1, 2... or by the model name.  In the above,
+assembly[0] and assembly['xray'] refer to the same model.
+Assemblies have insert and append functions for adding new
+models, and "del model[idx]" for removing them.
+
+Once the assembly is created computing the values for the system
+is a matter of calling::
+
+    assembly.eval()
+    print "Chi**2",assembly.chisq
+    print "Reduced chi**2",assembly.chisq/assembly.degrees_of_freedom
+    plot(arange(len(assembly.residuals)), assembly.residuals)
+
+This defines the attributes residuals, degrees_of_freedom and chisq,
+which is what the optimizer uses as the cost function to minimize.
+
+assembly.eval uses the current values for the parameters in the
+individual models.  These parameters can be changed directly
+in the model.  In the reflectometry example above, you could
+set the gold thickness using xray.layer.Au.depth=156, or
+something similar (the details are model specific).  Parameters
+can also be changed through the assembly parameter set.  In the same
+example, this would be assembly.parameterset['xray']['Au']['depth'].
+See parameter set for details.
+
+In the process of modeling data, particularly with multiple
+datasets, you will sometimes want to temporarily ignore
+how well one of the datasets matches so that you
+can more quickly refine the model for the other datasets,
+or see how particular models are influencing the fit.  To
+temporarily ignore the xray data in the example above use::
+
+    assembly.parts[0].isfitted = False
+
+The model itself isn't ignored since its parameters may be
+needed to compute the parameters for other models.  To
+reenable checking against the xray data, you would assign
+a True value instead.  More subtle weighting of the models
+can be controlled using assembly.parts[idx].weight, but
+see below for a note on model weighting.
+
+A note on model weighting
+-------------------------
+
+Changing the weight is equivalent to scaling the error bars
+on the given model by a factor of weight/n where n is the
+number of data points.  It is better to set the correct error
+bars on your data in the first place than to adjust the weights.
+If you have the correct error bars, then you should expect
+roughly 2/3 of the data points to lie within one error bar of
+the theory curve.  If consecutive data points have largely
+overlapping errorbars, then your uncertainty is overestimated.
+
+Another case where weights are adjusted (abused?) is to
+compensate for systematic errors in the data by forcing the
+errorbars to be large enough to cover the systematic bias.
+This is a poor approach to the problem.  A better strategy
+is to capture the systematic effects in the model, and treat
+the measurement of the independent variable as an additional
+data point in the fit.  This is still not statistically sound
+as there is likely to be a large correlation between the
+uncertainty of the measurement and the values of all the
+other variables.
+
+That said, adjusting the weight on a dataset is a quick way
+of reducing its influence on the entire fit.  Please use it
+with care.
+"""
+
+__all__ = ['Assembly', 'Fitness']
+import numpy
+
+import park
+from park.parameter import Parameter,ParameterSet
+from park.fitresult import FitParameter
+import park.expression
+
+
+
+class Fitness(object):
+    """
+    Container for theory and data.
+
+    The fit object compares theory with data.
+
+    TODO: what to do with fittable metadata (e.g., footprint correction)?
+    """
+    data = None
+    model = None
+    def __init__(self, model=None,data=None):
+        self.data,self.model = data,model
+    def _parameterset(self):
+        return self.model.parameterset
+    parameterset = property(_parameterset)
+    def residuals(self):
+        return self.data.residuals(self.model.eval)
+    def residuals_deriv(self, pars=[]):
+        return self.data.residuals_deriv(self.model.eval_derivs,pars=pars)
+    def set(self, **kw):
+        """
+        Set parameters in the model.
+
+        User convenience function.  This allows a user with an assembly
+        of models in a script to for example set the fit range for
+        parameter 'a' of the model::
+            assembly[0].set(a=[5,6])
+
+        Raises KeyError if the parameter is not in parameterset.
+        """
+        self.model.set(**kw)
+    def abort(self):
+        if hasattr(self.model,'abort'): self.model.abort()
+
+class Part(object):
+    """
+    Part of a fitting assembly.  Part holds the model itself and
+    associated data.  The part can be initialized with a fitness
+    object or with a pair (model,data) for the default fitness function.
+
+    fitness (Fitness)
+        object implementing the `park.assembly.Fitness` interface.  In
+        particular, fitness should provide a parameterset attribute
+        containing a ParameterSet and a residuals method returning a vector
+        of residuals.
+    weight (dimensionless)
+        weight for the model.  See comments in assembly.py for details.
+    isfitted (boolean)
+        True if the model residuals should be included in the fit.
+        The model parameters may still be used in parameter
+        expressions, but there will be no comparison to the data.
+    residuals (vector)
+        Residuals for the model if they have been calculated, or None
+    degrees_of_freedom
+        Number of residuals minus number of fitted parameters.
+        Degrees of freedom for individual models does not make
+        sense in the presence of expressions combining models,
+        particularly in the case where a model has many parameters
+        but no data or many computed parameters.  The degrees of
+        freedom for the model is set to be at least one.
+    chisq
+        sum(residuals**2); use chisq/degrees_of_freedom to
+        get the reduced chisq value.
+
+        Get/set the weight on the given model.
+
+        assembly.weight(3) returns the weight on model 3 (0-origin)
+        assembly.weight(3,0.5) sets the weight on model 3 (0-origin)
+    """
+
+    def __init__(self, fitness, weight=1., isfitted=True):
+        if isinstance(fitness, tuple):
+            fitness = park.Fitness(*fitness)
+        self.fitness = fitness
+        self.weight = weight
+        self.isfitted = isfitted
+        self.residuals = None
+        self.chisq = numpy.Inf
+        self.degrees_of_freedom = 1
+
+
+class Assembly(object):
+    """
+    Collection of fit models.
+
+    Assembly implements the `park.fit.Objective` interface.
+
+    See `park.assembly` for usage.
+
+    Instance variables:
+
+    residuals : array
+        a vector of residuals spanning all models, with model
+        weights applied as appropriate.
+    degrees_of_freedom : integer
+        length of the residuals - number of fitted parameters
+    chisq : float
+        sum squared residuals; this is not the reduced chisq, which
+        you can get using chisq/degrees_of_freedom
+
+    These fields are defined for the individual models as well, with
+    degrees of freedom adjusted to the length of the individual data
+    set.  If the model is not fitted or the weight is zero, the residual
+    will not be calculated.
+
+    The residuals fields are available only after the model has been
+    evaluated.
+    """
+
+    def __init__(self, models=[]):
+        """Build an assembly from a list of models."""
+        self.parts = []
+        for m in models:
+            self.parts.append(Part(m))
+        self._reset()
+
+    def __iter__(self):
+        """Iterate through the models in order"""
+        for m in self.parts: yield m
+
+    def __getitem__(self, n):
+        """Return the nth model"""
+        return self.parts[n].fitness
+
+    def __setitem__(self, n, fitness):
+        """Replace the nth model"""
+        self.parts[n].fitness = fitness
+        self._reset()
+
+    def __delitem__(self, n):
+        """Delete the nth model"""
+        del self.parts[n]
+        self._reset()
+
+    def weight(self, idx, value=None):
+        """
+        Query the weight on a particular model.
+
+        Set weight to value if value is supplied.
+
+        :Parameters:
+         idx : integer
+           model number
+         value : float
+           model weight
+        :return: model weight
+        """
+        if value is not None:
+            self.parts[idx].weight = value
+        return self.parts[idx].weight
+
+    def isfitted(self, idx, value=None):
+        """
+        Query if a particular model is fitted.
+
+        Set isfitted to value if value is supplied.
+
+        :param idx: model number
+        :type idx: integer
+        :param value:
+        """
+        if value is not None:
+            self.parts[idx].isfitted = value
+        return self.parts[idx].isfitted
+
+    def append(self, fitness, weight=1.0, isfitted=True):
+        """
+        Add a model to the end of set.
+
+        :param fitness: the fitting model
+        The fitting model can be an instance of `park.assembly.Fitness`,
+        or a tuple of (`park.model.Model`,`park.data.Data1D`)
+        :param weight: model weighting (usually 1.0)
+        :param isfitted: whether model should be fit (equivalent to weight 0.)
+        """
+        self.parts.append(Part(fitness,weight,isfitted))
+        self._reset()
+
+    def insert(self, idx, fitness, weight=1.0, isfitted=True):
+        """Add a model to a particular position in the set."""
+        self.parts.insert(idx,Part(fitness,weight,isfitted))
+        self._reset()
+
+    def _reset(self):
+        """Adjust the parameter set after the addition of a new model."""
+        subsets = [m.fitness.parameterset for m in self]
+        self.parameterset = ParameterSet('root',subsets)
+        self.parameterset.setprefix()
+        #print [p.path for p in self.parameterset.flatten()]
+
+    def eval(self):
+        """
+        Recalculate the theory functions, and from them, the
+        residuals and chisq.
+
+        :note: Call this after the parameters have been updated.
+        """
+        # Handle abort from a separate thread.
+        self._cancel = False
+
+        # Evaluate the computed parameters
+        self._fitexpression()
+
+        # Check that the resulting parameters are in a feasible region.
+        if not self.isfeasible(): return numpy.inf
+
+        resid = []
+        k = len(self._fitparameters)
+        for m in self.parts:
+            # In order to support abort, need to be able to propagate an
+            # external abort signal from self.abort() into an abort signal
+            # for the particular model.  Can't see a way to do this which
+            # doesn't involve setting a state variable.
+            self._current_model = m
+            if self._cancel: return numpy.inf
+            if m.isfitted and m.weight != 0:
+                m.residuals = m.fitness.residuals()
+                N = len(m.residuals)
+                m.degrees_of_freedom = N-k if N>k else 1
+                m.chisq = numpy.sum(m.residuals**2)
+                resid.append(m.weight*m.residuals)
+        self.residuals = numpy.hstack(resid)
+        N = len(self.residuals)
+        self.degrees_of_freedom = N-k if N>k else 1
+        self.chisq = numpy.sum(self.residuals**2)
+        return self.chisq
+
+    def jacobian(self, pvec, step=1e-8):
+        """
+        Returns the derivative wrt the fit parameters at point p.
+
+        Numeric derivatives are calculated based on step, where step is
+        the portion of the total range for parameter j, or the portion of
+        point value p_j if the range on parameter j is infinite.
+        """
+        # Make sure the input vector is an array
+        pvec = numpy.asarray(pvec)
+        # We are being lazy here.  We can precompute the bounds, we can
+        # use the residuals_deriv from the sub-models which have analytic
+        # derivatives and we need only recompute the models which depend
+        # on the varying parameters.
+        # Meanwhile, let's compute the numeric derivative using the
+        # three point formula.
+        # We are not checking that the varied parameter in numeric
+        # differentiation is indeed feasible in the interval of interest.
+        range = zip(*[p.range for p in self._fitparameters])
+        lo,hi = [numpy.asarray(v) for v in range]
+        delta = (hi-lo)*step
+        # For infinite ranges, use p*1e-8 for the step size
+        idx = numpy.isinf(delta)
+        #print "J",idx,delta,pvec,type(idx),type(delta),type(pvec)
+        delta[idx] = pvec[idx]*step
+        delta[delta==0] = step
+
+        # Set the initial value
+        for k,v in enumerate(pvec):
+            self._fitparameters[k].value = v
+        # Gather the residuals
+        r = []
+        for k,v in enumerate(pvec):
+            # Center point formula:
+            #     df/dv = lim_{h->0} ( f(v+h)-f(v-h) ) / ( 2h )
+            h = delta[k]
+            self._fitparameters[k].value = v + h
+            self.eval()
+            rk = self.residuals
+            self._fitparameters[k].value = v - h
+            self.eval()
+            rk -= self.residuals
+            self._fitparameters[k].value = v
+            r.append(rk/(2*h))
+        # return the jacobian
+        return numpy.vstack(r).T
+
+
+    def cov(self, pvec):
+        """
+        Return the covariance matrix inv(J'J) at point p.
+        """
+
+        # Find cov of f at p
+        #     cov(f,p) = inv(J'J)
+        # Use SVD
+        #     J = U S V'
+        #     J'J = (U S V')' (U S V')
+        #         = V S' U' U S V'
+        #         = V S S V'
+        #     inv(J'J) = inv(V S S V')
+        #              = inv(V') inv(S S) inv(V)
+        #              = V inv (S S) V'
+        J = self.jacobian(pvec)
+        u,s,vh = numpy.linalg.svd(J,0)
+        JTJinv = numpy.dot(vh.T.conj()/s**2,vh)
+        return JTJinv
+
+    def stderr(self, pvec):
+        """
+        Return parameter uncertainty.
+
+        This is just the sqrt diagonal of covariance matrix inv(J'J) at point p.
+        """
+        return numpy.sqrt(numpy.diag(self.cov(pvec)))
+
+    def isfeasible(self):
+        """
+        Returns true if the parameter set is in a feasible region of the
+        modeling space.
+        """
+        return True
+
+    # Fitting service interface
+    def fit_parameters(self):
+        """
+        Return an alphabetical list of the fitting parameters.
+
+        This function is called once at the beginning of a fit,
+        and serves as a convenient place to precalculate what
+        can be precalculated such as the set of fitting parameters
+        and the parameter expressions evaluator.
+        """
+        self.parameterset.setprefix()
+        self._fitparameters = self.parameterset.fitted
+        self._restraints = self.parameterset.restrained
+        pars = self.parameterset.flatten()
+        context = self.parameterset.gather_context()
+        self._fitexpression = park.expression.build_eval(pars,context)
+        #print "constraints",self._fitexpression.__doc__
+
+        self._fitparameters.sort(lambda a,b: cmp(a.path,b.path))
+        # Convert to fitparameter a object
+        fitpars = [FitParameter(p.path,p.range,p.value)
+                   for p in self._fitparameters]
+        return fitpars
+
+    def set_result(self, result):
+        """
+        Set the parameters resulting from the fit into the parameter set,
+        and update the calculated expression.
+
+        The parameter values may be retrieved by walking the assembly.parameterset
+        tree, checking each parameter for isfitted, iscomputed, or isfixed.
+        For example::
+
+            assembly.set_result(result)
+            for p in assembly.parameterset.flatten():
+                if p.isfitted():
+                    print "%s %g in [%g,%g]"%(p.path,p.value,p.range[0],p.range[1])
+                elif p.iscomputed():
+                    print "%s computed as %g"%(p.path.p.value)
+
+        This does not calculate the function or the residuals for these parameters.
+        You can call assembly.eval() to do this.  The residuals will be set in
+        assembly[i].residuals.  The theory and data are model specific, and can
+        be found in assembly[i].fitness.data.
+        """
+        for n,p in enumerate(result.parameters):
+            self._fitparameters[n] = p.value
+        self._fitexpression()
+
+    def all_results(self, result):
+        """
+        Extend result from the fit with the calculated parameters.
+        """
+        calcpars = [FitParameter(p.path,p.range,p.value)
+                    for p in self.parameterset.computed]
+        result.parameters += calcpars
+
+    def result(self, status='step'):
+        """
+        Details to send back to the fitting client on an improved fit.
+
+        status is 'start', 'step' or 'end' depending if this is the
+        first result to return, an improved result, or the final result.
+
+        [Not implemented]
+        """
+        return None
+
+    def fresiduals(self, pvec):
+        chisq = self.__call__(pvec)
+        return self.residuals
+
+    def __call__(self, pvec):
+        """
+        Cost function.
+
+        Evaluate the system for the parameter vector pvec, returning chisq
+        as the cost function to be minimized.
+
+        Raises a runtime error if the number of fit parameters is
+        different than the length of the vector.
+        """
+        # Plug fit parameters into model
+        #print "Trying",pvec
+        pars = self._fitparameters
+        if len(pvec) != len(pars):
+            raise RuntimeError("Unexpected number of parameters")
+        for n,value in enumerate(pvec):
+            pars[n].value = value
+        # Evaluate model
+        chisq = self.eval()
+        # Evaluate additional restraints based on parameter value
+        # likelihood
+        restraints_penalty = 0
+        for p in self._restraints:
+            restraints_penalty += p.likelihood(p.value)
+        # Return total cost function
+        return self.chisq + restraints_penalty
+
+    def abort(self):
+        """
+        Interrupt the current function evaluation.
+
+        Forward this to the currently executing model if possible.
+        """
+        self._cancel = True
+        if hasattr(self._current_model,'abort'):
+            self._current_model.abort()
+
+class _Exp(park.Model):
+    """
+    Sample model for testing assembly.
+    """
+    parameters = ['a','c']
+    def eval(self,x):
+        return self.a*numpy.exp(self.c*x)
+class _Linear(park.Model):
+    parameters = ['a','c']
+    def eval(self,x):
+        #print "eval",self.a,self.c,x,self.a*x+self.c
+        return self.a*x+self.c
+def example():
+    """
+    Return an example assembly consisting of a pair of functions,
+        M1.a*exp(M1.c*x), M2.a*exp(2*M1.c*x)
+    and ideal data for
+        M1.a=1, M1.c=1.5, M2.a=2.5
+    """
+    import numpy
+    import park
+    from numpy import inf
+    # Make some fake data
+    x1 = numpy.linspace(0,1,11)
+    x2 = numpy.linspace(0,1,12)
+    # Define a shared model
+    if True: # Exp model
+        y1,y2 = numpy.exp(1.5*x1),2.5*numpy.exp(3*x2)
+        M1 = _Exp('M1',a=[1,3],c=[1,3])
+        M2 = _Exp('M2',a=[1,3],c='2*M1.c')
+        #M2 = _Exp('M2',a=[1,3],c=3)
+    else:  # Linear model
+        y1,y2 = x1+1.5, 2.5*x2+3
+        M1 = _Linear('M1',a=[1,3],c=[1,3])
+        M2 = _Linear('M2',a=[1,3],c='2*M1.c')
+    if False: # Unbounded
+        M1.a = [-inf,inf]
+        M1.c = [-inf,inf]
+        M2.a = [-inf,inf]
+    D1 = park.Data1D(x=x1, y=y1)
+    D2 = park.Data1D(x=x2, y=y2)
+    # Construct the assembly
+    assembly = park.Assembly([(M1,D1),(M2,D2)])
+    return assembly
+
+class _Sphere(park.Model):
+    parameters = ['a','b','c','d','e']
+    def eval(self,x):
+        return self.a*x**2+self.b*x+self.c + exp(self.d) - 3*sin(self.e)
+
+def example5():
+    import numpy
+    import park
+    from numpy import inf
+    # Make some fake data
+    x = numpy.linspace(0,1,11)
+    # Define a shared model
+    S = _Sphere(a=1,b=2,c=3,d=4,e=5)
+    y = S.eval(x1)
+    Sfit = _Sphere(a=[-inf,inf],b=[-inf,inf],c=[-inf,inf],d=[-inf,inf],e=[-inf,inf])
+    D = park.Data1D(x=x, y=y)
+    # Construct the assembly
+    assembly = park.Assembly([(Sfit,D)])
+    return assembly
+
+def test():
+    assembly = example()
+    assert assembly[0].parameterset.name == 'M1'
+
+    # extract the fitting parameters
+    pars = [p.name for p in assembly.fit_parameters()]
+    assert set(pars) == set(['M1.a','M1.c','M2.a'])
+    # Compute chisq and verify constraints are updated properly
+    assert assembly([1,1.5,2.5]) == 0
+    assert assembly[0].model.c == 1.5 and assembly[1].model.c == 3
+
+    # Try without constraints
+    assembly[1].set(c=3)
+    assembly.fit_parameters()  # Fit parameters have changed
+    assert assembly([1,1.5,2.5]) == 0
+
+    # Check that assembly.cov runs ... still need to check that it is correct!
+    C = assembly.cov(numpy.array([1,1.5,2.5]))
+
+if __name__ == "__main__": test()
Index: park-1.2.1/park/data.py
===================================================================
--- park-1.2.1/park/data.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/data.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,392 @@
+# This program is public domain
+"""
+Park 1-D data objects.
+
+The class Data1D represents simple 1-D data objects, along with
+an ascii file loader.  This format will work well for many
+uses, but it is likely that more specialized models will have
+their own data file formats.
+
+The minimal data format for park must supply the following
+methods:
+
+    residuals(fn)
+        returns the residuals vector for the model function.
+    residuals_deriv(fn_deriv,par)
+        returns the residuals vector for the model function, and
+        for the derivatives with respect to the given parameters.
+
+The function passed is going to be model.eval or in the case
+where derivatives are available, model.eval_deriv.  Normally
+this will take a vector of dependent variables and return the
+theory function for that vector but this is only convention.
+The fitting service only uses the parameter set and the residuals
+method from the model.
+
+The park GUI will make more demands on the interface, but the
+details are not yet resolved.
+"""
+from __future__ import with_statement  # Used only in test()
+
+__all__ = ['Data1D']
+
+import numpy
+
+try:
+    from park._modeling import convolve as _convolve
+    def convolve(xi,yi,x,dx):
+        """
+        Return convolution y of width dx at points x based on the
+        sampled input function yi = f(xi).
+        """
+        y = numpy.empty(x.shape,'d')
+        xi,yi,x,dx = [numpy.ascontiguousarray(v,'d') for v in xi,yi,x,dx]
+        _convolve(xi,yi,x,dx,y)
+        return y
+except:
+    def convolve(*args,**kw):
+        """
+        Return convolution y of width dx at points x based on the
+        sampled input function yi = f(xi).
+        
+        Note: C version is not available in this build, and python 
+        version is not implemented.
+        """
+        raise NotImplementedError("convolve is a compiled function")
+    
+
+__all__ = ['Data1D']
+    
+class Data1D(object):
+    """
+    Data representation for 1-D fitting.
+
+    Attributes
+    ==========
+    
+    filename
+        The source of the data.  This may be the empty string if the
+        data is simulation data.
+    x,y,dy
+        The data values.
+        x is the measurement points of data to be fitted. x must be sorted.
+        y is the measured value
+        dy is the measurement uncertainty.
+    dx
+        Resolution at the the measured points.  The resolution may be 0, 
+        constant, or defined for each data point.  dx is the 1-sigma
+        width of the Gaussian resolution function at point x.  Note that 
+        dx_FWHM = sqrt(8 ln 2) dx_sigma, so scale dx appropriately.
+
+
+    fit_x,fit_dx,fit_y,fit_dy
+        The points used in evaluating the residuals.
+    calc_x
+        The points at which to evaluate the theory function.  This may be 
+        different from the measured points for a number of reasons, such 
+        as a resolution function which suggests over or under sampling of 
+        the points (see below).  By default calc_x is x, but it can be set
+        explicitly by the user.
+    calc_y, fx
+        The value of the function at the theory points, and the value of
+        the function after resolution has been applied.  These values are
+        computed by a call to residuals.
+
+    Notes on calc_x
+    ===============   
+    
+    The contribution of Q to a resolution of width dQo at point Qo is::
+
+       p(Q) = 1/sqrt(2 pi dQo**2) exp ( (Q-Qo)**2/(2 dQo**2) )
+
+    We are approximating the convolution at Qo using a numerical
+    approximation to the integral over the measured points, with the 
+    integral is limited to p(Q_i)/p(0)>=0.001.  
+
+    Sometimes the function we are convoluting is rapidly changing.
+    That means the correct convolution should uniformly sample across
+    the entire width of the Gaussian.  This is not possible at the
+    end points unless you calculate the theory function beyond what is
+    strictly needed for the data. For a given dQ and step size,
+    you need enough steps that the following is true::
+
+        (n*step)**2 > -2 dQ**2 * ln 0.001
+
+    The choice of sampling density is particularly important near
+    critical points where the shape of the function changes.  In
+    reflectometry, the function goes from flat below the critical
+    edge to O(Q**4) above.  In one particular model, calculating every 
+    0.005 rather than every 0.02 changed a value above the critical 
+    edge by 15%.  In a fitting program, this would lead to a somewhat
+    larger estimate of the critical edge for this sample.
+
+    Sometimes the theory function is oscillating more rapidly than
+    the instrument can resolve.  This happens for example in reflectivity
+    measurements involving thick layers.  In these systems, the theory
+    function should be oversampled around the measured points Q.  With 
+    a single thick layer, oversampling can be limited to just one 
+    period 2 pi/d.  With multiple thick layers, oscillations will 
+    show interference patterns and it will be necessary to oversample 
+    uniformly through the entire width of the resolution.  If this is
+    not accommodated, then aliasing effects make it difficult to
+    compute the correct model.
+    """
+    filename = ""
+    x,y = None,None
+    dx,dy = 0,1
+    calc_x,calc_y = None,None
+    fit_x,fit_y = None,None
+    fit_dx,fit_dy = 0,1
+    fx = None
+    
+    def __init__(self,filename="",x=None,y=None,dx=0,dy=1):
+        """
+        Define the fitting data.
+        
+        Data can be loaded from a file using filename or
+        specified directly using x,y,dx,dy.  File loading
+        happens after assignment of x,y,dx,dy.
+        """
+        self.x,self.y,self.dx,self.dy = x,y,dx,dy
+        self.select(None)
+        if filename: self.load(filename)
+    
+    def resample(self,minstep=None):
+        """
+        Over/under sampling support.
+
+        Compute the calc_x points required to adequately sample
+        the function y=f(x) so that the value reported for each
+        measured point is supported by the resolution.  minstep
+        is the minimum allowed sampling density that should be
+        used.
+        """
+        self.calc_x = resample(self.x,self.dx,minstep)
+
+    def load(self, filename, **kw):
+        """
+        Load a multicolumn datafile.
+        
+        Data should be in columns, with the following defaults::
+        
+            x,y or x,y,dy or x,dx,y,dy
+            
+        Note that this resets the selected fitting points calc_x and the
+        computed results calc_y and fx.
+
+        Data is sorted after loading.
+        
+        Any extra keyword arguments are passed to the numpy loadtxt
+        function.  This allows you to select the columns you want,
+        skip rows, set the column separator, change the comment character,
+        amongst other things.
+        """
+        self.dx,self.dy = 0,1
+        self.calc_x,self.calc_y,self.fx = None,None,None
+        if filename:
+            columns = numpy.loadtxt(filename, **kw)
+            self.x = columns[:,0]
+            if columns.shape[1]==4:
+                self.dx = columns[:,1]
+                self.y = columns[:,2]
+                self.dy = columns[:,3]
+            elif columns.shape[1]==3:
+                self.dx = 0
+                self.y = columns[:,1]
+                self.dy = columns[:,2]
+            elif columns.shape[1]==2:
+                self.dx = 0
+                self.y = columns[:,1]
+                self.dy = 1
+            else:
+                raise IOError,"Unexpected number of columns in "+filename
+            idx = numpy.argsort(self.x)
+            self.x,self.y = self.x[idx],self.y[idx]
+            if not numpy.isscalar(self.dx): self.dx = self.dx[idx]
+            if not numpy.isscalar(self.dy): self.dy = self.dy[idx]
+        else:
+            self.x,self.dx,self.y,self.dy = None,None,None,None
+        
+        self.filename = filename
+        self.select(None)
+
+    def select(self, idx):
+        """
+        A selection vector for points to use in the evaluation of the 
+        residuals, or None if all points are to be used.
+        """
+        if idx is not None:
+            self.fit_x = self.x[idx]
+            self.fit_y = self.y[idx]
+            if not numpy.isscalar(self.dx): self.fit_dx = self.dx[idx]
+            if not numpy.isscalar(self.dy): self.fit_dy = self.dy[idx]
+        else:
+            self.fit_x = self.x
+            self.fit_dx = self.dx
+            self.fit_y = self.y
+            self.fit_dy = self.dy
+ 
+    def residuals(self, fn):
+        """
+        Compute the residuals of the data wrt to the given function.
+
+        y = fn(x) should be a callable accepting a list of points at which
+        to calculate the function, returning the values at those
+        points.
+
+        Any resolution function will be applied after the theory points
+        are calculated.
+        """
+        calc_x = self.calc_x if self.calc_x else self.fit_x
+        self.calc_y = fn(calc_x)
+        self.fx = resolution(calc_x,self.calc_y,self.fit_x,self.fit_dx)
+        return (self.fit_y - self.fx)/self.fit_dy
+
+    def residuals_deriv(self, fn, pars=[]):
+        """
+        Compute residuals and derivatives wrt the given parameters.
+
+        fdf = fn(x,pars=pars) should be a callable accepting a list 
+        of points at which to calculate the function and a keyword 
+        argument listing the parameters for which the derivative will
+        be calculated.
+
+        Returns a list of the residuals and the derivative wrt the
+        parameter for each parameter.
+
+        Any resolution function will be applied after the theory points
+        and derivatives are calculated.
+        """
+        calc_x = self.calc_x if self.calc_x else self.fit_x
+
+        # Compute f and derivatives
+        fdf = fn(calc_x,pars=pars)
+        self.calc_y = fdf[0]
+
+        # Apply resolution
+        fdf = [resolution(calc_x,y,self.fit_x,self.fit_dx) for y in fdf]
+        self.fx = fdf[0]
+        delta = (self.fx-self.fit_x)/self.fit_dy
+        raise RuntimeError('check whether we want df/dp or dR/dp where R=residuals^2')
+
+        # R = (F(x;p)-y)/sigma => dR/dp  = 1/sigma dF(x;p)/dp
+        # dR^2/dp = 2 R /sigma dF(x;p)/dp 
+        df = [ v/self.fit_dy for v in fdf_res[1:] ]
+
+        return [delta]+df
+
+def equal_spaced(x,tol=1e-5):
+    """
+    Return true if data is regularly spaced within tolerance.  Tolerance
+    uses relative error.
+    """
+    step = numpy.diff(x)
+    step_bar = numpy.mean(step)
+    return (abs(step-step_bar) < tol*step_bar).all()
+
+def resample(x,dx,minstep):
+    """
+    Defining the minimum support basis.
+
+    Compute the calc_x points required to adequately sample
+    the function y=f(x) so that the value reported for each
+    measured point is supported by the resolution.  minstep
+    is the minimum allowed sampling density that should be used.
+    """
+    raise NotImplementedError
+
+def resolution(calcx,calcy,fitx,fitdx):
+    """
+    Apply resolution function.  If there is no resolution function, then
+    interpolate from the calculated points to the desired theory points.
+    If the data are irregular, use a brute force convolution function.
+    
+    If the data are regular and the resolution is fixed, then you can
+    deconvolute the data before fitting, saving time and complexity.
+    """
+    if numpy.any(fitdx != 0):
+        if numpy.isscalar(fitdx):
+            fitdx = numpy.ones(fitx.shape)*fitdx
+        fx = convolve(calcx, calcy, fitx, fitdx)
+    elif calcx is fitx:
+        fx = calcy
+    else:
+        fx = numpy.interp(fitx,calcx,calcy)
+    return fx
+
+
+class TempData:
+    """
+    Create a temporary file with a given data set and remove it when done.
+
+    Example::
+
+        from __future__ import with_statement
+        ...
+        with TempData("1 2 3\n4 5 6") as filename:
+            # run tests of loading from filename
+
+    This class is useful for testing data file readers.
+    """
+    def __init__(self,data,suffix='.dat',prefix='',text=True):
+        import os,tempfile
+        fid,self.filename = tempfile.mkstemp('.dat',prefix,text=True)
+        os.write(fid,data)
+        os.close(fid)
+    def __enter__(self):
+        return self.filename
+    def __exit__(self, exc_type, exc_value, traceback):
+        import os
+        os.unlink(self.filename)
+
+D2 = "# x y\n1 1\n2 2\n3 4\n2 5\n"
+"""x,y dataset for TempData"""
+D3 = "# x y dy\n1 1 .1\n2 2 .2\n3 4 .4\n2 5 .5\n"
+"""x,y,dy dataset for TempData"""
+D4 = "# x dx y dy\n1 .1 1 .1\n2 .2 2 .2\n3 .3 4 .4\n2 .3 5 .5\n"
+"""x,dx,y,dy dataset for TempData"""
+    
+def test():
+    import os
+    import numpy
+
+    # Check that two column data loading works
+    with TempData(D2) as filename:
+        data = Data1D(filename=filename)
+    assert numpy.all(data.x == [1,2,2,3])
+    assert numpy.all(data.y == [1,2,5,4])
+    assert data.dx == 0
+    assert data.dy == 1
+    assert data.calc_x is None
+    assert data.residuals(lambda x: x)[3] == 1
+    
+    # Check that interpolation works
+    data.calc_x = [1,1.5,3]
+    assert data.residuals(lambda x: x)[1] == 0
+    assert numpy.all(data.calc_y == data.calc_x) 
+    # Note: calc_y is updated by data.residuals, so be careful with this test
+
+    # Check that three column data loading works
+    with TempData(D3) as filename:
+        data = Data1D(filename=filename)
+    assert numpy.all(data.x == [1,2,2,3])
+    assert numpy.all(data.y == [1,2,5,4])
+    assert numpy.all(data.dy == [.1,.2,.5,.4])
+    assert data.dx == 0
+    assert data.calc_x is None
+    assert data.residuals(lambda x: x)[3] == 1/.4
+
+    # Check that four column data loading works
+    with TempData(D4) as filename:
+        data = Data1D(filename=filename)
+    assert numpy.all(data.x == [1,2,2,3])
+    assert numpy.all(data.dx == [.1,.2,.3,.3])
+    assert numpy.all(data.y == [1,2,5,4])
+    assert numpy.all(data.dy == [.1,.2,.5,.4])
+
+    # Test residuals
+    print "Fix the convolution function!"
+    print data.residuals(lambda x: x)
+    assert data.calc_x is None
+
+if __name__ == "__main__": test()
Index: park-1.2.1/park/deps.py
===================================================================
--- park-1.2.1/park/deps.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/deps.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,94 @@
+# This program is public domain
+# Author: Paul Kienzle
+"""
+Functions for manipulating dependencies.
+
+Parameter values must be updated in the correct order.  If parameter A
+depends on parameter B, then parameter B must be evaluated first.  
+"""
+
+def order_dependencies(pairs):
+    """
+    Order elements from pairs so that b comes before a in the
+    ordered list for all pairs (a,b).
+    """
+    #print "order_dependencies",pairs
+    emptyset = set()
+    order = []
+
+    # Break pairs into left set and right set
+    left,right = [set(s) for s in zip(*pairs)] if pairs != [] else ([],[])
+    while pairs != []:
+        #print "within",pairs
+        # Find which items only occur on the right
+        independent = right - left
+        if independent == emptyset:
+            cycleset = ", ".join(str(s) for s in left)
+            raise ValueError,"Cyclic dependencies amongst %s"%cycleset
+
+        # The possibly resolvable items are those that depend on the independents
+        dependent = set([a for a,b in pairs if b in independent])
+        pairs = [(a,b) for a,b in pairs if b not in independent]
+        if pairs == []:
+            resolved = dependent
+        else:
+            left,right = [set(s) for s in zip(*pairs)]
+            resolved = dependent - left
+        #print "independent",independent,"dependent",dependent,"resolvable",resolved
+        order += resolved
+        #print "new order",order
+    order.reverse()
+    return order
+
+# ========= Test code ========
+def _check(msg,pairs):
+    """
+    Verify that the list n contains the given items, and that the list
+    satisfies the partial ordering given by the pairs in partial order.
+    """
+    left,right = zip(*pairs) if pairs != [] else ([],[])
+    items = set(left)
+    n = order_dependencies(pairs)
+    if set(n) != items or len(n) != len(items):
+        n.sort()
+        items = list(items); items.sort()
+        raise Exception,"%s expect %s to contain %s for %s"%(msg,n,items,pairs)
+    for lo,hi in pairs:
+        if lo in n and hi in n and n.index(lo) >= n.index(hi):
+            raise Exception,"%s expect %s before %s in %s for %s"%(msg,lo,hi,n,pairs)
+
+def test():
+    import numpy
+
+    # Null case
+    _check("test empty",[])
+
+    # Some dependencies
+    _check("test1",[(2,7),(1,5),(1,4),(2,1),(3,1),(5,6)])
+    _check("test1 renumbered",[(6,1),(7,3),(7,4),(6,7),(5,7),(3,2)])
+    _check("test1 numpy",numpy.array([(2,7),(1,5),(1,4),(2,1),(3,1),(5,6)]))
+
+    # No dependencies
+    _check("test2",[(4,1),(3,2),(8,4)])
+
+    # Cycle test
+    pairs = [(1,4),(4,3),(4,5),(5,1)] 
+    try: n = order_dependencies(pairs)
+    except ValueError: pass
+    else: raise Exception,"test3 expect ValueError exception for %s"%(pairs,)
+
+    # large test for gross speed check
+    A = numpy.random.randint(4000,size=(1000,2))
+    A[:,1] += 4000  # Avoid cycles
+    _check("test-large",A)
+
+    # depth tests
+    k = 200
+    A = numpy.array([range(0,k),range(1,k+1)]).T
+    _check("depth-1",A)
+
+    A = numpy.array([range(1,k+1),range(0,k)]).T
+    _check("depth-2",A)
+
+if __name__ == "__main__":
+    test()
Index: park-1.2.1/park/diffev.py
===================================================================
--- park-1.2.1/park/diffev.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/diffev.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,221 @@
+from random import randint, random
+import numpy
+
+class DE:
+    """
+    Differential evolution test implementation
+    Implements the Scheme_DE_rand_1 variant
+
+       F: float
+           weighting factor which determines the magnitude of perturbation
+           occurring in each generation.
+
+       crossover_rate:  float
+           In general, 0 < crossover_rate < 1.  Usually
+           considerably lower than 1.  Convergence slows as the value
+           increases, but higher crossover rates may be necessary when
+
+       func: w = f(p)
+           The function to be minimized, of the form w = f(p), where p
+           is a vector of length func_dim and w is a scalar
+
+       func_dim: int
+           The dimension of the objective function
+
+       pop: array
+           The initial population.  This should be an iterable composed of
+           Rank-1 numpy arrays.  The population size must be at least 4,
+           preferably much greater.
+
+       l_bound, u_bound: vector
+           arrays of size func_dim representing the upper and lower bounds
+           for the parameters in the solution vectors
+
+       tol: float
+           if (max(pop_values) - min(pop_values) <= conv), the population
+           has converged and the evolution will stop
+       """
+
+    def __init__(self, F, crossover_rate, func, func_dim, pop, l_bound,
+                 u_bound, tol=1e-7):
+        self.pop_size          = len(pop)
+        self.dimension         = func_dim
+        self.F_orig            = F
+        self.F                 = F
+        self.crossover_rate    = crossover_rate
+        self.func              = func
+        self.tol               = tol
+        self.l_bound           = l_bound
+        self.u_bound           = u_bound
+        self.population        = pop
+        self.generations       = 0
+        self.pop_values        = [self.func(p) for p in self.population]
+        self.best_gene,self.best_value = self.get_best_gene()
+        self.best_gene_history = [self.best_gene]
+        self.best_val_history  = [self.best_value]
+        self.ncalls            = 0
+
+    #//////////////////////////////////////////////////
+    def evolve(self, numgens=1000, monitor=None):
+        '''Evolve the population for numgens generations, or until it converges.
+            Returns the best vector from the run'''
+        try:
+            import psyco
+            psyco.full()
+        except ImportError:
+            pass
+
+        start_gen = self.generations
+        for gen in xrange(self.generations + 1, self.generations + numgens + 1):
+            candidate_pop = []
+
+            for i in xrange(self.pop_size):
+                (r1, r2, r3) = self.select_participants(i)
+
+                #perturbation/mutation
+                candidate = self.population[r1] + self.F*(self.population[r2]
+                                                          - self.population[r3])
+
+                #crossover
+                candidate = self.crossover(candidate, i)
+
+                #check bound constraints
+                if not self.is_within_bounds(candidate):
+                    candidate = self.get_random_gene()
+
+                #test fitness to determine membership in next gen
+                candidate_val = self.func(candidate)
+                if candidate_val < self.pop_values[i]:
+                    candidate_pop.append(candidate)
+                    self.pop_values[i] = float(candidate_val)
+                else:
+                    candidate_pop.append(self.population[i])
+            self.ncalls += self.pop_size
+
+            self.population  = candidate_pop
+            x,fx = self.get_best_gene()
+            if fx < self.best_value:
+                self.best_gene, self.best_value = x,fx
+                if monitor is not None:
+                    monitor.improvement(x,fx,self.ncalls)
+            self.best_val_history.append(self.best_value)
+            self.best_gene_history.append(self.best_gene)
+            self.generations = gen
+
+            if monitor is not None:
+                monitor.progress(gen-start_gen,numgens)
+
+            if self.is_converged():
+                break
+
+        return self.best_gene
+
+
+    #////////////////////////////////////////////////
+    def get_random_gene(self):
+        '''Generate a random gene within the bounds constraints.
+           Used to replace out-of-bounds genes'''
+        result = [numpy.random.uniform(low, high)
+                  for low, high in zip(self.l_bound, self.u_bound)]
+        return numpy.array(result)
+
+    #////////////////////////////////////////////////
+    def is_within_bounds(self, gene):
+        '''Determine whether the gene meets the bounds constraints'''
+        return numpy.all( (self.l_bound < gene) & (self.u_bound > gene) )
+
+    #////////////////////////////////////////////////
+    #This appears to be a total failure.  I'll leave the code in,
+    #but it's not called from the main evolution loop anymore.
+    def adjust_F(self):
+        '''Adjust F to account for stagnation of the population '''
+        #has the best vector improved this generation?
+        idx = len(self.best_val_history)
+        if self.best_val_history[idx - 1] == self.best_val_history[idx - 2]:
+            self.stagnant_gens += 1
+        else:
+            self.stagnant_gens  = 0
+
+        #adapt F to account for stagnation
+        self.F = min( (self.F_orig + (self.stagnant_gens * 0.01)), 2.0)
+
+    #////////////////////////////////////////////////
+    def get_best_gene(self):
+        '''find the most fit gene in the current population'''
+        #print "pop", self.pop_values
+
+        best_index = numpy.argmin(self.pop_values)
+        return self.population[best_index],self.pop_values[best_index],
+
+    #////////////////////////////////////////////////
+    def select_participants(self, i):
+        '''generate r1, r2, and r3 randomly from the range [0, NP-1]
+            such that they are distinct values not equal to i'''
+        choices = [i]
+        for choice in xrange(3):
+            while True:
+                j = numpy.random.randint(0, self.pop_size-1)
+                if j not in choices:
+                    break
+            choices.append(j)
+        return choices[1:]
+    #////////////////////////////////////////////////
+    def crossover(self, candidate, i):
+        '''Perform a crossover between the candidate and the ith member of
+            the previous generation.'''
+        result = []
+
+        #generate lower bound of crossover (this bound is inclusive)
+        low = randint(0, self.dimension-1)
+
+        #determine the size of the crossover
+        L = 0
+        while random() < self.crossover_rate and L < (self.dimension - low):
+            L += 1
+
+        #calculate the upper bound of crossover (this bound is exclusive)
+        high = low + L
+
+        high  = numpy.repeat(high, self.dimension)
+        low   = numpy.repeat(low,  self.dimension)
+        seq   = numpy.arange(0, self.dimension)
+
+        result = numpy.choose( ((seq >= low)&(seq < high)),
+                             [candidate, self.population[i]] )
+        return result
+
+    def is_converged(self):
+        '''check for convergence'''
+        return max(self.pop_values) - min(self.pop_values) <= self.tol
+
+def diffev(fn, bounds, x0=None, pop_scale=4, crossover_rate=0.8,
+           Fscale=1, tolerance=1e-5, maxiter=1000, monitor=None):
+    """
+    Run differential evolution, returning x,fx,ncalls
+    """
+    lo, hi = numpy.asarray(bounds[0]), numpy.asarray(bounds[1])
+    dim = len(lo)
+    pop = gen_pop(dim*pop_scale, lo, hi, dim)
+    if x0 is not None:
+        pop[0] = numpy.asarray(x0)
+    evolver = DE(Fscale, crossover_rate, fn, dim, pop, lo, hi, tolerance)
+    evolver.evolve(maxiter, monitor=monitor)
+    return evolver.best_gene, evolver.best_value, evolver.ncalls
+
+
+
+
+########################################################################
+#end DE def
+########################################################################
+def gen_pop(size, l_bound, u_bound, dimension):
+    '''generate a random population of vectors within the given bounds.  dimension
+       indicates the length of the vectors.  l_bound and u_bound should be vectors
+       of length dimension (any iterable container should work)'''
+    result = []
+    for i in xrange(size):
+        entry = []
+        for j in xrange(dimension):
+            entry.append( ((u_bound[j] - l_bound[j]) * random()) + l_bound[j] )
+        result.append(numpy.array(entry))
+    return result
Index: park-1.2.1/park/expression.py
===================================================================
--- park-1.2.1/park/expression.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/expression.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,260 @@
+# This program is public domain
+"""
+Functions for manipulating expressions.   
+"""
+import math
+import re
+from deps import order_dependencies
+
+# simple pattern which matches symbols.  Note that it will also match
+# invalid substrings such as a3...9, but given syntactically correct
+# input it will only match symbols.
+_symbol_pattern = re.compile('([a-zA-Z][a-zA-Z_0-9.]*)')
+
+def symbols(expr,symtab):
+    """
+    Given an expression string and a symbol table, return the set of symbols
+    used in the expression.  Symbols are only returned once even if they
+    occur multiple times.  The return value is a set with the elements in
+    no particular order.
+    
+    This is the first step in computing a dependency graph.
+    """
+    matches = [m.group(0) for m in _symbol_pattern.finditer(expr)]
+    return set([symtab[m] for m in matches if m in symtab])
+
+def substitute(expr,mapping):
+    """
+    Replace all occurrences of symbol s with mapping[s] for s in mapping.
+    """
+    # Find the symbols and the mapping
+    matches = [(m.start(),m.end(),mapping[m.group(1)])
+               for m in _symbol_pattern.finditer(expr)
+               if m.group(1) in mapping]
+    
+    # Split the expression in to pieces, with new symbols replacing old
+    pieces = []
+    offset = 0
+    for start,end,text in matches:
+        pieces += [expr[offset:start],text]
+        offset = end
+    pieces.append(expr[offset:])
+    
+    # Join the pieces and return them
+    return "".join(pieces)
+
+def find_dependencies(pars):
+    """
+    Returns a list of pair-wise dependencies from the parameter expressions.
+    
+    For example, if p3 = p1+p2, then find_dependencies([p1,p2,p3]) will
+    return [(p3,p1),(p3,p2)].  For base expressions without dependencies,
+    such as p4 = 2*pi, this should return [(p4, None)]
+    """
+    symtab = dict([(p.path, p) for p in pars])
+    # Hack to deal with expressions without dependencies --- return a fake
+    # dependency of None.  
+    # The better solution is fix order_dependencies so that it takes a 
+    # dictionary of {symbol: dependency_list}, for which no dependencies 
+    # is simply []; fix in parameter_mapping as well
+    def symbols_or_none(expr,symtab):
+        syms = symbols(expr,symtab)
+        return syms if len(syms) else [None]
+    deps = [(p,dep) 
+            for p in pars if p.iscomputed()
+            for dep in symbols_or_none(p.expression,symtab)]
+    return deps
+
+def parameter_mapping(pairs):
+    """
+    Find the parameter substitution we need so that expressions can
+    be evaluated without having to traverse a chain of 
+    model.layer.parameter.value
+    """
+    left,right = zip(*pairs)
+    pars = set(left+right)
+    symtab = dict( ('P%d'%i,p) for i,p in enumerate(pars) )
+    # p is None when there is an expression with no dependencies
+    mapping = dict( (p.path,'P%d.value'%i) 
+                    for i,p in enumerate(pars) 
+                    if p is not None)
+    return symtab,mapping
+
+def no_constraints(): 
+    """
+    This parameter set has no constraints between the parameters.
+    """
+    pass
+
+def build_eval(pars, context={}):
+    """
+    Build and return a function to evaluate all parameter expressions in
+    the proper order.
+    
+    Inputs:
+        pars is a list of parameters
+        context is a dictionary of additional symbols for the expressions
+
+    Output:
+        updater function
+
+    Raises:
+       AssertionError - model, parameter or function is missing
+       SyntaxError - improper expression syntax
+       ValueError - expressions have circular dependencies
+
+    This function is not terribly sophisticated, and it would be easy to
+    trick.  However it handles the common cases cleanly and generates
+    reasonable messages for the common errors.
+
+    This code has not been fully audited for security.  While we have
+    removed the builtins and the ability to import modules, there may
+    be other vectors for users to perform more than simple function
+    evaluations.  Unauthenticated users should not be running this code.
+
+    Parameter names are assumed to contain only _.a-zA-Z0-9#[]
+    
+    The list of parameters is probably something like::
+    
+        parset.setprefix()
+        pars = parset.flatten()
+    
+    Note that math uses acos while numpy uses arccos.  To avoid confusion
+    we allow both.
+    
+    Should try running the function to identify syntax errors before
+    running it in a fit.
+    
+    Use help(fn) to see the code generated for the returned function fn.
+    dis.dis(fn) will show the corresponding python vm instructions.
+    """
+
+    # Initialize dictionary with available functions
+    globals = {}
+    globals.update(math.__dict__)
+    globals.update(dict(arcsin=math.asin,arccos=math.acos,
+                        arctan=math.atan,arctan2=math.atan2))
+    globals.update(context)
+
+    # Sort the parameters in the order they need to be evaluated
+    deps = find_dependencies(pars)
+    if deps == []: return no_constraints
+    par_table,par_mapping = parameter_mapping(deps)
+    order = order_dependencies(deps)
+    
+    # Finish setting up the global and local namespace
+    globals.update(par_table)
+    locals = {}
+
+    # Define the function body
+    exprs = [p.path+"="+p.expression for p in order]
+    code = [substitute(s,par_mapping) for s in exprs]
+        
+    # Define the constraints function
+    functiondef = """
+def eval_expressions():
+    '''
+    %s
+    '''
+    %s
+"""%("\n    ".join(exprs),"\n    ".join(code))
+
+    #print "Function:",function
+    exec functiondef in globals,locals
+    retfn = locals['eval_expressions']
+
+    # Remove garbage added to globals by exec
+    globals.pop('__doc__',None)
+    globals.pop('__name__',None)
+    globals.pop('__file__',None)
+    globals.pop('__builtins__')
+    #print globals.keys()
+
+    return retfn
+
+def test():
+    import inspect, dis
+    import math
+    
+    symtab = {'a.b.x':1, 'a.c':2, 'a.b':3, 'b.x':4}
+    expr = 'a.b.x + sin(4*pi*a.c) + a.b.x/a.b'
+    
+    # Check symbol lookup
+    assert symbols(expr, symtab) == set([1,2,3])
+
+    # Check symbol rename
+    assert substitute(expr,{'a.b.x':'Q'}) == 'Q + sin(4*pi*a.c) + Q/a.b'
+
+
+    # Check dependency builder
+    # Fake parameter class
+    class Parameter:
+        def __init__(self, name, value=0, expression=''):
+            self.path = name
+            self.value = value
+            self.expression = expression
+        def iscomputed(self): return (self.expression != '')
+        def __repr__(self): return self.path
+    p1 = Parameter('G0.sigma',5)
+    p2 = Parameter('other',expression='2*pi*sin(G0.sigma/.1875) + M1.G1')
+    p3 = Parameter('M1.G1',6)
+    p4 = Parameter('constant',expression='2*pi*35')
+    # Simple chain
+    assert set(find_dependencies([p1,p2,p3])) == set([(p2,p1),(p2,p3)])
+    # Constant expression
+    assert set(find_dependencies([p1,p4])) == set([(p4,None)])
+    # No dependencies
+    assert set(find_dependencies([p1,p3])) == set([])
+
+    # Check function builder
+    fn = build_eval([p1,p2,p3])
+
+    # Inspect the resulting function
+    if False:
+        print inspect.getdoc(fn)
+        print dis.dis(fn)
+
+    # Evaluate the function and see if it updates the
+    # target value as expected
+    fn()
+    expected = 2*math.pi*math.sin(5/.1875) + 6
+    assert p2.value == expected,"Value was %s, not %s"%(p2.value,expected)
+    
+    # Check empty dependency set doesn't crash
+    fn = build_eval([p1,p3])
+    fn()
+
+    # Check that constants are evaluated properly
+    fn = build_eval([p4])
+    fn()
+    assert p4.value == 2*math.pi*35
+
+    # Check additional context example; this also tests multiple
+    # expressions
+    class Table:
+        Si = 2.09
+        values = {'Si': 2.07}
+    tbl = Table()
+    p5 = Parameter('lookup',expression="tbl.Si")
+    fn = build_eval([p1,p2,p3,p5],context=dict(tbl=tbl))
+    fn()
+    assert p5.value == 2.09,"Value for %s was %s"%(p5.expression,p5.value)
+    p5.expression = "tbl.values['Si']"
+    fn = build_eval([p1,p2,p3,p5],context=dict(tbl=tbl))
+    fn()
+    assert p5.value == 2.07,"Value for %s was %s"%(p5.expression,p5.value)
+    
+
+    # Verify that we capture invalid expressions
+    for expr in ['G4.cage', 'M0.cage', 'M1.G1 + *2', 
+                 'piddle',
+                 'import sys; print "p0wned"',
+                 '__import__("sys").argv']:
+        try:
+            p6 = Parameter('broken',expression=expr)
+            fn = build_eval([p6])
+            fn()
+        except Exception,msg: pass
+        else:  raise "Failed to raise error for %s"%expr
+
+if __name__ == "__main__": test()
Index: park-1.2.1/park/fit.py
===================================================================
--- park-1.2.1/park/fit.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/fit.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,332 @@
+# This program is public domain
+"""
+Fitting service interface.
+
+A fit consists of a set of models and a fitting engine.  The models are
+collected in an assembly, which manages the parameter set and the
+constraints between them.  The models themselves are tightly coupled
+to the data that they are modeling and the data is invisible to the fit.
+
+The fitting engine can use a variety of methods depending on model.
+
+
+Usage
+=====
+
+The fitter can be run directly on the local machine::
+
+    import park
+    M1 = park.models.Peaks(datafile=park.sampledata('peak.dat'))
+    M1.add_peak('P1', 'gaussian', A=[4,6], mu=[0.2, 0.5], sigma=0.1)
+    result = park.fit(models=[M1])
+    print result
+
+The default settings print results every time the fit improves, and
+print a global result when the fit is complete.  This is a suitable
+interface for a fitting script.
+
+For larger fit jobs you will want to run the fit on a remote server.
+The model setup is identical, but the fit call is different::
+
+    service = park.FitService('server:port')
+    result = park.fit(models=[M1], service=service)
+    print result
+
+Again, the default settings print results every time the fit improves,
+and print a global result when the fit is complete.
+
+For long running fit jobs, you want to be able to disconnect from
+the server after submitting the job, and later reconnect to fetch
+the results.  An additional email field will send notification by
+email when the fit starts and ends, and daily updates on the status
+of all fits::
+
+    service = park.FitService('server:port')
+    service.notify(email='me@my.email.address',update='daily')
+    fit = park.Fit(models=[M1])
+    id = service.submit_job(fit, jobname='peaks')
+    print id
+
+The results can be retrieved either by id returned from the server,
+or by the given jobname::
+
+    import park
+    service = park.FitService('server:port',user='userid')
+    fitlist = service.retrieve('peaks')
+    for fit in fitlist:
+        print fit.summary()
+
+The fit itself is a complicated object, including the model, the
+optimizer, and the type of uncertainty analysis to perform.
+
+GUI Usage
+=========
+
+When used from a graphical user interface, a different programming
+interface is needed.  In this case, the user may want to watch
+the progress of the fit and perhaps stop it.  Also, as fits can
+take some time to complete, the user would like to be able to
+set up additional fits and run them at the same time, switching
+between them as necessary to monitor progress.
+
+"""
+import time, thread
+
+import numpy
+
+import assembly, fitresult
+
+class Objective(object):
+    """
+    Abstract interface to the fitness function for the park minimizer
+    classes.
+
+    Park provides a specific implementation `park.assembly.Assembly`.
+
+    TODO: add a results() method to return model specific info to the
+    TODO: fit handler.
+    """
+    def residuals(self, p):
+        """
+        Some fitters, notably Levenberg-Marquardt, operate directly on the
+        residuals vector.  If the individual residuals are not available,
+        then LM cannot be used.
+
+        This method is optional.
+        """
+        raise NotImplementedError
+
+    def residuals_deriv(self, p):
+        """
+        Returns residuals and derivatives with respect to the given
+        parameters.
+
+        If these are unavailable in the model, then they can be approximated
+        by numerical derivatives, though it is generally better to use a
+        derivative free optimizer such as coliny or cobyla which can use the
+        function evaluations more efficiently.  In any case, your objective
+        function is responsible for calculating these.
+
+        This method is optional.
+        """
+        raise NotImplementedError
+
+    def fit_parameters(self):
+        """
+        Returns a list of fit parameters.  Each parameter has a name,
+        an initial value and a range.
+
+        See `park.fitresult.FitParameter` for an example.
+
+        On each function evaluation a new parameter set will be passed
+        to the fitter, with values in the same order as the list of
+        parameters.
+        """
+        raise NotImplementedError
+
+    def __call__(self, p):
+        """
+        Returns the objective value for parameter set p .
+        """
+        raise NotImplementedError
+
+    def abort(self):
+        """
+        Halts the current function evaluation, and has it return inf.
+        This will be called from a separate thread.  If the function
+        contains an expensive calculation, it should reset an abort
+        flag before each evaluation and test it periodically.
+
+        This method is optional.
+        """
+
+class Fitter(object):
+    """Abstract interface for a fitness optimizer.
+
+    A fitter has a single method, fit, which takes an objective
+    function (`park.fit.Objective`) and a handler.
+
+    For a concrete instance see `park.fitmc.FitMC`.
+    """
+    def __init__(self, **kw):
+        for k,v in kw.items():
+            if hasattr(self,k):
+                setattr(self,k,v)
+            else:
+                raise AttributeError(k+" is not an attribute of "+self.__class__.__name__)
+
+    def _threaded(self, fn, *args, **kw):
+        thread.start_new_thread(fn,args,kw)
+
+
+    def _fit(self, objective, x0, bounds):
+        """
+        Run the actual fit in a separate thread
+
+        Each cycle k of n:
+            self.handler.progress(k,n)
+        Each improvement:
+            self.handler.result.update(x,fx,ncalls)
+            self.handler.improvement()
+        On completion (if not already performed):
+            self.hander.result.update(x,fx,ncalls)
+            self.handler.done
+            self.handler.finalize()
+        """
+        raise NotImplementedError
+
+    def fit(self, fitness, handler):
+        """
+        Global optimizer.
+
+        This function should return immediately
+        """
+        # Determine initial value and bounds
+        pars = fitness.fit_parameters()
+        bounds = numpy.array([p.range for p in pars]).T
+        x0 = [p.value for p in pars]
+
+        # Initialize the monitor and results.
+        # Need to make our own copy of the fit results so that the
+        # values don't get stomped on by the next fit iteration.
+        handler.done = False
+        self.handler = handler
+        fitpars = [fitresult.FitParameter(pars[i].name,pars[i].range,v)
+                   for i,v in enumerate(x0)]
+        handler.result = fitresult.FitResult(fitpars, 0, numpy.NaN)
+
+        # Run the fit (fit should perform _progress and _improvement updates)
+        # This function may return before the fit is complete.
+        self._fit(fitness, x0, bounds)
+
+class FitJob(object):
+    """
+    Fit job.
+
+    This implements `park.job.Job`.
+    """
+    def __init__(self, objective=None, fitter=None, handler=None):
+        self.fitter = fitter
+        self.objective = objective
+        self.handler = handler
+    def run(self):
+        self.fitter.fit(self.objective, self.handler)
+
+class LocalQueue(object):
+    """
+    Simple interface to the local job queue.  Currently supports start and
+    wait.  Needs to support stop and status.  Also, needs to be a proper queue,
+    and needs to allow jobs to run in separate processes according to priority,
+    etc.  All the essentials of the remote queuing system without the remote
+    calls.
+
+    Unlike the remote queue, the local queue need not maintain persistence.
+    """
+    running = False
+    def start(self, job):
+        self.job = job
+        job.run()
+        return id(job)
+
+    def wait(self, interval=.1):
+        """
+        Wait for the job to complete.  This is used in scripts to impose
+        a synchronous interface to the fitting service.
+        """
+        while not self.job.handler.done:
+            time.sleep(interval)
+        return self.job.handler.result
+
+def fit(models=None, fitter=None, service=None, handler=None):
+    """
+    Start a fit with a set of models.  The model set must be
+    in a form accepted by `park.assembly.Assembly`.
+
+    This is a convenience function which sets up the default
+    optimizer and uses the local fitting engine to do the work.
+    Progress reports are printed as they are received.
+
+    The choice of fitter, service and handler can be specified
+    by the caller.
+
+    The default fitter is FitMC, which is a monte carlo Nelder-Mead
+    simplex local optimizer with 100 random start points.
+
+    The default handler does nothing.  Instead, ConsoleUpdate could
+    be used to report progress during the fit.
+
+    The default service is to run in a separate thread with FitThread.
+    Note that this will change soon to run in a separate process on
+    the local machine so that python's global interpreter lock does
+    not interfere with parallelism.
+    """
+    if models is None: raise RuntimeError('fit expected a list of models')
+    if service is None: service = LocalQueue()
+    if fitter is None:
+        import fitmc
+        fitter = fitmc.FitMC()
+    if handler is None: handler = fitresult.FitHandler()
+
+    objective = assembly.Assembly(models) if isinstance(models,list) else models
+    job = FitJob(objective,fitter,handler)
+    service.start(job)
+    return service.wait()
+
+
+def assembly_example():
+    import park, time
+    problem = park.assembly.example()
+    #result = fit(problem)
+    #result.print_summary()
+    handler=fitresult.ConsoleUpdate(improvement_delta=0.1,progress_delta=1)
+    #result = fit(problem, handler=handler)
+    result = fit(problem)
+    print "=== Fit complete ==="
+    result.print_summary()
+    print "=== Target values ==="
+    print "M1: a=1, c=1.5"
+    print "M2: a=2.5, c=3"
+
+    if False:  # Detailed results
+        print "parameter vector",result.pvec
+        problem(result.pvec)
+        print "residuals",problem.residuals
+        for k,m in enumerate(problem.parts):
+            print "Model",k,"chisq",m.chisq,"weight",m.weight
+            print "pars",m.fitness.model.a,m.fitness.model.c
+            print "x",m.fitness.data.fit_x
+            print "y",m.fitness.data.fit_y
+            print "f(x)",m.fitness.data.fx
+            print "(y-f(x))/dy",m.residuals
+
+
+def demo(fitter=None):
+    """Multiple minima example"""
+    import time, math
+    class MultiMin(object):
+        def fit_parameters(self):
+            return [fitresult.FitParameter('x1',[-5,5],1)]
+        def __call__(self, p):
+            x=p[0]
+            fx = x**2 + math.sin(2*math.pi*x+3*math.pi/2)
+            return fx
+    handler = fitresult.ConsoleUpdate() # Show updates on the console
+    handler.progress_delta = 1          # Update user every second
+    handler.improvement_delta = 0.1     # Show improvements almost immediately
+    fitter.fit(MultiMin(), handler)
+    while not handler.done: time.sleep(1)
+
+def demo2(fitter=None):
+    import park, time
+    problem = park.assembly.example()
+    handler = fitresult.ConsoleUpdate() # Show updates on the console
+    handler.progress_delta = 1          # Update user every second
+    handler.improvement_delta = 1       # Show improvements at the same rate
+    fitter.fit(problem, handler)
+    while not handler.done: time.sleep(1)
+
+
+
+if __name__ == "__main__":
+    #main()
+    assembly_example()
Index: park-1.2.1/park/fitmc.py
===================================================================
--- park-1.2.1/park/fitmc.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/fitmc.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,180 @@
+
+# The job queue needs to be in a transaction/rollback protected database.
+# If the server is rebooted, long running jobs should continue to work.
+#
+from __future__ import division
+import numpy
+
+import simplex
+
+import fitresult, pmap, fit
+
+__all__ = ['fitmc', 'FitMCJob']
+
+class LocalFit(object):
+    """
+    Abstract interface for a local minimizer
+
+    See `park.fitmc.FitSimplex` for a concrete implementation.
+    """
+    def fit(self, objective, x0):
+        """
+        Minimize the value of a fitness function.
+
+        See `park.fitmc.Fitness` for the definition of the goodness of fit
+        object.  x0 is a vector containing the initial value for the fit.
+        """
+    def abort(self):
+        """
+        Cancel the fit.  This will be called from a separate execution
+        thread.  The fit should terminate as soon as possible after this
+        function has been called.  Ideally this would interrupt the
+        cur
+        """
+
+class FitSimplex(LocalFit):
+    """
+    Local minimizer using Nelder-Mead simplex algorithm.
+
+    Simplex is robust and derivative free, though not very efficient.
+
+    This class wraps the bounds contrained Nelder-Mead simplex
+    implementation for `park.simplex.simplex`.
+    """
+    radius = 0.05
+    """Size of the initial simplex; this is a portion between 0 and 1"""
+    xtol = 1
+    #xtol = 1e-4
+    """Stop when simplex vertices are within xtol of each other"""
+    ftol = 1e-4
+    """Stop when vertex values are within ftol of each other"""
+    maxiter = None
+    """Maximum number of iterations before fit terminates"""
+    def fit(self, fitness, x0):
+        """Run the fit"""
+        self.cancel = False
+        pars = fitness.fit_parameters()
+        bounds = numpy.array([p.range for p in pars]).T
+        result = simplex.simplex(fitness, x0, bounds=bounds,
+                                 radius=self.radius, xtol=self.xtol,
+                                 ftol=self.ftol, maxiter=self.maxiter,
+                                 abort_test=self._iscancelled)
+        #print "calls:",result.calls
+        #print "simplex returned",result.x,result.fx
+        # Need to make our own copy of the fit results so that the
+        # values don't get stomped on by the next fit iteration.
+        fitpars = [fitresult.FitParameter(pars[i].name,pars[i].range,v)
+                   for i,v in enumerate(result.x)]
+        res = fitresult.FitResult(fitpars, result.calls, result.fx)
+        # Compute the parameter uncertainties from the jacobian
+        res.calc_cov(fitness)
+        return res
+
+    def abort(self):
+        """Cancel the fit in progress from another thread of execution"""
+        # Effectively an atomic operation; rely on GIL to protect it.
+        self.cancel = True
+        # Abort the current function evaluation if possible.
+        if hasattr(fitness,'abort'): self.fitness.abort()
+
+    def _iscancelled(self): return self.cancel
+
+class MapMC(object):
+    """
+    Evaluate a local fit at a particular start point.
+
+    This is the function to be mapped across a set of start points for the
+    monte carlo map-reduce implementation.
+
+    See `park.pmap.Mapper` for required interface.
+    """
+    def __init__(self, minimizer, fitness):
+        self.minimizer, self.fitness = minimizer, fitness
+    def __call__(self, x0):
+        return self.minimizer.fit(self.fitness,x0)
+    def abort(self):
+        self.minimizer.abort()
+
+class CollectMC(object):
+    """
+    Collect the results from multiple start points in a Monte Carlo fit engine.
+
+    See `park.pmap.Collector` for required interface.
+    """
+    def __init__(self, number_expected, handler):
+        self.number_expected = number_expected
+        """Number of starting points to check with local optimizer"""
+        self.iterations = 0
+        self.best = None
+        self.calls = 0
+        self.handler = handler
+        self.handler.done = False
+    def __call__(self, result):
+        # Keep track of the amount of work done
+        self.iterations += 1
+        self.calls += result.calls
+        if self.best is None or result.fitness < self.best.fitness:
+            self.best = result
+            self.handler.result = result
+            self.handler.improvement()
+        # Progress should go after improvement in case the fit handler
+        # wants to suppress intermediate improvements
+        self.handler.progress(self.iterations, self.number_expected)
+    def abort(self):
+        self.handler.done = True
+        self.handler.abort()
+    def finalize(self):
+        self.handler.result.calls = self.calls
+        self.handler.done = True
+        self.handler.finalize()
+    def error(self, msg):
+        self.handler.done = True
+        self.handler.error(msg)
+
+def fitmc(fitness, x0, bounds, localfit, n, handler):
+    """
+    Run a monte carlo fit.
+
+    This procedure maps a local optimizer across a set of n initial points.
+    The initial parameter value defined by the fitness parameters defines
+    one initial point.  The remainder are randomly generated within the
+    bounds of the problem.
+
+    localfit is the local optimizer to use.  It should be a bounded
+    optimizer following the `park.fitmc.LocalFit` interface.
+
+    handler accepts updates to the current best set of fit parameters.
+    See `park.fitresult.FitHandler` for details.
+    """
+    lo,hi = bounds
+    delta = hi-lo
+    delta[numpy.isinf(delta)] = 1
+    lo[numpy.isinf(lo)] = hi-1
+    lo[numpy.isinf(lo)] = 0
+    P = numpy.random.rand(n,len(x0))*delta+lo
+    #print "Population",P
+    P[0] = x0
+    pmap.pmapreduce(MapMC(localfit,fitness),
+                    CollectMC(n,handler),
+                    P)
+
+class FitMC(fit.Fitter):
+    """
+    Monte Carlo optimizer.
+
+    This implements `park.fit.Fitter`.
+    """
+    localfit = FitSimplex()
+    start_points = 10
+
+    def _fit(self, objective, x0, bounds):
+        """
+        Run a monte carlo fit.
+
+        This procedure maps a local optimizer across a set of initial points.
+        """
+        fitmc(objective, x0, bounds, self.localfit,
+              self.start_points, self.handler)
+
+if __name__ == "__main__":
+    fit.demo2(FitMC(localfit=FitSimplex(),start_points=10))
Index: park-1.2.1/park/fitresult.py
===================================================================
--- park-1.2.1/park/fitresult.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/fitresult.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,237 @@
+import sys, math, time
+import numpy
+
+from formatnum import format_uncertainty
+
+class FitHandler(object):
+    """
+    Abstract interface for fit thread handler.
+
+    The methods in this class are called by the optimizer as the fit
+    progresses.
+
+    Note that it is up to the optimizer to call the fit handler correctly,
+    reporting all status changes and maintaining the 'done' flag.
+    """
+    done = False
+    """True when the fit job is complete"""
+    result = None
+    """The current best result of the fit"""
+
+    def improvement(self):
+        """
+        Called when a result is observed which is better than previous
+        results from the fit.
+
+        result is a FitResult object, with parameters, #calls and fitness.
+        """
+    def error(self, msg):
+        """
+        Model had an error; print traceback
+        """
+    def progress(self, current, expected):
+        """
+        Called each cycle of the fit, reporting the current and the
+        expected amount of work.   The meaning of these values is
+        optimizer dependent, but they can be converted into a percent
+        complete using (100*current)//expected.
+
+        Progress is updated each iteration of the fit, whatever that
+        means for the particular optimization algorithm.  It is called
+        after any calls to improvement for the iteration so that the
+        update handler can control I/O bandwidth by suppressing
+        intermediate improvements until the fit is complete.
+        """
+    def finalize(self):
+        """
+        Fit is complete; best results are reported
+        """
+    def abort(self):
+        """
+        Fit was aborted.
+        """
+
+class ConsoleUpdate(FitHandler):
+    """
+    Print progress to the console.
+    """
+    isbetter = False
+    """Record whether results improved since last update"""
+    progress_delta =  60
+    """Number of seconds between progress updates"""
+    improvement_delta = 5
+    """Number of seconds between improvement updates"""
+    def __init__(self,quiet=False,progress_delta=60,improvement_delta=5):
+        """
+        If quiet is true, only print out final summary, not progress and
+        improvements.
+        """
+        self.progress_time = time.time()
+        self.progress_percent = 0
+        self.improvement_time = self.progress_time
+        self.isbetter = False
+        self.quiet = quiet
+        self.progress_delta = progress_delta
+        self.improvement_delta = improvement_delta
+
+    def progress(self, k, n):
+        """
+        Report on progress.
+        """
+        if self.quiet: return
+        t = time.time()
+        p = int((100*k)//n)
+
+        # Show improvements if there are any
+        dt = t - self.improvement_time
+        if self.isbetter and dt > self.improvement_delta:
+            self.result.print_summary()
+            self.isbetter = False
+            self.improvement_time = t
+
+        # Update percent complete
+        dp = p-self.progress_percent
+        if dp < 1: return
+        dt = t - self.progress_time
+        if dt > self.progress_delta:
+            if 1 <= dp <= 2:
+                print "%d%% complete"%p
+                self.progress_percent = p
+                self.progress_time = t
+            elif 2 < dp <= 5:
+                if p//5 != self.progress_percent//5:
+                    print "%d%% complete"%(5*(p//5))
+                    self.progress_percent = p
+                    self.progress_time = t
+            else:
+                if p//10 != self.progress_percent//10:
+                    print "%d%% complete"%(10*(p//10))
+                    self.progress_percent = p
+                    self.progress_time = t
+
+    def improvement(self):
+        """
+        Called when a result is observed which is better than previous
+        results from the fit.
+        """
+        self.isbetter = True
+
+    def error(self, msg):
+        """
+        Model had an error; print traceback
+        """
+        if self.isbetter:
+            self.result.print_summary()
+        print msg
+
+    def finalize(self):
+        if self.isbetter:
+            self.result.print_summary()
+        print "Total function calls:",self.result.calls
+
+    def abort(self):
+        if self.isbetter:
+            self.result.print_summary()
+
+
+class FitParameter(object):
+    """
+    Fit result for an individual parameter.
+    """
+    def __init__(self, name, range, value):
+        self.name = name
+        self.range = range
+        self.value = value
+        self.stderr = None
+    def summarize(self):
+        """
+        Return parameter range string.
+
+        E.g.,  "       Gold .....|.... 5.2043 in [2,7]"
+        """
+        bar = ['.']*10
+        lo,hi = self.range
+        if numpy.isfinite(lo)and numpy.isfinite(hi):
+            portion = (self.value-lo)/(hi-lo)
+            if portion < 0: portion = 0.
+            elif portion >= 1: portion = 0.99999999
+            barpos = int(math.floor(portion*len(bar)))
+            bar[barpos] = '|'
+        bar = "".join(bar)
+        lostr = "[%g"%lo if numpy.isfinite(lo) else "(-inf"
+        histr = "%g]"%hi if numpy.isfinite(hi) else "inf)"
+        valstr = format_uncertainty(self.value, self.stderr)
+        return "%25s %s %s in %s,%s"  % (self.name,bar,valstr,lostr,histr)
+    def __repr__(self):
+        return "FitParameter('%s')"%self.name
+
+class FitResult(object):
+    """
+    Container for reporting fit results.
+    """
+    def __init__(self, parameters, calls, fitness):
+        self.parameters = parameters
+        """Fit parameter list, each with name, range and value attributes."""
+        self.calls = calls
+        """Number of function calls"""
+        self.fitness = fitness
+        """Value of the goodness of fit metric"""
+        self.pvec = numpy.array([p.value for p in self.parameters])
+        """Parameter vector"""
+        self.stderr = None
+        """Parameter uncertainties"""
+        self.cov = None
+        """Covariance matrix"""
+
+    def update(self, pvec, fitness, calls):
+        self.calls = calls
+        self.fitness = fitness
+        self.pvec = pvec.copy()
+        for i,p in enumerate(self.parameters):
+            p.value = pvec[i]
+
+    def calc_cov(self, fn):
+        """
+        Return the covariance matrix inv(J'J) at point p.
+        """
+        if hasattr(fn, 'jacobian'):
+            # Find cov of f at p
+            #     cov(f,p) = inv(J'J)
+            # Use SVD
+            #     J = U S V'
+            #     J'J = (U S V')' (U S V')
+            #         = V S' U' U S V'
+            #         = V S S V'
+            #     inv(J'J) = inv(V S S V')
+            #              = inv(V') inv(S S) inv(V)
+            #              = V inv (S S) V'
+            J = fn.jacobian(self.pvec)
+            u,s,vh = numpy.linalg.svd(J,0)
+            JTJinv = numpy.dot(vh.T.conj()/s**2,vh)
+            self.set_cov(JTJinv)
+
+    def set_cov(self, cov):
+        """
+        Return the covariance matrix inv(J'J) at point p.
+        """
+        self.cov = cov
+        if cov is not None:
+            self.stderr = numpy.sqrt(numpy.diag(self.cov))
+            # Set the uncertainties on the individual parameters
+            for k,p in enumerate(self.parameters):
+                p.stderr = self.stderr[k]
+        else:
+            self.stderr = None
+            # Reset the uncertainties on the individual parameters
+            for k,p in enumerate(self.parameters):
+                p.stderr = None
+
+    def __str__(self):
+        if self.parameters == None: return "No results"
+        L = ["P%-3d %s"%(n+1,p.summarize()) for n,p in enumerate(self.parameters)]
+        L.append("=== goodness of fit: %g"%(self.fitness))
+        return "\n".join(L)
+
+    def print_summary(self, fid=sys.stdout):
+        print >>fid, self
+
Index: park-1.2.1/park/fitservice.py
===================================================================
--- park-1.2.1/park/fitservice.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/fitservice.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,126 @@
+# This program is public domain
+"""
+Interface to the PARK fitting service.
+
+*** WARNING *** This is a design sketch and is not presently used.
+"""
+
+class FitClient(object):
+    """
+    Client-side view of the fitting service.
+    """
+    def __init__(self, url=DEFAULT_URL,user=DEFAULT_USER):
+        self.url = url
+        self.user = user
+        self.server = xmlrpclib.Server(self.server)
+        self.context = self.server.welcome(self.user)
+
+    # User management
+    def notify(self,email=None,rate=None):
+        """Set the email address and update frequency for user notification"""
+        if email is not None:
+            self.server.setemail(self.context, email)
+        if rate is not None:
+            self.server.setrate(self.context, rate)
+
+    # Job management
+    def find_by_name(self, name):
+        return self.server.find_by_name(self.context,status)
+    def find_jobs(self,status='active'):
+        """List active and recently completed jobs"""
+        return self.server.find_jobs(self.context,status)
+    def stop(self, jobid):
+        """Delete job"""
+        self.server.delete(self.context, jobid)
+    def start(self, job):
+        s = encode(job)
+        return self.server.submit(self.context, jobname, job)
+    def fetch_job(self, jobid):
+        """Retrieve job"""
+        text = self.server.retrieve(self.context, jobid)
+        return encode(text)
+    def status(self, jobid):
+        """Return the current job status"""
+        return self.server.status(self.context, jobid)
+    
+    # Message Queue
+    def _listener(self, callback):
+        while True:
+            message = self.server.next_message(self.user)
+            event,result = pickle.loads(text)
+            callback(event, self, result)
+
+    def listen(self, callback=on_message):
+        """
+        Listen to the message queue for information about running jobs.
+
+        The listener runs in a separate thread, allowing the listen
+        call to return immediately.  The callback is called each time
+        there is a new message on the queue.  If no callback is
+        supplied, then a simple print is used.
+        """
+        thread.start_new_thread(self._listener,(callback,))
+
+
+
+class Fit(object): 
+    """ 
+    Coordinate the fit.
+    """
+    update = on_fit_update
+    complete = on_fit_complete
+    def __init__(self, models = [],
+                 optimizer = None,
+                 service = None,
+                 update = None,
+                 complete = None):
+        """
+        Set up the fit.
+
+        models []
+           List of models which make up the assembly.
+        optimizer
+           Optimizer to use for the fit
+        """
+        self.assembly = assembly.Assembly(models)
+        self.optimizer = optimizer
+        self.complete = complete
+        self.update = update
+        self.isrunning = False
+        self._stopping = False
+
+    def stop(self):
+        """
+        Interrupt a running fit.
+        """
+        self._stopping = True
+        self.service.stop()
+
+    def start(self):
+        """
+        Start the fit.
+        """
+        self.service.start()
+
+    def _update(self, result):
+        self.result = result
+        if self.update: self.update(self,result)
+
+    def _complete(self, result):
+        self.result = result
+        if self.complete: self.complete(self,result)
+
+    def wait(self):
+        """
+        Wait for a fit to complete.
+        """
+        while self.isrunning: time.sleep(0.1)
+
+    def fit(self):
+        """ 
+        Start the fit and wait for the result.
+        Returns the result.
+        """
+        self.start()
+        self.wait()
+        return self.result
Index: park-1.2.1/park/formatnum.py
===================================================================
--- park-1.2.1/park/formatnum.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/formatnum.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,312 @@
+# This program is public domain
+# Author Paul Kienzle
+"""
+Format numbers nicely for printing.
+
+Usage::
+
+   >> from danse.common.util.formatnum import *
+   >> v,dv = 757.2356,0.01032
+   >> print format_uncertainty_pm(v,dv)
+   757.235 +/- 0.010
+   >> format_uncertainty_compact(v,dv)
+   757.235(10)
+   >> format_uncertainty(v,dv)
+   757.235(10)
+
+Set format_uncertainty.compact to False to use the +/-
+format by default, otherwise leave it at True for compact
+value(##) format.
+"""
+from __future__ import division
+
+import math
+import numpy
+__all__ = ['format_uncertainty', 'format_uncertainty_pm', 
+           'format_uncertainty_compact']
+
+# These routines need work for +/- formatting::
+# - pm formats are not rigorously tested
+# - pm formats do not try to align the scale to multiples of 1000
+# - pm formats do not try to align the value scale to uncertainty scale
+
+# Coordinating scales across a set of numbers is not supported.  For easy
+# comparison a set of numbers should be shown in the same scale.  One could
+# force this from the outside by adding scale parameter (either 10**n, n, or
+# a string representing the desired SI prefix) and having a separate routine
+# which computes the scale given a set of values.
+
+# Coordinating scales with units offers its own problems.  Again, the user
+# may want to force particular units.  This can be done by outside of the
+# formatting routines by scaling the numbers to the appropriate units then
+# forcing them to print with scale 10**0.  If this is a common operation,
+# however, it may want to happen inside.
+
+# The value e<n> is currently formatted into the number.  Alternatively this
+# scale factor could be returned so that the user can choose the appriate
+# SI prefix when printing the units.  This gets tricky when talking about
+# composite units such as 2.3e-3 m**2 -> 2300 mm**2, and with volumes
+# such as 1 g/cm**3 -> 1 kg/L.
+
+def format_uncertainty_pm(value, uncertainty=None):
+    """
+    Given *value* v and *uncertainty* dv, return a string v +/- dv.
+    
+    The returned string uses only the number of digits warranted by
+    the uncertainty in the measurement.
+
+    If the uncertainty is 0 or not otherwise provided, the simple
+    %g floating point format option is used.
+
+    Infinite and indefinite numbers are represented as inf and NaN.
+    """
+    return _format_uncertainty(value, uncertainty, compact=False)
+
+def format_uncertainty_compact(value, uncertainty=None):
+    """
+    Given *value* v and *uncertainty* dv, return the compact
+    representation v(##), where ## are the first two digits of 
+    the uncertainty.
+    
+    The returned string uses only the number of digits warranted by
+    the uncertainty in the measurement.
+
+    If the uncertainty is 0 or not otherwise provided, the simple
+    %g floating point format option is used.
+
+    Infinite and indefinite numbers are represented as inf and NaN.
+    """
+    return _format_uncertainty(value, uncertainty, compact=True)
+
+class FormatUncertainty:
+    """
+    Given *value* and *uncertainty*, return a concise string representation.
+
+    This will either by the +/- form of :func:`format_uncertainty_pm` or
+    the compact form of :func:`format_uncertainty_compact` depending on 
+    whether *compact* is specified or whether *format_uncertainty.compact* 
+    is True or False.
+    """
+    compact = True
+    def __call__(self, value, uncertainty):
+        return _format_uncertainty(value, uncertainty, self.compact)
+format_uncertainty = FormatUncertainty()
+
+def _format_uncertainty(value,uncertainty,compact):
+    """
+    Implementation of both the compact and the +/- formats.
+    """
+    if numpy.isinf(value):
+        return "inf" if value > 0 else "-inf"
+    
+    if numpy.isnan(value):
+        return "NaN"
+
+    # Uncertainty check must come after indefinite check since the %g
+    # format string doesn't handle indefinite numbers consistently
+    # across platforms.
+    if uncertainty == None or uncertainty <= 0 or numpy.isnan(uncertainty):
+        return "%g"%value
+    if numpy.isinf(uncertainty):
+        if compact:
+            return "%g(inf)"%value
+        else:
+            return "%g +/- inf"%value
+    
+    # Process sign
+    sign = "-" if value < 0 else ""
+    value = abs(value)
+    
+    # Determine the number of digits in the value and the error
+    # Note that uncertainty <= 0 is handled above
+    err_place = int(math.floor(math.log10(uncertainty)))
+    if value == 0:
+        val_place = err_place-1
+    else:
+        val_place = int(math.floor(math.log10(value)))
+
+    # If pm, return a simple v +/- dv
+    if not compact:
+        scale = 10**(err_place-1)
+        val_digits = val_place-err_place+2
+        return "%s%.*g +/- %.2g"%(sign,val_digits,value,uncertainty)
+
+    # The remainder is for the v(dv) case 
+    err_str = "(%2d)"%int(uncertainty/10.**(err_place-1)+0.5)
+    if err_place > val_place:
+        # Degenerate case: error bigger than value
+        # The mantissa is 0.#(##)e#, 0.0#(##)e# or 0.00#(##)e#
+        if err_place - val_place > 2: value = 0
+        val_place = int(math.floor((err_place+2)/3.))*3
+        digits_after_decimal = val_place - err_place + 1
+        val_str = "%.*f%s"%(digits_after_decimal,value/10.**val_place,err_str)
+        if val_place != 0: val_str += "e%d"%val_place
+    elif err_place == val_place:
+        # Degenerate case: error and value the same order of magnitude
+        # The value is ##(##)e#, #.#(##)e# or 0.##(##)e#
+        val_place = int(math.floor((err_place+1)/3.))*3
+        digits_after_decimal = val_place - err_place + 1
+        val_str = "%.*f%s"%(digits_after_decimal,value/10.**val_place,err_str)
+        if val_place != 0: val_str += "e%d"%val_place
+    elif err_place <= 1 and val_place >= -3:
+        # Normal case: nice numbers and errors
+        # The value is ###.###(##)
+        digits_after_decimal = abs(err_place-1)
+        val_str = "%.*f%s"%(digits_after_decimal,value,err_str)
+    else:
+        # Extreme cases: zeros before value or after error
+        # The value is ###.###(##)e#, ##.####(##)e# or #.#####(##)e#
+        total_digits = val_place - err_place + 2
+        val_place = int(math.floor(val_place/3.))*3
+        val_str = "%.*g%se%d"%(total_digits,
+                               value/10.**val_place,
+                               err_str,val_place)
+
+    return sign+val_str
+
+
+
+def test():
+    # Oops... renamed function after writing tests
+    value_str = format_uncertainty_compact
+    
+    # val_place > err_place
+    assert value_str(1235670,766000) == "1.24(77)e6"
+    assert value_str(123567.,76600) == "124(77)e3"
+    assert value_str(12356.7,7660) == "12.4(77)e3"
+    assert value_str(1235.67,766) == "1.24(77)e3"
+    assert value_str(123.567,76.6) == "124(77)"
+    assert value_str(12.3567,7.66) == "12.4(77)"
+    assert value_str(1.23567,.766) == "1.24(77)"
+    assert value_str(.123567,.0766) == "0.124(77)"
+    assert value_str(.0123567,.00766) == "0.0124(77)"
+    assert value_str(.00123567,.000766) == "0.00124(77)"
+    assert value_str(.000123567,.0000766) == "124(77)e-6"
+    assert value_str(.0000123567,.00000766) == "12.4(77)e-6"
+    assert value_str(.00000123567,.000000766) == "1.24(77)e-6"
+    assert value_str(.000000123567,.0000000766) == "124(77)e-9"
+    assert value_str(.00000123567,.0000000766) == "1.236(77)e-6"
+    assert value_str(.0000123567,.0000000766) == "12.357(77)e-6"
+    assert value_str(.000123567,.0000000766) == "123.567(77)e-6"
+    assert value_str(.00123567,.000000766) == "0.00123567(77)"
+    assert value_str(.0123567,.00000766) == "0.0123567(77)"
+    assert value_str(.123567,.0000766) == "0.123567(77)"
+    assert value_str(1.23567,.000766) == "1.23567(77)"
+    assert value_str(12.3567,.00766) == "12.3567(77)"
+    assert value_str(123.567,.0764) == "123.567(76)"
+    assert value_str(1235.67,.764) == "1235.67(76)"
+    assert value_str(12356.7,7.64) == "12356.7(76)"
+    assert value_str(123567,76.4) == "123567(76)"
+    assert value_str(1235670,764) == "1.23567(76)e6"
+    assert value_str(12356700,764) == "12.3567(76)e6"
+    assert value_str(123567000,7640) == "123.567(76)e6"
+    assert value_str(1235670000,76400) == "1.23567(76)e9"
+    
+    # val_place == err_place
+    assert value_str(123567,764000) == "0.12(76)e6"
+    assert value_str(12356.7,76400) == "12(76)e3"
+    assert value_str(1235.67,7640) == "1.2(76)e3"
+    assert value_str(123.567,764) == "0.12(76)e3"
+    assert value_str(12.3567,76.4) == "12(76)"
+    assert value_str(1.23567,7.64) == "1.2(76)"
+    assert value_str(.123567,.764) == "0.12(76)"
+    assert value_str(.0123567,.0764) == "12(76)e-3"
+    assert value_str(.00123567,.00764) == "1.2(76)e-3"
+    assert value_str(.000123567,.000764) == "0.12(76)e-3"
+
+    # val_place == err_place-1
+    assert value_str(123567,7640000) == "0.1(76)e6"
+    assert value_str(12356.7,764000) == "0.01(76)e6"
+    assert value_str(1235.67,76400) == "0.001(76)e6"
+    assert value_str(123.567,7640) == "0.1(76)e3"
+    assert value_str(12.3567,764) == "0.01(76)e3"
+    assert value_str(1.23567,76.4) == "0.001(76)e3"
+    assert value_str(.123567,7.64) == "0.1(76)"
+    assert value_str(.0123567,.764) == "0.01(76)"
+    assert value_str(.00123567,.0764) == "0.001(76)"
+    assert value_str(.000123567,.00764) == "0.1(76)e-3"
+
+    # val_place == err_place-2
+    assert value_str(12356700,7640000000) == "0.0(76)e9"
+    assert value_str(1235670,764000000) == "0.00(76)e9"
+    assert value_str(123567,76400000) == "0.000(76)e9"
+    assert value_str(12356,7640000) == "0.0(76)e6"
+    assert value_str(1235,764000) == "0.00(76)e6"
+    assert value_str(123,76400) == "0.000(76)e6"
+    assert value_str(12,7640) == "0.0(76)e3"
+    assert value_str(1,764) == "0.00(76)e3"
+    assert value_str(0.1,76.4) == "0.000(76)e3"
+    assert value_str(0.01,7.64) == "0.0(76)"
+    assert value_str(0.001,0.764) == "0.00(76)"
+    assert value_str(0.0001,0.0764) == "0.000(76)"
+    assert value_str(0.00001,0.00764) == "0.0(76)e-3"
+
+    # val_place == err_place-3
+    assert value_str(12356700,76400000000) == "0.000(76)e12"
+    assert value_str(1235670,7640000000) == "0.0(76)e9"
+    assert value_str(123567,764000000) == "0.00(76)e9"
+    assert value_str(12356,76400000) == "0.000(76)e9"
+    assert value_str(1235,7640000) == "0.0(76)e6"
+    assert value_str(123,764000) == "0.00(76)e6"
+    assert value_str(12,76400) == "0.000(76)e6"
+    assert value_str(1,7640) == "0.0(76)e3"
+    assert value_str(0.1,764) == "0.00(76)e3"
+    assert value_str(0.01,76.4) == "0.000(76)e3"
+    assert value_str(0.001,7.64) == "0.0(76)"
+    assert value_str(0.0001,0.764) == "0.00(76)"
+    assert value_str(0.00001,0.0764) == "0.000(76)"
+    assert value_str(0.000001,0.00764) == "0.0(76)e-3"
+
+    # negative values
+    assert value_str(-1235670,765000) == "-1.24(77)e6"
+    assert value_str(-1.23567,.765) == "-1.24(77)"
+    assert value_str(-.00000123567,.0000000765) == "-1.236(77)e-6"
+    assert value_str(-12356.7,7.64) == "-12356.7(76)"
+    assert value_str(-123.567,764) == "-0.12(76)e3"
+    assert value_str(-1235.67,76400) == "-0.001(76)e6"
+    assert value_str(-.000123567,.00764) == "-0.1(76)e-3"
+    assert value_str(-12356,7640000) == "-0.0(76)e6"
+    assert value_str(-12,76400) == "-0.000(76)e6"
+    assert value_str(-0.0001,0.764) == "-0.00(76)"
+
+    # zero values
+    assert value_str(0,0) == "0"
+    assert value_str(-1.23567,0) == "-1.23567"
+    assert value_str(0,765000) == "0.00(77)e6"
+    assert value_str(0,.765) == "0.00(77)"
+    assert value_str(0,.000000765) == "0.00(77)e-6"
+    assert value_str(0,76400) == "0.000(76)e6"
+    assert value_str(0,7640) == "0.0(76)e3"
+
+    # marked values
+    assert value_str(-numpy.inf,None) == "-inf"
+    assert value_str(numpy.inf,None) == "inf"
+    assert value_str(numpy.NaN,None) == "NaN"
+    
+    # plus/minus form
+    assert format_uncertainty_pm(-1.23567,0.765) == "-1.24 +/- 0.77"
+    assert format_uncertainty_compact(-1.23567,0.765) == "-1.24(77)"
+    assert format_uncertainty_pm(752.3567,0.01) == "752.357 +/- 0.01"
+    assert format_uncertainty(-1.23567,0.765) == "-1.24(77)"
+
+    # bad uncertainty
+    assert format_uncertainty_pm(-1.23567,numpy.NaN) == "-1.23567"
+    assert format_uncertainty_pm(-1.23567,-numpy.inf) == "-1.23567"
+    assert format_uncertainty_pm(-1.23567,-0.1) == "-1.23567"
+    assert format_uncertainty_compact(-1.23567,numpy.NaN) == "-1.23567"
+    assert format_uncertainty_compact(-1.23567,-numpy.inf) == "-1.23567"
+    assert format_uncertainty_compact(-1.23567,-0.1) == "-1.23567"
+
+    # no uncertainty
+    assert format_uncertainty_pm(-1.23567,0) == "-1.23567"
+    assert format_uncertainty_compact(-1.23567,0) == "-1.23567"
+    assert format_uncertainty_pm(-1.23567,None) == "-1.23567"
+    assert format_uncertainty_compact(-1.23567,None) == "-1.23567"
+
+    # inf uncertainty
+    assert format_uncertainty_pm(-1.23567,numpy.inf) == "-1.23567 +/- inf"
+    assert format_uncertainty_compact(-1.23567,numpy.inf) == "-1.23567(inf)"
+
+
+if __name__ == "__main__": test()
+    
Index: park-1.2.1/park/lib/modeling.cc
===================================================================
--- park-1.2.1/park/lib/modeling.cc	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/lib/modeling.cc	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,76 @@
+/* This program is public domain. */
+
+#include <iostream>
+#include <iomanip>
+#include <Python.h>
+
+extern "C" void
+resolution(int Nin, const double Qin[], const double Rin[],
+           int N, const double Q[], const double dQ[], double R[]);
+
+#define INVECTOR(obj,buf,len)										\
+    do { \
+        int err = PyObject_AsReadBuffer(obj, (const void **)(&buf), &len); \
+        if (err < 0) return NULL; \
+        len /= sizeof(*buf); \
+    } while (0)
+    
+#define OUTVECTOR(obj,buf,len) \
+    do { \
+        int err = PyObject_AsWriteBuffer(obj, (void **)(&buf), &len); \
+        if (err < 0) return NULL; \
+        len /= sizeof(*buf); \
+    } while (0)
+
+PyObject* Pconvolve(PyObject *obj, PyObject *args)
+{
+  PyObject *Qi_obj,*Ri_obj,*Q_obj,*dQ_obj,*R_obj;
+  const double *Qi, *Ri, *Q, *dQ;
+  double *R;
+  Py_ssize_t nQi, nRi, nQ, ndQ, nR;
+  
+  if (!PyArg_ParseTuple(args, "OOOOO:resolution", 
+                        &Qi_obj,&Ri_obj,&Q_obj,&dQ_obj,&R_obj)) return NULL;
+  INVECTOR(Qi_obj,Qi,nQi);
+  INVECTOR(Ri_obj,Ri,nRi);
+  INVECTOR(Q_obj,Q,nQ);
+  INVECTOR(dQ_obj,dQ,ndQ);
+  OUTVECTOR(R_obj,R,nR);
+  if (nQi != nRi) {
+#ifndef BROKEN_EXCEPTIONS
+    PyErr_SetString(PyExc_ValueError, "_librefl.convolve: Qi and Ri have different lengths");
+#endif
+    return NULL;
+  }
+  if (nQ != ndQ || nQ != nR) {
+#ifndef BROKEN_EXCEPTIONS
+    PyErr_SetString(PyExc_ValueError, "_librefl.convolve: Q, dQ and R have different lengths");
+#endif
+    return NULL;
+  }
+  resolution(nQi,Qi,Ri,nQ,Q,dQ,R);
+  return Py_BuildValue("");
+}
+
+static PyMethodDef methods[] = {
+   {"convolve",
+    Pconvolve,
+    METH_VARARGS,
+    "convolve(Qi,Ri,Q,dQ,R): compute convolution of width dQ[k] at points Q[k], returned in R[k]"},
+    {0}
+} ;
+
+
+#if defined(WIN32) && !defined(__MINGW32__)
+__declspec(dllexport)
+#endif
+
+	
+extern "C" void init_modeling(void) {
+  Py_InitModule4("_modeling",
+  		methods,
+		"Modeling C Library",
+		0,
+		PYTHON_API_VERSION
+		);
+}
Index: park-1.2.1/park/lib/resolution.c
===================================================================
--- park-1.2.1/park/lib/resolution.c	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/lib/resolution.c	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,252 @@
+/* This program is public domain. */
+
+#include <stdio.h>
+#include <assert.h>
+#include <string.h>
+#include <math.h>
+#ifdef SGI
+#include <ieeefp.h>
+#endif
+
+/* What to do at the endpoints --- USE_TRUNCATED_NORMALIZATION will
+ * avoid the assumption that the data is zero where it hasn't been
+ * measured.
+ */
+#define USE_TRUNCATED_NORMALIZATION
+
+/* Computed using extended precision with Octave's symbolic toolbox. */
+#define PI4          12.56637061435917295385
+#define PI_180        0.01745329251994329576
+#define LN256         5.54517744447956247533
+#define SQRT2         1.41421356237309504880
+#define SQRT2PI       2.50662827463100050241
+
+/* Choose the resolution limit based on relative contribution to resolution
+ * rather than absolute contribution.  So if we ignore everything below,
+ * e.g. 0.1% of the peak, that limit occurs when G(x)/G(0) = 0.001 for
+ * gaussian G of width sima, or when x = sqrt(-2 sigma^2 log(0.001)). */
+#define LOG_RESLIMIT -6.90775527898213703123
+
+/** \file 
+The resolution function returns the convolution of the reflectometry
+curve with a Q-dependent gaussian.
+
+We provide the following function:
+   resolution(Nin, Qin, Rin, N, Q, dQ, R)  returns convolution
+
+where
+   Nin is the number of theory points
+   Qin,Rin are the computed theory points
+   N is the number of Q points to calculate
+   Q are the locations of the measured data points
+   dQ are the width (sigma) of the convolution at each measured point
+   R is the returned convolution.
+
+Note that FWHM = sqrt(8 ln 2) dQ, so scale dQ appropriately.
+
+The contribution of Q to a resolution of width dQo at point Qo is:
+
+   p(Q) = 1/sqrt(2 pi dQo^2) exp ( (Q-Qo)^2/(2 dQo^2) )
+
+We are approximating the convolution at Qo using a numerical
+approximation to the integral over the measured points.  For 
+efficiency, the integral is limited to p(Q_i)/p(0)>=0.001.  
+
+Note that the function we are convoluting is falling off as Q^4.
+That means the correct convolution should uniformly sample across
+the entire width of the Gaussian.  This is not possible at the
+end points unless you calculate the reflectivity beyond what is
+strictly needed for the data. For a given dQ and step size,
+you need enough steps that the following is true:
+
+    (n*step)^2 > -2 dQ^2 * ln 0.001
+
+The choice of sampling density is particularly important near the 
+critical edge.  This is where the resolution calculation has the 
+largest effect on the reflectivity curve. In one particular model, 
+calculating every 0.001 rather than every 0.02 changed one value 
+above the critical edge by 15%.  This is likely to be a problem for 
+any system with a well defined critical edge.  The solution is to
+compute the theory function over a finer mesh where the derivative
+is changing rapidly.  For the critical edge, I have found a sampling
+density of 0.005 to be good enough.
+
+For systems involving thick layers, the theory function oscillates 
+rapidly around the measured points.  This is a problem when the
+period of the oscillation, 2 pi/d for total sample depth d, is on
+the order of the width of the resolution function. This is true even 
+for gradually changing profiles in materials with very high roughness
+values.  In these systems, the theory function should be oversampled
+around the measured points Q.  With a single thick layer, oversampling
+can be limited to just one period.  With multiple thick layers,
+oscillations will show interference patterns and it will be necessary 
+to oversample uniformly between the measured points.  When oversampled
+spacing is less than about 2 pi/7 d, it is possible to see aliasing
+effects.  
+
+FIXME is it better to use random sampling or strictly
+regular spacing when you are undersampling?
+
+===============================================================
+*/
+
+#undef USE_TRAPEZOID_RULE
+#ifdef USE_TRAPEZOID_RULE
+#warning This code does strange things with small sigma and large spacing
+/* FIXME trapezoid performs very badly for large spacing unless
+   we normalize to the unit width.  For very small sigma, the gaussian
+   is a spike, but we are approximating it by a triangle so it is
+   not surprising it works so poorly.  A slightly better solution is
+   to use the inner limits rather than the outer limits, but this will
+   still break down if the Q spacing is approximately equal to limit. 
+   Best is to not use trapezoid.
+*/
+
+/* Trapezoid rule for numerical integration of convolution */
+double 
+convolve_point(const double Qin[], const double Rin[], int k, int n,
+        double Qo, double limit, double sigma)
+{
+  const double two_sigma_sq = 2. * sigma * sigma;
+  double z, Glo, RGlo, R, norm;
+  
+  z = Qo - Qin[k];
+  Glo = exp(-z*z/two_sigma_sq);
+  RGlo = Rin[k]*Glo;
+  norm = R = 0.;
+  while (++k < n) {
+    /* Compute the next endpoint */
+    const double zhi = Qo - Qin[k];
+    const double Ghi = exp(-zhi*zhi/two_sigma_sq);
+    const double RGhi = Rin[k] * Ghi;
+    const double halfstep = 0.5*(Qin[k] - Qin[k-1]);
+    
+    /* Add the trapezoidal area. */
+    norm += halfstep * (Ghi + Glo);
+    R += halfstep * (RGhi + RGlo);
+    
+    /* Save the endpoint for next trapezoid. */
+    Glo = Ghi;
+    RGlo = RGhi;
+    
+    /* Check if we've calculated far enough */
+    if (Qin[k] >= Qo+limit) break;
+  }
+
+  /* Scale to area of the linear spline distribution we actually used. */
+  return R / norm;
+
+  /* Scale to gaussian of unit area */
+  /* Fails badly for small sigma or large Q steps---do not use. */
+  /* return R / sigma*SQRT2PI; */
+}
+
+#else /* !USE_TRAPEZOID_RULE */
+
+/* Analytic convolution of gaussian with linear spline */
+/* More expensive but more reliable */
+double 
+convolve_point(const double Qin[], const double Rin[], int k, int n,
+               double Qo, double limit, double sigma)
+{
+  const double two_sigma_sq = 2. * sigma * sigma;
+  double z, Glo, erflo, erfmin, R;
+  
+  z = Qo - Qin[k];
+  Glo = exp(-z*z/two_sigma_sq);
+  erfmin = erflo = erf(-z/(SQRT2*sigma));
+  R = 0.;
+  /* printf("%5.3f: (%5.3f,%11.5g)",Qo,Qin[k],Rin[k]); */
+  while (++k < n) {
+  	/* No additional contribution from duplicate points. */
+  	if (Qin[k] == Qin[k-1]) continue;
+ 
+    /* Compute the next endpoint */
+    const double zhi = Qo - Qin[k];
+    const double Ghi = exp(-zhi*zhi/two_sigma_sq);
+    const double erfhi = erf(-zhi/(SQRT2*sigma));
+    const double m = (Rin[k]-Rin[k-1])/(Qin[k]-Qin[k-1]);
+    const double b = Rin[k] - m * Qin[k];
+
+    /* Add the integrals. */
+    R += 0.5*(m*Qo+b)*(erfhi-erflo) - sigma/SQRT2PI*m*(Ghi-Glo);
+
+    /* Debug computation failures. */
+    /*
+    if isnan(R) {
+    	printf("NaN from %d: zhi=%g, Ghi=%g, erfhi=%g, m=%g, b=%g\n",
+    	       k,zhi,Ghi,erfhi,m,b);
+    }
+    */
+    
+    /* Save the endpoint for next trapezoid. */
+    Glo = Ghi;
+    erflo = erfhi;
+    
+    /* Check if we've calculated far enough */
+    if (Qin[k] >= Qo+limit) break;
+  }
+  /* printf(" (%5.3f,%11.5g)",Qin[k<n?k:n-1],Rin[k<n?k:n-1]); */
+
+#ifdef USE_TRUNCATED_NORMALIZATION
+  /* Normalize by the area of the truncated gaussian */
+  /* At this point erflo = erfmax */
+  /* printf ("---> %11.5g\n",2*R/(erflo-erfmin)); */
+  return 2 * R / (erflo - erfmin);
+#else
+  /* Return unnormalized (we used a gaussian of unit area) */
+  /* printf ("---> %11.5g\n",R); */
+  return R;
+#endif
+}
+
+#endif /* !USE_TRAPEZOID_RULE */
+
+void
+resolution(int Nin, const double Qin[], const double Rin[],
+           int N, const double Q[], const double dQ[], double R[])
+{
+  int lo,out;
+
+  /* FIXME fails if Qin are not sorted; slow if Q not sorted */
+  assert(Nin>1);
+
+  /* Scan through all Q values to be calculated */
+  lo = 0;
+  for (out=0; out < N; out++) {
+    /* width of resolution window for Q is w = 2 dQ^2. */
+    const double sigma = dQ[out];
+    const double Qo = Q[out];
+    const double limit = sqrt(-2.*sigma*sigma* LOG_RESLIMIT);
+
+    /* if (out%20==0) printf("%d: Q,dQ = %g,%g\n",out,Qo,sigma); */
+
+    /* Line up the left edge of the convolution window */
+    /* It is probably forward from the current position, */
+    /* but if the next dQ is a lot higher than the current */
+    /* dQ or if the Q are not sorted, then it may be before */
+    /* the current position. */
+    /* FIXME verify that the convolution window is just right */
+    while (lo < Nin-1 && Qin[lo] < Qo-limit) lo++;
+    while (lo > 0 && Qin[lo] > Qo-limit) lo--;
+
+    /* Special handling to avoid 0/0 for w=0. */
+    if (sigma > 0.) {
+      R[out] = convolve_point(Qin,Rin,lo,Nin,Qo,limit,sigma);
+    } else if (lo < Nin-1) {
+      /* Linear interpolation */
+      double m = (Rin[lo+1]-Rin[lo])/(Qin[lo+1]-Qin[lo]);
+      double b = Rin[lo] - m*Qin[lo];
+      R[out] = m*Qo + b;
+    } else if (lo > 0) {
+      /* Linear extrapolation */
+      double m = (Rin[lo]-Rin[lo-1])/(Qin[lo]-Qin[lo-1]);
+      double b = Rin[lo] - m*Qin[lo];
+      R[out] = m*Qo + b;
+    } else {
+      /* Can't happen because there is more than one point in Qin. */
+      assert(Nin>1);
+    }
+  }
+
+}
Index: park-1.2.1/park/message.py
===================================================================
--- park-1.2.1/park/message.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/message.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,171 @@
+# This program is public domain
+"""
+Asynchronous message streams.
+
+When you have multiple listeners to a process, some of which can connect and
+disconnect at different times, you need to dispatch the incoming messages to 
+all the listeners.  When a listener joins a running stream, they need to get 
+an immediate status update for the computation.  This is stored as the 
+stream header. Then, without dropping any messages, the listeners should 
+receive all subsequent messages in the stream in order.
+
+The message stream is multi-channelled, and indexed by a key.  Within 
+the service framework the key is likely to be the jobid associated with 
+the message stream.
+
+The contents of the message stream are expected to be monitor messages
+(see `park.monitor` for details), with the stream header being a
+`park.monitor.Join` message.  When a new listener is registered,
+the header is immediately put on the queue (if there is one), then all
+subsequent message are sent until the listener calls hangup().
+
+To attach to a message stream you need an object which accepts put().
+An asynchronous queue object is a good choice since it allows you to
+buffer the messages in one thread and pull them off as needed in another.
+You can also use a `park.monitor.Monitor` object to process the messages
+directly.
+"""
+from __future__ import with_statement
+
+__all__ = ['message_stream']
+
+import thread
+
+
+# Design note: we cannot buffer the messages in the stream.  The only
+# logical buffer size in this case would be the entire stream history
+# and that may large.  Furthermore, leading messages in the stream are
+# not of value to the computation.  The final result will be available
+# for as long as the job is stored on the server, which is a lifetime
+# much longer than the message queue.
+
+class MessageStream:
+    """
+    Message streams.
+
+    For each active job on the system there is a message stream
+    containing the job header and a list of listening objects.
+    Listeners should accept a put() message to process the next
+    item in the stream.
+    """
+    def __init__(self):
+        self._lock = thread.allocate_lock()
+        self.listeners = {}
+        self.headers = {}
+
+    def __getitem__(self, stream):
+        """Get the header for message stream 'stream'."""
+        with self._lock:
+            return self.headers.get(stream, None)
+    def __setitem__(self, stream, message):
+        """Set the header for message stream 'stream'."""
+        with self._lock:
+            self.headers[stream] = message
+    def __delitem__(self, stream):
+        """Delete the message stream 'stream'."""
+        self.headers.pop(stream,None)
+        self.listeners.pop(stream,None)
+
+    def put(self, stream, message):
+        """
+        Put a message on stream 'stream', transfering it to all listening queues.
+        """
+        #print "post message on",stream,message
+        with self._lock:
+            #print "post",stream,message,"->",self.listeners.get(stream,[])
+            for queue in self.listeners.get(stream,[]):
+                queue.put(message)
+
+    def listen(self, stream, queue):
+        """
+        Listen to message stream 'stream', adding all new messages to queue.
+
+        The stream header will be the first message queued.
+        """
+        #print "listening on",stream
+        with self._lock:
+            queues = self.listeners.setdefault(stream,[])
+            queues.append(queue)
+            #print "listen",stream,"->",self.listeners.get(stream,[])
+            # Make sure that the Join header is the first item in the queue.
+            header = self.headers.get(stream, None)
+            if header is not None:
+                queue.put(header)
+
+    def hangup(self, stream, queue):
+        """
+        Stop listening to message stream 'stream' with queue.
+
+        If stream is None then remove queue from all message streams.
+        """
+        with self._lock:
+            if stream is not None:
+                try:
+                    queuelist = self.listeners[stream]
+                    queuelist.remove(queue)
+                    if queuelist == []: del self.listeners[stream]
+                except KeyError:
+                    pass
+                except ValueError:
+                    pass
+            else:
+                purge = []
+                for stream,queuelist in self.listeners.iteritems():
+                    try:
+                        queuelist.remove(queue)
+                        if queuelist == []: purge.append(stream)
+                    except ValueError:
+                        pass
+                for stream in purge: del self.listeners[stream]
+
+message_stream = MessageStream()
+stream = message_stream # so message.stream works
+
+
+def demo():
+    import Queue
+    import time
+    t0 = time.time()
+    class NamedQueue(Queue.Queue):
+        def __init__(self, name):
+            Queue.Queue.__init__(self)
+            self.name = name
+        def __str__(self): return self.name
+        def __repr__(self): return "Queue('%s')"%self.name
+    def process_queue(queue):
+        print ">>>",queue
+        while True:
+            value = queue.get()
+            print "recv",queue,":",value[2],"<-",value[:2],"time",time.time()-t0
+    def post(id, stream, deltas):
+        for k,t in enumerate(deltas):
+            time.sleep(t)
+            print "       send ",id,":",stream,"<-",(id,k)
+            message_stream.put(stream,(id,k,stream))
+    def listen(stream,queue):
+        print "+++",queue,":",stream
+        message_stream.listen(stream,queue)
+    def hangup(stream,queue):
+        print "---",queue,":",stream
+        message_stream.hangup(stream,queue)
+
+    q1 = NamedQueue('Q1')
+    q2 = NamedQueue('Q2')
+    q3 = NamedQueue('Q3')
+    for q in [q1,q2,q3]: thread.start_new_thread(process_queue,(q,))
+    message_queue['M2'] = ('H2',0,'M2')
+    thread.start_new_thread(post,('S1','M1',[1,2,1,2,3,2]))
+    thread.start_new_thread(post,('S2','M1',[1,2,1,2,3,2]))
+    thread.start_new_thread(post,('S3','M2',[1,2,1,2,3,2]))
+    thread.start_new_thread(post,('S4','M2',[1,2,1,2,3,2]))
+    listen('M1',q1)
+    listen('M1',q3)
+    listen('M2',q3)
+    time.sleep(5)
+    hangup(None,q3)
+    hangup('M1',q1)
+    listen('M2',q2)
+    time.sleep(5)
+    time.sleep(15)
+
+if __name__ == "__main__": demo()
Index: park-1.2.1/park/model.py
===================================================================
--- park-1.2.1/park/model.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/model.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,511 @@
+# This program is public domain
+"""
+Define a park fitting model.
+
+Usage
+-----
+
+The simplest sort of fitting model is something like the following::
+
+    import numpy
+    import park
+    def G(x,mu,sigma):
+        return numpy.exp(-0.5*(x-mu)**2/sigma**2)
+
+    class Gauss(park.Model):
+        parameters = ['center','width','scale']
+        def eval(self,x):
+            return self.scale * G(x,self.center,self.width)
+
+It has a function which is evaluated at a series of x values and
+a set of adjustable parameters controlling the shape of f(x).
+
+You can check your module with something like the following::
+
+    $ ipython -pylab
+
+    from gauss import Gauss
+
+    g = Gauss(center=5,width=1,scale=10)
+    x = asarray([1,2,3,4,5])
+    y = g(x)
+    plot(x,y)
+
+This should produce a plot of the Gaussian peak.
+
+You will then want to try your model with some data.  Create a file
+with some dummy data, such as gauss.dat::
+
+    # x y
+    4 2
+    5 6
+    6 7
+    7 3
+
+In order to match the model to data, you need to define a fitness
+object.  This shows the difference between the model and the data,
+which you can then plot or sum to create a weighted chisq value::
+
+    f = park.Fitness(g, 'gauss.dat')
+    plot(f.data.fit_x, f.residuals())
+
+The data file can have up to four columns, either x,y or x,y,dy
+or x,dx,y,dy where x,y are the measurement points and values,
+dx is the instrument resolution in x and dy is the uncertainty
+in the measurement value y.  You can pass any keyword arguments
+to data.load that are accepted by numpy.loadtxt.  For example,
+you can reorder columns or skip rows.  Additionally, you can modify
+data attributes directly x,y,dx,dy and calc_x.  See help on park.Data1D
+for details.
+
+Once you have your model debugged you can use it in a fit::
+
+    g = Gauss(center=[3,5],scale=7.2,width=[1,3])
+    result = park.fit((g, 'gauss.dat'))
+    result.plot()
+
+In this example, center and width are allowed to vary but scale is fixed.
+
+Existing models can be readily adapted to Park::
+
+    class Gauss(object):
+        "Existing model"
+        center,width,scale = 0,1,1
+        def __init__(self,**kw):
+            for k,v in kw.items(): setattr(self,k,v)
+        def eval(self,x):
+            return self.scale *G(x,self.center,self.width)
+
+    class GaussAdaptor(Gauss,Model):
+        "PARK adaptor"
+        parameters = ['center','width','scale']
+        def __init__(self,*args,**kw):
+            Model.__init__(self,*args)
+            Gauss.__init__(self,**kw)
+
+    g = GaussAdaptor(center=[3,5],scale=7.2,width=[1,3])
+    result = park.fit((g, 'gauss.dat'))
+    result.plot()
+
+Models can become much more complex than the ones described above,
+including multilevel models where fitting parameters can be added
+and removed dynamically.
+
+In many cases the park optimizer will need an adaptor for pre-existing
+models.  The adaptor above relies on python properties to translate
+model.par access into model._par.get() and model._par.set() where _par
+is the internal name for par.  This technique works for simple static
+models, but will not work well for sophisticated models which have,
+for example, a dynamic parameter set where the model parameters cannot
+be set as properties.  A solution to this problem is to subclass the
+park.Parameter and override the value attribute as a property.
+
+Once models are defined they can be used in a variety of contexts, such
+as simultaneous fitting with constraints between the parameters.  With
+some care in defining the model, computationally intensive fits can
+be distributed across multiple processors.  We provide a simple user
+interface for interacting with the model parameters and managing fits.
+This can be extended with model specialized model editors which format
+the parameters in a sensible way for the model, or allow direct manipulation
+of the model structure.  The underlying fitting engine can also be
+used directly from your own user interface.
+
+"""
+
+__all__ = ['Model']
+
+import park
+from park.parameter import Parameter, ParameterSet
+from copy import copy, deepcopy
+
+
+class ParameterProperty(object):
+    """
+    Attribute accessor for a named parameter.
+
+    Objects of class ParameterProperty act similarly to normal property
+    objects in that assignment to object.attr triggers the __set__
+    method of ParameterProperty and queries of the value object.attr
+    triggers the __get__ method.  These methods look up the actual
+    parameter in the dictionary model.parameterset, delegating the
+    the set/get methods of the underlying parameter object.
+
+    For example::
+
+        model.P1 = 5
+        print model.P1
+
+    is equivalent to::
+
+        model.parameterset['P1'].set(5)
+        print model.parameterset['P1'].get()
+
+    To enable this behaviour, the model developer must assign a
+    ParameterProperty object to a class attribute of the model.  It
+    must be a class attribute.  Properties assigned as an instance
+    attribute in the __init__ of the class will not be recognized
+    as properties.
+
+    An example model will look something like the following::
+
+        class MyModel:
+            # A property must be assigned as a class attribute!
+            P1 = ParameterProperty('P1')
+            P2 = ParameterProperty('P2')
+            def __init__(self, P1=None, P2=None):
+                parP1 = Parameter('P1')
+                if P1 is not None: parP1.set(P1)
+                parP2 = Parameter('P2')
+                if P2 is not None: parP2.set(P2)
+                self.parameterset = { 'P1': parP1, 'P2': parP2 }
+
+    This work is done implicitly by MetaModel, and by subclassing
+    the class Model, the model developer does not ever need to
+    use ParameterProperty.
+    """
+    def __init__(self,name,**kw):
+        self.name = name
+    def __get__(self,instance,cls):
+        return instance.parameterset[self.name].get()
+    def __set__(self,instance,value):
+        instance.parameterset[self.name].set(value)
+
+
+class MetaModel(type):
+    """
+    Interpret parameters.
+
+    This is a meta class, and the usual meta class rules apply.  That is,
+    the Model base class should be defined like::
+
+        class Model(object):
+            __metaclass__ = MetaModel
+            ...
+
+    The MetaModel.__new__ method is called whenever a new Model
+    class is created.  The name of the model, its superclasses and
+    its attributes are passed to MetaModel.__new__ which creates the
+    actual class.  It is not called for new instances of the model.
+
+    MetaModel is designed to simplify the definition of parameters for
+    the model.  When processing the class, MetaModel defines
+    parameters which is a list of names of all the parameters in the
+    model, and parameterset, which is a dictionary of the actual
+    parameters and any parameter sets in the model.
+
+    Parameters can be defined using a list of names in the parameter
+    attribute, or by defining the individual parameters as separate
+    attributes of class Parameter.
+
+    For example, the following defines parameters P1, P2, and P3::
+
+        class MyModel(Model):
+            parameters = ['P1','P2']
+            P2 = Parameter(range=[0,inf])
+            P3 = Parameter()
+
+    For each parameter, MetaModel will define a parameter accessor,
+    and add the parameter definition to the parameter set. The accessor
+    delegates query and assignment to the Parameter get/set methods. The
+    attributes of the particular parameter instance can be
+    adjusted using model.parameterset['name'].attribute = value.
+    """
+    def __new__(cls, name, bases, vars):
+        #print 'calling model meta',vars
+
+        # Find all the parameters for the model.  The listed parameters
+        # are defined using::
+        #    parameters = ['P1', 'P2', ...]
+        # The remaining parameters are defined individually::
+        #    P3 = Parameter()
+        #    P4 = Parameter()
+        # The undefined parameters are those that are listed but not defined.
+        listed = vars.get('parameters',[])
+        defined = [k for k,v in vars.items() if isinstance(v,Parameter)]
+        undefined = [k for k in listed if k not in defined]
+        #print listed,defined,undefined
+
+        # Define a getter/setter for every parameter so that the user
+        # can say model.name to access parameter name.
+        #
+        # Create a parameter object for every undefined parameter.
+        # Check if the base class defines a default value, and use
+        # that for the initial value.  We don't want to do this for
+        # parameters explicitly defined since the user may have
+        # given them a default value already.
+        #
+        # For predefined parameters we must set the name explicitly.
+        # This saves the user from having to use, e.g.::
+        #     P1 = Parameter('P1')
+        pars = []
+        for p in undefined:
+            # Check if the attribute value is initialized in a base class
+            for b in bases:
+                if hasattr(b,p):
+                    value = getattr(b,p)
+                    break
+            else:
+                value = 0.
+            #print "looking up value in base as",value
+            pars.append(Parameter(name=p,value=value))
+            vars[p] = ParameterProperty(p)
+        for p in defined:
+            # Make sure parameter name is correct.  Note that we are using
+            # _name rather than .name to access the name attribute since
+            # name is a read-only parameter.
+            vars[p]._name = p
+            pars.append(vars[p])
+            vars[p] = ParameterProperty(p)
+
+        # Construct the class attributes 'parameters' and 'parameterset'.
+        # Listed parameters are given first, with unlisted parameters
+        # following alphabetically.  For hierarchical parameter sets,
+        # we also need to include the defined sets.
+        unlisted = list(set(defined+undefined) - set(listed))
+        unlisted.sort()
+        parameters = listed + unlisted
+        parsets = [k for k,v in vars.items() if isinstance(v,ParameterSet)]
+        vars['parameters'] = parameters
+        vars['parameterset'] = ParameterSet(pars=pars+parsets)
+        #print 'pars',vars['parameters']
+        #print 'parset',vars['parameterset']
+
+        # Return a new specialized instance of the class with parameters
+        # and parameterset made explicit.
+        return type.__new__(cls, name, bases, vars)
+
+class Model(object):
+    """
+    Model definition.
+
+    The model manages attribute access to the fitting parameters and
+    also manages the dataset.
+
+    derivatives ['p1','p2',...]
+        List of parameters for which the model can calculate
+        derivatives.  The derivs
+        The model function can compute the derivative with respect
+        to this parameter.  The function model.derivs(x,[p1,p2,...])
+        will return (f(x),df/dp1(x), ...).  The parameters and their
+        order are determined by the fitting engine.
+
+        Note: This is a property of the model, not the fit.
+        Numerical derivatives will be used if the parameter is
+        used in an expression or if no analytic derivative is
+        available for the parameter.  Automatic differentiation
+        on parameter expressions is possible, but beyond the scope
+        of this project.
+
+    eval(x)
+        Evaluate the model at x.  This must be defined by the subclass.
+
+    eval_deriv(x,pars=[])
+        Evaluate the model and the derivatives at x.  This must be
+        defined by the subclass.
+
+    parameters
+        The names of the model parameters.  If this is not provided, then
+        the model will search the subclass for park.Parameter attributes
+        and construct the list of names from that.  Any parameters in the
+        list not already defined as park.Parameter attributes will be
+        defined as parameters with a default of 0.
+
+    parameterset
+        The set of parameters defined by the model.  These are the
+        parameters themselves, gathered into a park.ParameterSet.
+
+    The minimum fittng model if you choose not to subclass park.Model
+    requires parameterset and a residuals() method.
+    """
+    __metaclass__ = MetaModel
+    derivatives = []
+    def __init__(self,*args,**kw):
+        """
+        Initialize the model.
+
+        Model('name',P1=value, P2=Value, ...)
+
+        When overriding __init__ in the subclass be sure to call
+        Model.__init__(self, *args, **kw).  This makes a private
+        copy of the parameterset for the model and initializes
+        any parameters set using keyword arguments.
+        """
+        #print 'calling model init on',id(self)
+        # To avoid trashing the Model.parname = Parameter() template
+        # when we create an instance and accessor for the parameter
+        # we need to make sure that we are using our own copy of the
+        # model dictionary stored in vars
+        if len(args)>1: raise TypeError("expect name as model parameter")
+        name = args[0] if len(args)>=1 else 'unknown'
+        self.parameterset = deepcopy(self.parameterset)
+        for k,v in kw.items():
+            self.parameterset[k].set(v)
+        self.parameterset._name = name
+        #print 'done',id(self)
+
+    def __call__(self, x, pars=[]):
+        return self.eval(x)
+
+    def eval(self, x):
+        """
+        Evaluate the model at x.
+
+        This method needs to be specialized in the model to evaluate the
+        model function.  Alternatively, the model can implement is own
+        version of residuals which calculates the residuals directly
+        instead of calling eval.
+        """
+        raise NotImplementedError('Model needs to implement eval')
+
+    def eval_derivs(self, x, pars=[]):
+        """
+        Evaluate the model and derivatives wrt pars at x.
+
+        pars is a list of the names of the parameters for which derivatives
+        are desired.
+
+        This method needs to be specialized in the model to evaluate the
+        model function.  Alternatively, the model can implement is own
+        version of residuals which calculates the residuals directly
+        instead of calling eval.
+        """
+        raise NotImplementedError('Model needs to implement eval_derivs')
+
+    def set(self, **kw):
+        """
+        Set the initial value for a set of parameters.
+
+        E.g., model.set(width=3,center=5)
+        """
+        # This is a convenience funciton for the user.
+        #
+        for k,v in kw.items():
+            self.parameterset[k].set(v)
+
+
+def add_parameter(model, name, **kw):
+    """
+    Add a parameter to an existing park model.
+
+    Note: this may not work if it is operating on a BaseModel.
+    """
+    setattr(model.__class__, name, park.model.ParameterProperty(name))
+    model.parameterset.append(park.Parameter(name=name, **kw))
+
+
+def test():
+    # Gauss theory function
+    import numpy
+    eps,inf = numpy.finfo('d').eps,numpy.inf
+    def G(x,mu,sigma):
+        mu,sigma = numpy.asarray(mu),numpy.asarray(sigma)
+        return numpy.exp(-0.5*(x-mu)**2/sigma**2)
+
+    # === Minimal model ===
+    # just list the fitting parameters and define the function.
+    class Gauss(Model):
+        parameters = ['center','width','scale']
+        def eval(self,x):
+            return self.scale * G(x,self.center,self.width)
+
+    # Create a couple of models and make sure they don't conflict
+    g1 = Gauss(center=5,width=1,scale=2)
+    assert g1.center == 5
+    g2 = Gauss(center=6,width=1)
+    assert g1.center == 5
+    assert g2.center == 6
+    assert g1(5) == 2
+    assert g2(6) == 0
+
+    # === explore parameters ===
+    # dedicated model using park parameters directly, and defining derivatives
+    class Gauss(Model):
+        center = Parameter(value=0,
+                           tip='Peak center')
+        width = Parameter(value=1,
+                          limits=(eps,inf),
+                          tip='Peak width (1-sigma)')
+        scale = Parameter(value=1,
+                          limits=(0,inf),
+                          tip='Peak scale (>0)')
+        def eval(self,x):
+            return self.scale * G(x,self.center,self.width)
+
+        # Derivatives
+        derivatives = ['center','width','scale']
+        def dscale(self, g, x):
+            return g/self.scale
+        def dcenter(self, g, x):
+            return 2*(x-self.center)/self.width*g
+        def dwidth(self, g, x):
+            return 2*(x-self.center)/self.width**3*g
+        dmap = dict(scale=dscale,center=dcenter,width=dwidth)
+        def eval_derivs(self,x,pars=[]):
+            """
+            Calculate function value and derivatives wrt to the parameters.
+
+            pars is a list of parameter names, possibly consisting of any
+            parameter with deriv=True.
+            """
+            g = self.eval(x)
+            dg = [self.dmap[p](self,g,x) for p in pars]
+            return [g]+dg
+
+    g1 = Gauss(center=5)
+    g1.parameterset['center'].tip = 'This is the center'
+    assert g1.center == 5
+    g2 = Gauss(center=6)
+    assert g1.center == 5
+    assert g2.center == 6
+    assert g1(5) == 1
+    assert g2(6) == 1
+
+    # ====== Test wrapper =======
+    # delegate to existing model via inheritence
+    class Gauss(object):
+        """Pre-existing model"""
+        center,width,scale = 0,1,1
+        def __init__(self,**kw):
+            #print "calling BaseGauss init on",id(self)
+            for k,v in kw.items(): setattr(self,k,v)
+            #print "done",id(self)
+        def eval(self,x):
+            return self.scale *G(x,self.center,self.width)
+
+    class GaussAdaptor(Gauss,Model):
+        """PARK wrapper"""
+        parameters = ['center','width','scale']
+        def __init__(self,*args,**kw):
+            #print "calling Gauss init on",id(self)
+            Model.__init__(self,*args)
+            Gauss.__init__(self,**kw)
+            #print "done",id(self)
+
+    g1 = GaussAdaptor(center=5)
+    g2 = GaussAdaptor(center=6)
+    g3 = GaussAdaptor()
+    assert g1.center == 5
+    assert g2.center == 6
+    assert g3.scale == 1
+    assert g1(5) == 1
+    assert g2(6) == 1
+
+    g4 = GaussAdaptor('M3',center=6)
+    assert g4.parameterset.name == 'M3'
+
+    # dedicated multilevel model using park parameters directly
+    class MultiGauss(Model):
+        def add(self,name,model):
+            pass
+
+    # wrapped model using park parameters indirectly
+    class BaseMultiGauss(object):
+        def __init__(self):
+            self.models = []
+        def add(self,**kw):
+            self.models.append(BaseGauss(**kw))
+    class WrapMultiGauss(BaseMultiGauss,Model):
+        def __init__(self):
+            pass
+
+if __name__ == "__main__": test()
Index: park-1.2.1/park/monitor.py
===================================================================
--- park-1.2.1/park/monitor.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/monitor.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,250 @@
+# This program is public domain
+"""
+Asychronous execution monitoring service.
+
+Long running computations need to convey status information to the user.
+This status can take multiple forms, such as output to the console or
+activity on a GUI, or even mail to your inbox.
+
+park.monitor defines several standard message types::
+
+    `Start` for job start
+    `Join` first message when joining an already running job
+    `Progress` for job activity
+    `Improvement` for partial results
+    `Complete` for final result
+    `Abort` when job is killed
+    `Error` when job has an error
+    `Log` for various debugging messages
+
+Individual services may have specialized message types.
+
+park.monitor also defines `Monitor` to process the various kinds of messages,
+and dispatch them to the various user defined handlers.
+
+For each message type, the Monitor dispatcher will look for a function
+named onMonitorQQQ where QQQ is the message type.  For example,
+onMonitorStart(self, message) will be called in response to a Start message.
+If onMonitorQQQ is not defined, then onMonitorMessage will be called.  The
+default behaviour of onMonitorMessage is to print the message on the console.
+
+Log messages are sent to the standard system logger.  See logging in the
+python standard library for details.
+
+The Monitor class has methods for onMonitorStart(message), etc.
+In panel, be sure to have methods for onMonitorStart(message), 
+onMonitorProgress(message), etc., for the kinds of monitor messages 
+the application will send.  The catch-all method is onMonitorMessage.
+
+See `park.monitor` for details on the message types.  Individual services
+may have additional message types.
+
+"""
+__all__ = ['Monitor']
+
+import sys
+import logging
+import traceback
+
+class Message(object):
+    """
+    Message type
+    """
+
+class Start(Message):
+    """
+    Start.
+
+    Sent when the job has started processing.
+    """
+    def __str__(self): return "Start"
+
+class Join(Message):
+    """
+    Join: k units of n with partial result
+    
+    Sent when the listener is attached to a running job.  This is
+    a combination of Progress and Improvement.
+    """
+    def __init__(self, k, n, partial):
+        self.total = n
+        """Total work to complete"""
+        self.complete = k
+        """Amount of work complete"""
+        self.result = partial
+        """The partial result completed; this is job specific"""
+    def __str__(self): return "Join: "+str(self.result)
+    
+class Progress(Message):
+    """
+    Progress: k units of n.
+
+    Sent when a certain amount of progress has happened.
+
+    Use the job controller to specify the reporting
+    frequency (time and/or percentage).
+    """
+    def __init__(self, k, n, units=None):
+        self.total = n
+        """Total work to complete"""
+        self.complete = k
+        """Amount of work complete"""
+        self.units = units
+        """Units of work, or None"""
+    def __str__(self):
+        if self.units is not None:
+            return "Progress: %s %s of %s"%(self.complete, self.units, self.total)
+        else:
+            return "Progress: %s of %s"%(self.complete, self.total)
+
+class Improvement(Message):
+    """
+    Improvement: partial result.
+
+    Use the job controller to specify the improvement frequency
+    (time and/or percentage).
+    """
+    def __init__(self, partial):
+        self.result = partial
+        """The partial result completed; this is job specific"""
+    def __str__(self):
+        return "Improvement: "+str(self.result)
+
+class Complete(Message):
+    """
+    Complete: final result.
+    """
+    def __init__(self, final):
+        self.result = final
+        """The final completed result; this is job specific"""
+    def __str__(self):
+        return "Complete: "+str(self.result)
+
+class Error(Message):
+    """
+    Traceback stack trace.
+    """
+    def __init__(self, trace=None):
+        if trace == None: trace = sys.exc_info()
+        self.trace = trace
+        """The stack trace returned from exc_info()"""
+    def __str__(self):
+        #print "traceback",traceback.format_exception(*self.trace)
+        try:
+            return "".join(traceback.format_exception(*self.trace))
+        except TypeError:
+            return "Error: "+str(self.trace)
+
+class Abort(Message):
+    """
+    Abort: partial result
+
+    Use the job controller to signal an abort.
+    """
+    def __init__(self, partial):
+        self.result = partial
+        """The partial result completed; this is job specific"""
+    def __str__(self):
+        return "Abort: "+str(self.result)
+
+class Log(Message):
+    """
+    Log module.function: log record
+    """
+    formatter = logging.Formatter("Log %(module)s.%(funcName)s: %(message)s")
+    def __init__(self, record):
+        self.record = record
+        """The partial result completed; this is job specific"""
+    def __str__(self):
+        return self.formatter.format(self.record)
+
+class Monitor(object):
+    """
+    Messages that are received during the processing of the job.
+
+    Standard message types::
+
+        `Start`, `Progress`, `Improvement`, `Complete`, `Error`, `Abort`, `Log`
+
+    Specific job types may have their own monitor messages.
+
+    The messages themselves should all produce nicely formatted results
+    in response to str(message).
+
+    The message dispatch calls on<Class>(message) if the on<Class>
+    method exists for the message type.  If not, then dispatch
+    calls otherwise(message).  By default onLog(message) submits the
+    log record to the logger.
+
+    Subclass Monitor to define your own behaviours.
+    """
+    def put(self, message):
+        """
+        Called from thread when new message has arrived.
+        """
+        fn = getattr(self, 
+                     "onMonitor"+message.__class__.__name__, 
+                     self.onMonitorMessage)
+        fn(message)
+
+    def onMonitorMessage(self, message):
+        """
+        What to do if the message handler is not found.
+
+        Default is to ignore the message.
+        """
+        print ">",str(message)
+
+    def onMonitorLog(self, message):
+        """
+        Called when the job sends a logging record.
+
+        The logging record contains a normal python logging record.
+
+        The default behaviour is to tie into the application logging
+        system using::
+
+            logger = logging.getLogger(message.record.name)
+            logger.handle(message.record)
+
+        Logging levels are set in the job controller.
+        """
+        logging.basicConfig()
+        logger = logging.getLogger(message.record.name)
+        logger.handle(message.record)
+
+
+def demo(rate=0):
+    import sys, time, thread, logging
+    import park.monitor
+
+    monitor = Monitor()
+    def messagestream(monitor,rate,stream):
+        for m in stream:
+            time.sleep(rate)
+            monitor.put(m)
+        time.sleep(rate)
+    R = logging.LogRecord('hi',60,'hello.py',3,'log message',(),None,'here')
+    try: raise Exception('Test exception')
+    except: trace = sys.exc_info()
+    stream=[park.monitor.Start(),
+            park.monitor.Progress(1,10),
+            park.monitor.Progress(2,10),
+            park.monitor.Progress(3,10),
+            park.monitor.Join('Good'),
+            park.monitor.Improvement('Better!'),
+            park.monitor.Abort('Abandoned'),
+            park.monitor.Start(),
+            park.monitor.Progress(1,10,'seconds'),
+            park.monitor.Improvement('Better!'),
+            park.monitor.Progress(8,10),
+            park.monitor.Complete('Best!'),
+            park.monitor.Start(),
+            park.monitor.Log(R),
+            park.monitor.Progress(6,10),
+            park.monitor.Error(trace)]
+    thread.start_new_thread(messagestream, (monitor,rate,stream))
+
+    time.sleep(20*(rate+0.01))
+
+if __name__ == "__main__": demo(rate=0.1)
Index: park-1.2.1/park/parameter.py
===================================================================
--- park-1.2.1/park/parameter.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/parameter.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,542 @@
+# This program is public domain
+"""
+Parameters and parameter sets.
+
+Parameter defines an individual parameter, and ParameterSet groups them
+into a hierarchy.
+
+Individual models need to provide a parameter set with the correct
+properties, either by using park.ParameterSet in their model definition,
+or by providing a wrapper which can translate assignment to parameter.value
+into the appropriate change in the wrapped model.  See wrapper.py for
+an example.
+"""
+
+__all__ = ['Parameter', 'ParameterSet']
+
+import numpy
+import expression        
+
+class Pnormal(object):
+    """
+    Negative log likelihood function for a parameter from a Gaussian
+    distribution.
+    
+    Given P(v;mu,sigma) = exp(-1/2 (mu-v)^2/sigma^2)/sqrt(2 pi sigma^2)
+    then -log(P) = 1/2 (mu-v)^2/sigma^2 + 1/2 log(2*pi*sigma^2)
+    
+    Assuming that parameter P is selected from a normal distribution,
+    then P.likelihood = Pnormal(mu,sigma)
+    """
+    def __init__(self, mean, std):
+        self.mean = mean
+        self.std = std
+        self.const = math.log(2*math.pi*std**2)/2
+    def __call__(self, value):
+        return ((value-self.mean)/self.std)**2/2 + self.const
+    
+inf = numpy.inf
+class Parameter(object):
+    """
+    A parameter is a box for communicating with the fitting service.
+    Parameters can have a number of properties, 
+    
+    Parameters have a number of properties:
+    
+    name "string"
+        name of the parameter within the parameter set.  
+        
+        The name is read only.  You can rename a parameter but only 
+        in the context of the parameter set which contains it, using
+        parameterset.rename(par,name).  This will change all expressions
+        containing the named parameter.
+    
+    path
+        dotted name of the parameter within the set of models.  The
+        dotted name is automatically generated by the parameter set
+        before expressions are parsed and evaluated.  There are
+        some operations on parameter sets (such as renaming the
+        layer containing a parameter) which will force an adjustment
+        of all the underlying parameter names, as well as any
+        expressions in which they are referenced.
+    
+    limits (low, high)
+        hard limits on the range of allowed parameter values, dictated
+        by physics rather than by knowledge of the particular system.
+        For example, thickness parameters would have limits (0,inf)
+        because negative thicknesses are unphysical.  These limits
+        are enforced when setting range for the fit.
+
+    units "string"
+        units for the parameter.  This should be a string, but
+        be parsable by whatever units package your application
+        supports.
+
+    tip "string"
+        parameter description, suitable for use in a tooltip
+    
+    value double
+        current value of the parameter, either fixed, fitted or computed
+    
+    range (low, high)
+        range of expected values for the parameter in the model
+    
+    expression "string"
+        expression for the parameter in the model.  This is a string
+        containing a formula for updating the parameter value based
+        on the values of other parameters in the system.  The expression
+        is ignored if 'calculated' is False.
+
+        Note: the list of symbols available to the expression evaluator
+        defaults to the contents of the math module.  The caller will be
+        able to override this within the fitting fitting class.
+    
+    status 'fixed'|'computed'|'fitted'
+        the parameter type.  Choose 'fixed' if the values is to
+        remain fixed throughout the fit, even if a range and an
+        expression have been specified.  Choose 'computed' if the
+        value of the parameter is to be computed each time the
+        parameters are updated.  Choose 'fitted' if an optimization
+        algorithm is supposed to search the parameter space.
+
+    likelihood
+        function to return the negative log likelihood of seeing a
+        particular parameter value.  2*likelihood(value) will be added
+        to the total cost function for the particular parameter set
+        during the fit.  This will be on par with the probabilty
+        of seeing the particular theory function given the observed
+        datapoints when performing the fit (the residual term is
+        closely related to the log likelihood of the normal distribution).
+        
+        Note: we are minimizing chi^2 = sum [ ((y-f(x;p))/dy)^2 ] rather 
+        than  -log P = sum [ ((y-f(x;p))/dy)^2/2 + log(2 pi dy^2) ],
+        where P is the probability of seeing f(x;p) given y,dy as the
+        mean and standard deviation of a normal distribution.  Because
+        chi^2_p = - 2 * log P_p + constant, the minima of p are the same
+        for chi^2 and negative log likelihood.  However, to weight the
+        likelihood properly when adding likelihood values to chisq, we
+        need the extra factor of 2 mentioned above.  The usual statistical
+        implications of normalized chi^2 will of course be suspect, both
+        because the assumption of independence between the points in
+        chi^2 (which definitely do not hold for the new 'point' p_k), and
+        because of the additional 2 log(2 pi dp_k^2) constant, but given
+        the uncertainty in the estimate of the distribution parameters,
+        this is likely a minor point.
+    """
+    # Protect parameter name from user modification
+    _name = "unknown"
+    def _getname(self): return self._name
+    name = property(_getname,doc="parameter name")
+    
+    # Do checking on the parameter range to make sure the model stays physical
+    _range = (-inf,inf)
+    def _setrange(self, r):
+        if self.limits[0]<=r[0]<=r[1] <=self.limits[1]:
+            self._range = r
+        else:
+            raise ValueError("invalid range %s for %s"%(r,self.name))
+    def _getrange(self): return self._range
+    range = property(_getrange,_setrange,doc="parameter range")
+
+    path = ""
+    value = 0.
+    limits = (-inf,inf)
+    expression = ""
+    status = 'fixed' # fixed, computed or fitted
+    likelihood = None
+    units = ""
+    tip = "Fitting parameter"
+    deriv = False
+    def __init__(self, name="unknown", **kw):
+        self._name = name
+        for k,v in kw.items():
+            if hasattr(self,k): setattr(self,k,v)
+            else: raise AttributeError("Unknown attribute %s"%k)
+    def __str__(self): return self.path if self.path != '' else self.name
+    def __repr__(self): return "Parameter('%s')"%self.name
+    def summarize(self):
+        """
+        Return parameter range string.
+        
+        E.g.,  "       Gold .....|.... 5.2043 in [2,7]"
+        """
+        range = ['.']*10
+        lo,hi = p.range
+        portion = (p.value-lo)/(hi-lo)
+        if portion < 0: portion = 0.
+        elif portion >= 1: portion = 0.99999999
+        bar = math.floor(portion*len(range))
+        range[bar] = '|'
+        range = "".join(range)
+        return "%25s %s %g in [%g,%g]"  % (p.name,range,p.value,lo,hi)
+
+    def isfitted(self): return self.status == 'fitted'
+    def iscomputed(self): return self.status == 'computed'
+    def isfixed(self): return self.status == 'fixed'
+    def isrestrained(self): return self.likelihood is not None
+    def isfeasible(self, value):
+        """Return true if the value is in the range"""
+        return self._range[0] <= value <= self._range[1]
+    def setprefix(self, prefix):
+        """
+        Set the full path to the parameter as used in expressions involving
+        the parameter name.
+        """
+        self.path = prefix+self.name
+    def get(self):
+        """
+        Return the current value for a parameter.
+        """
+        return self.value
+    def set(self, value):
+        """
+        Set a parameter to a value, a range or an expression.  If it is a value,
+        the parameter will be fixed for the fit.  If it is a range, the value
+        will be varying for the fit.  If it is an expression, the parameter will
+        be calculated from the values of other parameters in the fit.
+
+        Raises ValueError if the value could not be interpreted.
+        """
+        
+        # Expression
+        if isinstance(value,basestring):
+            self.expression = value
+            self.status = 'computed'
+            return
+        
+        # Fixed value
+        if numpy.isscalar(value):
+            self.value = value
+            self.status = 'fixed'
+            return
+        
+        # Likelihood
+        if hasattr(value,'__call__'):
+            self.range = value.range
+            self.likelihood = value
+            self.status = 'fitted'
+            return
+        
+        # Range
+        try:
+            lo,hi = numpy.asarray(value)
+            if not numpy.isscalar(lo) or not numpy.isscalar(hi):
+                raise Exception
+            self.range = (lo,hi)
+            self.status = 'fitted'
+            return
+        except:
+            pass
+        
+        raise ValueError,\
+            "parameter %s expects value, expression or range: %s"%(self.name,value)
+
+class ParameterSet(list):
+    """
+    The set of parameters used to fit theory to data.
+    
+    ParameterSet forms a hierarchy of parameters.  The parameters themselves
+    are referred to by the path through the hierarchy, usually::
+    
+         fitname.component.parameter
+         
+    Though more or fewer levels are permitted.  Parameters are assumed to
+    have a unique label throughout the fit.  This is required so that
+    expressions tying the results of one fit to another can uniquely
+    reference a parameter.
+    
+    Attributes:
+    
+    name
+        the name of the parameter set
+    path
+        the full dotted name of the parameter set
+    context
+        a dictionary providing additional context for evaluating parameters;
+        Note that this namespace is shared with other theory functions, so
+        populate it carefully.
+    """
+    # Protect parameter set name from user modification
+    def _getname(self): return self._name
+    def _setname(self, name): 
+        raise NotImplementedError("parameter.name is protected; use fit.rename_parameter() to change it")
+    name = property(_getname,doc="parameter name")
+    path = ""
+    def __init__(self, name="unknown", pars=[]):
+        super(ParameterSet,self).__init__(pars)
+        self._name = name
+        self.context = {}
+
+    def _byname(self, parts):
+        """byname recursion function"""
+        if len(parts) == 1: return self
+        for p in self:
+            if parts[1] == p.name:
+                if len(pars) == 2:
+                    return p
+                elif isinstance(p, ParameterSet):
+                    return p._byname(parts[1:])
+        return None
+
+    def byname(self, name):
+        """Lookup parameter from dotted path"""
+        parts = name.split('.')
+        if parts[0] == self.name:
+            p =  _byname(self, name.split('.'))
+            if p: return p
+        raise KeyError("parameter %s not in parameter set"%name)
+    
+    def __getitem__(self, idx):
+        """Allow indexing by name or by number"""
+        if isinstance(idx,basestring):
+            for p in self:
+                if p.name == idx: 
+                    return p
+            raise KeyError("%s is not in the parameter set"%idx)
+        else:
+            return super(ParameterSet,self).__getitem__(idx)
+
+    def flatten(self):
+        """
+        Iterate over the elements in depth first order.
+        """
+        import park
+        L = []
+        for p in self:
+            # Yuck! I really only want to try flattening if p is a
+            # ParameterSet but it seems we have parameter.ParameterSet
+            # and park.parameter.ParameterSet as separate entities,
+            # depending on how things were included.  The solution is
+            # probably to force absolute include paths always.
+            try:
+                L += p.flatten()
+            except:
+                L.append(p)
+        return L
+        """
+        # Iterators are cute but too hard to use since you can
+        # only use them in a [p for p in generator()] once.
+        for p in self:
+            if isinstance(p, ParameterSet):
+                for subp in p.flatten(): 
+                    yield subp
+            else: 
+                yield p
+        """
+
+    def _fixed(self):
+        """
+        Return the subset of the parameters which are fixed
+        """
+        return [p for p in self.flatten() if p.isfixed()]
+    fixed = property(_fixed,doc=_fixed.__doc__)
+
+    def _fitted(self):
+        """
+        Return the subset of the paramters which are varying
+        """
+        return [p for p in self.flatten() if p.isfitted()]
+    fitted = property(_fitted,doc=_fitted.__doc__)
+
+    def _computed(self):
+        """
+        Return the subset of the parameters which are calculated
+        """
+        return [p for p in self.flatten() if p.iscomputed()]
+    computed = property(_computed,doc=_computed.__doc__)
+
+    def _restrained(self):
+        """
+        Return the subset of the parameters which have a likelihood
+        function associated with them.
+        """
+        return [p for p in self.flatten() if p.isrestrained()]
+    restrained = property(_restrained,doc=_restrained.__doc__)
+
+    def setprefix(self,prefix=None):
+        """
+        Fill in the full path name for all the parameters in the tree.
+        
+        Note: this function must be called from the root parameter set
+        to build proper path names.
+        
+        This is required before converting parameter expressions into
+        function calls.
+        """
+        if prefix == None:
+            # We are called from root, so we don't have a path
+            self.path = ""
+            prefix = ""
+        else:
+            self.path = prefix+self.name
+            prefix = self.path+'.'
+        for p in self:
+            #print "setting prefix for",p,prefix
+            p.setprefix(prefix)
+
+    def rename(self, par, name):
+        """
+        Rename the parameter to something new.
+        Called from root of the parameter hierarchy, rename the particular
+        parameter object to something else.
+
+        This changes the internal name of the parameter, as well as all
+        expressions in which it occurs.  If the parameter is actually
+        a parameter set, then it renames all parameters in the set.
+        """
+        
+        # Must run from root for self.setprefix and self.computed to work
+        if self.path != "":  
+            raise RuntimeError,"rename must be called from root parameter set"
+
+        # Change the name of the node
+        par._name = name
+        
+        # Determine which nodes (self and children) are affected
+        if isinstance(par,ParameterSet):
+            changeset = par.flatten()
+        else:
+            changeset = [par]
+        
+        #  Map the old names into the new names
+        old = [p.path for p in changeset]
+        self.setprefix() # Reset the path names of all parameters
+        new = [p.path for p in changeset]
+        mapping = dict(zip(old,new))
+
+        # Perform the substitution into all of the expressions
+        exprs = self.computed
+        for p in exprs:
+            p.expression = expression.substitute(p.expression, mapping)
+
+    def gather_context(self):
+        """
+        Gather all additional symbols that can be used in expressions.
+        
+        For example, if reflectometry provides a volume fraction
+        function volfrac(rho1,rho2,frac) to compute densities, then
+        this function can be added as a context dictionary to the
+        reflectometry parameter set.  Note that there is no guarantee
+        which function will be used if the same function exists in
+        two separate contexts.
+        """
+        context = {} # Create a new dictionary
+        context.update(self.context)
+        for p in self:
+            if hasattr(p,'gather_context'):
+                context.update(p.gather_context())
+        return context
+
+
+def test():
+    # Check single parameter
+    a = Parameter('a')
+    assert a.name == 'a'
+    a.value = 5
+    assert a.value == 5
+    # Check the setters
+    a.set(7)
+    assert a.value == 7 and a.status == 'fixed' and a.isfixed()
+    a.set([3,5])
+    assert a.value == 7 and a.range[0] == 3 and a.range[1]==5 and a.isfitted()
+    a.set('3*M.b')
+    assert a.iscomputed() and a.expression == '3*M.b'
+
+    # Check limits
+    a.limits = (0,inf)
+    try: 
+        a.range = (-1,1)
+        raise Exception,"Failed to check range in limits"
+    except ValueError: pass # Correct failure
+    
+    # Check that we can't change name directly
+    try:
+        a.name = 'Q'
+        raise Exception,"Failed to protect name"
+    except AttributeError: pass # Correct failure
+
+    assert str(a) == 'a'  # Before setpath, just print name
+    a.setprefix('M.')
+    assert a.path == 'M.a'
+    assert str(a) == 'M.a'  # After setpath, print path
+    assert a.units == ''
+    a.units = 'kg'
+    assert a.units == 'kg'
+    assert repr(a) == "Parameter('a')"
+
+    # Check parameter set
+    b,c = Parameter('b'),Parameter('c')
+    M1 = ParameterSet('M1',[a,b,c])
+    assert M1[0] is a
+    assert M1['a'] is a
+    a.set(5) # one value
+    b.set([3,5])
+    c.set('3*M1.a')
+    assert M1.computed == [c]
+    assert M1.fitted == [b]
+    assert M1.fixed == [a]
+    d = Parameter('d')
+    M1.insert(0,d)
+    assert M1.fixed == [d,a]
+    e = Parameter('e')
+    M1.append(e)
+    assert M1.fixed == [d,a,e]
+
+    a2,b2,c2 = [Parameter(s) for s in ('a','b','c')]
+    a2.set(15)
+    M2 = ParameterSet('M2',[a2,b2,c2])
+    # Adjust parameter in set
+    b2.set([3,5])
+    assert M2.fitted == [b2]
+    
+    # Hierarchical parameter sets
+    r = Parameter('r')
+    root = ParameterSet('root',[M1,r,M2])
+    assert root.fixed == [d,a,e,r,a2,c2]
+    assert root.fitted == [b,b2]
+    assert root.computed == [c]
+    root.setprefix()
+    assert a2.path == "M2.a"
+    # Rename individual parameter
+    root.rename(a,'a1')
+    assert a.path == "M1.a1"
+    assert c.expression == "3*M1.a1"
+    # Rename parameter set
+    root.rename(M1,'m1')
+    assert c.path == "m1.c"
+    assert c.expression == "3*m1.a1"
+
+    # Integration test: parameter and expression working together.
+    fn = expression.build_eval(root.flatten(), root.gather_context())
+    #import dis; dis.dis(fn)
+    fn()
+    assert c.value == 3*a.value
+
+    # Test context
+    M2.context['plus2'] = lambda x: x+2
+    c2.set('plus2(M2.a)')
+    assert set(root.computed) == set([c,c2])
+    fn = expression.build_eval(root.flatten(), root.gather_context())
+    #print dis.dis(fn)
+    fn()
+    assert c2.value == a2.value+2
+
+    # Multilevel hierarchy
+    # Forming M3.a.x, M3.a.y, M3.b with M3.a.y = 2*M3.b+M3.a.x
+    x = Parameter('x'); x.set(15)
+    y = Parameter('y'); y.set('2*M3.b+M3.a.x')
+    b = Parameter('b'); b.set(10)
+    a = ParameterSet('a',[x,y])
+    M3 = ParameterSet('M3',[a,b])
+    root = ParameterSet('root',[M3])
+    root.setprefix()
+    fn = expression.build_eval(root.flatten(), root.gather_context())
+    #import dis; dis.dis(fn)
+    fn()
+    #print "y-value:",y.value
+    assert y.value == 35
+    
+
+
+
+if __name__ == "__main__": test()
Index: park-1.2.1/park/parkde.py
===================================================================
--- park-1.2.1/park/parkde.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/parkde.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,56 @@
+"""
+"""
+from __future__ import division
+import thread
+import numpy
+
+import diffev
+
+import fit,fitresult
+
+__all__ = ['DiffEv']
+
+class DiffEv(fit.Fitter):
+    """
+    Differential evolution optimizer
+
+    This implements `park.fit.Fitter`.
+    """
+
+    maxiter=1000
+    """Maximum number of iterations"""
+    pop_scale=4
+    """Number of active points per dimension"""
+    crossover_rate=0.9
+    """Amount of mixing in population"""
+    Fscale=0.5
+    """Step size along difference vector"""
+    tolerance=1e-5
+    """Fit tolerance"""
+
+    def progress(self, k, n):
+        self.handler.progress(k,n)
+    def improvement(self, x, fx, ncalls=-1):
+        self.handler.result.update(x,fx,-1)
+        self.handler.improvement()
+
+    def _call(self, objective, x0, bounds):
+        x,fx,calls = diffev.diffev(objective, bounds,
+                                   maxiter=self.maxiter,
+                                   pop_scale=self.pop_scale,
+                                   crossover_rate=self.crossover_rate,
+                                   Fscale=self.Fscale,
+                                   tolerance=self.tolerance,
+                                   x0=x0, monitor=self)
+
+        # Post the result
+        self.handler.result.update(x, fx, calls)
+        self.handler.result.calc_cov(objective)
+        self.handler.done = True
+        self.handler.finalize()
+
+    def _fit(self, objective, x0, bounds):
+        self._threaded(self._call, objective, x0, bounds)
+
+if __name__ == "__main__":
+    fit.demo2(fitter=DiffEv())
Index: park-1.2.1/park/parksnob.py
===================================================================
--- park-1.2.1/park/parksnob.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/parksnob.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,51 @@
+"""
+"""
+from __future__ import division
+import thread
+import numpy
+
+from snobfit.snobfit import snobfit
+
+import fit,fitresult
+
+__all__ = ['Snobfit']
+
+class Snobfit(fit.Fitter):
+    """
+    Response surface optimizer
+
+    This implements `park.fit.Fitter`.
+    """
+
+    p=0.5
+    """Locality: 1.0 for completely global, 0. for local"""
+    dn=5
+    """Number of surplus points in quadratic fit; increase for noisy functions"""
+    maxiter=1000
+    """Maximum number of iterations"""
+    maxfun=1000
+    """Maximum number of function evaluations"""
+    nstop=50
+    """Number of times no improvement is tolerated"""
+
+    def _monitor(self, k, x, fx, improved):
+        self.handler.progress(k,self.maxiter)
+        if improved:
+            self.handler.result.update(x,fx,-1)
+            self.handler.improvement()
+
+    def _call_snobfit(self, objective, x0, bounds):
+        x,fx,calls = snobfit(objective, x0, bounds, fglob=0,
+                             callback=self._monitor)
+
+        # Post the result
+        self.handler.result.update(x, fx, calls)
+        self.handler.result.calc_cov(objective)
+        self.handler.done = True
+        self.handler.finalize()
+
+    def _fit(self, objective, x0, bounds):
+        self._threaded(self._call_snobfit, objective, x0, bounds)
+
+if __name__ == "__main__":
+    fit.demo2(fitter=Snobfit())
Index: park-1.2.1/park/peaks.py
===================================================================
--- park-1.2.1/park/peaks.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/peaks.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,252 @@
+# This program is public domain
+"""  
+A park model implementing a multipeak fitter.
+
+*** WARNING *** this example was used to inform the design process,
+and has not yet been updated to correspond to the current implementation.
+Do not use this as a tutorial.
+
+This is an example model showing how to put together a multi-part
+fit objective function.
+
+:group main: Peaks
+:group peaks: Gaussian, Lorentzian, Voigt
+:group backgrounds: Constant, Quadratic, Linear
+"""
+import numpy
+import park
+
+class Peaks(park.Model): 
+    """Peak fitter"""
+    # Increment the minor version when adding new functions
+    # Increment the major version if refactoring the class
+    __version__ = "0.9"
+    peak_types = {}
+
+    def __init__(self, datafile=None):
+        """
+        Create a new multipeak model.  If datafile is present it will
+        be loaded directly.  If not, then it can be set later with
+        self.data.load(datafile), or replaced with a specialized
+        data loader for the particular instrument.
+        """
+        # The initializer creates self.parameterset which will hold
+        # the set of fitting parameters for all models.  It will also
+        # define self.data = park.data.Data1D() which computes the
+        # residuals given an eval function.
+        park.Model.__init__(self,filename=datafile)
+        self.parameterset = park.ParameterSet()
+        self.peaks = {}
+        self.peaknum = 0
+
+    def eval(self, x):
+        """
+        Returns the raw theory value at x, unconvoluted and unweighted.
+        
+        The residuals will be calculated by park.Model.residuals by
+        calling 
+        """
+        x = self.data[0]
+        y = numpy.empty(x.shape)
+        for name,peak in self.peaks.itervalues():
+            y += peak(x)
+        return y
+        
+    def add_peak(self, type, name=None, **kw):
+        """ 
+        Add a peak to the model.
+        
+        The name of the peak is used to distinguish it from other peaks 
+        in the model, and appears as part of the parameter name in
+        constraint expressions.  For example, if name is 'P3' and this
+        is part of model 'M1' then the sigma parameter of a gaussian
+        peak would be 'M1.P3.sigma'.
+
+        The peak can be specified either by peak type and initial arguments
+        or by the peak itself, which is a function of x with a parameters 
+        attribute.
+        
+        Available peak types include:
+            gaussian(scale=1, center=0, sigma=1) 
+                = scale exp ( -0.5 (x-center)**2 / sigma**2 )
+            lorentzian(scale=1, center=0, gamma=1) 
+                = scale/pi gamma / ((x-center)**2 + gamma**2))
+            voigt(scale=1, center=0, sigma=1, gamma=1)
+                = scale (lorentzian(mu,gamma) * gaussian(mu,sigma))(x-center)
+                where * represents convolution
+        
+        Available background functions include:
+            constant(C=0)
+                = C
+            slope(C=0, B=1)
+                = B x + C
+            quadratic(A=1, B=0, C=0)
+                = A x**2 + B x + C
+
+        Additional peak types can be registered.  The complete dictionary
+        of available types is in Peaks.peak_types.
+        """
+        self.peaknum
+        if name == None: name = 'P'+str(self.peaknum)
+        if isinstance(type,basestring):
+            peak = self.peak_types[type](**kw)
+        else:
+            type = peak
+        self.peaks[name] = peak
+        self.parameters[name] = peak.parameters
+        setattr(self,name,peak)
+
+    def remove_peak(self, name):
+        """
+        Remove the named peak from the model.
+        """
+        del self.peaks[name]
+        del self.parameterset[name]
+
+    @classmethod
+    def register_peak_generator(cls, type, factory):
+        """
+        Register a new peak generator.  This will return a callable
+        object with a parameter set that can be used in the peak
+        fitting class.
+
+        The individual peak functions a number of attributes:
+        
+        __name__
+            the name to be displayed in the list of peaks; the
+            default is __class__.__name__
+        __factory__
+            the name of a factory function in the python
+            namespace which returns a new peak.  These parameters
+            will be set from the stored parameter list when the
+            peak model is reloaded.  This can be omitted if it
+            is simply the module and name of the class.
+        __doc__
+            the tool tip to be displayed with the peak
+        __call__(x)
+            evaluate the peak function at x
+        """
+        cls.peak_types[type]=factory
+
+
+
+# ================== Peaks and backgrounds ===================
+class Gaussian(object):
+    """Gaussian peak: scale exp ( -0.5 (x-center)**2 / sigma**2 )"""
+    def __init__(self, name, scale=1., center=0, sigma=1):
+        park.define_parameters(name, ['scale','center','sigma'],
+                               scale=scale,center=center,sigma=sigma)
+        self.parameters['scale'].deriv = self.dscale
+        self.parameters['center'].deriv = self.dcenter
+        self.parameters['sigma'].deriv = self.dsigma
+    def __call__(self, x):
+        return self.scale * numpy.exp(-((x-self.center)/self.sigma)**2)
+    def dscale(self, x):
+        return self(x)/self.scale
+    def dcenter(self, x):
+        return 2*(x-self.center)/self.sigma*self(x)
+    def dsigma(self, x):
+        return 2*(x-self.center)/self.sigma**3*self(x)
+
+class Lorentzian(object):
+    """Lorentzian peak (HWHM): scale/pi  gamma/((x-center)**2 + gamma**2)"""
+    def __init__(self, name, scale=1., center=0, gamma=1):
+        park.define_parameters(name, ['scale','center','gamma'],
+                               scale=scale,center=center,gamma=gamma)
+        self.parameters['scale'].deriv = self.dscale
+        self.parameters['center'].deriv = self.dcenter
+        self.parameters['gamma'].deriv = self.dgamma
+    def __call__(self, x):
+        return self.scale/numpy.pi * self.gamma/((x-self.center)**2 + self.gamma**2)
+    def dscale(self, x):
+        return self(x)/self.scale
+    def dcenter(self, x):
+        return 2*(x-self.center)/((x-self.center)**2+self.gamma**2)*self(x)
+    def dgamma(self, x):
+        return self(x)*(1/self.gamma - self.gamma/((x-self.center)**2 + self.gamma**2))
+
+class Voigt(object):
+    """Voigt peak (HWHM,sigma): A [G(sigma) * L(gamma)](x-center)"""
+    def __init__(self, name, scale=1, center=0, sigma=1, gamma=1):
+        park.define_parameters(name, ['scale','center','sigma','gamma'],
+                               scale=scale, center=center, sigma=sigma,
+                               gamma=gamma)
+    def __call__(self, x):
+        return self.scale*voigt(x-self.center, 
+                                sigma=self.sigma, gamma=self.gamma)
+
+class Quadratic(object):
+    """Quadratic background: A x**2 + B x + C"""
+    name = "Background: quadratic"
+    def __init__(self, name, A=1, B=0, C=0):
+        park.define_parameters(name, ['A', 'B', 'C'], A=A, B=B, C=C)
+        self.parameters['A'].deriv = self.dA
+        self.parameters['B'].deriv = self.dB
+        self.parameters['C'].deriv = self.dC
+    def __call__(self, x):
+        return numpy.polyval([self.a.value, self.b.value, self.c.value],x)
+    def dA(self, x): return x**2
+    def dB(self, x): return x
+    def dC(self, x): return numpy.ones(x.shape)
+
+class Linear(object):
+    """Linear background: B x**2 + C"""
+    name = "Background: linear"
+    def __init__(self, name, B=0, C=0):
+        park.define_parameters(name, ['B', 'C'], B=B, C=C)
+        self.parameters['B'].deriv = self.dB
+        self.parameters['C'].deriv = self.dC
+    def __call__(self, x):
+        return numpy.polyval([self.a.value, self.b.value, self.c.value],x)
+    def dB(self, x): return x
+    def dC(self, x): return numpy.ones(x.shape)
+
+class Constant(object):
+    """Constant background: C"""
+    name = "Background: constant"
+    def __init__(self, name, C=0):
+        park.define_parameters(name, ['C'], C=C)
+        self.parameters['C'].deriv = self.dC
+    def __call__(self, x):
+        return numpy.polyval([self.a.value, self.b.value, self.c.value],x)
+    def dC(self, x): return numpy.ones(x.shape)
+
+
+def voigt(x, sigma, gamma):
+    """
+    Return the voigt function, which is the convolution of a Lorentz
+    function with a Gaussian.
+    
+    :Parameters:
+     gamma : real
+      The half-width half-maximum of the Lorentzian
+     sigma : real
+      The 1-sigma width of the Gaussian, which is one standard deviation.
+
+    Ref: W.I.F. David, J. Appl. Cryst. (1986). 19, 63-64
+    
+    Note: adjusted to use stddev and HWHM rather than FWHM parameters
+    """
+    # wofz function = w(z) = Fad[d][e][y]eva function = exp(-z**2)erfc(-iz)
+    from scipy.special import wofz
+    z = (x+1j*gamma)/(sigma*sqrt(2))
+    V = wofz(z)/(sqrt(2*pi)*sigma)
+    return V.real
+
+
+def init():
+    """
+    Register peak types with the Peaks model.
+    """
+    for cls in [Gaussian, Lorentzian, Voigt, Quadratic, Linear, Constant]:
+        Peaks.register_peak_generator(cls.__name__,cls)
+init()
+
+# ========================================
+
+def test():
+    x = numpy.linspace(-5,5,10)
+    
+    pass
+
+if __name__ == "__main__": test()
Index: park-1.2.1/park/pmap.py
===================================================================
--- park-1.2.1/park/pmap.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/pmap.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,227 @@
+"""
+Parallel map-reduce implementation using threads.
+"""
+
+import traceback
+import thread
+
+class Collector(object):
+    """
+    Abstract interface to map-reduce accumulator function.
+    """
+    def __call__(self, part):
+        """Receive next part, storing it in the accumulated result"""
+    def finalize(self):
+        """Called when all parts have been accumulated"""
+    def error(self, part, msg):
+        """
+        Exception seen on executing map or reduce.  The collector
+        can adjust the accumulated result appropriately to reflect
+        the error.
+        """
+
+class Mapper(object):
+    """
+    Abstract interface to map-reduce mapper function.
+    """
+    def __call__(self, value):
+        """Evaluate part"""
+    def abort(self):
+        """Stop the mapper"""
+
+def pmap(mapper, inputs):
+    """
+    Apply function mapper to all inputs.
+    
+    This is the serial version of a parallel iterator, yielding the next 
+    sequence value as soon as it is available.  There is no guarantee 
+    that the order of the inputs will be preserved in the parallel
+    version, so don't depend on it!
+    """
+    for item in inputs:
+        yield mapper(item)
+
+def preduce(collector, outputs):
+    """
+    Collect all outputs, calling collector(item) for each item in the sequence.
+    """
+    for item in outputs:
+        collector(item)
+    return collector
+
+def _pmapreduce_thread(fn, collector, inputs):
+    try:
+        preduce(collector, pmap(fn,inputs))
+        collector.finalize()
+    except KeyboardInterrupt:
+        fn.abort()
+        thread.interrupt_main()
+    #except:
+    #    raise
+    #    msg = traceback.format_exc()
+    #    collector.error(msg)
+
+def _pmapreduce_profile(fn, collector, inputs):
+    import cProfile, pstats, os
+    def mapr():
+        _pmapreduce_thread(fn, collector, inputs)
+    cProfile.runctx('mapr()', dict(mapr=mapr), {}, 'mapr.out')
+    stats = pstats.Stats('mapr.out')
+    #stats.sort_stats('time')
+    stats.sort_stats('calls')
+    stats.print_stats()
+    os.unlink('mapr.out')
+
+profile_mapper = False
+"""True if the mapper cost should be profiled."""
+
+def pmapreduce(mapper, collector, inputs):
+    """
+    Apply function mapper to inputs, accumulating the results in collector.
+    
+    Collector is a function which accepts the result of mapper(item) for 
+    each item of inputs.  There is no guarantee that the outputs will be
+    received in order.
+    
+    The map is executed in a separate thread so the function returns
+    to the caller immediately.
+    """
+    global profile_mapper
+    fn = _pmapreduce_profile if profile_mapper else _pmapreduce_thread        
+    thread.start_new_thread(fn,(mapper,collector, inputs))
+
+
+def main():
+    import time,numpy
+    class TestCollector(object):
+        def __init__(self):
+            self.done = False
+        def __call__(self, part):
+            print "collecting",part,'for',id(self)
+        def finalize(self):
+            self.done = True
+            print "finalizing"
+        def abort(self):
+            self.done = True
+            print "aborting"
+        def error(self,msg):
+            print "error",msg
+    
+    class TestMapper(object):
+        def __call__(self, x): 
+            print "mapping",x,'for',id(self)
+            if x == 8: raise Exception('x is 8')
+            time.sleep(4*numpy.random.rand())
+            return x
+
+    collector1 = TestCollector()
+    collector2 = TestCollector()
+    pmapreduce(TestMapper(), collector1, [1,2,3,4,5])
+    pmapreduce(TestMapper(), collector2, [1,2,3,8])
+    while not collector1.done and not collector2.done:
+        time.sleep(1)
+if __name__ == "__main__": main()
+
+    
+_ = '''
+# The choice of job to do next is complicated.
+# 1. Strongly prefer a job of the same type as is already running
+# on the node.  If this job requires significant resources (e.g.,
+# a large file transfer) increase that preference.
+# 2. Strongly prefer sending a user's own job to their own machine.
+# That way at least they can make progress even if the cluster is busy.
+# 3. Try to start each job as soon as possible.  That way if there are
+# errors, then the user gets feedback early in the process.
+# 4. Try to balance the load across users.  Rather than first come
+# first serve, jobs use round robin amongst users.
+# 5. Prefer high priority jobs.
+
+
+def map(apply,collect,items,priority=1):
+    mapper = MapJob(apply, items, collect, priority)
+    
+class MapJob(object):
+    """
+    Keep track of which jobs have been submitted and which are complete
+    """
+    def __init__(self, workfn, worklist, manager, priority):
+        self.workfn = workfn
+        self.worklist = worklist
+        self.manager = manager
+        self.priority = priority
+        self._priority_edge = 0
+    def next_work(self):
+        
+
+class MapServer(object):
+    """
+    Keep track of work units.
+    """
+    def __init__(self):
+        self.workingset = {}
+        
+    def add_work(self, workfn, worklist, manager, priority):
+        """
+        Add a new work list to distributed to worker objects.  The manager
+        gathers the results of the work.  Work is assigned from the queue
+        based on priority.
+        """
+        job = MapJob(workfn, worklist, manager, priority)
+
+        # add work to the queue in priority order
+        for i,job in enumerate(self.jobs):
+            if job.priority < priority: break
+        self.jobs.insert(i,job)
+            
+        # Create an entry in a persistent store for each job to
+        # capture completed work units and to recover from server
+        # reboot.
+
+        # Assign _priority_edge to cumsum(priority)/total_priority.
+        # This allows us to select the next job according to priority
+        # with a random number generator.
+        # NOTE: scalability is a bit of a concern --- the lookup
+        # operation is linear in the number of active jobs.  This
+        # can be mitigated by limiting the number of active jobs.
+        total_priority = 0.
+        for job in self.jobs: total_priority += job.priority
+        edge = 0.
+        for job in self.jobs:
+            edge += job.priority/total_priority
+            self.job._priority_edge = edge
+        
+    
+    def register(self, pool=None):
+        """
+        Called by a worker when they are registering for work.
+        
+        Returns the next piece of work.
+        """
+        P = numpy.random.rand()
+        for job in self.jobs:
+            if P < j._priority_edge:
+                return job.new_work()
+
+        return NoWork
+
+    def update(self, jobid, result):
+        """
+        Called by a worker when the work unit is complete.
+        
+        Returns the next piece of work.
+        """
+        current_job = self.lookup(jobid)
+        current_job.reduce(result)
+        if numpy.random.rand() < current_job.switch_probability:
+            return current_job.next_work()
+        
+        P = numpy.random.rand()
+        for job in self.jobs:
+            if P < job._priority_edge:
+                if job == current_job:
+                    return curent_job.next_work()
+                else:
+                    return job.new_work()
+        
+        return NoWork
+'''
Index: park-1.2.1/park/rangemap.py
===================================================================
--- park-1.2.1/park/rangemap.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/rangemap.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,359 @@
+"""
+Defines transformations between the a fit space and the parameter space.
+
+The parameter space maps a set of possibly bounded dimensions into a
+real-valued fitness::
+
+    f: [a,b]**n |-> R
+
+The fit space is constrained to map values in the unit box::
+
+    f': [0,1]**n |-> R 
+
+Here we can define f' using the zero-one mapper::
+
+    f' = ZeroOneMapper(f,a,b)
+    f'(v) = f(x)
+
+An asymptote function maps x to [0,1], preserving at least 12 digits of
+precision.
+
+Using this mapping, the optimizer can operate in a well known space
+independent of parameter precision.  This allows the use of reasonable
+constants for items such as a step size and initial value.
+
+Note: we do not yet support fitness functions with analytic derivatives.
+Given df/dp and mapping function f'(v) = f(M(v)) then df'/dv = df/dM dM/dv.
+So the derivatives need to be multiplied by the derivative of the
+parameter mapper.  Not difficult computationally, but the definition of
+the fitness function does not currently support this organization.
+"""
+
+
+class FitCobyla(LocalFit):
+    """
+    Local minimizer using COBYLA, Constrained Optimization by Linear
+    Approximation.
+    
+    COBYLA minimizes an objective function F(X) subject to M inequality 
+    constraints on X, where X is a vector of variables that has N components. 
+    The algorithm employs linear approximations to the objective and 
+    constraint functions, the approximations being formed by linear 
+    interpolation at N+1 points in the space of the variables.  These 
+    interpolation points can be regarded as vertices of a simplex. The
+    parameter rho controls the size of the simplex and it is reduced
+    automatically from radius to xtol. For each RHO the subroutine tries
+    to achieve a good vector of variables for the current size, and then
+    RHO is reduced until the value xtol is reached. Therefore radius and
+    xtol should be set to reasonable initial changes to and the required   
+    accuracy in the variables respectively, but this accuracy should be
+    viewed as a subject for experimentation because it is not guaranteed.
+
+    Basic bounds constraints are built into the interface.  These bounds
+    map the space into the n-D unit cube, so values values for radius
+    and xtol should be chosen accordingly.  Regardless of xtol, the
+    fit is stopped after maxiter function evaluations.
+    """
+    radius = 0.05
+    """Size of the initial simplex; this is a portion between 0 and 1"""
+    xtol = 1e-4
+    """Stop when simplex vertices are within xtol of each other"""
+    maxiter = 1000
+    def __call__(self):
+        from scipy.optimize.cobyla import fmin_cobyla as cobyla
+
+class ArctanAsymptote(object):
+    """
+    Arctan asymptote function for mapping (-inf,inf) to [0,1].
+
+    With all attributes constant, either the class or an instance can be used.
+
+    An asymptote function is monotonically increasing and finite everywhere.
+    Assuming the function is a at -inf and b at +inf, scale is 1/(b-a) and
+    offset is a.  Together with an inverse function, this allows us to
+    transform between an infinite range and the range [0,1].  Semidefinite
+    map to (-inf,0] and [0,inf).
+
+    ArctanAsymptote is approximately linear in [-0.5,0.5], and is useful 
+    out to the range [-1e15,1e15],  though with much reduced precision 
+    for larger values.
+    """
+    scale = 1/numpy.pi
+    offset = -numpy.pi/2
+    forward = numpy.arctan
+    inverse = numpy.tan
+
+def fp_forward(x):
+    """
+    Linearize floating point values using an exponential scale.
+    
+    Convert sign*m*2^e to sign*(e+1023+m), yielding a value in [-2047,2047]
+    """
+    (m,e) = numpy.frexp(x)
+    s = numpy.sign(x)
+    v = (e+1023+m*s)*s
+    return v
+def fp_inverse(v):
+    """
+    Restore floating point value from linear form on an exponential scale.
+    
+    Convert sign*(e+1023+m) to sign*m*2^e, yielding a value in [-1e308,1e308]
+    """
+    s = sign(v)
+    v *= s
+    m = floor(v)
+    e = v-m
+    x = ldexp(s*m,e)
+    return x
+fp_min, fp_max = -(1024+1023), (1024+1023)
+class Asymptote:
+    """
+    Logarithmic asymptote function.
+    
+    This is an asymptote function which follows the floating point 
+    representation, using the following mapping::
+    
+        [-1e308,-1e-308] -> [0,0.5)
+        0 -> 0.5
+        [1e-308,1e308] -> (0.5+eps,1]
+
+    With 11 bits of the 53 bits of available precision for the exponent, 
+    this leaves a tolerance of 1e-12 on the fitting parameters, which is 
+    easily good enough for any real unbounded fitting problem.  The user 
+    can increased the precision on the bounded parameters by using tighter 
+    bounds.
+    """
+    scale = 1/(fp_max-fp_min)
+    offset = fp_min
+    forward = fp_forward
+    inverse = fp_inverse
+
+class ZeroOneMapper(object):
+    """
+    Map function range into [0,1]**n.
+
+    Given f: [a,b] -> R
+    Defines f': [0,1] -> R
+    
+    The method encode(x) returns a value in [0,1]
+    The method decode(v) returns a value in [a,b]
+    The method __call__(v) returns f(decode(v))
+    
+    For indefinite and semidefinte ranges, an asymptote transformation
+    is needed.  Normally this is `park.rangemap.Asymptote` which supports
+    the full floating point range of values, but is limited to about
+    12 digits of precision.  The arctan function can also be used
+    (see `park.rangemap.ArctanAsymptote`).
+    
+    Range is determined by the bounds low and high.  Each dimension may be 
+    unbounded, semi-definite or bounded. Bounded functions use a linear 
+    mapping between low and high.  Unbounded functions use an asymptote 
+    function, which should be -1 at -inf, 0 at 0 and 1 at +inf, and 
+    approximate the identity function between [-0.5,0.5].
+    
+    This can be used to turn any unconstrained optimizer into a [0,1]
+    bounded optimizer by sending infeasible points to infinity.
+
+    Note: Newton-style optimizers will not work in this regime, but 
+    instead require a hint about the direction of the unconstrained
+    region.
+    """
+    def __init__(self, base_function, low, high, asymptote=Asymptote,
+                 range_check=False):
+        """
+        Function range mapper.  See `park.rangemap.ZeroOneMapper` for details.
+
+        base_function is the function being wrapped.
+        
+        low[k],high[k] is the range of values for each fitting parameter k.
+        
+        range_check is True if calls to the mapper may include out of range
+        values.  In this case, f(v) returns inf for out of range values and
+        encode/decode raise exceptions.
+
+        asymptote is the asymptote function to use when linearizing
+        indefinite ranges.
+        """
+
+        unbounded_low = numpy.isinf(low)
+        unbounded_high = numpy.isinf(high)
+        unbounded = unbounded_low|unbounded_high
+    
+        # Determine linear scale and offset for v = (x-offset)/scale
+        # For semi-infinite ranges, use the known bound as the offset.  This 
+        # transforms [a,inf) to [0.5,1] and (-inf,b] to [0,0.5].
+        # Infinite ranges should use offset 0.
+        scale = (high-low)
+        offset = -low
+        scale[unbounded_low | unbounded_high] = 0.
+        offset[unbounded_low] = high[unbounded_low]
+        offset[unbounded_high] = low[unbounded_high]
+        offset[unbounded_low & unbounded_high] = 0. # low&high must be last
+
+        # Determine transformed scale and offset
+        # If the value is unbounded, then use asymptote directly.
+        # If bounded below, then transform [0.5,1] to [0,1] using (v-0.5)/0.5
+        # If bounded above, then transform [0,0.5] to [0,1] using v/0.5
+        a_scale = numpy.ones(scale.shape)*asymptote.scale
+        a_scale[unbounded_low ^ unbounded_high] *= 0.5
+        a_offset = numpy.zeros(offset.shape)+asymptote.offset
+        a_offset[unbounded_high & ~unbounded_low] += 0.5*asymptote.scale
+        a_forward = asymptote.forward
+        a_inverse = asymptote.inverse
+
+        # Preselect the relevant elements
+        a_scale = a_scale[unbounded]
+        a_offset = a_offset[unbounded]
+
+        # Remember the transformation parameters
+        self.low,self.high = low,high
+        self.base_function = base_function
+        self.scale,self.offset = scale,offset
+        self.a_scale,self.a_offset = a_scale,a_offset
+        self.a_forward,self.a_inverse = a_forward,a_inverse
+        
+        # To minimize function call overhead, use only the linear bounds
+        # tranform if the problem is completely bounded.
+        if numpy.any(unbounded):
+            self.__call__ = self._indefinite
+        else:
+            self.__call__ = self._definite    
+
+    def _indefinite(self, v):
+        """Use this function if some bounds are indefinite"""
+        x = v.copy()
+        x[unbounded] = self.a_inverse(x[unbounded])*self.a_scale + self.a_offset
+        x = x*self.scale+self.offset
+        return self.base_function(x)
+    def _definite(self, v):
+        """Use this function if all bounds are definite"""
+        return self.base_function(v*self.scale+self.offset)
+    def _indefinite_checked(self, v):
+        """Use this function if some bounds are indefinite and ranges need
+        to be checked"""
+        if numpy.any(v<0|v>1): return numpy.inf
+        x = v.copy()
+        x[unbounded] = self.a_inverse(x[unbounded])*self.a_scale + self.a_offset
+        x = x*self.scale+self.offset
+        return self.base_function(x)
+    def _definite_checked(self, v):
+        """Use this function if all bounds are definite, but ranges need
+        to be checked"""
+        if numpy.any(v<0|v>1): return numpy.inf
+        return self.base_function(v*self.scale+self.offset)
+
+    def decode(v):
+        """Transform from range [0,1] to [a,b]."""
+        if numpy.any(v<0|v>1): 
+            raise RuntimeError("value out of range [0,1]")
+        x = v.copy()
+        x[unbounded] = self.a_forward(x[unbounded])*a_scale + a_offset
+        x = x*scale+offset
+        return x
+    def encode(x):
+        """Transform from range [a,b] to [0,1]."""
+        if numpy.any(x<self.low|x>self.high): 
+            raise RuntimeError("value out of range [a,b]")
+        v = (x-offset)/scale
+        v[unbounded] = (self.a_inverse(v[unbounded])-a_offset)/a_scale
+        return v
+
+# Functional implementation: return transformation functions.
+# This is faster but less pythonic
+def zero_one_mapper(base_function, low, high, asymptote=ArctanAsymptote):
+    """
+    Map function range into [0,1]**n.   
+    
+    Returns a pair of functions f and inv.  The function f takes an
+    encoded value in 
+     encode, which takes x and returns a value 
+    in [0,1] and decode, which takes a value in [0,1] and returns x.
+    
+    Range is determined by the bounds low and high.  Each dimension may be 
+    unbounded, semi-definite or bounded. Bounded functions use a linear 
+    mapping between low and high.  Unbounded functions use an asymptote 
+    function, which should be -1 at -inf, 0 at 0 and 1 at +inf, and 
+    approximate the identity function between [-0.5,0.5].
+    
+    This can be used to turn any unconstrained optimizer into a [0,1]
+    bounded optimizer by sending infeasible points to infinity.
+
+    Note: Newton-style optimizers will not work in this regime, but 
+    instead require a hint about the direction of the unconstrained
+    region.
+    """
+
+    unbounded_low = numpy.isinf(low)
+    unbounded_high = numpy.isinf(high)
+    unbounded = unbounded_low|unbounded_high
+    
+    # Determine linear scale and offset for v = (x-offset)/scale
+    # For semi-infinite ranges, use the known bound as the offset.  This 
+    # transforms [a,inf) to [0.5,1] and (-inf,b] to [0,0.5].
+    # Infinite ranges should use offset 0.
+    scale = (high-low)
+    offset = -low
+    scale[unbounded_low | unbounded_high] = 0.
+    offset[unbounded_low] = high[unbounded_low]
+    offset[unbounded_high] = low[unbounded_high]
+    offset[unbounded_low & unbounded_high] = 0. # low&high must be last
+
+    # Determine transformed scale and offset
+    # If the value is unbounded, then use asymptote directly.
+    # If bounded below, then transform [0.5,1] to [0,1] using (v-0.5)/0.5
+    # If bounded above, then transform [0,0.5] to [0,1] using v/0.5
+    a_scale = numpy.ones(scale.shape)*asymptote.scale
+    a_scale[unbounded_low ^ unbounded_high] *= 0.5
+    a_offset = numpy.zeros(offset.shape)+asymptote.offset
+    a_offset[unbounded_high & ~unbounded_low] += 0.5*asymptote.scale
+
+    # Preselect the relevant elements
+    a_scale = a_scale[unbounded]
+    a_offset = a_offset[unbounded]
+    
+
+    if numpy.any(unbounded):
+        def function(v):
+            x = v.copy()
+            x[unbounded] = asymptote.inverse(x[unbounded])*a_scale + a_offset
+            x = x*scale+offset
+            return base_function(x)
+        def decode(v):
+            x = v.copy()
+            x[unbounded] = asymptote.inverse(x[unbounded])*a_scale + a_offset
+            x = x*scale+offset
+            return x
+        def encode(x):
+            v = (x-offset)/scale
+            v[unbounded] = (asymptote.function(v[unbounded])-a_offset)/a_scale
+            return v
+    else:
+        def function(v):  return base_function( (x-offset)/scale )
+        def decode(v):  return (v*scale+offset)
+        def function(x):  return base_function( (x-offset)/scale )
+        def decode(v):  return (v*scale+offset)
+
+    return function, encode, decode
+
+
+def bounded(function, low, high):
+    """
+    Evaluate a function with bounds checking.
+    
+    This can be used to turn any unconstrained optimizer into a
+    constrained optimizer by sending infeasible points to infinity.
+
+    Note: Newton-style optimizers will not work in this regime, but 
+    instead require a hint about the direction of the unconstrained
+    region.
+    
+    Note: unused function which may be removed in future versions.
+    """
+    def function_wrapper(x):
+        if  numpy.any((x<low) | (x>high)):
+            return numpy.Inf
+        else:
+            return function(x)
+    function_wrapper.__doc__ = function.__doc__
+    return function_wrapper
+
Index: park-1.2.1/park/serial.py
===================================================================
--- park-1.2.1/park/serial.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/serial.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,457 @@
+# This program is public domain
+"""
+Service to handle serialization and deserialization of fit objects.
+
+Object serialization is useful for long term storage, interlanguage
+communication and network transmission.  In all cases, the process
+involves an initial encode() followed by a later decode().
+
+We need the following properties for serialization/deserialization:
+
+1. human readable so that disaster recovery is possible
+2. readable/writable by other languages and environments
+3. support for numerics: complex, nan, inf, arrays, full precision
+4. version support: load object into newer versions of the
+   program even if the class structure has changed
+5. refactoring support: load object into newer versions of the
+   program even if the classes have been moved or renamed
+
+A complete solution would also support self referential data structures,
+but that is beyond our needs.
+
+Python's builtin serialization, pickle/cPickle, cannot meet these
+needs.  It is python specific, and not friendly to human readers
+or readers from other environments such as IDL which may want to
+load or receive data from a python program.  Pickle inf/nan doesn't
+work on windows --- some of our models may use inf data, and some of
+our results may be nan.  pickle has minimal support for versioning:
+users can write __setstate__ which accepts a dictionary and adjusts
+it accordingly.  Beware though that version must be an instance
+variable rather than a class variable, since class variables are not
+seen by pickle.  If the class is renamed, then pickle can do nothing
+to recover it.
+
+Instead of pickle, we break the problem into steps: structucture and
+encoding.  A pair of functions deconstruct() and reconstruct() work
+directly with the structure.  Deconstruct extracts the state of the 
+python object defined using a limited set of python primitives.  
+Reconstruct takes an extracted state and rebuilds the complete python 
+object.  See documentation on the individual functions for details.
+
+For serial encoding we will use json.  The json format is human 
+readable and easily parsed. json iteself does not define support 
+of Inf/NaN, though some json tools support it using the native
+javascript values of Infinity and Nan. Various xml encodings are 
+also possible, though somewhat more difficult to work with.
+
+Object persistence for long term storage places particular burdens
+on the serialization protocol.  In particular, the class may have
+changed since the instance was serialized.  To aid the process of
+maintaining classes over the long term, the class definition can
+contain the following magic names:
+
+__version__
+    Strict version number of the class.  See isnewer() for
+    details, or distutils.version.StrictVersion.
+__factory__
+    Name of a factory function to return a new instance of
+    the class.  This will be stored as the class name, and
+    should include the complete path so that it can be
+    imported by python.
+__reconstruct__
+    Method which takes a structure tree and rebuilds the object.
+    This is different from __setstate__ in that __setstate__
+    assumes its children have already been reconstructed.  This
+    is the difference between top-down and bottom-up
+    interpretation.  Bottom-up is usually easiear and sufficient,
+    but top-down is required for radical restructuring of the
+    object representation.
+
+
+Example
+=======
+
+The following example shows how to use reconstruct and factory to get
+maximum flexibility when restoring an object::
+
+    from danse.common.serial import isnewer, reconstruct, setstate
+    def data():
+        from data import Data
+        return Data()
+    class Data(object):
+        __version__ = '1.2'
+        __factory__ = 'danse.builder.data'
+        def __reconstruct__(self,instance):
+            '''
+            Reconstruct the state from 
+            '''
+            if isnewer('1.0',instance['version']):
+                raise RuntimeError('pre-1.0 data objects no longer supported')
+            if isnewer('1.1',instance['version']):
+                # Version 1.1 added uncertainty; default it to zero
+                instance['state']['uncertainty'] = 0
+            setstate(self,reconstruct(instance['state']))
+"""
+
+import types
+import sys
+import demjson
+
+def encode(obj):
+    """
+    Convert structure to a string.
+   
+    Basic python types (list, string, dictionary, numbers, boolean, None) 
+    are converted directly to the corresponding string representation.
+    tuples and sets are converted to lists, and str is converted to unicode.
+   
+    Python objects are represented by::
+    
+        { 
+            '.class': 'module.classname', 
+            '.version': 'versionstring',
+            '.state': { object state }
+        }
+    
+    where state comes from the object __getstate__, the object __dict__ or
+    the object __slots__.  See the pickle documentation for details.
+
+    Python functions are represented by::
+    
+        {
+            '.function': 'module.functionname'
+        }
+        
+    """
+    return demjson.encode(deconstruct(obj))
+
+def decode(string):
+    """
+    Convert string to structure, reconstructing classes as needed.  See
+    pickle documentation for details.  This function will fail with a
+    RuntimeError if the version of the class in the string is newer
+    than the version of the class in the python path.
+    """
+    return reconstruct(demjson.decode(string))
+
+
+def deconstruct(obj):
+    """
+    Convert an object hierarchy into python primitives.
+    
+    The primitives used are int, float, str, unicode, bool, None,
+    list, tuple, and dict.
+    
+    Classes are encoded as a dict with keys '.class', '.version', and '.state'.  
+    Version is copied from the attribute __version__ if it exists.
+    
+    Functions are encoded as a dict with key '.function'.
+    
+    Raises RuntimeError if object cannot be deconstructed.  For example,
+    deconstruct on deconstruct will cause problems since '.class' will
+    be in the dictionary of a deconstructed object.
+    """
+    if type(obj) in [int, float, str, unicode, bool] or obj is None:
+        return obj
+    elif type(obj) in [list, tuple, set]:
+        return type(obj)(deconstruct(el) for el in obj)
+    elif type(obj) == dict:
+        # Check for errors
+        for name in ['.class', '.function']:
+            if name in obj:
+                raise RuntimeError("Cannot deconstruct dict containing "+name)
+        return dict((k,deconstruct(v)) for k,v in obj.items())
+    elif type(obj) == types.FunctionType:
+        return {
+            '.function'  : obj.__module__+'.'+obj.__name__
+        }
+    else:
+        cls = _getclass(obj)
+        version = _getversion(obj)
+        return { 
+            '.class'   : _getclass(obj), 
+            '.version' : _getversion(obj), 
+            '.state'   : deconstruct(_getstate(obj))
+        }
+ 
+def reconstruct(tree):
+    """
+    Reconstruct an object hierarchy from a tree of primitives.
+    
+    The tree is generated by deconstruct from python primitives 
+    (list, dict, string, number, boolean, None) with classes 
+    encoded as a particular kind of dict.
+    
+    Unlike pickle, we do not make an exact copy of the original
+    object.  In particular, the serialization format may not
+    distinguish between list and tuples, or str and unicode.  We
+    also have no support for self-referential structures.
+    
+    Raises RuntimeError if could not reconstruct
+    """
+    if type(tree) in [int, float, str, unicode, bool] or tree is None:
+        return tree
+    elif type(tree) in [list, tuple, set]:
+        return type(tree)(reconstruct(el) for el in tree)
+    elif type(tree) == dict:
+        if '.class' in tree:
+            # Chain if program version is newer than stored version (too cold)
+            fn = _lookup_refactor(tree['.class'],tree['.version'])
+            if fn is not None: return fn(tree)
+
+            # Fail if program version is older than stored version (too hot)
+            obj = _createobj(tree['.class'])
+            if isnewer(tree['.version'],_getversion(obj)):
+                raise RuntimeError('Version of %s is out of date'%tree['.class'])
+            # Reconstruct if program version matches stored version (just right)
+            if hasattr(obj, '__reconstruct__'):
+                obj.__reconstruct__(tree['.state'])
+            else:
+                _setstate(obj,reconstruct(tree['.state']))
+            return obj
+        elif '.function' in tree:
+            return _import_symbol(tree['.function'])
+        else:
+            return dict((k,reconstruct(v)) for k,v in tree.items())
+    else:
+        raise RuntimeError('Could not reconstruct '+type(obj).__name__)
+
+def _getversion(obj):
+    version = getattr(obj,'__version__','0.0')
+    try:
+        # Force parsing of version number to check format
+        isnewer(version,'0.0')
+    except ValueError,msg:
+        raise ValueError("%s for class %s"%(msg,obj.__class__.__name__))
+    return version
+
+def _getclass(obj):
+    if hasattr(obj,'__factory__'): return obj.__factory__
+    return obj.__class__.__module__+'.'+obj.__class__.__name__
+
+def _getstate(obj):
+    if hasattr(obj,'__getinitargs__') or hasattr(obj,'__getnewargs__'):
+        # Laziness: we could fetch the initargs and store them, but until
+        # we need to do so, I'm not going to add the complexity.
+        raise RuntimeError('Cannot serialize a class with initialization arguments')
+    elif hasattr(obj,'__getstate__'):
+        state = obj.__getstate__()
+    elif hasattr(obj,'__slots__'):
+        state = dict((s,getattr(obj,s)) for s in obj.__slots__ if hasattr(obj,s))
+    elif hasattr(obj,'__dict__'):
+        state = obj.__dict__
+    else:
+        state = {}
+    return state
+
+def _setstate(obj,kw):
+    if hasattr(obj,'__setstate__'):
+        obj.__setstate__(kw)
+    elif hasattr(obj,'__slots__'):
+        for k,v in kw.items(): setattr(obj,k,v)
+    elif hasattr(obj,'__dict__'):
+        obj.__dict__ = kw
+    else:
+        pass
+    return obj
+
+def _lookup_refactor(cls,ver):
+    return None
+
+class _EmptyClass: pass
+def _import_symbol(path):
+    """
+    Recover symbol from path.
+    """
+    parts = path.split('.')
+    module_name = ".".join(parts[:-1])
+    symbol_name = parts[-1]
+    __import__(module_name)
+    module = sys.modules[module_name]
+    symbol = getattr(module,symbol_name)
+    return symbol
+
+def _createobj(path):
+    """
+    Create an empty object which we can update with __setstate__
+    """
+    factory = _import_symbol(path)
+    if type(factory) is types.FunctionType:
+        # Factory method to return an empty class instance
+        obj = factory()
+    elif type(factory) is types.ClassType:
+        # Old-style class: create an empty class and override its __class__
+        obj = _EmptyClass()
+        obj.__class__ = factory
+    elif type(factory) is types.TypeType:
+        obj = factory.__new__(factory)
+    else:
+        raise RuntimeError('%s should be a function, class or type'%path)
+    return obj
+
+def isnewer(version,target):
+    """
+    Version comparison function.  Returns true if version is at least
+    as new as the target version.
+
+    A version number consists of two or three dot-separated numeric 
+    components, with an optional "pre-release" tag on the end.  The 
+    pre-release tag consists of the letter 'a' or 'b' followed by 
+    a number.  If the numeric components of two version numbers 
+    are equal, then one with a pre-release tag will always
+    be deemed earlier (lesser) than one without.
+
+    The following will be true for version numbers::
+
+      8.2 < 8.19a1 < 8.19 == 8.19.0
+
+    
+    You should follow the rule of incrementing the minor version number
+    if you add attributes to your models, and the major version number
+    if you remove attributes.  Then assuming you are working with
+    e.g., version 2.2, your model loading code will look like::
+    
+        if isnewer(version, Model.__version__):
+            raise IOError('software is older than model')
+        elif isnewer(xml.version, '2.0'):
+            instantiate current model from xml
+        elif isnewer(xml.version, '1.0'):
+            instantiate old model from xml
+            copy old model format to new model format
+        else:
+            raise IOError('pre-1.0 models not supported')
+
+    Based on distutils.version.StrictVersion
+    """
+    from distutils.version import StrictVersion as Version
+    return Version(version) > Version(target)
+
+class _RefactoringRegistry(object):
+    """
+    Directory of renamed classes.
+    
+    """
+    registry = {}
+    
+    @classmethod
+    def register(cls,oldname,newname,asof_version):
+        """
+        As of the target version, references to the old name are no
+        longer valid (e.g., when reconstructing stored objects), and
+        should be resolved by the new name (or None if they should
+        just raise an error.)  The old name can then be reused for
+        new objects or abandoned.
+        """
+        # Insert (asof_version,newname) in the right place in the
+        # list of rename targets for the object.  This list will
+        # be empty unless the name is reused.
+        if name not in cls.registry: cls.registry[name] = []
+        for idx,(version,name) in cls.registry[name]:
+            if isnewer(asof_version, version):
+                cls.registry[name].insert(idx,(asof_version, newname))
+                break
+        else:
+            cls.registry[name].append((asof_version, newname))
+    
+    @classmethod
+    def redirect(cls, oldname, newname, version):
+        if oldname not in cls.registry[oldname]: return None
+        for idx,(target_version,newname) in cls.registry[name]:
+            if isnewer(target_version, version):
+                return target_version
+        # error conditions at this point
+
+def refactor(oldname,newname,asof_version):
+    """
+    Register the renaming of a class.  
+   
+    As code is developed and maintained over time, it is sometimes 
+    beneficial to restructure the source to support new features.  
+    However, the structure and location of particular objects is 
+    encoded in the saved file format.
+
+    When you move a class that may be stored in a model,
+    be sure to put an entry into the registry saying where
+    the model was moved, or None if the model is no longer
+    supported.
+    
+    reconstructor as a function to build a python object from
+    a particular class/version, presumably older than the current
+    version.  This is necessary, e.g., to set default values for new
+    fields or to modify components of the model which are now
+    represented differently.
+    
+    The reconstructor function takes the structure above as 
+    its argument and returns a python instance.  You are free
+    to restructure the state and version fields as needed to
+    bring the object in line with the next version, then call
+    setstate(tree) to build the return object.  Indeed this
+    technique will chain, and you can morph an ancient version
+    of your models into the latest version.
+    """
+
+    return _RefactoringRegistry.redirect(oldname, newname, asof_version)
+
+# Test classes need to be at the top level for reconstruct to find them
+class _Simple: x = 5
+class _SimpleNew(object): x = 5
+class _Slotted(object): __slots__ = ['a','b']
+class _Controlled:
+    def __getstate__(self): return ["mystate",self.__dict__]
+    def __setstate__(self, state):
+        if state[0] != "mystate": raise RuntimeError("didn't get back my state")
+        self.__dict__ = state[1]
+class _Factory: __factory__ = __name__ + "._factory"
+def _factory():
+    obj =  _Factory()
+    # Note: can't modify obj because state will be overridden
+    _Factory.fromfactory = True
+    return obj
+class _VersionError:
+    __version__ = "3.5."
+def _hello():
+    return 'hello'
+def test():
+    primitives = ['list',1,{'of':'dict',2:'really'},True,None]
+    assert deconstruct(primitives) == primitives
+    # Hmmm... dicts with non-string keys are not permitted by strict json
+    # I'm not sure we care for our purposes, but it would be best to avoid
+    # them and instead have a list of tuples which can be converted to and
+    # from a dict if the need arises
+    assert encode(primitives) == '["list",1,{"of":"dict",2:"really"},true,null]'
+
+    h = _Simple()
+    h.a = 2
+    #print encode(deconstruct(h))
+    assert decode(encode(h)).a == h.a
+    
+    assert decode(encode(_hello))() == 'hello'
+
+    h = _SimpleNew()
+    h.a = 2
+    #print encode(deconstruct(h))
+    assert decode(encode(h)).a == h.a
+
+    h = _Slotted()
+    h.a = 2
+    #print encode(deconstruct(h))
+    assert decode(encode(h)).a == h.a
+
+    h = _Controlled()
+    h.a = 2
+    #print encode(deconstruct(h))
+    assert decode(encode(h)).a == h.a
+
+    h = _Factory()
+    h.a = 2
+    #print encode(deconstruct(h))
+    assert decode(encode(h)).a == h.a
+    assert hasattr(h,'fromfactory')
+    
+    try:
+        deconstruct(_VersionError())
+        raise RuntimeError("should have raised a version error")
+    except ValueError,msg:
+        assert "_VersionError" in str(msg)
+
+if __name__ == "__main__": test()
Index: park-1.2.1/park/setup.py
===================================================================
--- park-1.2.1/park/setup.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/setup.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,27 @@
+#!/usr/bin/env python
+# This program is public domain
+
+import os.path
+from numpy.distutils.misc_util import Configuration
+from numpy.distutils.core      import setup
+
+
+def configuration(parent_package='',
+                  top_path=None
+                  ):
+    config = Configuration('park', parent_package, top_path)
+
+    # Extension reflmodule
+    sources = [os.path.join(config.package_path,'lib',s)
+               for s in ('modeling.cc','resolution.c')]
+
+    config.add_extension('_modeling',
+                         include_dirs=[config.package_path],
+                         sources=sources,
+                         )
+
+
+    return config
+
+if __name__ == '__main__':
+    setup(**configuration(top_path='').todict())
Index: park-1.2.1/park/simplex.py
===================================================================
--- park-1.2.1/park/simplex.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/simplex.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,333 @@
+#__docformat__ = "restructuredtext en"
+# ******NOTICE***************
+# from optimize.py module by Travis E. Oliphant
+#
+# You may copy and use this module as you see fit with no
+# guarantee implied provided you keep this notice in all copies.
+# *****END NOTICE************
+#
+# Modified by Paul Kienzle to support bounded minimization
+"""
+Downhill simplex optimizer.
+"""
+
+__all__ = ['simplex']
+
+__docformat__ = "restructuredtext en"
+
+import numpy
+__version__="0.7"
+
+def wrap_function(function, bounds):
+    ncalls = [0]
+    if bounds is not None:
+        lo, hi = [numpy.asarray(v) for v in bounds]
+        def function_wrapper(x):
+            ncalls[0] += 1
+            if numpy.any((x<lo)|(x>hi)):
+                return numpy.inf
+            else:
+                return function(x)
+    else:
+        def function_wrapper(x):
+            ncalls[0] += 1
+            return function(x)
+    return ncalls, function_wrapper
+
+class Result:
+    """
+    Results from the fit.
+    
+    x : ndarray
+        Best parameter set
+    fx : float
+        Best value
+    iters : int
+        Number of iterations
+    calls : int
+        Number of function calls
+    success : boolean
+        True if the fit completed successful, false if terminated early
+        because of too many iterations.
+    """
+    def __init__(self, x, fx, iters, calls, status):
+        self.x,self.fx,self.iters,self.calls=x,fx,iters,calls
+        self.status = status
+    def __str__(self):
+        return "Minimum %g at %s after %d calls"%(self.fx,self.x,self.calls)
+
+def dont_abort(): return False
+
+def simplex(f, x0=None, bounds=None, radius=0.05,
+            xtol=1e-4, ftol=1e-4, maxiter=None,
+            update_handler=None, abort_test=dont_abort):
+    """
+    Minimize a function using Nelder-Mead downhill simplex algorithm.
+
+    This optimizer is also known as Amoeba (from Numerical Recipes) and
+    the Nealder-Mead simplex algorithm.  This is not the simplex algorithm
+    for solving constrained linear systems.
+
+    Downhill simplex is a robust derivative free algorithm for finding
+    minima.  It proceeds by choosing a set of points (the simplex) forming
+    an n-dimensional triangle, and transforming that triangle so that the
+    worst vertex is improved, either by stretching, shrinking or reflecting
+    it about the center of the triangle.  This algorithm is not known for
+    its speed, but for its simplicity and robustness, and is a good algorithm
+    to start your problem with.
+
+    *Parameters*:
+
+        f : callable f(x,*args)
+            The objective function to be minimized.
+        x0 : ndarray
+            Initial guess.
+        bounds : (ndarray,ndarray) or None
+            Bounds on the parameter values for the function.
+        radius: float
+            Size of the initial simplex.  For bounded parameters (those
+            which have finite lower and upper bounds), radius is clipped
+            to a value in (0,0.5] representing the portion of the 
+            range to use as the size of the initial simplex.
+
+    *Returns*: Result (`park.simplex.Result`)
+
+        x : ndarray
+            Parameter that minimizes function.
+        fx : float
+            Value of function at minimum: ``fopt = func(xopt)``.
+        iters : int
+            Number of iterations performed.
+        calls : int
+            Number of function calls made.
+        success : boolean
+            True if fit completed successfully.
+
+    *Other Parameters*:
+
+        xtol : float
+            Relative error in xopt acceptable for convergence.
+        ftol : number
+            Relative error in func(xopt) acceptable for convergence.
+        maxiter : int=200*N
+            Maximum number of iterations to perform.  Defaults
+        update_handler : callable
+            Called after each iteration, as callback(xk,fxk), where xk
+            is the current parameter vector and fxk is the function value.
+            Returns True if the fit should continue.
+
+
+    *Notes*
+
+        Uses a Nelder-Mead simplex algorithm to find the minimum of
+        function of one or more variables.
+
+    """
+    fcalls, func = wrap_function(f, bounds)
+    x0 = numpy.asfarray(x0).flatten()
+    #print "x0",x0
+    N = len(x0)
+    rank = len(x0.shape)
+    if not -1 < rank < 2:
+        raise ValueError, "Initial guess must be a scalar or rank-1 sequence."
+
+    if maxiter is None:
+        maxiter = N * 200
+
+    rho = 1; chi = 2; psi = 0.5; sigma = 0.5;
+
+    if rank == 0:
+        sim = numpy.zeros((N+1,), dtype=x0.dtype)
+    else:
+        sim = numpy.zeros((N+1,N), dtype=x0.dtype)
+    fsim = numpy.zeros((N+1,), float)
+    sim[0] = x0
+    fsim[0] = func(x0)
+
+    # Metropolitan simplex: simplex has vertices at x0 and at
+    # x0 + j*radius for each unit vector j.  Radius is a percentage 
+    # change from the initial value, or just the radius if the initial
+    # value is 0.  For bounded problems, the radius is a percentage of
+    # the bounded range in dimension j.
+    val = x0*(1+radius)
+    val[val == 0] = radius
+    if bounds is not None:
+        radius = numpy.clip(radius,0,0.5)
+        lo,hi = [numpy.asarray(v) for v in bounds]
+
+        # Keep the initial simplex inside the bounds
+        x0[x0<lo] = lo
+        x0[x0>hi] = hi
+        bounded = ~numpy.isinf(lo) & ~numpy.isinf(hi)
+        val[bounded] = x0[bounded] + (hi[bounded]-lo[bounded])*radius
+        val[val<lo] = lo
+        val[val>hi] = hi
+        
+        # If the initial point was at or beyond an upper bound, then bounds
+        # projection will put x0 and x0+j*radius at the same point.  We
+        # need to detect these collisions and reverse the radius step
+        # direction when such collisions occur.  The only time the collision
+        # can occur at the lower bound is when upper and lower bound are
+        # identical.  In that case, we are already done.
+        collision = val==x0
+        if numpy.any(collision):
+            reverse = x0*(1-radius)
+            reverse[reverse==0] = -radius
+            reverse[bounded] = x0[bounded] - (hi[bounded]-lo[bounded])*radius
+            val[collision] = reverse[collision]
+        
+        # Make tolerance relative for bounded parameters
+        tol = numpy.ones(x0.shape)*xtol
+        tol[bounded] = (hi[bounded]-lo[bounded])*xtol
+        xtol = tol
+
+    # Compute values at the simplex vertices
+    for k in range(0,N):
+        y = x0+0
+        y[k] = val[k]
+        sim[k+1] = y
+        fsim[k+1] = func(y)
+
+    #print sim
+    ind = numpy.argsort(fsim)
+    fsim = numpy.take(fsim,ind,0)
+    # sort so sim[0,:] has the lowest function value
+    sim = numpy.take(sim,ind,0)
+    #print sim
+
+    iterations = 1
+    while iterations < maxiter:
+        if numpy.all(abs(sim[1:]-sim[0]) <= xtol) \
+            and max(abs(fsim[0]-fsim[1:])) <= ftol:
+            #print abs(sim[1:]-sim[0])
+            break
+
+        xbar = numpy.sum(sim[:-1],0) / N
+        xr = (1+rho)*xbar - rho*sim[-1]
+        #print "xbar" ,xbar,rho,sim[-1],N
+        #break
+        fxr = func(xr)
+        doshrink = 0
+
+        if fxr < fsim[0]:
+            xe = (1+rho*chi)*xbar - rho*chi*sim[-1]
+            fxe = func(xe)
+
+            if fxe < fxr:
+                sim[-1] = xe
+                fsim[-1] = fxe
+            else:
+                sim[-1] = xr
+                fsim[-1] = fxr
+        else: # fsim[0] <= fxr
+            if fxr < fsim[-2]:
+                sim[-1] = xr
+                fsim[-1] = fxr
+            else: # fxr >= fsim[-2]
+                # Perform contraction
+                if fxr < fsim[-1]:
+                    xc = (1+psi*rho)*xbar - psi*rho*sim[-1]
+                    fxc = func(xc)
+
+                    if fxc <= fxr:
+                        sim[-1] = xc
+                        fsim[-1] = fxc
+                    else:
+                        doshrink=1
+                else:
+                    # Perform an inside contraction
+                    xcc = (1-psi)*xbar + psi*sim[-1]
+                    fxcc = func(xcc)
+
+                    if fxcc < fsim[-1]:
+                        sim[-1] = xcc
+                        fsim[-1] = fxcc
+                    else:
+                        doshrink = 1
+
+                if doshrink:
+                    for j in xrange(1,N+1):
+                        sim[j] = sim[0] + sigma*(sim[j] - sim[0])
+                        fsim[j] = func(sim[j])
+
+        ind = numpy.argsort(fsim)
+        sim = numpy.take(sim,ind,0)
+        fsim = numpy.take(fsim,ind,0)
+        if update_handler is not None:
+            update_handler(sim[0],fsim[0])
+        iterations += 1
+        if abort_test(): break
+
+    status = 0 if iterations < maxiter else 1
+    res = Result(sim[0],fsim[0],iterations,fcalls[0], status)
+    return res
+
+def main():
+    import time
+    def rosen(x):  # The Rosenbrock function
+        return numpy.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0,axis=0)
+
+
+    x0 = [0.8,1.2,0.7]
+    print "Nelder-Mead Simplex"
+    print "==================="
+    start = time.time()
+    x = simplex(rosen,x0)
+    print x
+    print "Time:",time.time() - start
+
+    x0 = [0]*3
+    print "Nelder-Mead Simplex"
+    print "==================="
+    print "starting at zero"
+    start = time.time()
+    x = simplex(rosen,x0)
+    print x
+    print "Time:",time.time() - start
+
+    x0 = [0.8,1.2,0.7]
+    lo,hi = [0]*3, [1]*3
+    print "Bounded Nelder-Mead Simplex"
+    print "==========================="
+    start = time.time()
+    x = simplex(rosen,x0,bounds=(lo,hi))
+    print x
+    print "Time:",time.time() - start
+
+
+    x0 = [0.8,1.2,0.7]
+    lo,hi = [0.999]*3, [1.001]*3
+    print "Bounded Nelder-Mead Simplex"
+    print "==========================="
+    print "tight bounds"
+    print "simplex is smaller than 1e-7 in every dimension, but you can't"
+    print "see this without uncommenting the print statement simplex function"
+    start = time.time()
+    x = simplex(rosen,x0,bounds=(lo,hi),xtol=1e-4)
+    print x
+    print "Time:",time.time() - start
+
+
+    x0 = [0]*3
+    hi,lo = [-0.999]*3, [-1.001]*3
+    print "Bounded Nelder-Mead Simplex"
+    print "==========================="
+    print "tight bounds, x0=0 outside bounds from above"
+    start = time.time()
+    x = simplex(lambda x:rosen(-x),x0,bounds=(lo,hi),xtol=1e-4)
+    print x
+    print "Time:",time.time() - start
+
+
+    x0 = [0.8,1.2,0.7]
+    lo,hi = [-numpy.inf]*3, [numpy.inf]*3
+    print "Bounded Nelder-Mead Simplex"
+    print "==========================="
+    print "infinite bounds"
+    start = time.time()
+    x = simplex(rosen,x0,bounds=(lo,hi),xtol=1e-4)
+    print x
+    print "Time:",time.time() - start
+
+if __name__ == "__main__":
+    main()
Index: park-1.2.1/park/version.py
===================================================================
--- park-1.2.1/park/version.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/version.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,35 @@
+# This program is public domain
+
+"""
+version and package information about PARK.
+"""
+
+# Version for park. 
+version = '1.2'
+
+# Version for park source. 
+VERSION_SRC = '%s.src' %(version)
+
+# Package name for park server. 
+PARK_NAME = 'park'
+
+# Package name for park server + fitting service. 
+PARK_FIT_NAME = 'fitPark'
+
+#Package name for park server + fitting service + GUI client. 
+PARK_GUI_NAME = 'guiPark'
+
+# Authors for the package 
+AUTHORS = ''
+
+# Email list of Authors
+EMAILS = ''
+
+# URL of PARK 
+URL = 'http://www.reflectometry.org/danse/park'
+
+# Download url 
+DOWNLOAD_URL = 'svn co svn://svn@danse.us/park/trunk'
+
+
+
Index: park-1.2.1/park/wxmonitor.py
===================================================================
--- park-1.2.1/park/wxmonitor.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
+++ park-1.2.1/park/wxmonitor.py	(revision 3570545bb48d1a639dd12458e3c717368477e67b)
@@ -0,0 +1,221 @@
+# This program is public domain
+"""
+Asynchronous monitoring service for wx applications.
+
+Define a monitor using park.wxmonitor.wxMonitor(panel) where panel is
+the window which will receive the monitor updates.
+
+In panel, be sure to have methods for onMonitorStart(message), 
+onMonitorProgress(message), etc., for the kinds of monitor messages 
+the application will send.  The catch-all method is onMonitorMessage,
+which by default will print the messages on the console.  If you 
+don't catch onMonitorLog messages then the log messages will be 
+sent to the standard python logger.
+
+See `park.monitor` for details on the message types.
+
+Example
+=======
+
+The following defines a panel which responds to monitor messages::
+
+    import wx
+
+    class Panel(wx.Panel):
+        def __init__(self, *args, **kw):
+            wx.Panel.__init__(self, *args, **kw)
+            self.text = wx.TextCtrl(self, size=(200,100), style=wx.TE_MULTILINE)
+            self.gauge = wx.Gauge(self, range=100)
+            sizer = wx.BoxSizer(wx.VERTICAL)
+            sizer.Add(self.text, 0, wx.LEFT | wx.EXPAND)
+            sizer.Add(self.gauge, 0, wx.LEFT | wx.EXPAND)
+            self.SetSizer(sizer)
+            self.text.SetValue('starting value')
+        def onMonitorMessage(self, message):
+            self.text.SetValue(str(message))
+        def onMonitorStart(self, message):
+            self.text.SetValue(str(message))
+            self.gauge.SetValue(0)
+        def onMonitorProgress(self, message):
+            self.text.SetValue(str(message))
+            self.gauge.SetValue(int(100*message.complete/message.total))
+        def onMonitorComplete(self, message):
+            self.text.SetValue(str(message))
+            self.gauge.SetValue(100)
+
+We can put this panel in a simple app::
+
+    app = wx.PySimpleApp()
+    frame = wx.Frame(None, -1, 'Test Monitor')
+    panel = Panel(frame)
+    frame.Show()
+
+Next we attach attach the monitor to this panel and feed some messages from
+another thread::
+
+    import time,thread
+    import park.wxmonitor, park.monitor
+    from park.monitor import Start, Progress, Improvement, Complete
+    monitor = park.wxmonitor.wxMonitor(panel)
+    msgs = [Start(), Progress(1,10), Progress(3,10), 
+            Improvement('Better!'), Progerss(6,10), Complete('Best!')]:
+    def message_stream(monitor,msgs):
+        time.sleep(1)
+        for message in msgs:
+          monitor.put(message)
+          time.sleep(1)
+    thread.start_new_thread(message_stream, (monitor,msgs))
+    app.MainLoop()
+    
+You should see the progress bar jump from 10% to 30% to 60% then all the way 
+to the end.
+"""
+import logging
+import time
+
+import wx
+import wx.lib.newevent
+
+import park.monitor
+
+(MonitorEvent, EVT_MONITOR) = wx.lib.newevent.NewEvent()
+
+# For wx on Mac OS X we need to sleep after posting a message from
+# a thread in order to give the GUI a chance to update itself.
+SLEEP_TIME = 0.01
+class wxMonitor(park.monitor.Monitor):
+    """
+    Attach a job monitor to a panel.
+    
+    The monitor will perform callbacks to onMonitorStart(message),
+    onMonitorProgress(message), etc. if the associated method is
+    defined.  If the type specific method is not defined, then the
+    monitor will call onMonitorMessage(message).  Otherwise the
+    message is dropped.
+    
+    See `park.monitor` for a description of the usual messages.
+    """
+    def __init__(self, win):
+        """
+        Window to receive the monitoring events.  This is running in the
+        GUI thread.
+        """
+        self.win = win
+        win.Bind(EVT_MONITOR, self.dispatch)
+
+    def put(self, message):
+        """
+        Intercept an event received from an asynchronous monitor.  This is
+        running in the asynchronous thread.
+        """
+        #print "dispatch",message
+        event = MonitorEvent(message=message)
+        wx.PostEvent(self.win, event)
+        time.sleep(SLEEP_TIME)
+
+    def dispatch(self, event):
+        """
+        Dispatch the event from the asynchronous monitor.  This is running
+        in the GUI thread.
+        """
+        message = event.message
+        #print "window dispatch",message
+ 
+        # First check for a handler in the monitor window
+        fn = getattr(self.win, 'onMonitor'+message.__class__.__name__, None)
+        # If none, then check in our class (we have a default onMonitorLog)
+        if fn is None:
+            fn = getattr(self, 'onMonitor'+message.__class__.__name__, None)
+        # If still none, then look for the generic handler
+        if fn is None:
+            fn = getattr(self.win, 'onMonitorMessage', self.onMonitorMessage)
+        # Process the message
+        fn(message)
+
+    def onMonitorMessage(self, message):
+        """
+        Generic message handler: do nothing.
+        """
+        print ">",str(message)
+
+    def onMonitorLog(self, message):
+        """
+        Called when the job sends a logging record.
+
+        The logging record contains a normal python logging record.
+
+        The default behaviour is to tie into the application logging
+        system using::
+
+            logger = logging.getLogger(message.record.name)
+            logger.handle(message.record)
+
+        Logging levels are set in the job controller.
+        """
+        logging.basicConfig()
+        logger = logging.getLogger(message.record.name)
+        logger.handle(message.record)
+
+
+def demo(rate=0):
+    import thread
+    import time
+    import sys
+    import logging
+    
+    class Panel(wx.Panel):
+        def __init__(self, *args, **kw):
+            wx.Panel.__init__(self, *args, **kw)
+            self.text = wx.TextCtrl(self, size=(200,100), style=wx.TE_MULTILINE)
+            self.gauge = wx.Gauge(self, range=100)
+            sizer = wx.BoxSizer(wx.VERTICAL)
+            sizer.Add(self.text, 0, wx.LEFT | wx.EXPAND)
+            sizer.Add(self.gauge, 0, wx.LEFT | wx.EXPAND)
+            self.SetSizer(sizer)
+            self.text.SetValue('starting value')
+        def onMonitorMessage(self, message):
+            self.text.SetValue(str(message))
+        def onMonitorStart(self, message):
+            self.text.SetValue(str(message))
+            self.gauge.SetValue(0)
+        def onMonitorProgress(self, message):
+            self.text.SetValue(str(message))
+            self.gauge.SetValue(int(100*message.complete/message.total))
+        def onMonitorComplete(self, message):
+            self.text.SetValue(str(message))
+            self.gauge.SetValue(100)
+
+    app = wx.PySimpleApp()
+    frame = wx.Frame(None, -1, 'Test Monitor')
+    panel = Panel(frame)
+    frame.Show()
+    monitor = wxMonitor(panel)
+
+    def messagestream(monitor,rate,stream):
+        for m in stream:
+            time.sleep(rate)
+            monitor.put(m)
+        time.sleep(rate)
+        wx.CallAfter(wx.Exit)
+    R = logging.LogRecord('hi',60,'hello.py',3,'log message',(),None,'here')
+    try: raise Exception('Test exception')
+    except: trace = sys.exc_info()
+    stream=[park.monitor.Start(),
+            park.monitor.Progress(1,10),
+            park.monitor.Progress(2,10),
+            park.monitor.Progress(3,10),
+            park.monitor.Improvement('Better!'),
+            park.monitor.Abort('Abandoned'),
+            park.monitor.Start(),
+            park.monitor.Progress(1,10,'seconds'),
+            park.monitor.Improvement('Better!'),
+            park.monitor.Progress(8,10),
+            park.monitor.Complete('Best!'),
+            park.monitor.Start(),
+            park.monitor.Log(R),
+            park.monitor.Progress(6,10),
+            park.monitor.Error(trace)]
+    thread.start_new_thread(messagestream, (monitor,rate,stream))
+    app.MainLoop()
+
+if __name__  == "__main__": demo(rate=1)
