Source code for PyMca5.PyMcaMath.linalg

#    Copyright (c) 2008-2014 V.A. Sole, ESRF
#
#   Permission to use and redistribute the source code or binary forms of
#   this software and its documentation, with or without modification is
#   hereby granted provided that the above notice of copyright, these
#   terms of use, and the disclaimer of warranty below appear in the
#   source code and documentation, and that none of the names of The
#   European Synchrotron Radiation Facility, or the authors
#   appear in advertising or endorsement of works derived from this
#   software without specific prior written permission from all parties.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
#   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
#   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
#   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
#   CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
#   TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
#   SOFTWARE OR THE USE OR OTHER DEALINGS IN THIS SOFTWARE.
#
import numpy
__license__ = "BSD"
__author__ = "V.A. Sole - ESRF Data Analysis"
__doc__ = """
Similar function to the numpy lstsq function with a more rigorous uncertainty
treatement besides other optimizations in view of simultaneously solving several
equations of the form `a x = b`. Hopefully licensed under the same terms as
numpy itself (BSD license).
"""

# Linear Least Squares

[docs]def lstsq(a, b, rcond=None, sigma_b=None, weight=False, uncertainties=True, covariances=False, digested_output=False, svd=True, last_svd=None): """ Return the least-squares solution to a linear matrix equation. Solves the equation `a x = b` by computing a vector `x` that minimizes the Euclidean 2-norm `|| b - a x ||^2`. The equation may be under-, well-, or over- determined (i.e., the number of linearly independent rows of `a` can be less than, equal to, or greater than its number of linearly independent columns). If `a` is square and of full rank, then `x` (but for round-off error) is the "exact" solution of the equation. Parameters ---------- a : array_like, shape (M, N) "Model" matrix. b : array_like, shape (M,) or (M, K) Ordinate or "dependent variable" values. If `b` is two-dimensional, the least-squares solution is calculated for each of the `K` columns of `b`. sigma_b : uncertainties on the b values or None. If sigma_b has shape (M,) or (M, 1) and b has dimension (M, K), the uncertainty will be the same for all spectra. weight: 0 - No data weighting. If required, uncertainties will be calculated using either the supplied experimental uncertainties or an experimental uncertainty of 1 for each data point. 1 - Statistical weight. Weighted fit using the supplied experimental uncertainties or the square root of the b values. svd: If not true, a simple matrix inversion will be used in case of weighting with unequal data weights. Ignored in any other cases. last_svd: Tuple containing U, s, V of the weighted model matrix or None. This is to prevent recalculation on repeated fits. uncertainties: If False, no uncertainties will be calculated unless the covariance matrix is requested. covariances: If True, an array of covariance matrix/matrices will be returned. digested_output: If True, returns a dictionnary with explicit keys Returns ------- x : ndarray, shape (N,) or (N, K) Least-squares solution. The shape of `x` depends on the shape of `b`. uncertainties: ndarray, shape (N,) or (N, K) covariances: ndarray, shape (N, N) or (K, N, N) Examples -------- Fit a line, ``y = mx + c``, through some noisy data-points: >>> x = np.array([0, 1, 2, 3]) >>> y = np.array([-1, 0.2, 0.9, 2.1]) By examining the coefficients, we see that the line should have a gradient of roughly 1 and cut the y-axis at, more or less, -1. We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]`` and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`: >>> A = np.vstack([x, np.ones(len(x))]).T >>> A array([[ 0., 1.], [ 1., 1.], [ 2., 1.], [ 3., 1.]]) >>> m, c = np.linalg.lstsq(A, y)[0] >>> print m, c 1.0 -0.95 Plot the data along with the fitted line: >>> import matplotlib.pyplot as plt >>> plt.plot(x, y, 'o', label='Original data', markersize=10) >>> plt.plot(x, m*x + c, 'r', label='Fitted line') >>> plt.legend() >>> plt.show() """ a = numpy.array(a, dtype=numpy.float, copy=False) b = numpy.array(b, dtype=numpy.float, copy=False) a_shape = a.shape b_shape = b.shape original = b_shape if len(a_shape) != 2: raise ValueError("Model matrix must be two dimensional") if len(b_shape) == 1: b.shape = b_shape[0], 1 b_shape = b.shape m = a.shape[0] n = a.shape[1] if m != b.shape[0]: raise ValueError('Incompatible dimensions between A and b matrices') fastest = False if weight: if sigma_b is not None: # experimental uncertainties provided these are the ones to use (if any) w = numpy.abs(numpy.array(sigma_b, dtype=numpy.float, copy=False)) w = w + numpy.equal(w, 0) if w.size == b_shape[0]: # same uncertainty for every spectrum fastest = True w.shape = b.shape[0] else: w.shape = b_shape else: # "statistical" weight # we are asked to somehow weight the data but no uncertainties provided # assume the uncertainties are the square root of the b values ... w = numpy.sqrt(numpy.abs(b)) w = w + numpy.equal(w, 0) else: # we have an unweighted fit with no uncertainties # assume all the uncertainties equal to 1 fastest = True w = numpy.ones(b.shape, numpy.float) if covariances: covarianceMatrix = numpy.zeros((b_shape[1], n, n), numpy.float) if not weight: # no weight is applied # get the SVD decomposition of the A matrix # One could avoid calculating U, s, V each time ... if last_svd is not None: U, s, V = last_svd else: U, s, V = numpy.linalg.svd(a, full_matrices=False) if rcond is None: s_cutoff = n * numpy.finfo(numpy.float).eps else: s_cutoff = rcond * s[0] s[s < s_cutoff] = numpy.inf # and get the parameters s.shape = -1 dummy = numpy.dot(V.T, numpy.eye(n)*(1./s)) parameters = numpy.dot(dummy, numpy.dot(U.T, b)) parameters.shape = n, b.shape[1] if uncertainties or covariances: # get the uncertainties #(in the no-weight case without experimental uncertainties, # the uncertainties on the data points are ignored and the # uncertainty on the fitted parameters are independent of the input data!!!!) if fastest: # This is correct for all weights equal to 1 _covariance = numpy.dot(dummy, dummy.T) sigmapar = numpy.sqrt(numpy.diag(_covariance)) sigmapar = numpy.outer(sigmapar, numpy.ones(b_shape[1])) sigmapar.shape = n, b_shape[1] if covariances: covarianceMatrix[:] = _covariance elif covariances: # loop in order not to use potentially big matrices # but calculates the covariance matrices # It only makes sense if the covariance matrix is requested sigmapar = numpy.zeros((n, b_shape[1]), numpy.float) for k in range(b_shape[1]): pseudoData = numpy.eye(b_shape[0]) * w[:, k] tmpTerm = numpy.dot(dummy, numpy.dot(U.T, pseudoData)) _covariance[:, :] = numpy.dot(tmpTerm, tmpTerm.T) sigmapar[:, k] = numpy.sqrt(numpy.diag(_covariance)) covarianceMatrix[k] = _covariance else: # loop in order not to use potentially big matrices # but not calculating the covariance matrix d = numpy.zeros(b.shape, numpy.float) sigmapar = numpy.zeros((n, b_shape[1])) for k in range(b_shape[0]): d[k] = w[k] sigmapar += (numpy.dot(dummy, numpy.dot(U.T, d))) ** 2 d[k] = 0.0 sigmapar[:, :] = numpy.sqrt(sigmapar) elif fastest: # same weight for all spectra # it could be made by the calling routine, because it is equivalent to supplying a # different model and different independent values ... # That way one could avoid calculating U, s, V each time A = a / weight b = b / weight # get the SVD decomposition of the A matrix if last_svd is not None: U, s, V = last_svd else: U, s, V = numpy.linalg.svd(A, full_matrices=False) if rcond is None: s_cutoff = n * numpy.finfo(numpy.float).eps else: s_cutoff = rcond * s[0] s[s < s_cutoff] = numpy.inf # and get the parameters s.shape = -1 dummy = numpy.dot(V.T, numpy.eye(n)*(1./s)) parameters = numpy.dot(dummy, numpy.dot(U.T, b)) parameters.shape = n, b.shape[1] if uncertainties or covariances: _covariance = numpy.dot(dummy, dummy.T) sigmapar = numpy.sqrt(numpy.diag(_covariance)) sigmapar = numpy.outer(sigmapar, numpy.ones(b_shape[1])) sigmapar.shape = n, b_shape[1] if covariances: covarianceMatrix[:] = _covariance else: parameters = numpy.zeros((n, b_shape[1]), numpy.float) sigmapar = numpy.zeros((n, b_shape[1]), numpy.float) if svd: # SVD - slower by a factor 2 for i in range(b_shape[1]): tmpWeight = w[:, i:i+1] tmpData = b[:, i:i+1] / tmpWeight A = a / tmpWeight U, s, V = numpy.linalg.svd(A, full_matrices=False) if rcond is None: s_cutoff = n * numpy.finfo(numpy.float).eps else: s_cutoff = rcond * s[0] s[s < s_cutoff] = numpy.inf s.shape = -1 dummy = numpy.dot(V.T, numpy.eye(n)*(1./s)) parameters[:, i:i+1] = numpy.dot(dummy, numpy.dot(U.T, tmpData)) if uncertainties or covariances: # get the uncertainties _covariance = numpy.dot(dummy, dummy.T) sigmapar[:, i] = numpy.sqrt(numpy.diag(_covariance)) if covariances: covarianceMatrix[i] = _covariance elif 1: # Pure matrix inversion (faster than SVD) # I do not seem to gain anything by re-using the storage #alpha = numpy.empty((n, n), numpy.float) #beta = numpy.empty((n, 1), numpy.float) for i in range(b_shape[1]): tmpWeight = w[:, i:i+1] A = a / tmpWeight tmpData = b[:, i:i+1] / tmpWeight #numpy.dot(A.T, A, alpha) #numpy.dot(A.T, tmpData, beta) alpha = numpy.dot(A.T, A) beta = numpy.dot(A.T, tmpData) try: _covariance = numpy.linalg.inv(alpha) except: print("Exception") print("Exception", sys.exc_info()[1]) continue parameters[:, i:i+1] = numpy.dot(_covariance, beta) if uncertainties: sigmapar[:, i] = numpy.sqrt(numpy.diag(_covariance)) if covariances: covarianceMatrix[i] = covariance else: # Matrix inversion with buffers does not improve bufferProduct = numpy.empty((n, n + 1), numpy.float) bufferAB = numpy.empty((b_shape[0], n+1), numpy.float) alpha = numpy.empty((n, n), numpy.float) for i in range(b_shape[1]): tmpWeight = w[:, i:i+1] A = a / tmpWeight tmpData = b[:, i:i+1] / tmpWeight bufferAB [:, :n] = A bufferAB [:, n:n+1] = tmpData numpy.dot(A.T, bufferAB, bufferProduct) alpha[:, :] = bufferProduct[:n, :n] beta = bufferProduct[:,n] try: _covariance = numpy.linalg.inv(alpha) except: print("Exception") print("Exception", sys.exc_inf()) continue parameters[:, i] = numpy.dot(_covariance, beta) if uncertainties: sigmapar[:, i] = numpy.sqrt(numpy.diag(_covariance)) if covariances: covarianceMatrix[i] = covariance if len(original) == 1: parameters.shape = -1 if covariances: sigmapar.shape = parameters.shape if len(original) == 1: covarianceMatrix.shape = parameters.shape[0], parameters.shape[0] result = [parameters, sigmapar, covarianceMatrix] elif uncertainties: sigmapar.shape = parameters.shape result = [parameters, sigmapar] else: result = [parameters] if digested_output: ddict = {} ddict['parameters'] = result[0] if len(result) > 1: ddict['uncertainties'] = result[1] elif covariances: ddict['covariances'] = result[2] if svd or fastest: ddict['svd'] = (U, s, V) return ddict else: return result
[docs]def getModelMatrixFromFunction(model_function, dummy_parameters, xdata, derivative=None): nPoints = xdata.size nParameters = len(dummy_parameters) modelMatrix = numpy.zeros((nPoints, nParameters) , numpy.float) pwork = dummy_parameters * 1 for i in range(len(dummy_parameters)): fitparam = dummy_parameters[i] if derivative is None: delta = (pwork[i] + numpy.equal(fitparam, 0.0)) * 0.00001 pwork[i] = fitparam + delta f1 = model_function(pwork, xdata) pwork[i] = fitparam - delta f2 = model_function(pwork, xdata) help0 = (f1-f2) / (2.0 * delta) pwork[i] = fitparam else: help0 = derivative(pwork, i, xdata) help0.shape = -1 modelMatrix[:, i] = help0 return modelMatrix
[docs]def modelFunction(p, x): return p[0] + (p[1] + p[2] * x) * x
[docs]def test1(): x = numpy.arange(10000.) x.shape = -1, 1 y = modelFunction([100., 50., 4.], x) A = getModelMatrixFromFunction(modelFunction, [0.0, 0.0, 0.0], x) parameters, uncertainties = lstsq(A, y, uncertainties=True, weight=False) print("Expected = 100., 50., 4.") print("Obtained = %f, %f, %f" % (parameters[0], parameters[1], parameters[2]))
[docs]def test2(): import time try: from PyMca5 import Gefit GEFIT = True def f(p, x): return p[1] * x + p[0] except: GEFIT = False data = "0 0.8214 0.1 1 2.8471 0.3 2 4.852 0.5 3 7.5347 0.7 4 10.2464 0.9 5 10.2707 1.1 6 12.8011 1.3 7 13.7108 1.5 8 17.8501 1.7 9 15.3667 1.9 10 19.3933 2.1" data = numpy.array([float(x) for x in data.split()]) data.shape = -1, 3 # the model matrix for a straight line A = numpy.ones((data.shape[0],2), numpy.float) A[:, 1] = data[:, 0] print("Unweighted results:") t0 = time.time() y = numpy.ones((data.shape[0], 1000), numpy.float) * data[:, 1:2] sigmay = numpy.ones((data.shape[0], 1000), numpy.float) * data[:, 2:3] parameters, uncertainties = lstsq(A, y, #sigma_b=sigmay, #sigma_b=numpy.ones(sigmay.shape), uncertainties=True, weight=False) print("Elapsed = %f" % (time.time() - t0)) print("Parameters = %f, %f" % (parameters[0,100], parameters[1, 100])) print("Uncertainties = %f, %f" % (uncertainties[0,100], uncertainties[1, 100])) if GEFIT: t0 = time.time() for i in range(y.shape[1]): parameters, chisq, uncertainties = Gefit.LeastSquaresFit(f, [0.0, 0.0], xdata=data[:,0], ydata=data[:,1], sigmadata=data[:,2], weightflag=0, linear=1) print("Elapsed = %f" % (time.time() - t0)) print("Gefit results:") print("Parameters = %f, %f" % (parameters[0], parameters[1])) print("Uncertainties = %f, %f" % (uncertainties[0], uncertainties[1])) print("Mathematica results:") print("Parameters = %f, %f" % (1.57043, 1.78945)) print("Uncertainties = %f, %f" % (0.68363, 0.11555)) print("Weighted results") t0 = time.time() #parameters, uncertainties = lstsq(A, data[:, 1], sigma_b=data[:,2], parameters, uncertainties = lstsq(A, y, sigma_b=numpy.outer(data[:,2], numpy.ones((1000, 1))), uncertainties=True, weight=True) print("Elapsed = %f" % (time.time() - t0)) print("Parameters = %f, %f" % (parameters[0, 100], parameters[1, 100])) print("Uncertainties = %f, %f" % (uncertainties[0, 100], uncertainties[1, 100])) if GEFIT: parameters, chisq, uncertainties = Gefit.LeastSquaresFit(f, [0.0, 0.0], xdata=data[:,0], ydata=data[:,1], sigmadata=data[:,2], weightflag=1, linear=1) print("Gefit results:") print("Parameters = %f, %f" % (parameters[0], parameters[1])) print("Uncertainties = %f, %f" % (uncertainties[0], uncertainties[1])) print("Mathematica results:") print("Parameters = %f, %f" % (0.843827, 1.97982)) print("Uncertainties = %f, %f" % (0.092449, 0.07262)) return data
if __name__ == "__main__": test1() test2()