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()
