feat: init LSC fitting

cuda and cpu-based LSC fitting using cupy and numpy

Signed-off-by: Dev Singh <dev@devksingh.com>

analysis-master/requirements.txt

@@ -3,4 +3,5 @@ numpy
 scipy
 scikit-learn
 six
-matplotlib
+matplotlib
+cupy

analysis-master/tra_analysis/regression.py

@@ -2,7 +2,7 @@
 # Written by Arthur Lu & Jacob Levine
 # Notes:
 #   this module has been automatically inegrated into analysis.py, and should be callable as a class from the package
-#   this module is cuda-optimized and vectorized (except for one small part)
+#   this module is cuda-optimized (as appropriate) and vectorized (except for one small part)
 # setup:
 
 __version__ = "0.0.4"
@@ -25,7 +25,8 @@ __changelog__ = """
 
 __author__ = (
 	"Jacob Levine <jlevine@imsa.edu>",
-	"Arthur Lu <learthurgo@gmail.com>"
+	"Arthur Lu <learthurgo@gmail.com>",
+	"Dev Singh <dev@devksingh.com>"
 )
 
 __all__ = [
@@ -44,6 +45,8 @@ __all__ = [
 ]
 
 import torch
+import cupy as cp
+import numpy as np
 
 global device
 
@@ -217,4 +220,62 @@ def CustomTrain(self, kernel, optim, data, ground, loss=torch.nn.MSELoss(), iter
 				ls=loss(pred,ground_cuda)
 				ls.backward()
 				optim.step()
-		return kernel
+		return kernel
+
+class CircleFit:
+	"""Class to fit data to a circle using both the Least Square Circle (LSC) method and the HyperFit method"""
+	# For more information on the LSC method, see: 
+	# http://www.dtcenter.org/sites/default/files/community-code/met/docs/write-ups/circle_fit.pdf
+	def __init__(self, x, y, xy=None):
+		if data != None: 
+			self.coords = data
+			self.ournp = np if device === "cpu" else cp # use the correct numpy implementation based on resources available
+		else: 
+			# following block combines x and y into one array if not already done
+			self.coords = self.ournp.vstack(([x_data.T], [y_data.T])).T
+			if device !== "cpu"
+				cp.cuda.Stream.null.synchronize() # ensure code finishes executing on GPU before continuing 
+	def calc_R(x, y, xc, yc):
+		"""Returns distance between center and point"""
+		return self.ournp.sqrt((x-xc)**2 + (y-yc)**2)
+	def f(c, x, y):
+    	"""Returns distance between point and circle at c"""
+    	Ri = calc_R(x, y, *c)
+    	return Ri - Ri.mean()
+	def LSC(self):
+		"""Fits given data to a circle and returns the center, radius, and variance"""
+		x = coords[:, 0]
+		y = coords[:, 1]
+		# guessing at a center
+		x_m = self.ournp.mean(x)
+		y_m = self.ournp.mean(y)
+
+		# calculation of the reduced coordinates
+		u = x - x_m
+		v = y - y_m
+
+		# linear system defining the center (uc, vc) in reduced coordinates:
+		#    Suu * uc +  Suv * vc = (Suuu + Suvv)/2
+		#    Suv * uc +  Svv * vc = (Suuv + Svvv)/2
+		Suv  = self.ournp.sum(u*v)
+		Suu  = self.ournp.sum(u**2)
+		Svv  = self.ournp.sum(v**2)
+		Suuv = self.ournp.sum(u**2 * v)
+		Suvv = self.ournp.sum(u * v**2)
+		Suuu = self.ournp.sum(u**3)
+		Svvv = self.ournp.sum(v**3)
+
+		# Solving the linear system
+		A = self.ournp.array([ [ Suu, Suv ], [Suv, Svv]])
+		B = self.ournp.array([ Suuu + Suvv, Svvv + Suuv ])/2.0
+		uc, vc = self.ournp.linalg.solve(A, B)
+
+		xc_1 = x_m + uc
+		yc_1 = y_m + vc
+
+		# Calculate the distances from center (xc_1, yc_1)
+		Ri_1     = self.ournp.sqrt((x-xc_1)**2 + (y-yc_1)**2)
+		R_1      = self.ournp.mean(Ri_1)
+		# calcualte residual error
+		residu_1 = self.ournp.sum((Ri_1-R_1)**2)
+		return xc_1, yc_1, R_1, residu_1