testing/common/fit.py

   1 # -*- coding: utf-8 -*-
   2
   3 """Fit histogram data coming in on stdin to a user-defined model.  You
   4 should write a python script sublcassing fit-hist.Model."""
   5
   6 import scipy.optimize, scipy
   7 import sys
   8
   9 kB=1.3806503e-23 # m²·kg/s²·K
  10
  11 class Model (object) :
  12
  13     """For minimizing fit error: returns the normalized rms difference
  14     between a dataset and a function y(X).  In order to use this
  15     class, you must sublcass it and overwrite Model.model(x,params).
  16     It is recommended that you also overwrite
  17     Model.guess_initial_params(), and then not pass initializing
  18     params to Model.model().  This helps ensure good automation for
  19     unsupervised testing.  Developing good heuristics for
  20     Model.guessInitialParams() is usually the hardest part writing a
  21     theoretical histogram test.
  22
  23     For bonus points, you can override Model.printModel() and
  24     Model.printParams(), to give your users an easy to understand idea
  25     of what's going on.  It's better to be verbose and clear.
  26     Remember, this is the testing branch.  It needs to be solid.  Save
  27     your wonderful, sneaky hacks for the main code ;).
  28
  29     """
  30     def __init__(self, numParams, data) :
  31         """Initialize the Model instance with a set of data, for
  32         example: data = [(0,4),(1,3),(2,3),(4,0)]"""
  33         self.data = data # an iterable filled with (x,y) pairs
  34         # N (# unfolding events) = the sum of the y values (# hits per bin)
  35         self.N = sum(zip(*data)[1]) # zip(* takes data to [[x0,x1,..],[y..]]
  36         # assumes constant binwidth
  37         self.binwidth = self.data[1][0] - self.data[0][0]
  38         self.numParams = numParams
  39     def nrms(self, params) :
  40         """Returns the normalized root mean squagred between the data
  41         and a model for a given set of parameters"""
  42         s = 0.0
  43         for x,y in self.data :
  44             yResidual = y-self.model(x, params)
  45             s += yResidual**2
  46         s /= float(len(self.data))
  47         return s
  48     def fit(self, x0=None) :
  49         """Minimize error given an initial guess vector x0, the heart of
  50         the class."""
  51         if x0 == None :
  52             x0 = self.guessInitialParams()
  53         xopt = scipy.optimize.fmin_powell(self.nrms, x0)
  54         return xopt
  55     def printModel(self) :
  56         "OVERRIDE: Give the user some idea of what model we're fitting with"
  57         print "model string not defined"
  58     def printParams(self) :
  59         "OVERRIDE: Give the user nicely named format for a parameter vector"
  60         for i,p in zip(range(len(params)), params) :
  61             print "param_%d\t=\t%g\n" % (i,p)
  62     def model(self, x, params) :
  63         "OVERRIDE: Implement your model"
  64         assert len(params) == self.numParams
  65         raise Exception, "model method not defined"
  66         return 1.0
  67     def guessInitialParams(self) :
  68         """OVERRIDE: Provide some intelligent heuristic to pick a
  69         starting point for minimizing your model"""
  70         return [1.0] * self.numParams
  71
  72 def readData(fileObj) :
  73     """Read histogram data from a whitespace-seperated, two-column
  74     ASCII text-file into a list of tuples"""
  75     data = []
  76     for line in fileObj :
  77         if line[0] == "#" : continue
  78         data.append([float(x) for x in line.split()])
  79     return data
  80
  81 def run(numParams, modelClass) :
  82     import sys
  83     data = readData(sys.stdin)
  84     mod = modelClass(numParams, data)
  85     popt = mod.fit()
  86     mod.printModel()
  87     mod.printParams(popt)
  88     print "N:\t%g" % mod.N
  89     print "bin:\t%g" % mod.binwidth