-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdiagonalShifter.py
328 lines (317 loc) · 17.4 KB
/
diagonalShifter.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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 19 02:45:56 2016
@author: Oddman
"""
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal
def shiftMat(sourceMat,angles=None,savgol=False):
'''This function takes a 2d array / matrix as its argument ("SourceMat")
and produces a list of matrices where the rows (shape[0]) are shifted
diagonally. The items in the list are pairs of the square matrix turned
into parallellograms with the triangles at the end cut off, like so:
______________ ______________ ______________
| | ( /! ! / \ ! !\ )
| | -> ( / ! ! / \ ! ! \ )
|_____________| ( /__!__________!/ , \!__________!__\ )
The index in the list of outputs is a function of the various possible
angles of shifting. The optional argument "angles" specifies the number of
angle-subdivisions from vertical to 45°. When None, that number is equal
to the number of rows in the source matrix.
'''
## Identify the number of rows in the source matrix.
rowNum = sourceMat.shape[0]
outList = []
## Assert that the source matrix is at least four times as long as it is
## wide - owing to the total loss (2 * number of rows), one really does not
## want less than that.
## (When applying to stripe plots, remember that rings is on the number of
## rows and that time is on columns, so this should not be a problem.
## If it is, run longer simulations.)
assert rowNum * 4 < sourceMat.shape[1], 'Source matrix too tall'
if angles == None:
## angles not specified? Use the number of rows.
angles = rowNum
else:
## Divisor must be an integer. (We need iteration steps.)
assert type(angles) == int, '"angles" must be an integer.'
## One more angle; the 0th angle is the base array, only truncated.
angleIndices = range(angles+1)
## Define the target shape.
targetShape = list(sourceMat.shape)
targetShape[1] = targetShape[1] - (2 * (rowNum - 1))
print 'target shape: {0}'.format(targetShape)
## Define what the expanded array will look like.
expandedTargetShape = targetShape[:]
expandedTargetShape[1] *= angles
## To do the slanting, we chop the cells of the source array into shares;
## every cell into the number of angles. This basically expands the array,
## making its rows angles times as long.
##
## Create the expanded array, which will be reused for the various
## angle-options:
sourceMatSplit = np.split(sourceMat,sourceMat.shape[1],axis=1)
concat = []
[concat.extend([mat/float(angles) for i in range(angles)]) for mat in sourceMatSplit]
expandedMat = np.concatenate(concat,axis=1)
if savgol == True:
savgolWindowLength = angles*2
savgolWindowLength += 1
print 'herpderp'
expandedMat = signal.savgol_filter(expandedMat,window_length = savgolWindowLength, polyorder = 3, axis = 1)
plt.imshow(expandedMat,interpolation='nearest',cmap=plt.cm.gray)
for angleIndex in angleIndices:
## For every possible angle, plus the case of no shift at all, do:
## create two matrices; left and right-slanting
rightSlantExpanded = np.zeros(expandedTargetShape,dtype=expandedMat.dtype)
leftSlantExpanded = np.zeros(expandedTargetShape,dtype=expandedMat.dtype)
## Then, for every row in the matrix, calculate the to-and-from
## indices for the expanded matrix to result in the correctly slanted
## output:
for rowIndex in range(rowNum):
shiftRightAmount = angleIndex*rowIndex
shiftLeftAmount = (angles*(rowNum-1))-(angleIndex*rowIndex)
## defensive check to see whether the size is correct (numpy won't
## give an error if the size is incorrect, it will 'broadcast')
assert np.shape(expandedMat[
rowIndex,
shiftRightAmount
:
expandedTargetShape[1]+shiftRightAmount
]) == np.shape(rightSlantExpanded[rowIndex,:])
rightSlantExpanded[rowIndex,:] = expandedMat[
rowIndex,
shiftRightAmount
:
expandedTargetShape[1]+shiftRightAmount
]
## Same defensive business for the right-hand case
assert np.shape(expandedMat[
rowIndex,
shiftLeftAmount
:
expandedTargetShape[1]+shiftLeftAmount
]) == np.shape(leftSlantExpanded[rowIndex,:])
leftSlantExpanded[rowIndex,:] = expandedMat[
rowIndex,
shiftLeftAmount
:
expandedTargetShape[1]+shiftLeftAmount
]
## The cutting has been done. Now, to condense the matrices.
leftSlant = np.zeros(shape=targetShape,dtype=expandedMat.dtype)
rightSlant = np.zeros(shape=targetShape,dtype=expandedMat.dtype)
for colNum in range(expandedTargetShape[1]):
## Add the expanded rows to the condensed matrix one by one
targetColNum = colNum // angles
leftSlant[:,targetColNum] += leftSlantExpanded[:,colNum]
rightSlant[:,targetColNum] += rightSlantExpanded[:,colNum]
## For every angle, put in a right-slanted and left-slanted matrix.
outList.append([rightSlant,leftSlant])
return outList
def angleSignificance(listOfAnglePairs,product=False):
angles = len(listOfAnglePairs)
## some defensive shape checking
## (should be no problem if shiftMat output is used)
baseShape = listOfAnglePairs[0][0].shape
base_dtype = listOfAnglePairs[0][0].dtype
for pair in listOfAnglePairs:
for mat in pair:
assert mat.shape == baseShape
allSamples = []
angleMatShape = np.copy(baseShape)
angleMatShape[1] *= 2
totalVec = np.zeros(shape = angleMatShape[1]*angles, dtype = base_dtype)
for i,pair in enumerate(listOfAnglePairs):
totalPair = np.concatenate(pair,axis=1)
multiMat = np.ones(totalPair.shape,dtype=np.float64)
multiMat = multiMat * np.mean(totalPair)
normed = totalPair / multiMat
normed -= np.mean(normed)
normed *= normed
# angleVec = np.zeros(shape = angleVecLen, dtype = base_dtype)
## For ease of computation, stick the opposed angles together in a
## single vector.
# angleVec[0:baseShape[1]] = np.sum(pair[0],axis = 0)
# angleVec[baseShape[1]:] = np.sum(pair[1],axis = 0)
if product:
angleVec = normed + .1
angleVec = np.product(angleVec)
else:
angleVec = np.sum(normed,axis=0)
totalVec[i*angleMatShape[1]:i*angleMatShape[1]+angleMatShape[1]] = angleVec
allSamples.append(angleVec)
overallMean = np.mean(totalVec)
# print totalVec
# print overallMean
BesselCorrectedNumerator = angleMatShape[1] - 1
sdPerAngle = []
for sample in allSamples:
deviance = (sample - np.mean(sample)) ** 2
# deviance = sample
sd = np.sqrt(np.sum(deviance) / float(BesselCorrectedNumerator))
sdPerAngle.append(sd)
# sdPerAngle.append(np.sum(deviance))
# thisAngleSummedPerFrame = [np.sum(mat,axis = 0) for mat in pair]
# mean =
return sdPerAngle
#def angleSquaredErrorPerRow(listOfAnglePairs):
# angles = len(listOfAnglePairs)
# ## some defensive shape checking
# ## (should be no problem if shiftMat output is used)
# baseShape = listOfAnglePairs[0][0].shape
# base_dtype = listOfAnglePairs[0][0].dtype
# for pair in listOfAnglePairs:
# for mat in pair:
# assert mat.shape == baseShape
# ## Set up target data containers
# allSamples = []
# angleMatShape = np.copy(baseShape)
# angleMatShape[1] *= 2
# totalVec = np.zeros(shape = angleMatShape[1]*angles, dtype = base_dtype)
# ## Fill the data containers
# for i,pair in enumerate(listOfAnglePairs):
#
# totalPair = np.concatenate(pair,axis=1)
# multiMat = np.ones(totalPair.shape,dtype=np.float64)
# multiMat = multiMat * np.mean(totalPair)
# normed = totalPair / multiMat
## prodded = np.product(normed,axis=0)
#
#
## print i
## angleVec = np.zeros(shape = angleVecLen, dtype = base_dtype)
# ## For ease of computation, stick the opposed angles together in a
# ## single vector.
## angleVec[0:baseShape[1]] = np.sum(pair[0],axis = 0)
## angleVec[baseShape[1]:] = np.sum(pair[1],axis = 0)
# totalVec[i*angleMatShape[1]:i*angleMatShape[1]+angleMatShape[1]] = normed
# allSamples.append(normed)
# overallMean = np.mean(totalVec)
## print totalVec
## print overallMean
# BesselCorrectedNumerator = angleMatShape[1] - 1
# sdPerAngle = []
# for sample in allSamples:
# deviance = (sample - np.mean(sample)) ** 2
# sdPerAngle.append(np.sqrt(np.sum(deviance) / float(BesselCorrectedNumerator)))
#
## thisAngleSummedPerFrame = [np.sum(mat,axis = 0) for mat in pair]
## mean =
# return sdPerAngle
#
#def angleSignificanceProduct(listOfAnglePairs):
# angles = len(listOfAnglePairs)
# ## some defensive shape checking
# ## (should be no problem if shiftMat output is used)
# baseShape = listOfAnglePairs[0][0].shape
# base_dtype = listOfAnglePairs[0][0].dtype
# for pair in listOfAnglePairs:
# for mat in pair:
# assert mat.shape == baseShape
# allSamples = []
# angleMatShape = np.copy(baseShape)
# angleMatShape[1] *= 2
# totalVec = np.zeros(shape = angleMatShape[1]*angles, dtype = base_dtype)
# for i,pair in enumerate(listOfAnglePairs):
## print pair
# totalPair = np.concatenate(pair,axis=1)
## print totalPair
## print np.max(totalPair)
# multiMat = np.ones(totalPair.shape,dtype=np.float64)
# multiMat = multiMat * np.mean(totalPair)
## print multiMat
# normed = totalPair / multiMat
# normed += .1
# prodded = np.product(normed,axis=0)
#
## print prodded
# totalVec[i*angleMatShape[1]:i*angleMatShape[1]+angleMatShape[1]] = prodded
# allSamples.append(prodded)
# for i,pair in enumerate(listOfAnglePairs):
## print i
# angleVec = np.zeros(shape = angleVecLen, dtype = base_dtype)
# ## For ease of computation, stick the opposed angles together in a
# ## single vector.
# angleVec[0:baseShape[1]] = np.sum(pair[0],axis = 0)
# angleVec[baseShape[1]:] = np.sum(pair[1],axis = 0)
# totalVec[i*angleVecLen:i*angleVecLen+angleVecLen] = angleVec
# allSamples.append(angleVec)
overallMean = np.mean(totalVec)
# print totalVec
# print overallMean
BesselCorrectedNumerator = angleMatShape[1] - 1
sdPerAngle = []
for sample in allSamples:
deviance = (sample - np.mean(sample)) ** 2
sdPerAngle.append(np.sqrt(np.sum(deviance) / float(BesselCorrectedNumerator)))
# thisAngleSummedPerFrame = [np.sum(mat,axis = 0) for mat in pair]
# mean =
return sdPerAngle
if __name__ == '__main__':
#==============================================================================
# Only runs when called as script. This is essentially unit testing code.
#==============================================================================
## Create a sample array.
testMat1 = np.array([[0,1,2,1],[1,2,1,0],[2,1,0,1],[1,0,1,2],[0,1,2,1],[1,2,1,0],[2,1,0,1],[1,0,1,2],[0,1,2,1],[1,2,1,0],[2,1,0,1],[1,0,1,2],[0,1,2,1],[1,2,1,0],[2,1,0,1],[1,0,1,2],[0,1,2,1],[1,2,1,0]],dtype=np.float64)
testMat = np.array([[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1],[1,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],[0,0,0,1,0,0,0,0,0,0],[0,0,0,0,1,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0],[0,0,0,0,0,0,1,0,0,0],[0,0,0,0,0,0,0,1,0,0],[0,0,0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,0,0,1]],dtype=np.float64)
testMat = np.transpose(testMat)
#==============================================================================
# Case 1: no Savitzky-Golay filtering.
#==============================================================================
## Apply the shifting.
shiftedMats = shiftMat(testMat,5)
mats = [[testMat,testMat]] + shiftedMats
matNo = len(mats)
## Display the result.
stripeFig = plt.figure(figsize = (8,8))
for i,mat in enumerate(mats):
mat = mats[i]
# ax1 = stripeFig.add_subplot('{0}2{1}'.format(matNo,i*2+1))
# ax2 = stripeFig.add_subplot('{0}2{1}'.format(matNo,i*2+2))
ax1 = stripeFig.add_subplot(matNo,2,i*2+1)
ax2 = stripeFig.add_subplot(matNo,2,i*2+2)
ax1.imshow(mat[0],interpolation='nearest',cmap=plt.cm.gray)
ax2.imshow(mat[1],interpolation='nearest',cmap=plt.cm.gray)
stripeFig.show()
## Apply the variance analysis.
print angleSignificance(shiftedMats,product=False)
# shiftMat(testMat,savgol=True)
#==============================================================================
# Case 2: Savitzky-Golay filtering on expanded matrix.
#==============================================================================
## Apply the shifting.
shiftedMats = shiftMat(testMat,5,savgol=True)
mats = [[testMat,testMat]] + shiftedMats
matNo = len(mats)
## Display the result.
stripeFig2 = plt.figure(figsize = (8,8))
for i,mat in enumerate(mats):
mat = mats[i]
# ax1 = stripeFig.add_subplot('{0}2{1}'.format(matNo,i*2+1))
# ax2 = stripeFig.add_subplot('{0}2{1}'.format(matNo,i*2+2))
ax1 = stripeFig2.add_subplot(matNo,2,i*2+1)
ax2 = stripeFig2.add_subplot(matNo,2,i*2+2)
ax1.imshow(mat[0],interpolation='nearest',cmap=plt.cm.gray)
ax2.imshow(mat[1],interpolation='nearest',cmap=plt.cm.gray)
stripeFig2.show()
#==============================================================================
# Case 2: Savitzky-Golay filtering on input.
#==============================================================================
## Apply the shifting.
testMat = signal.savgol_filter(testMat,window_length = 5, polyorder = 3, axis = 1)
shiftedMats = shiftMat(testMat,5,savgol=False)
mats = [[testMat,testMat]] + shiftedMats
matNo = len(mats)
## Display the result.
stripeFig3 = plt.figure(figsize = (8,8))
for i,mat in enumerate(mats):
mat = mats[i]
# ax1 = stripeFig.add_subplot('{0}2{1}'.format(matNo,i*2+1))
# ax2 = stripeFig.add_subplot('{0}2{1}'.format(matNo,i*2+2))
ax1 = stripeFig3.add_subplot(matNo,2,i*2+1)
ax2 = stripeFig3.add_subplot(matNo,2,i*2+2)
ax1.imshow(mat[0],interpolation='nearest',cmap=plt.cm.gray)
ax2.imshow(mat[1],interpolation='nearest',cmap=plt.cm.gray)
stripeFig3.show()