Skip to content

Commit 6cc7421

Browse files
azukovazukov
authored
Merge pull request #71 from shishlo/main
The calculation of coefficient errors in the Polynomial fitting class
2 parents f99113c + e23176f commit 6cc7421

File tree

3 files changed

+139
-58
lines changed

3 files changed

+139
-58
lines changed

‎examples/OrbitUtils_Tests/Polynomials/polynomial_fitting_test.py‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
#----------------------------------------------
2+
# This is an example of a polymonial fitting
3+
#----------------------------------------------
4+
5+
import sys
6+
import math
7+
import random
8+
9+
from orbit.core.orbit_utils import Matrix
10+
from orbit.core.orbit_utils import Function
11+
from orbit.core.orbit_utils import SplineCH
12+
13+
from orbit.utils.fitting import PolynomialFit
14+
15+
random_err_range = 0.001
16+
f = Function()
17+
for i in range(15):
18+
x = 0.1*i
19+
y = -1+1.0*x**2+2.0*x**3 + 0.01*x**4
20+
y += random.uniform(-random_err_range,+random_err_range)
21+
f.add(x,y)
22+
23+
print ("polynomial for fitting = y = -1+1.0*x**2+2.0*x**3 + 0.01*x**4")
24+
print ("================Function====================================")
25+
poly = PolynomialFit(4)
26+
poly.fitFunction(f)
27+
[coef_arr,err_arr] = poly.getCoefficientsAndErr()
28+
for i in range(len(coef_arr)):
29+
print ("i=",i," coef = %8.5f +- %8.5f"%(coef_arr[i],err_arr[i]))
30+
31+
print ("================SplineCH====================================")
32+
spline = SplineCH()
33+
spline.compile(f)
34+
poly.fitSpline(f)
35+
[coef_arr,err_arr] = poly.getCoefficientsAndErr()
36+
for i in range(len(coef_arr)):
37+
print ("i=",i," coef = %8.5f +- %8.5f"%(coef_arr[i],err_arr[i]))
38+
39+
print ("Stop.")
40+
sys.exit(0)

‎examples/OrbitUtils_Tests/Polynomials/polynomial_test.py‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
#-------------------------------------------------------------
2+
# This is an example of Polynomial class
3+
# The methods "derivativeTo","derivative", "copyTo" and "value" are included.
4+
#-------------------------------------------------------------
5+
6+
import sys
7+
import math
8+
9+
from orbit.core.orbit_utils import Polynomial
10+
11+
poly = Polynomial(4)
12+
print ("initial order = ",poly.order())
13+
14+
poly.coefficient(2,2.0)
15+
16+
print ("========Polynomial=======")
17+
order = poly.order()
18+
for i in range(order+1):
19+
print ("i=",i," coef=",poly.coefficient(i))
20+
21+
x = 10.
22+
y = poly.value(x)
23+
print ("x=",x," y=",y)
24+
25+
poly1 = Polynomial()
26+
poly.derivativeTo(poly1)
27+
28+
print ("========Derivative polynomial=======")
29+
order = poly1.order()
30+
for i in range(order+1):
31+
print ("i=",i," coef=",poly1.coefficient(i))
32+
33+
x = 10.
34+
y = poly1.value(x)
35+
yp = poly.derivative(x)
36+
print ("x=",x," y_prime=",y," yp =",yp)
37+
38+
print ("========Copy of initial polynomial=======")
39+
poly2 = Polynomial()
40+
poly.copyTo(poly2)
41+
order = poly2.order()
42+
for i in range(order+1):
43+
print ("i=",i," coef=",poly2.coefficient(i))
44+
45+
x = 10.
46+
y = poly2.value(x)
47+
yp = poly2.derivative(x)
48+
print ("x=",x," y=",y," y_prime=",yp)
49+
50+
"""
51+
# Memory leak test
52+
count = 1
53+
while(1 < 2):
54+
poly1 = Polynomial()
55+
poly2 = Polynomial()
56+
poly.derivativeTo(poly1)
57+
poly.copyTo(poly2)
58+
count += 1
59+
if(count % 100000 == 0):
60+
print ("count=",count)
61+
"""
62+
63+
print ("Stop.")
64+
sys.exit(0)
65+

‎py/orbit/utils/fitting/PolynomialFit.py‎

Lines changed: 34 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,9 @@
77
import sys
88
import math
99

10-
from orbit.core.orbit_utils import Function, SplineCH, Matrix, Polynomial
10+
from orbit.core.orbit_utils import Function, SplineCH
11+
from orbit.core.orbit_utils import PhaseVector, Matrix
12+
from orbit.core.orbit_utils import Polynomial
1113

1214

1315
class PolynomialFit:
@@ -60,68 +62,42 @@ def fitSpline(self, spline):
6062

