The calculation of coefficient errors in the Polynomial fitting class

examples/OrbitUtils_Tests/Polynomials/polynomial_fitting_test.py

@@ -0,0 +1,40 @@
+#----------------------------------------------
+# This is an example of a polymonial fitting 
+#----------------------------------------------
+
+import sys
+import math
+import random
+
+from orbit.core.orbit_utils import Matrix
+from orbit.core.orbit_utils import Function
+from orbit.core.orbit_utils import SplineCH
+
+from orbit.utils.fitting import PolynomialFit
+
+random_err_range = 0.001
+f = Function()
+for i in range(15):
+	x = 0.1*i
+	y = -1+1.0*x**2+2.0*x**3 + 0.01*x**4
+	y += random.uniform(-random_err_range,+random_err_range)
+	f.add(x,y)
+
+print ("polynomial for fitting = y = -1+1.0*x**2+2.0*x**3 + 0.01*x**4")
+print ("================Function====================================")
+poly = PolynomialFit(4)
+poly.fitFunction(f)
+[coef_arr,err_arr] = poly.getCoefficientsAndErr()
+for i in range(len(coef_arr)):
+	print ("i=",i," coef = %8.5f +- %8.5f"%(coef_arr[i],err_arr[i]))
+
+print ("================SplineCH====================================")
+spline = SplineCH()
+spline.compile(f)
+poly.fitSpline(f)	
+[coef_arr,err_arr] = poly.getCoefficientsAndErr()
+for i in range(len(coef_arr)):
+	print ("i=",i," coef = %8.5f +- %8.5f"%(coef_arr[i],err_arr[i]))
+
+print ("Stop.")
+sys.exit(0)

examples/OrbitUtils_Tests/Polynomials/polynomial_test.py

@@ -0,0 +1,65 @@
+#-------------------------------------------------------------
+# This is an example of Polynomial class
+# The methods "derivativeTo","derivative", "copyTo" and "value" are included.
+#-------------------------------------------------------------
+
+import sys
+import math
+
+from orbit.core.orbit_utils import Polynomial
+
+poly = Polynomial(4)
+print ("initial order = ",poly.order())
+
+poly.coefficient(2,2.0)
+
+print ("========Polynomial=======")
+order = poly.order()
+for i in range(order+1):
+	print ("i=",i," coef=",poly.coefficient(i))
+
+x = 10.
+y = poly.value(x)
+print ("x=",x," y=",y)
+
+poly1 = Polynomial()
+poly.derivativeTo(poly1)
+
+print ("========Derivative polynomial=======")
+order = poly1.order()
+for i in range(order+1):
+	print ("i=",i," coef=",poly1.coefficient(i))
+
+x = 10.
+y = poly1.value(x)
+yp = poly.derivative(x)
+print ("x=",x," y_prime=",y,"  yp =",yp)
+
+print ("========Copy of initial polynomial=======")
+poly2 = Polynomial()
+poly.copyTo(poly2)
+order = poly2.order()
+for i in range(order+1):
+	print ("i=",i," coef=",poly2.coefficient(i))
+
+x = 10.
+y = poly2.value(x)
+yp = poly2.derivative(x)
+print ("x=",x," y=",y," y_prime=",yp)
+
+"""
+# Memory leak test
+count = 1
+while(1 < 2):
+ poly1 = Polynomial()
+ poly2 = Polynomial()
+ poly.derivativeTo(poly1)
+ poly.copyTo(poly2)
+ count += 1
+ if(count % 100000 == 0):
+	 print ("count=",count)
+"""
+
+print ("Stop.")
+sys.exit(0)
+

py/orbit/utils/fitting/PolynomialFit.py

@@ -7,7 +7,9 @@
 import sys
 import math
 
-from orbit.core.orbit_utils import Function, SplineCH, Matrix, Polynomial
+from orbit.core.orbit_utils import Function, SplineCH
+from orbit.core.orbit_utils import PhaseVector, Matrix
+from orbit.core.orbit_utils import Polynomial
 
 
 class PolynomialFit:
@@ -60,68 +62,42 @@ def fitSpline(self, spline):
 
     def _makePolynomial(self):
         nPoints = len(self.x_y_err_arr)
-        if nPoints < (self.order + 1):
+        if(nPoints < (self.order + 1)):
             self.order = nPoints - 1
