source: sasview/park-1.2.1/park/fitresult.py @ 8d27cac

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