6163
def _makePolynomial(self):
6264
nPoints = len(self.x_y_err_arr)
63-
if nPoints < (self.order + 1):
65+
if(nPoints < (self.order + 1)):
6466
self.order = nPoints - 1
65-
# check if just one of errors is zero
66-
infoZeroErr = 1.0
67-
for [x, y, err] in self.x_y_err_arr:
68-
infoZeroErr *= err
69-
for i in range(nPoints):
70-
[x, y, err] = self.x_y_err_arr[i]
71-
sigma = 1.0
72-
if infoZeroErr != 0.0:
73-
sigma = 1.0 / (err * err)
74-
self.x_y_err_arr[i][2] = sigma
75-
# now make A matrix
76-
aMatr = Matrix(nPoints, self.order + 1)
67+
#---- now make A and A^T matrices and y-vector
68+
aMtr = Matrix(nPoints, self.order + 1)
69+
yVct = PhaseVector(nPoints)
7770
for i in range(nPoints):
71+
y = self.x_y_err_arr[i][1]
72+
yVct.set(i,y)
7873
for j in range(self.order + 1):
7974
x = self.x_y_err_arr[i][0]
80-
aMatr.set(i, j, math.pow(x, j))
81-
aTCa = Matrix(self.order + 1, self.order + 1)
82-
for i in range(self.order + 1):
83-
for j in range(self.order + 1):
84-
a = 0.0
85-
for k in range(nPoints):
86-
sigma = self.x_y_err_arr[k][2]
87-
a += aMatr.get(k, i) * sigma * aMatr.get(k, j)
88-
aTCa.set(i, j, a)
75+
aMtr.set(i, j, math.pow(x, j))
76+
aMtrT = aMtr.transpose()
8977
# now the resuting coefficients and errors
90-
aTCaI = aTCa.invert()
91-
if aTCaI == None:
92-
print("python PolynomialFit: Problem with data.")
93-
for i in range(nPoints):
94-
x = self.x_y_err_arr[i][0]
95-
y = self.x_y_err_arr[i][1]
96-
err = self.x_y_err_arr[i][2]
97-
print(" x,y,err = %12.5g %12.5g %12.5g " % (x, y, err))
98-
print("Stop.")
99-
sys.exit(1)
100-
e = aTCaI.mult(aTCa)
101-
coef_arr = [0.0] * (self.order + 1)
102-
err_arr = [0.0] * (self.order + 1)
78+
ATA_inv = (aMtrT.mult(aMtr)).invert()
79+
coeffVct = (ATA_inv.mult(aMtrT)).mult(yVct)
80+
# polinimial coefficients have been found -> coeffVct
81+
coef_arr = [0.0] * (self.order + 1)
82+
self.polynomial.order(self.order)
10383
for i in range(self.order + 1):
104-
err_arr[i] = math.sqrt(math.fabs(aTCaI.get(i, i)))
84+
coef_arr[i] = coeffVct.get(i)
85+
self.polynomial.coefficient(i,coef_arr[i])
86+
#---- here we estimate errors by deviation y_fit from y_init only
87+
#---- It means that we do not pay attention to different weights
88+
#---- of the initial Y data errors
89+
coeff_err_arr = [0.0] * (self.order + 1)
10590
for i in range(self.order + 1):
106-
coef_arr[i] = 0.0
107-
for j in range(self.order + 1):
108-
for k in range(nPoints):
109-
sigma = self.x_y_err_arr[k][2]
110-
y = self.x_y_err_arr[k][1]
111-
coef_arr[i] += aTCaI.get(i, j) * aMatr.get(k, j) * sigma * y
112-
# polinimial coefficients have been found
113-
self.polynomial.order(self.order)
114-
for i in range(len(coef_arr)):
115-
self.polynomial.coefficient(i, coef_arr[i])
116-
# now let's calculate errors
117-
if infoZeroErr == 0.0:
118-
total_sigma = 0.0
119-
for k in range(nPoints):
120-
x = self.x_y_err_arr[k][0]
121-
y = self.x_y_err_arr[k][1]
122-
total_sigma += (self.polynomial.value(x) - y) ** 2
123-
total_sigma = math.sqrt(total_sigma / (nPoints - 2))
124-
for i in range(len(err_arr)):
125-
err_arr[i] *= total_sigma
91+
coeff_err_arr[i] = math.sqrt(ATA_inv.get(i,i))
92+
total_sigma = 0.0
93+
for k in range(nPoints):
94+
x = self.x_y_err_arr[k][0]
95+
y = self.x_y_err_arr[k][1]
96+
total_sigma += (self.polynomial.value(x) - y) ** 2
97+
degrees_of_freedom = int(abs(nPoints - (self.order + 1)))
98+
if(degrees_of_freedom == 0): degrees_of_freedom = 1
99+
total_sigma = math.sqrt(total_sigma / degrees_of_freedom)
100+
for i in range(self.order + 1):
101+
coeff_err_arr[i] *= total_sigma
126102
# set the resulting coefficients and errors array
127-
self.coef_err_arr = [coef_arr, err_arr]
103+
self.coef_err_arr = [coef_arr, coeff_err_arr]

0 commit comments

Comments
 (0)