-        # check if just one of errors is zero
-        infoZeroErr = 1.0
-        for [x, y, err] in self.x_y_err_arr:
-            infoZeroErr *= err
-        for i in range(nPoints):
-            [x, y, err] = self.x_y_err_arr[i]
-            sigma = 1.0
-            if infoZeroErr != 0.0:
-                sigma = 1.0 / (err * err)
-            self.x_y_err_arr[i][2] = sigma
-        # now make A matrix
-        aMatr = Matrix(nPoints, self.order + 1)
+        #---- now make A and A^T matrices and y-vector
+        aMtr = Matrix(nPoints, self.order + 1)
+        yVct = PhaseVector(nPoints)
         for i in range(nPoints):
+            y = self.x_y_err_arr[i][1]
+            yVct.set(i,y)
             for j in range(self.order + 1):
                 x = self.x_y_err_arr[i][0]
-                aMatr.set(i, j, math.pow(x, j))
-        aTCa = Matrix(self.order + 1, self.order + 1)
-        for i in range(self.order + 1):
-            for j in range(self.order + 1):
-                a = 0.0
-                for k in range(nPoints):
-                    sigma = self.x_y_err_arr[k][2]
-                    a += aMatr.get(k, i) * sigma * aMatr.get(k, j)
-                aTCa.set(i, j, a)
+                aMtr.set(i, j, math.pow(x, j))
+        aMtrT = aMtr.transpose()
         # now the resuting coefficients and errors
-        aTCaI = aTCa.invert()
-        if aTCaI == None:
-            print("python PolynomialFit: Problem with data.")
-            for i in range(nPoints):
-                x = self.x_y_err_arr[i][0]
-                y = self.x_y_err_arr[i][1]
-                err = self.x_y_err_arr[i][2]
-                print(" x,y,err = %12.5g %12.5g %12.5g " % (x, y, err))
-            print("Stop.")
-            sys.exit(1)
-        e = aTCaI.mult(aTCa)
-        coef_arr = [0.0] * (self.order + 1)
-        err_arr = [0.0] * (self.order + 1)
+        ATA_inv = (aMtrT.mult(aMtr)).invert()
+        coeffVct = (ATA_inv.mult(aMtrT)).mult(yVct)
+        # polinimial coefficients have been found -> coeffVct
+        coef_arr = [0.0] * (self.order + 1)   
+        self.polynomial.order(self.order)
         for i in range(self.order + 1):
-            err_arr[i] = math.sqrt(math.fabs(aTCaI.get(i, i)))
+            coef_arr[i] = coeffVct.get(i)
+            self.polynomial.coefficient(i,coef_arr[i])
+        #---- here we estimate errors by deviation y_fit from y_init only
+        #---- It means that we do not pay attention to different weights
+        #---- of the initial Y data errors
+        coeff_err_arr = [0.0] * (self.order + 1)
         for i in range(self.order + 1):
-            coef_arr[i] = 0.0
-            for j in range(self.order + 1):
-                for k in range(nPoints):
-                    sigma = self.x_y_err_arr[k][2]
-                    y = self.x_y_err_arr[k][1]
-                    coef_arr[i] += aTCaI.get(i, j) * aMatr.get(k, j) * sigma * y
-        # polinimial coefficients have been found
-        self.polynomial.order(self.order)
-        for i in range(len(coef_arr)):
-            self.polynomial.coefficient(i, coef_arr[i])
-        # now let's calculate errors
-        if infoZeroErr == 0.0:
-            total_sigma = 0.0
-            for k in range(nPoints):
-                x = self.x_y_err_arr[k][0]
-                y = self.x_y_err_arr[k][1]
-                total_sigma += (self.polynomial.value(x) - y) ** 2
-            total_sigma = math.sqrt(total_sigma / (nPoints - 2))
-            for i in range(len(err_arr)):
-                err_arr[i] *= total_sigma
+           coeff_err_arr[i] = math.sqrt(ATA_inv.get(i,i))
+        total_sigma = 0.0
+        for k in range(nPoints):
+            x = self.x_y_err_arr[k][0]
+            y = self.x_y_err_arr[k][1]
+            total_sigma += (self.polynomial.value(x) - y) ** 2
+        degrees_of_freedom = int(abs(nPoints - (self.order + 1)))
+        if(degrees_of_freedom == 0): degrees_of_freedom = 1
+        total_sigma = math.sqrt(total_sigma / degrees_of_freedom)
+        for i in range(self.order + 1):
+            coeff_err_arr[i] *= total_sigma
         # set the resulting coefficients and errors array
-        self.coef_err_arr = [coef_arr, err_arr]
+        self.coef_err_arr = [coef_arr, coeff_err_arr]