-
Notifications
You must be signed in to change notification settings - Fork 0
/
TV5.py
370 lines (308 loc) · 14.7 KB
/
TV5.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
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
### To run this code on the AO Lab Linuxes:
## 1) Move to the directory with TV5.py in a terminal window.
## 2) type "ur_setup"
## 3) type "python"
## 4) type "from TV5 import telviz" at the python prompt
## 5a) then telviz('subdirectory','filenum') at the python prompt.
## 5b) then telviz('subdirectory','filenum',True) at the python prompt to send
## plots to .pngs.
##
## Both the subdirectory and the filenum should have single quotes around them
## to make them a string. So as it is currently in the Dropbox, you'd type
## telviz('TelemetryOpenLoop', '20140621_000545') as an example.
##
##
## To run on the lab Mac, skip step 2.
## Modified by Sarah Hale
## can email at [email protected]
# default values for subdirectory and filenum match Sarah's setup
def telviz(subdirectory='runs',filenum='20140701_214748',normalize_recon=True,normalize_dm=True,normalize_intensity=True,run_summary=False,save=False):
'''
Visualizes slopes, intensities, and newpos data in a heatmap grid, and
plots tip/tilt and pinned actuators as a function of time.
Inputs: subdirectory & filenum = strings, save = optional boolean
Outputs: 1 figure with first and last time step visualizations of each tel
file type, and 2 plots as a function of time
'''
from matplotlib import pyplot as plt
from matplotlib.widgets import Slider
from matplotlib.font_manager import FontProperties
import numpy as np
from FTR import FourierTransformReconstructor as FTRecon
from FTR.utils import circle_aperture, remove_piston, remove_tiptilt
import os.path
import imp
from kapaolibplus import (subaps_to_grid, overlay_indices, Slopes, IntensityMap, DMPositions, newpos_to_grid, overlay_indices_newpos, slope_to_grid, overlay_indices_slope)
# If we are told to run the ConfigSummary, do so
if run_summary:
# Non-sketch way does not work :(
#summary = imp.load_source('ConfigSummary', '/home/student/Pomona_v3/Config/ConfigSummary.py')
#summary.main()
# Use sketch way
execfile('/home/student/Pomona_v3/Config/ConfigSummary.py')
filenum1, filenum2 = filenum.split('_')
SLOPE_X_FILE = subdirectory + '/' + 'slope_x_' + filenum + '.tel'
SLOPE_Y_FILE = subdirectory + '/' + 'slope_y_' + filenum + '.tel'
INTENSITY_MAP = subdirectory + '/' + 'intensity_map_' + filenum + '.tel'
NEWPOS_FILE = subdirectory + '/' + 'new_pos_' + filenum + '.tel'
if os.path.isfile(SLOPE_X_FILE)==False:
SLOPE_X_FILE = subdirectory + '/' + 'slope_x_' + filenum1 + '_' + str(int(filenum2)-1) + '.tel'
if os.path.isfile(SLOPE_X_FILE) == False:
SLOPE_X_FILE = subdirectory + '/' + 'slope_x_' + filenum1 + '_' + str(int(filenum2)+1) + '.tel'
if os.path.isfile(SLOPE_Y_FILE)==False:
SLOPE_Y_FILE = subdirectory + '/' + 'slope_y_' + filenum1 + '_' + str(int(filenum2)-1) + '.tel'
if os.path.isfile(SLOPE_Y_FILE) == False:
SLOPE_Y_FILE = subdirectory + '/' + 'slope_y_' + filenum1 + '_' + str(int(filenum2)+1) + '.tel'
if os.path.isfile(INTENSITY_MAP)==False:
INTENSITY_MAP = subdirectory + '/' + 'intensity_map_' + filenum1 + '_' + str(int(filenum2)-1) + '.tel'
if os.path.isfile(INTENSITY_MAP) == False:
INTENSITY_MAP = subdirectory + '/' + 'intensity_map_' + filenum1 + '_' + str(int(filenum2)+1) + '.tel'
if os.path.isfile(NEWPOS_FILE)==False:
NEWPOS_FILE = subdirectory + '/' + 'new_pos_' + filenum1 + '_' + str(int(filenum2)-1) + '.tel'
if os.path.isfile(NEWPOS_FILE) == False:
NEWPOS_FILE = subdirectory + '/' + 'new_pos_' + filenum1 + '_' + str(int(filenum2)+1) + '.tel'
slope_x, slope_y, intensity_map = Slopes(SLOPE_X_FILE), Slopes(SLOPE_Y_FILE), IntensityMap(INTENSITY_MAP)
new_pos = DMPositions(NEWPOS_FILE)
### Helper functions here
## Copied from Sarah's modified phase5.py
shape = (11, 11)
r = 5.5
ap = circle_aperture(shape, r)
def slope_to_recon(slope_frame=slope_x.data[0]):
"""
Turn a 1D list of slope values to an 11x11 grid with zeros in empty spots
"""
grid = np.zeros((11,11))
grid[0][3:8] = slope_frame[0:5]
grid[1][2:9] = slope_frame[5:12]
grid[2][1:10] = slope_frame[12:21]
grid[3] = slope_frame[21:32]
grid[4] = slope_frame[32:43]
grid[5] = slope_frame[43:54]
grid[6] = slope_frame[54:65]
grid[7] = slope_frame[65:76]
grid[8][1:10] = slope_frame[76:85]
grid[9][2:9] = slope_frame[85:92]
grid[10][3:8] = slope_frame[92:97]
# transpose the grid
grid = grid.T # we're filling in row-by-row from the top, but numbering starts
# in the bottom left with zero and proceeds column-by-column
return grid
def recon(timestep):
'''
Reconstructs the wavefront at a given timestamp using x and y
slope data.
'''
xs = slope_to_recon(slope_x.data[timestep])
ys = slope_to_recon(slope_y.data[timestep])
recon = FTRecon(ap, filter='mod_hud', suppress_tt=True)
phi = recon(xs, ys)
return phi
def remove_ttp(data):
pr, p = remove_piston(ap, data)
tr, tx, ty = remove_tiptilt(ap, pr)
return tr
def recon2(timestep):
xs = slope_to_recon(slope_x.data[timestep])
ys = slope_to_recon(slope_y.data[timestep])
recon = FTRecon(ap, filter='mod_hud', suppress_tt=True)
phi = recon(xs, ys)
reconflat = remove_ttp(phi).flatten()
return reconflat
lenx = len(slope_x.data)
leny = len(slope_y.data)
leni = len(intensity_map.data)
lenn = len(new_pos.data)
# Go through the slopes data and recon for each timestep (and store the recon)
wave_recon = []
intensity = []
new_dm = []
max_time = len(slope_x.data) - 1 # the maximum index that exists for the time
for t in range(max_time):
wave_recon.append(recon(t))
# Go through the intensity data and new pos data and store the processed version too
intensity.append(subaps_to_grid(intensity_map.data[t]))
new_dm.append(newpos_to_grid(new_pos.data[t]))
# Find the min and max values for each heatmap plot
# do this in order to get the correct normalization for the colorbars
#print wave_recon[4][1] # prints an array
#print intensity[10][1] # prints an array
#print new_dm[10][1] # prints an array
# the wave reconstruction is an array of 2d arrays
# the intensity is an array of arrays
min_recon_val = 1000000 # big number
max_recon_val = -1000000 # small number
max_slope_x_val = -1000000
min_dm_val = 1000000
max_dm_val = -1000000
min_intensity_val = 1000000
max_intensity_val = -1000000
print 'max time', max_time
#print 'len of intensity', len(intensity)
# loop over all the times
for t in range(max_time):
if normalize_recon:
# loop over all the individual data values in the wave reconstruction
wave_max_i = len(wave_recon[0])
for i in range(wave_max_i):
wave_max_j = len(wave_recon[0][0])
for j in range(wave_max_j):
current_recon_val = wave_recon[t][i][j]
if min_recon_val > current_recon_val:
min_recon_val = current_recon_val
if max_recon_val < current_recon_val:
max_recon_val = current_recon_val
#
#
if normalize_intensity:
# loop over all the individual data values in the intensity
inten_max_i = len(intensity[0])
for i in range(inten_max_i):
inten_max_j = len(intensity[0][0])
for j in range(inten_max_j):
current_inten_val = intensity[t][i][j]
if min_intensity_val > current_inten_val:
min_intensity_val = current_inten_val
if max_intensity_val < current_inten_val:
max_intensity_val = current_inten_val
#
#
if normalize_dm:
# loop over all the individual data values in the new dm position
dm_max_i = len(new_dm[0])
for i in range(dm_max_i):
dm_max_j = len(new_dm[0][0])
for j in range(dm_max_j):
current_dm_val = new_dm[t][i][j]
if min_dm_val > current_dm_val:
min_dm_val = current_dm_val
if max_dm_val < current_dm_val:
max_dm_val = current_dm_val
#
#
#print "wave recon slice", wave_recon[0]
print "max recon val", max_recon_val
print "min recon val", min_recon_val
#print "intensity slice", intensity[0]
print "max intensity val", max_intensity_val
print "min intensity val", min_intensity_val
#print "new dm slice", new_dm[0]
print "max dm val", max_dm_val
print "min dm val", min_dm_val
# Make initial data so that the colorbars are ok
recon_init = wave_recon[0]
if normalize_recon:
recon_init[0][0] = max_recon_val
recon_init[0][1] = min_recon_val
dm_init = new_dm[0]
if normalize_dm:
dm_init[0][0] = max_dm_val
dm_init[0][1] = min_dm_val
inten_init = intensity[0]
if normalize_intensity:
inten_init[0][0] = max_intensity_val
inten_init[0][1] = min_intensity_val
# Initialize our data variables
# the variables will be updated when the time slider changes
#slope_x_data = slope_x.data[0]
#slope_y_data = slope_y.data[0]
recon_data = recon_init#wave_recon[0]
new_pos_data = dm_init#newpos_to_grid(new_pos.data[0])
intensity_map_data = inten_init#subaps_to_grid(intensity_map.data[0])
# Plots
# Summary Panel
figall = plt.figure(figsize=(10,8))
# Set font size to be smaller
plt.rcParams.update({'font.size': 10})
# Set the subplots to have spacing between them
plt.subplots_adjust(hspace=0.4, wspace=0.4)
# Plot the three 'square' plots vertically on the left side
plt.subplot2grid((3,2),(0,0))
intensity_map_im = plt.imshow(intensity_map_data, origin='lower', interpolation='none')
plt.colorbar()
plt.title('Intensity')
plt.subplot2grid((3,2),(1,0))
new_pos_im = plt.imshow(new_pos_data, origin='lower', interpolation='none')
plt.colorbar()
plt.title('DM Position')
plt.subplot2grid((3,2),(2,0))
recon_im = plt.imshow(recon_data, origin='lower', interpolation='none')
plt.colorbar()
plt.title('Wave Reconstruction')
# determine the time step between each time index
timestep = slope_x.timestamps[1] - slope_x.timestamps[0]
# Add axes for time slider
axes = figall.add_axes([0.25, 0.02, 0.5, 0.02])
timeslider = Slider(axes, 'Timestep', 0, max_time, valinit=0, valfmt='%i')
#timeslider = Slider(axes, 'Time', 0, int(slope_x.timestamps[max_time] - slope_x.timestamps[0]), valinit=0, valfmt='%i')
# TODO: Make the time slider have units of time instead of timesteps
# Currently this is not working because of the conversion from datetime object to integer does not wokr
def update(val):
# Update the data
time_index = int(val)
#time_index = int(val/timestep)
recon_data = recon(time_index)#wave_recon[time_index]
new_pos_data = newpos_to_grid(new_pos.data[time_index])
intensity_map_data = subaps_to_grid(intensity_map.data[time_index])
# Set the image array to this
recon_im.set_array(recon_data)
new_pos_im.set_array(new_pos_data)
intensity_map_im.set_array(intensity_map_data)
# Redraw the plot
figall.canvas.draw()
# Whe the slider is slid, update the plot
timeslider.on_changed(update)
# Plot the three 'rectangular' plots vertically on the right side
## RMS as a funtion of time plot
rms_x_array = np.zeros(0)
for i in xrange(len(slope_x.data)):
rec = recon2(i)
rms = np.sqrt(np.mean(np.square(rec)))
rms_x_array = np.append(rms_x_array,rms)
rms_y_array = np.zeros(0)
for i in xrange(len(slope_y.data)):
rec = recon2(i)
rms = np.sqrt(np.mean(np.square(rec)))
rms_y_array = np.append(rms_y_array,rms)
plt.subplot2grid((3,2),(0,1))
#plt.plot(range(len(rms_x_array)),rms_x_array, 'b')
#plt.plot(range(len(rms_y_array)),rms_y_array, 'g')
plt.plot(slope_x.timestamps - slope_x.timestamps[0], rms_x_array, 'b')
plt.plot(slope_y.timestamps - slope_y.timestamps[0], rms_y_array, 'g')
#plt.xlabel('Timestep')
plt.xlabel('Time (ms)')
plt.ylabel('RMS')
plt.title('RMS a Function of Time')
## Tip/tilt as a function of time plot
# pulling out the 120th and 122nd index for each time step
tt_1 = np.zeros(lenn)
tt_2 = np.zeros(lenn)
for i in range(0,lenn):
tt_1[i] = new_pos.data[i][120] # Channel 2 (Left/Right on PDVShow & Andor)
tt_2[i] = new_pos.data[i][122] # Channel 1 (Up/Down on PDVShow & Andor)
plt.subplot2grid((3,2),(1,1))
plt.plot(new_pos.timestamps - new_pos.timestamps[0],tt_1,'.', new_pos.timestamps - new_pos.timestamps[0],tt_2,'.')
plt.ylabel("Tip/tilt")
plt.title('Tip/Tilt as Function of Time')
plt.xlabel("Time (ms)")
fontP = FontProperties()
fontP.set_size('x-small')
#plt.legend(['120 (L/R? PDVShow)', '122 (U/D? PDVShow)'],'best', prop=fontP)
lgn = plt.legend(['120 (L/R?)', '122 (U/D?)'],'best', prop=fontP)
lgn.draggable()
## Pinned actuators as a function of time plot
pinned = np.zeros(lenn)
for i in range(0,lenn):
for j in range(0,120):
if new_pos.data[i][j] <= 100 or new_pos.data[i][j] >= 64900:
pinned[i] = pinned[i] + 1
plt.subplot2grid((3,2),(2,1))
plt.plot(new_pos.timestamps - new_pos.timestamps[0],pinned,'.')
plt.ylabel("Number of pinned actuators")
plt.xlabel("Time (ms)")
plt.title("Pinned Actuators as a Function of Time")
plt.suptitle('All, Unfixed ('+filenum+')',fontsize=16)
if save == True:
figall.savefig('./' + subdirectory + '/' + 'fig_all_' + filenum + '.png', dpi=300)
update(0)
plt.show()