Source Code for Module mvpa.misc.fx

  1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 
  2  # vi: set ft=python sts=4 ts=4 sw=4 et: 
  3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  4  # 
  5  #   See COPYING file distributed along with the PyMVPA package for the 
  6  #   copyright and license terms. 
  7  # 
  8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  9  """Misc. functions (in the mathematical sense)""" 
 10   
 11  __docformat__ = 'restructuredtext' 
 12   
 13  import numpy as N 
 14   
 15   
 16 -def singleGammaHRF(t, A=5.4, W=5.2, K=1.0): 
 17      """Hemodynamic response function model. 
 18   
 19      The version consists of a single gamma function (also see 
 20      doubleGammaHRF()). 
 21   
 22      :Parameters: 
 23        t: float 
 24          Time. 
 25        A: float 
 26          Time to peak. 
 27        W: float 
 28          Full-width at half-maximum. 
 29        K: float 
 30          Scaling factor. 
 31      """ 
 32      A = float(A) 
 33      W = float(W) 
 34      K = float(K) 
 35      return K * (t / A) ** ((A ** 2) / (W ** 2) * 8.0 * N.log(2.0)) \ 
 36             * N.e ** ((t - A) / -((W ** 2) / A / 8.0 / N.log(2.0))) 
 37   
 38   
 39 -def doubleGammaHRF(t, A1=5.4, W1=5.2, K1=1.0, A2=10.8, W2=7.35, K2=0.35): 
 40      """Hemodynamic response function model. 
 41   
 42      The version is using two gamma functions (also see singleGammaHRF()). 
 43   
 44      :Parameters: 
 45        t: float 
 46          Time. 
 47        A: float 
 48          Time to peak. 
 49        W: float 
 50          Full-width at half-maximum. 
 51        K: float 
 52          Scaling factor. 
 53   
 54      Parameters A, W and K exists individually for each of both gamma 
 55      functions. 
 56      """ 
 57      A1 = float(A1) 
 58      W1 = float(W1) 
 59      K1 = float(K1) 
 60      A2 = float(A2) 
 61      W2 = float(W2) 
 62      K2 = float(K2) 
 63      return singleGammaHRF(t, A1, W1, K1) - singleGammaHRF(t, A2, W2, K2) 
 64   
 65   
 66 -def leastSqFit(fx, params, y, x=None, **kwargs): 
 67      """Simple convenience wrapper around SciPy's optimize.leastsq. 
 68   
 69      The advantage of using this wrapper instead of optimize.leastsq directly 
 70      is, that it automatically constructs an appropriate error function and 
 71      easily deals with 2d data arrays, i.e. each column with multiple values for 
 72      the same function argument (`x`-value). 
 73   
 74      :Parameters: 
 75        fx: functor 
 76          Function to be fitted to the data. It has to take a vector with 
 77          function arguments (`x`-values) as the first argument, followed by 
 78          an arbitrary number of (to be fitted) parameters. 
 79        params: sequence 
 80          Sequence of start values for all to be fitted parameters. During 
 81          fitting all parameters in this sequences are passed to the function 
 82          in the order in which they appear in this sequence. 
 83        y: 1d or 2d array 
 84          The data the function is fitted to. In the case of a 2d array, each 
 85          column in the array is considered to be multiple observations or 
 86          measurements of function values for the same `x`-value. 
 87        x: Corresponding function arguments (`x`-values) for each datapoint, i.e. 
 88          element in `y` or columns in `y', in the case of `y` being a 2d array. 
 89          If `x` is not provided it will be generated by `N.arange(m)`, where 
 90          `m` is either the length of `y` or the number of columns in `y`, if 
 91          `y` is a 2d array. 
 92        **kwargs: 
 93          All additonal keyword arguments are passed to `fx`. 
 94   
 95      :Returns: 
 96        tuple: as returned by scipy.optimize.leastsq 
 97          i.e. 2-tuple with list of final (fitted) parameters of `fx` and an 
 98          integer value indicating success or failure of the fitting procedure 
 99          (see leastsq docs for more information). 
100      """ 
101      # import here to not let the whole module depend on SciPy 
102      from scipy.optimize import leastsq 
103   
104      y = N.asanyarray(y) 
105   
106      if len(y.shape) > 1: 
107          nsamp, ylen = y.shape 
108      else: 
109          nsamp, ylen = (1, len(y)) 
110   
111      # contruct matching x-values if necessary 
112      if x is None: 
113          x = N.arange(ylen) 
114   
115      # transform x and y into 1d arrays 
116      if nsamp > 1: 
117          x = N.array([x] * nsamp).ravel() 
118          y = y.ravel() 
119   
120      # define error function 
121      def efx(p): 
122          err = y - fx(x, *p, **kwargs) 
123          return err 
124   
125      # do fit 
126      return leastsq(efx, params) 
127