COMPASS  5.4.4
End-to-end AO simulation tool using GPU acceleration
tools.py
1 
37 
38 import numpy as np
39 
40 import shlex
41 from subprocess import Popen, PIPE # , call
42 from sys import stdout
43 from time import sleep
44 
45 # from mpl_toolkits.mplot3d import Axes3D
46 from matplotlib import cm
47 import matplotlib.pyplot as plt
48 
49 
50 def clr(*figs):
51  """
52  THE Fab function
53 
54  clears the current figure (no arg) or specified window
55 
56  """
57  if (figs):
58  for fig in figs:
59  # fig = fig[i]
60  plt.figure(num=fig)
61  plt.clf()
62  else:
63  plt.clf()
64 
65 
66 def system(cmd, output=False):
67  """
68  Execute the external command
69  system("ls")
70 
71 
72  out = system("ls", out=True)
73  out = system("ls -l", out=True)
74 
75 
76 
77 
78  and get its stdout exitcode and stderr.
79  """
80  args = shlex.split(cmd)
81  proc = Popen(args, stdout=PIPE, stderr=PIPE)
82  out, err = proc.communicate()
83  exitcode = proc.returncode
84  #
85  if ('\n' in out):
86  out = out.split('\n')[:-1]
87 
88  for i in range(len(out)):
89  print((out[i]))
90 
91  if (output):
92  # print("here")
93  return out, exitcode, err
94 
95 
96 def pli(data, color='gist_earth', cmin=9998, cmax=9998, win=1, origin=None,
97  aspect='equal'):
98  """
99  plots the transpose of the data
100 
101  color maps keywords can be found in
102  http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
103 
104  """
105  options = ''
106  if (cmin != 9998):
107  exec('options += ",vmin=cmin"')
108 
109  if (cmax != 9998):
110  exec('options += ",vmax=cmax"')
111 
112  if (color == b'yorick'):
113  color = 'gist_earth'
114  if (origin is None):
115  origin = ""
116  if (aspect != 'auto'):
117  aspect = "\'" + aspect + "\'"
118  else:
119  aspect = "\'auto\'"
120 
121  exec('plt.matshow(data, aspect=' + aspect + ', fignum=win, cmap=color' + options +
122  origin + ")")
123 
124 
125 def binning(w, footprint):
126 
127  # the averaging block
128  # prelocate memory
129  binned = np.zeros(w.shape[0] * w.shape[1]).reshape(w.shape[0], w.shape[1])
130  # print(w)
131  for i in range(w.shape[0]):
132  for j in range(w.shape[1]):
133  binned[i, j] = w[i, j].sum() / (footprint * footprint + 0.0)
134 
135  return binned
136 
137 
138 def minmax(tab):
139  tabVect = np.reshape(tab, tab.size)
140  return [np.min(tabVect), np.max(tabVect)]
141 
142 
143 def plg(
144  data,
145  x="",
146  win=1,
147  xlog=0,
148  ylog=0,
149  color="black",
150 ):
151  """
152 
153 
154  color = "green"
155  color = "0.71" [0-1] gray scale
156  color = '#eeefff'
157  See also:
158 
159  http://matplotlib.org/api/colors_api.html
160 
161  """
162  fig = plt.figure(win)
163  ax = fig.add_subplot(1, 1, 1)
164  try:
165  data.ndim
166  if (data.ndim > 1):
167  print(("Warning %dD dimensions. Cannot plot data. Use pli instead. " %
168  data.ndim))
169  except:
170  return
171  if (x == ""):
172  ax.plot(data, color=color)
173  else:
174  ax.plot(x, data, color=color)
175 
176  if (xlog == 1):
177  ax.set_xscale('log')
178  else:
179  ax.set_xscale('linear')
180 
181  if (ylog == 1):
182  ax.set_yscale('log')
183  else:
184  ax.set_yscale('linear')
185  fig.show()
186  return fig, ax
187 
188 
189 def zcen(data):
190  data = np.array(data)
191  if (len(data.shape) > 1):
192  print("oups zcen with dims > 1 not coded yet...")
193  return 0
194  tmp = tmp2 = []
195  for i in range(len(data) - 1):
196  tmp = (float(data[i]) + float(data[i + 1])) / 2.
197  tmp2 = np.append(tmp2, tmp)
198  return tmp2
199 
200 
201 def getValidSubapArray(nssp, rext, rint, return2d=False):
202  # The Grata case, tip-tilt sensor only.
203  if (nssp == 1):
204  return [1]
205  # to avoid some bug that eliminates useful central subapertures when
206  # obs=0.286
207  if ((nssp == 7) and (rint > 0.285 and rint < 0.29)):
208  rint = 0.285
209  print("cas particulier")
210  x = zcen(np.linspace(-1, 1, num=nssp + 1))
211  xx = []
212  for i in range(nssp):
213  xx = np.hstack((xx, x))
214  x = np.reshape(xx, (nssp, nssp))
215  y = np.transpose(x)
216  r = np.sqrt(x * x + y * y)
217  valid2dext = ((r < rext)) * 1
218  valid2dint = ((r >= rint)) * 1
219  valid2d = valid2dint * valid2dext
220 
221  if (return2d):
222  return valid2d
223  else:
224  valid = np.reshape(valid2d, [nssp * nssp])
225 
226  return valid.tolist()
227 
228 
229 """
230 def plsh(slopesvector, nssp=14, rmax=0.98, obs=0, win=1, invertxy=False):
231 
232  tmp = getValidSubapArray( nssp, rmax, obs);
233  X,Y = meshgrid(np.linspace(-1, 1, nssp), np.linspace(-1, 1, nssp))
234  vx = np.zeros([nssp*nssp])
235  vy = np.zeros([nssp*nssp])
236  hart = where(tmp)[0]
237  vx.flat[hart] = slopesvector.flat[0:len(slopesvector)/2]
238  vy.flat[hart] = slopesvector.flat[len(slopesvector)/2+1:]
239  vx = vx.reshape([nssp, nssp])
240  vy = vy.reshape([nssp, nssp])
241 
242  figure(num=win)
243  if(invertxy):
244  Q = quiver(X,Y, vy, vx)
245  else:
246  Q = quiver(X,Y, vx, vy)
247  #qk = quiverkey(Q, 0.5, 0.92, 2, r'$2 \frac{m}{s}$', labelpos='W', fontproperties={'weight': 'bold'})
248  l,r,b,t = axis()
249  dx, dy = r-l, t-b
250  axis([l-0.05*dx, r+0.05*dx, b-0.05*dy, t+0.05*dy]) # MUST DO OTHERWISE THE AUTOSCALE CAN MISS SOME ARROWS
251  #title('Minimal arguments, no kwargs')
252 """
253 
254 
255 def plpyr(slopesvector, validArray):
256  """
257  wao.config.p_wfss[0]._isvalid
258  """
259  nslopes = slopesvector.shape[0] / 2
260  x, y = np.where(validArray.T)
261  plt.quiver(x, y, slopesvector[0:nslopes], slopesvector[nslopes:])
262 
263 
264 def plsh(slopesvector, nssp, validint, sparta=False, invertxy=False, returnquiver=False):
265  """
266  <slopesvector> is the input vector of slopes
267  <nssp> is the number of subapertures in the diameter of the pupil
268  <validint> is the normalized diameter of central obscuration (between 0 and 1.00)
269  <sparta> when==1, slopes are ordered xyxyxyxy...
270  when==0, slopes are xxxxxxyyyyyyy
271  <xy> when==1, swap x and y. Does nothing special when xy==0.
272 
273  The routine plots a field vector of subaperture gradients defined in
274  vector <slopesvector>.
275  The routine automatically adjusts/finds what are the valid subapertures
276  for plotting, depending on the number of elements in <slopesvector>. Only the
277  devalidated subapertures inside the central obscuration cannot be
278  known, that’s why <validint> has to be passed in the argument list.
279 
280  """
281  nsub = slopesvector.shape[0] // 2
282  x = np.linspace(-1, 1, nssp)
283  x, y = np.meshgrid(x, x)
284  r = np.sqrt(x * x + y * y)
285  # defines outer and inner radiuses that will decide of validity of subapertures
286  # inner radius <validint> is passed as an argument.
287  # outer one will be computed so that it will match the number of
288  # subapertures in slopesvector
289  rorder = np.sort(r.reshape(nssp * nssp))
290  # number of subapertures not valid due to central obscuration
291  ncentral = nssp * nssp - np.sum(r >= validint, dtype=np.int32)
292  # determine value of external radius so that the test (validint < r < validext)
293  # leads to the correct number of subapertures
294  validext = rorder[ncentral + nsub]
295  # get the indexes of valid subapertures in the nsspxnssp map
296  valid = (r < validext) & (r >= validint)
297  ivalid = np.where(valid)
298  # feeding data <slopesvector> into <vv>
299  vx = np.zeros([nssp, nssp])
300  vy = np.zeros([nssp, nssp])
301  if (sparta is False):
302  # Canary, compass, etc.. slopes ordered xxxxxxxyyyyyyy
303  vy[ivalid] = slopesvector[0:nsub]
304  vx[ivalid] = slopesvector[nsub:]
305  else:
306  # SPARTA case, slopes ordered xyxyxyxyxyxy...
307  vx[ivalid] = slopesvector[0::2]
308  vy[ivalid] = slopesvector[1::2]
309  if (invertxy is True):
310  # swaps X and Y
311  tmp = vx
312  vx = vy
313  vy = tmp
314  if (returnquiver):
315  return x, y, vx, vy
316  else:
317  plt.quiver(x, y, vx, vy, pivot='mid')
318 
319 
320 def pl3d(im):
321  """
322  ir = pyfits.get_data("/home/fvidal/data/Run2015/June2015_27_onsky/ir/ir_2015-06-28_06h27m40s_script44_gain.fits")
323 
324  JAMAIS TESTEE !!!!!!!!!!!!!!
325 
326  """
327  X = np.arange(-5, 5, 0.25)
328  Y = np.arange(-5, 5, 0.25)
329  X, Y = np.meshgrid(X, Y)
330  Z = im
331  plt.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet)
332  plt.show()
333 
334 
335 def FFThz(signal, fe, freq=0):
336  """ PSD = FFThz( signal, fe ) OU f = FFThz( 1024, fe, freq=1 )
337  On the first form, returns the power spectral density of signal.
338  If signal has units 'u', the PSD has units 'u^2/Hz'.
339  The frequency axis can be get by using the keyword freq=1."""
340  if freq == 1:
341  n = signal.size
342  d = np.linspace(0, fe, n + 1)[0:n / 2 + 1]
343  return d[1:]
344  else:
345  n = signal.size
346  d = np.abs(np.fft.fft(signal))[0:n / 2 + 1]
347  d = d**2 / (fe * n / 2)
348  d[n / 2] /= 2
349  return d[1:]
350 
351 
352 def computePSD(zerall, fe, izerNum, wfsNum):
353  if np.isscalar(wfsNum):
354  wfsNum = [wfsNum]
355 
356  for ii in wfsNum:
357  PSD = FFThz(zerall[ii][izerNum, :], fe)
358 
359  PSD /= len(wfsNum)
360  if (len(wfsNum) > 1):
361  ff = FFThz(zerall[wfsNum][izerNum, :], fe, freq=1)
362  else:
363  ff = FFThz(zerall[wfsNum[0]][izerNum, :], fe, freq=1)
364 
365  return PSD, ff
366 
367 
368 def countExample(seconds):
369  for i in range(1, int(seconds)):
370  stdout.write("\r%d" % i)
371  stdout.flush()
372  sleep(1)
373  stdout.write("\n")
374 
375 
376 def plotSubapRectangles(pup, isvalid, istart, jstart):
377  fig = plt.matshow(pup)
378  pdiam = istart[1] - istart[0]
379  for i in istart:
380  for j in jstart:
381  if (isvalid[i // pdiam, j // pdiam]):
382  color = "green"
383  else:
384  color = "red"
385  fig.axes.add_patch(
386  plt.Rectangle((i - 0.5, j - 0.5), pdiam, pdiam, fill=False,
387  color=color))
def getValidSubapArray(nssp, rext, rint, return2d=False)
Definition: tools.py:207
def pl3d(im)
Definition: tools.py:336
def system(cmd, output=False)
Execute the external command.
Definition: tools.py:81
def binning(w, footprint)
Definition: tools.py:127
def plsh(slopesvector, nssp, validint, sparta=False, invertxy=False, returnquiver=False)
<slopesvector> is the input vector of slopes <nssp> is the number of subapertures in the diameter of ...
Definition: tools.py:288
def clr(*figs)
THE Fab function.
Definition: tools.py:56
def minmax(tab)
Definition: tools.py:140
def pli(data, color='gist_earth', cmin=9998, cmax=9998, win=1, origin=None, aspect='equal')
plots the transpose of the data
Definition: tools.py:106
def plpyr(slopesvector, validArray)
Definition: tools.py:266
def plg(data, x="", win=1, xlog=0, ylog=0, color="black")
Definition: tools.py:167
def plotSubapRectangles(pup, isvalid, istart, jstart)
Definition: tools.py:386
def FFThz(signal, fe, freq=0)
PSD = FFThz( signal, fe ) OU f = FFThz( 1024, fe, freq=1 ) On the first form, returns the power spect...
Definition: tools.py:349
def zcen(data)
Definition: tools.py:195
def computePSD(zerall, fe, izerNum, wfsNum)
Definition: tools.py:362
def countExample(seconds)
Definition: tools.py:378