-
Notifications
You must be signed in to change notification settings - Fork 1
/
rotate.py
124 lines (99 loc) · 2.86 KB
/
rotate.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
import default
import kernel
import f
import J
import chi
import entropy
def readFiles(Greal, Gimag):
import os
import sys
import numpy as np
import sys
omega_n = []
G = []
G_real = []
G_imag = []
try:
ifile = open(Greal, "r")
except:
sys.exit(Greal + " does not exist. ")
for index, string in enumerate(ifile):
a = string.split()
omega_n.append(float(a[0]))
G_real.append(float(a[1]))
ifile.close()
try:
ifile = open(Gimag, "r")
except:
sys.exit(Gimag + " does not exist. ")
for index, string in enumerate(ifile):
a = string.split()
G_imag.append(float(a[1]))
ifile.close()
for i in range(len(G_real)):
G.append(G_real[i] + G_imag[i]*1j)
return omega_n, G
def printFile(x, y, fileName):
ofile = open(fileName, "w")
for i in range(len(x)):
ofile.write(str(x[i]) + " " + str(y[i]) + "\n")
ofile.close()
def isOrthogonal(matrix):
import sys
import numpy as np
t = np.shape(matrix)
if (t[0] != t[1]):
sys.exit("Matrix is not square. ")
dimension = t[0]
result = np.transpose(matrix).dot(matrix)
identity = np.zeros((dimension, dimension))
for i in range(dimension):
identity[i,i] = 1.0
diff = result - identity
norm = np.sqrt(np.trace(np.transpose(diff).dot(diff)))
return norm == 0
def main():
import os
import sys
import numpy as np
import numpy.linalg
Greal = "G_cc_real.txt"
Gimag = "G_cc_imag.txt"
omega_n, G = readFiles(Greal, Gimag)
Niom = len(omega_n)
C_real = np.zeros((Niom, Niom))
C_imag = np.zeros((Niom, Niom))
ifile = open("CM_cc_real.txt", "r")
for (index, string) in enumerate(ifile):
a = string.split()
rowIndex = int(a[0])-1
colIndex = int(a[1])-1
if (True):
C_real[rowIndex, colIndex] = float(a[2])
ifile.close()
ifile = open("CM_cc_imag.txt", "r")
for (index, string) in enumerate(ifile):
a = string.split()
rowIndex = int(a[0])-1
colIndex = int(a[1])-1
if (True):
C_imag[rowIndex, colIndex] = float(a[2])
ifile.close()
eigenvalues_real, eigenvectors_real = np.linalg.eig(C_real)
eigenvalues_imag, eigenvectors_imag = np.linalg.eig(C_imag)
O_real = numpy.zeros((Niom, Niom))
O_imag = numpy.zeros((Niom, Niom))
for i in range(Niom):
for j in range(Niom):
O_real[i, j] = eigenvectors_real[i][j]
O_imag[i, j] = eigenvectors_imag[i][j]
if (not isOrthogonal(O_real)):
print "O_real is not an orthogonal matrix. "
return -1
if (not isOrthogonal(O_imag)):
print "O_imag is not an orthogonal matrix. "
return -1
C_real_inv = numpy.linalg.inv(C_real)
C_imag_inv = numpy.linalg.inv(C_imag)
return 0
main()