""" FitHypersphere.py fit_hypersphere(collection of tuples or lists of real numbers) will return a hypersphere of the same dimension as the tuples: (radius, (center)) using the Hyper (hyperaccurate) algorithm of Ali Al-Sharadqah and Nikolai Chernov Error analysis for circle fitting algorithms Electronic Journal of Statistics Vol. 3 (2009) 886-911 DOI: 10.1214/09-EJS419 generalized to n dimensions Sun Apr 22 16:45:44 PDT 2012 Kevin Karplus Note: this version using numpy's eig() works with Hyper and Pratt (whose constraint matrices are non-singular), but not with Taubin and Kasa, which have singular constraint matrices. It fails for too-low noise (a singular M matrix). Creative Commons Attribution-ShareAlike 3.0 Unported License. http://creativecommons.org/licenses/by-sa/3.0/ """ import numpy as np from numpy import linalg from sys import stderr from math import sqrt def fit_hypersphere(data, method="Hyper"): """returns a hypersphere of the same dimension as the collection of input tuples (radius, (center)) Methods available for fitting are "algebraic" fitting methods Hyper Al-Sharadqah and Chernov's Hyperfit algorithm Pratt Vaughn Pratt's algorithm The following methods, though very similar, are not implemented yet, because the contraint matrix N would be singular, and so the N.I*M computation is not doable. Taubin G. Taubin's algorithm Kasa Kasa's algorithm """ num_points = len(data) # print >>stderr, "DEBUG: num_points=", num_points if num_points==0: return (0,None) if num_points==1: return (0,data[0]) dimen = len(data[0]) # dimensionality of hypersphere # print >>stderr, "DEBUG: dimen=", dimen if num_points>stderr, "DEBUG: square_mag=", square_mag # central n columns of matrix center = np.asmatrix(data, dtype=float) # print >>stderr, "DEBUG: center=", center # matrix of data data_M = np.bmat( [[square_mag, center, np.ones((num_points,1))]]) # print >> stderr, "DEBUG: data_M=",data_M # matrix of 2nd moments: M = data_M.transpose() * data_M Means = np.mean(data_M, axis=0) # print >> stderr, "DEBUG: M=",M # print >> stderr, "DEBUG: Means=",Means # construct constraint matrix N = np.asmatrix(np.identity(dimen+2, dtype=float)) N[0,0] = 0 N[-1,-1]=0 if method=="Hyper": row = 4*Means row[0,-1]=2 N[0,:] = row N[:,0] += row.transpose() elif method=="Pratt": # constraint is sum of squares for coordinates minus # 4 * N[0,-1]= -2 N[-1,0]=-2 elif method=="Taubin": row = 2*Means row[0,-1]=0 N[0,:] = row N[:,0] += row.transpose() elif method=="Kasa": N[0,0]=1 for i in range(1,dimen+1): N[i,i]=0 else: raise ValueError("Error: unknown method: {} should be one of 'Hyper', 'Pratt', or 'Taubin'") # print >> stderr, "DEBUG:", method, "N=",N # print >> stderr, "DEBUG: N.I=", N.I # print >> stderr, "DEBUG: N.I*N=",N.I*N NM = N.I*M eigen_vals,eigen_vects = linalg.eig(NM) # I'm assuming that the eigenvalues are all real. # print >> stderr, "DEBUG: eigen_vals=", eigen_vals # print >> stderr, "DEBUG: eigen_vects=", eigen_vects # print >> stderr, "DEBUG: list(eigen_vals)=", list(eigen_vals) positives = [x for x in eigen_vals if x>0] if len(positives)+1 != len(eigen_vals): # raise ValueError("Error: for method {} exactly one eigenvalue should be negative: {}".format(method,eigen_vals)) print>>stderr, "Warning: for method {} exactly one eigenvalue should be negative: {}".format(method,eigen_vals) smallest_positive = min(positives) # print >> stderr, "DEBUG: smallest_positive=", smallest_positive param_vect =eigen_vects[:,list(eigen_vals).index(smallest_positive)].transpose() params = np.asarray(param_vect)[0] # print >> stderr, "DEBUG: params=", repr(params) radius = 0.5* sqrt( sum(a*a for a in params[1:-1])- 4*params[0]*params[-1])/abs(params[0]) center = tuple(-0.5*params[1:-1]/params[0]) return (radius,center)