source: sasview/park-1.2.1/park/fitmc.py @ d8c4019

# 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.
    """
    # Generate random number within bounds.  If bounds are indefinite, use [0,1]
    # If bounds are semi-definite, use [low,low+1] or [high-1,high], depending
    # on which limit is unbounded.
    lo,hi = bounds
    inf_lo = numpy.isinf(lo)
    inf_hi = numpy.isinf(hi)
    delta = hi-lo
    delta[inf_lo|inf_hi] = 1.0
    lo[inf_lo] = hi[inf_lo] - 1.0
    lo[inf_lo&inf_hi] = 0.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))