""" 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:44:58 PDT 2012 Kevin Karplus Note: the version using scipy's eigh() code works with Taubin and Kasa algorithms as well as Hyper and Pratt, but still 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 scipy 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 try: eigen_vals,eigen_vects = linalg.eigh(N,M) # M is positive definite, so needs to go in the "b" position # but this means that our eigenvalues are inverted # so we want the largest eigenvalue, not the smallest positive one except np.linalg.linalg.LinAlgError: # probably M singular (noise so low that we have perfect hypersphere) # How do I get an eigenvector corresponding to a zero eigenvalue? raise ValueError("OOPS: eigh error in {}, haven't figured out what the workaround is yet".format(method)) # I'm assuming that the eigenvalues are all real, as they should be for a Hermitian matrix # print >> stderr, "DEBUG: eigen_vals=", eigen_vals # print >> stderr, "DEBUG: eigen_vects=", eigen_vects # print >> stderr, "DEBUG: list(eigen_vals)=", list(eigen_vals) param_vect =eigen_vects[:,list(eigen_vals).index(max(eigen_vals))].transpose() # print >> stderr, "DEBUG: param_vect=", param_vect # params = np.asarray(param_vect)[0] # numpy gets 1xn matrix for param_vect, scipy gets array of n params = np.asarray(param_vect) # print >> stderr, "DEBUG: params=", 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)