forked from gacevedobolton/myVTKPythonLibrary
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcomputePCA.py
159 lines (115 loc) · 4.22 KB
/
computePCA.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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
#coding=utf8
########################################################################
### ###
### Created by Martin Genet, 2012-2016 ###
### ###
### University of California at San Francisco (UCSF), USA ###
### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
### École Polytechnique, Palaiseau, France ###
### ###
########################################################################
import math
import numpy
import myVTKPythonLibrary as myVTK
########################################################################
def discretizeData(
farray_rr,
farray_cc,
farray_ll,
farray_h,
n_r,
n_c,
n_l,
verbose=True):
myVTK.myPrint(verbose, "*** discretizeData ***")
n_tuples = farray_ll.GetNumberOfTuples()
h = numpy.empty((n_r, n_c, n_l))
for k_l in xrange(n_l):
#print "k_l = "+str(k_l)
sel_l = [k_tuple for k_tuple in xrange(n_tuples) if (math.floor(n_l*farray_ll.GetTuple1(k_tuple)) == k_l)]
#print len(sel_l)
for k_c in xrange(n_c):
#print "k_c = "+str(k_c)
sel_c = [k_tuple for k_tuple in sel_l if (math.floor(n_c*farray_cc.GetTuple1(k_tuple)) == k_c)]
#print len(sel_c)
for k_r in xrange(n_r):
#print "k_r = "+str(k_r)
sel_r = [k_tuple for k_tuple in sel_c if (math.floor(n_r*farray_rr.GetTuple1(k_tuple)) == k_r)]
#print len(sel_r)
(m, s) = computeMeanStddevAngles([farray_h.GetTuple1(k_tuple) for k_tuple in sel_r], verbose=False)
h[k_r, k_c, k_l] = m
if (verbose >= 2): print h
return h
def computeSVD(
h,
verbose=True):
myVTK.myPrint(verbose, "*** computeSVD ***")
(n_r, n_c, n_l) = h.shape
M_r = numpy.empty((n_r, n_r))
for i in xrange(n_r):
for j in xrange(n_r):
M_r[i,j] = numpy.sum([h[i,k_c,k_l]*h[j,k_c,k_l] for k_c in xrange(n_c) for k_l in xrange(n_l)])
#print M_r
U_r, D_r, V_r = numpy.linalg.svd(M_r)
M_c = numpy.empty((n_c, n_c))
for i in xrange(n_c):
for j in xrange(n_c):
M_c[i,j] = numpy.sum([h[k_r,i,k_l]*h[k_r,i,k_l] for k_r in xrange(n_r) for k_l in xrange(n_l)])
#print M_c
U_c, D_c, V_c = numpy.linalg.svd(M_c)
M_l = numpy.empty((n_l, n_l))
for i in xrange(n_l):
for j in xrange(n_l):
M_l[i,j] = numpy.sum([h[k_r,k_c,i]*h[k_r,k_c,j] for k_r in xrange(n_r) for k_c in xrange(n_c)])
#print M_l
U_l, D_l, V_l = numpy.linalg.svd(M_l)
return U_r, D_r, V_r, U_c, D_c, V_c, U_l, D_l, V_l
def computePOD(
h,
verbose=True):
myVTK.myPrint(verbose, "*** computePOD ***")
(n_r, n_c, n_l) = h.shape
M_r = numpy.empty((n_r, n_r))
for i in xrange(n_r):
for j in xrange(n_r):
M_r[i,j] = numpy.sum([h[i,k_c,k_l]*h[j,k_c,k_l] for k_c in xrange(n_c) for k_l in xrange(n_l)])
#print M_r
w_r,v_r = numpy.linalg.eig(M_r)
#print w_r
#print v_r
w_r = w_r.real
v_r = v_r.real
#print w_r
#print v_r
idx = w_r.argsort()[::-1]
w_r = w_r[idx]
v_r = v_r[:,idx]
#print w_r
#print v_r[:,0]
M_c = numpy.empty((n_c, n_c))
for i in xrange(n_c):
for j in xrange(n_c):
M_c[i,j] = numpy.sum([h[k_r,i,k_l]*h[k_r,i,k_l] for k_r in xrange(n_r) for k_l in xrange(n_l)])
#print M_c
w_c,v_c = numpy.linalg.eig(M_c)
w_c = w_c.real
v_c = v_c.real
idx = w_c.argsort()[::-1]
w_c = w_c[idx]
v_c = v_c[:,idx]
#print w_c
#print v_c[:,0]
M_l = numpy.empty((n_l, n_l))
for i in xrange(n_l):
for j in xrange(n_l):
M_l[i,j] = numpy.sum([h[k_r,k_c,i]*h[k_r,k_c,j] for k_r in xrange(n_r) for k_c in xrange(n_c)])
#print M_l
w_l,v_l = numpy.linalg.eig(M_l)
w_l = w_l.real
v_l = v_l.real
idx = w_l.argsort()[::-1]
w_l = w_l[idx]
v_l = v_l[:,idx]
#print w_l
#print v_l[:,0]
return w_r,v_r,w_c,v_c,w_l,v_l