Skip to content

Commit 757ae68

Browse files
authored
Merge pull request #23 from shishlo/main
Matrix and Function classes changed
2 parents 4d826b4 + 41c7d36 commit 757ae68

File tree

8 files changed

+647
-62
lines changed

8 files changed

+647
-62
lines changed
Lines changed: 166 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,166 @@
1+
#------------------------------------------------------
2+
# This is an example of Function and SplineCH
3+
# They are containers for (x,y) (or (x,y,err) ) points
4+
# Function provides a linear interpolation, and SplineCH
5+
# uses 3-rd order polynomials. SplineCH can be used for
6+
# derivatives calculations.
7+
#-------------------------------------------------------
8+
9+
import sys
10+
import math
11+
12+
from orbit.core.orbit_utils import Function
13+
from orbit.core.orbit_utils import SplineCH
14+
15+
f = Function()
16+
17+
def FF(x):
18+
return math.sin(x)
19+
20+
def FFP(x):
21+
return math.cos(x)
22+
23+
#------------------------------------------
24+
# Let's create PyORBIT Function = sin(x)
25+
#------------------------------------------
26+
n = 50
27+
step = 2*math.pi/n
28+
for i in range(n):
29+
x = step*i+0.1*((1.0*i)/n)**2;
30+
y = FF(x)
31+
f.add(x,y)
32+
33+
"""
34+
#---- Test of creating of x,y,err lists
35+
(x_arr,y_arr,err_arr) = f.getXYErrLists()
36+
print ("debug x_arr =",x_arr)
37+
print ("debug y_arr =",y_arr)
38+
print ("debug err_arr=",err_arr)
39+
sys.exit(0)
40+
"""
41+
42+
"""
43+
#----- Test of initialization form PyLists
44+
x_arr = [0.,1.,2.,3.]
45+
y_arr = [1.,2.,3.,4.]
46+
err_arr = [0.1,0.2,0.3,0.4]
47+
f_test = Function()
48+
#---------------------------------------
49+
f_test.initFromLists(x_arr,y_arr)
50+
f_test.dump()
51+
#-----------------------------------------
52+
f_test.initFromLists(x_arr,y_arr,err_arr)
53+
f_test.dump()
54+
sys.exit(0)
55+
"""
56+
57+
#--------------------------------
58+
# This will print out Function
59+
#--------------------------------
60+
#f.dump()
61+
#f.dump("funcion_x_y_err.dat")
62+
63+
#--------------------------------------
64+
# This will create spline from Function
65+
#--------------------------------------
66+
spline = SplineCH()
67+
spline.compile(f)
68+
69+
#--------------------------------
70+
# This will print out Spline
71+
#--------------------------------
72+
#spline.dump()
73+
#spline.dump("spline_x_y_err.dat")
74+
75+
print ("====================================")
76+
n = 30
77+
step = 0.9*(f.getMaxX() - f.getMinX())/(n-1)
78+
y_dev_max = 0.
79+
yp_dev_max = 0.
80+
for j in range(n):
81+
x = f.getMinX() + j*step
82+
y = FF(x)
83+
yp = FFP(x)
84+
ys = spline.getY(x)
85+
yps = spline.getYP(x)
86+
ys = spline.getY(x)
87+
yps = spline.getYP(x)
88+
dy = abs(ys - y)
89+
dyp = abs(yps - yp)
90+
if(y_dev_max < dy): y_dev_max = dy
91+
if(yp_dev_max < dyp): yp_dev_max = dyp
92+
#------------------------
93+
st = " x= %+8.6f "%x
94+
st += " y = %+8.6f ys = %+8.6f "%(y,ys)
95+
st += " err = %+8.6f "%dy
96+
st += " yp = %+8.6f yps = %+8.6f "%(yp,yps)
97+
st += " err = %+8.6f "%dyp
98+
print (st)
99+
100+
print ("====================================")
101+
102+
print (" max deviation y =",y_dev_max)
103+
print (" max deviation yp =",yp_dev_max)
104+
105+
print ("Stop.")
106+
107+
#sys.exit(0)
108+
109+
#----------------------------------
110+
# Test removePoint and updatePoint
111+
#----------------------------------
112+
print ("====================================")
113+
f = Function()
114+
f.add(1.,1.)
115+
f.add(2.,2.)
116+
f.add(3.,3.)
117+
118+
f_test = Function()
119+
120+
iteration = 0
121+
122+
while(1 < 2):
123+
#f.dump()
124+
#print ("====================================")
125+
f.removePoint(0)
126+
#f.dump()
127+
#print ("====================================")
128+
f.add(1.,1.)
129+
#f.dump()
130+
#print ("====================================")
131+
f.removePoint(1)
132+
#f.dump()
133+
#print ("====================================")
134+
f.add(2.,2.)
135+
#f.dump()
136+
#print ("====================================")
137+
f.removePoint(2)
138+
#f.dump()
139+
#print ("====================================")
140+
f.add(3.,3.)
141+
#f.dump()
142+
#print ("====================================")
143+
f.updatePoint(0,1.)
144+
f.updatePoint(1,2.)
145+
f.updatePoint(2,2.)
146+
#print ("====================================")
147+
(x_arr,y_arr,err_arr) = f.getXYErrLists()
148+
#print ("====================================")
149+
f_test.initFromLists(x_arr,y_arr)
150+
#-----------------------------------------
151+
f_test.initFromLists(x_arr,y_arr,err_arr)
152+
#print ("====================================")
153+
if(iteration%100000 == 0):
154+
print ("iter=",iteration)
155+
iteration += 1
156+
157+
print ("====================================")
158+
print ("====================================")
159+
print ("====================================")
160+
161+
f.updatePoint(0,5.)
162+
f.dump()
163+
print ("====================================")
164+
f.updatePoint(2,7.,3.)
165+
f.dump()
166+
print ("====================================")
Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
#==================================
2+
# This script is a test for Matrix
3+
# 1. Determinant
4+
# 2. Transpose
5+
# 3. Inversion
6+
# calculations.
7+
#==================================
8+
9+
import sys
10+
import time
11+
12+
from orbit.core.orbit_utils import Matrix
13+
14+
def printM(m):
15+
print ("----matrix--- size=",m.size())
16+
for i in range(m.size()[0]):
17+
for j in range(m.size()[1]):
18+
print ("m(" + str(i) + "," + str(j)+")="+str(m.get(i,j)) + " ",end="")
19+
print("")
20+
21+
def setM(arr):
22+
n_row = len(arr)
23+
n_col = len(arr[0])
24+
m = Matrix(n_row,n_col)
25+
for i in range(n_row):
26+
for j in range(n_col):
27+
m.set(i,j,arr[i][j])
28+
return m
29+
30+
print ("Start.")
31+
32+
"""
33+
34+
arr = [[],[]]
35+
arr[0] = [ 1., 0.]
36+
arr[1] = [ 0., 3.]
37+
38+
arr = [[],[],[]]
39+
arr[0] = [ 1., 0., 0.]
40+
arr[1] = [ 0., 2., 0.]
41+
arr[2] = [ 0., 0., 3.]
42+
43+
arr = [[],[],[],[]]
44+
arr[0] = [ 1., 0., 0., 0.]
45+
arr[1] = [ 0., 2., 0., 0.]
46+
arr[2] = [ 0., 0., 3., 0.]
47+
arr[3] = [ 0., 0., 0., 4.]
48+
49+
"""
50+
51+
arr = [[],[],[],[],[]]
52+
arr[0] = [ 1., 1., 2., 3., 4.]
53+
arr[1] = [ 0., 2., 3., 4., 5.]
54+
arr[2] = [ 0., 0., 3., 6., 7.]
55+
arr[3] = [ 0., 0., 0., 4., 8.]
56+
arr[4] = [ 0., 0., 0., 0., 5.]
57+
58+
mtrx = setM(arr)
59+
print ("===================================")
60+
printM(mtrx)
61+
print ("===================================")
62+
print ("Det(mtrx)=",mtrx.det())
63+
print ("===========Inverted================")
64+
mtrx_invert = mtrx.invert()
65+
printM(mtrx_invert)
66+
print ("===========Caclulate M*M^-1 =======")
67+
mtrx_res = mtrx_invert.mult(mtrx)
68+
printM(mtrx_res)
69+
print ("===================================")
70+
71+
arr = [[],[],[]]
72+
arr[0] = [ 1., 1., 2., 3., 4., 5.]
73+
arr[1] = [ 0., 2., 3., 4., 5., 6.]
74+
arr[2] = [ 0., 0., 3., 6., 7., 8.]
75+
mtrx_new = setM(arr)
76+
printM(mtrx_new)
77+
print ("===============Transpose===========")
78+
mtrx_new_transp = mtrx_new.transpose()
79+
printM(mtrx_new_transp)
80+
print ("===============Initial=============")
81+
printM(mtrx_new)
82+
print ("===================================")
83+
84+
print ("Stop.")
85+
86+
#----------------------------------------------------------
87+
#---- Comment sys.exit(0) if you want to check memeory leak
88+
#---- in Matrix operations.
89+
#----------------------------------------------------------
90+
sys.exit(0)
91+
92+
count = 1
93+
start_time = time.time()
94+
while(1<2):
95+
det = mtrx.det()
96+
mtrx_invert = mtrx.invert()
97+
mtrx_res = mtrx_invert.mult(mtrx)
98+
mtrx_new_transp = mtrx_new.transpose()
99+
mtrx_invert = mtrx.invert()
100+
if(count % 100000 == 0):
101+
print ("count =",count," speed [calc/sec] = ",count/(time.time()-start_time))
102+
count += 1
103+

0 commit comments

Comments
 (0)