COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
closed_loop_fake_wfs.py
1 # -*- coding: utf-8 -*-
2 """
3 Created on Wed Oct 9 14:03:29 2017
4 
5 @author: sdurand
6 """
7 # import cProfile
8 # import pstats as ps
9 #@profile
10 import sys
11 import os
12 # import numpy as np
13 import carmaWrap as ch
14 import shesha as ao
15 import time
16 import matplotlib.pyplot as plt
17 import hdf5_util as h5u
18 import numpy as np
19 plt.ion()
20 sys.path.append('/home/sdurand/hracode/codes/PYRCADO/Python')
21 import PYRCADOCALIB as pyrcalib
22 from astropy.io import fits
23 from numba import autojit
24 #import gnumpy as gpu
25 
26 print("TEST SHESHA\n closed loop: call loop(int niter)")
27 
28 if (len(sys.argv) != 2):
29  error = 'command line should be:"python -i test.py parameters_filename"\n with "parameters_filename" the path to the parameters file'
30  raise StandardError(error)
31 
32 # get parameters from file
33 param_file = sys.argv[1]
34 if (param_file.split('.')[-1] == "py"):
35  filename = param_file.split('/')[-1]
36  param_path = param_file.split(filename)[0]
37  sys.path.insert(0, param_path)
38  exec("import %s as config" % filename.split(".py")[0])
39  sys.path.remove(param_path)
40 # elif (param_file.split('.')[-1] == "h5"):
41 # sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/data/par/par4bench/")
42 # import scao_16x16_8pix as config
43 # sys.path.remove(os.environ["SHESHA_ROOT"] + "/data/par/par4bench/")
44 # h5u.configFromH5(param_file, config)
45 else:
46  raise ValueError("Parameter file extension must be .py or .h5")
47 
48 print("param_file is", param_file)
49 
50 if (hasattr(config, "simul_name")):
51  if (config.simul_name is None):
52  simul_name = ""
53  else:
54  simul_name = config.simul_name
55  print("simul name is", simul_name)
56 else:
57  simul_name = ""
58 
59 clean = 1
60 matricesToLoad = {}
61 if (simul_name != ""):
62  clean = 0
63  param_dict = h5u.params_dictionary(config)
64  matricesToLoad = h5u.checkMatricesDataBase(os.environ["SHESHA_ROOT"] + "/data/",
65  config, param_dict)
66 
67 # initialisation:
68 
69 # context
70 # c = ch.carmaWrap_context(0)
71 # c = ch.carmaWrap_context(devices=np.array([0,1], dtype=np.int32))
72 # c.set_active_device(0) #useful only if you use ch.carmaWrap_context()
73 c = ch.carmaWrap_context(devices=config.p_loop.devices)
74 # wfs
75 print("->wfs")
76 wfs, tel = ao.wfs_init(config.p_wfss, config.p_atmos, config.p_tel, config.p_geom,
77  config.p_target, config.p_loop, config.p_dms)
78 
79 # atmos
80 print("->atmos")
81 atm = ao.atmos_init(c, config.p_atmos, config.p_tel, config.p_geom, config.p_loop,
82  config.p_wfss, wfs, config.p_target, clean=clean,
83  load=matricesToLoad)
84 
85 # dm
86 print("->dm")
87 dms = ao.dm_init(config.p_dms, config.p_wfss, wfs, config.p_geom, config.p_tel)
88 
89 # target
90 print("->target")
91 tar = ao.target_init(c, tel, config.p_target, config.p_atmos, config.p_geom,
92  config.p_tel, config.p_dms, config.p_wfss)
93 
94 print("->rtc")
95 # rtc
96 rtc = ao.rtc_init(tel, wfs, config.p_wfss, dms, config.p_dms, config.p_geom,
97  config.p_rtc, config.p_atmos, atm, config.p_tel, config.p_loop,
98  clean=clean, simul_name=simul_name, load=matricesToLoad)
99 
100 if not clean:
101  h5u.validDataBase(os.environ["SHESHA_ROOT"] + "/data/", matricesToLoad)
102 
103 print("====================")
104 print("init done")
105 print("====================")
106 print("objects initialzed on GPU:")
107 print("--------------------------------------------------------")
108 print(atm)
109 print(wfs)
110 print(dms)
111 print(tar)
112 print(rtc)
113 
114 
115 def import_im(nb_im, path):
116  im = fits.open(path)
117  size = im[0].data.shape[0]
118  pyr_im_cube = np.zeros((nb_im, size, size), dtype=np.float32)
119  for i in range(nb_im):
120  pyr_im_cube[i] = im[i].data
121  im.close()
122  return pyr_im_cube
123 
124 
125 def create_P(bin_factor, size):
126  return np.repeat(
127  np.identity(size / bin_factor, dtype=np.float32), bin_factor, axis=0)
128 
129 
130 def calib_pyr(centers, wfs_numbers, bin_factor=1, crop_factor=0):
131 
132  #initialisation
133  #offset 4 roi :
134  offset = np.zeros((2, 4))
135  j = [2, 1, 0, 3]
136  npup = config.p_wfss[wfs_numbers]._validsubsx.shape[0]
137  #decoupage 4 roi
138  for i in range(4):
139  #decoupage coordonnee
140  #x :
141  subx = config.p_wfss[wfs_numbers]._validsubsx[npup * (i) / 4:npup * (i + 1) / 4]
142  #y :
143  suby = config.p_wfss[wfs_numbers]._validsubsy[npup * (i) / 4:npup * (i + 1) / 4]
144  # calcul des 4 centres
145  center_compass = [((np.max(subx) - np.min(subx)) / 2.) + np.min(subx),
146  ((np.max(suby) - np.min(suby)) / 2.) + np.min(suby)]
147  # calcul des offsets
148  offset[:, i] = [
149  np.int32((centers[j[i]][0] - crop_factor / 2.) / bin_factor) -
150  center_compass[0],
151  np.int32((centers[j[i]][1] - crop_factor / 2.) / bin_factor) -
152  center_compass[1]
153  ]
154 
155  return offset
156 
157 
158 def pyr_aquisition(n=0):
159 
160  #fonction d'aquisition d'image pour la pyramide
161  #lib sesame python
162  # cam 10gbit.py
163  # ten gb class
164  # get_image(1, num_cam) -->
165  #im_path = ['pyrimgforSeb1.fits','pyrimgforSeb2.fits','pyrimgforSeb3.fits','pyrimgforSeb4.fits','pyrimgforSeb5.fits','pyrimgforSeb6.fits']
166  #im = fits.open('/home/sdurand/im_pyr_banc/'+ im_path[n])
167  #pyr_im = im[0].data
168  path = '/home/sdurand/RecordPyrImages_2017_06_06_07h49/pyrImageCube.fits'
169  im = fits.open(path)
170  pyr_im = im[n].data
171  im.close()
172  return pyr_im
173 
174 
175 def get_slope_pyrhr(npup, valid_pixel):
176 
177  pup = np.zeros((npup / 4, 4))
178  j = [0, 2, 3, 1]
179  for i in range(4):
180  pup[:, i] = valid_pixel[(npup / 4) * j[i]:(npup / 4) * (j[i] + 1)]
181  tot = np.sum(pup, axis=1)
182  t = np.average(tot)
183 
184  gx = (pup[:, 0] + pup[:, 2] - (pup[:, 1] + pup[:, 3])) / t
185  gy = (pup[:, 0] + pup[:, 1] - (pup[:, 2] + pup[:, 3])) / t
186  #gz = (pup[:,0] - pup[:,1] - pup[:,2] + pup[:,3]) / t
187 
188  slope = np.append(gx, gy) * (
189  (config.p_wfss[0].pyr_ampl * config.p_wfss[0].Lambda * 1e-6) /
190  config.p_tel.diam) * (180 / np.pi) * 3600
191 
192  return slope
193 
194 
195 def crop_im(im, taille_sortie):
196 
197  #im_crop = np.zeros((taille_sortie,taille_sortie),dtype=np.float32)
198  size = im.shape[0]
199  im_crop = im[np.int32((size / 2.) - (taille_sortie / 2.)):np.int32(
200  (size / 2.) + (taille_sortie / 2)),
201  np.int32((size / 2.) - (taille_sortie / 2.)):np.int32(
202  (size / 2.) + (taille_sortie / 2.))]
203 
204  return im_crop
205 
206 
207 @autojit
208 def binning_im(im, bin_factor):
209 
210  bin_factor = np.int32(bin_factor)
211  size = im.shape[0]
212  size_bin = size / bin_factor
213  binimage = np.zeros((size_bin, size_bin), dtype=np.float32)
214 
215  a = np.arange(size)
216  xx, yy = np.meshgrid(a, a)
217  xx = xx / bin_factor
218  yy = yy / bin_factor
219  for i in range(size):
220  for j in range(size):
221  binimage[xx[i, j], yy[i, j]] += im[i, j]
222  return binimage / (bin_factor**2)
223 
224 
225 #@autojit
226 def binning_im_2(im, bin_factor):
227  size = im.shape[0]
228  bin_factor = np.int32(bin_factor)
229  P = create_P(bin_factor, size) #
230  #GP = gpu.garray(P)
231  #Gim = gpu.garray(im)
232  binimage = ((P.T).dot(im)).dot(P)
233  #binimage = ((GP.T).dot(Gim)).dot(GP)
234 
235  return binimage / (bin_factor**2)
236 
237 
238 def loop(n, d_valid_pix=[], d_P=[], offset=[],
239  bool_fake_wfs=np.zeros(len(config.p_wfss)), bin_factor=[], crop_factor=[],
240  cube_im=[]):
241  print("----------------------------------------------------")
242  print("iter# | S.E. SR | L.E. SR | Est. Rem. | framerate")
243  print("----------------------------------------------------")
244  t0 = time.time()
245  #fake_pos = np.where(bool_fake_wfs==1)
246 
247  for i in range(n):
248  atm.move_atmos()
249  if (config.p_controllers[0].type_control == "geo"):
250  for t in range(config.p_target.ntargets):
251  tar.atmos_trace(t, atm, tel)
252  rtc.docontrol_geo(0, dms, tar, 0)
253  rtc.applycontrol(0, dms)
254  tar.dmtrace(0, dms)
255 
256  else:
257  for t in range(config.p_target.ntargets):
258  tar.atmos_trace(t, atm, tel)
259  tar.dmtrace(t, dms)
260 
261  fake_it = 0
262  for w in range(len(config.p_wfss)):
263  wfs.sensors_trace(w, "all", tel, atm, dms)
264  if bool_fake_wfs[w]: #verif fake_wfs
265  if (config.p_wfss[w].type_wfs == 'pyrhr'): # verif type_wfs = pyrhr
266  if (bin_factor[fake_it] > 1): # verif bining
267  if (cube_im == []): # verif bincube not here
268  pyr_im = pyr_aquisition(i) # aquistion image
269  pyr_im_crop = crop_im(
270  pyr_im,
271  pyr_im.shape[0] - crop_factor[w]) # crop image
272 
273  else:
274  pyr_im_crop = crop_im(cube_im[i],
275  cube_im[i].shape[0] - 2).astype(
276  np.float32) # crop image
277 
278  d_imhr = ch.carmaWrap_obj_Float2D(
279  ch.carmaWrap_context(), data=pyr_im_crop /
280  (bin_factor[fake_it]**2)) # inject pyr_image in GPU
281  d_imlr = d_P[fake_it].gemm(d_imhr, 't', 'n').gemm(
282  d_P[fake_it]) # bining GPU
283  else:
284  if (cube_im == []):
285  pyr_im = pyr_aquisition(i) # aquistion image
286  d_imlr = ch.carmaWrap_obj_Float2D(
287  ch.carmaWrap_context(),
288  data=pyr_im) # inject pyr_image in GPU
289  else:
290  d_imlr = ch.carmaWrap_obj_Float2D(
291  ch.carmaWrap_context(),
292  data=cube_im[i]) # inject pyr_image in GPU
293  # valable seulmement pour wf0 :
294  wfs.copy_pyrimg(
295  w, d_imlr, d_valid_pix[fake_it][0],
296  d_valid_pix[fake_it][1]) # envoie de l image pyramide
297 
298  elif (config.p_wfss[w].type_wfs == 'sh'): # verif type_wfs = pyrhr
299  print("TODO SH")
300  else:
301  print("error")
302  fake_it += 1 # increment for fake_wfs
303  else:
304  wfs.sensors_compimg(w) # normal wfs
305 
306  rtc.docentroids(0)
307  #slope_compass_0[:,i] = rtc.get_centroids(0)
308  rtc.docontrol(0)
309 
310  rtc.applycontrol(0, dms)
311 
312  if ((i + 1) % 100 == 0):
313  strehltmp = tar.get_strehl(0)
314  print(i + 1, "\t", strehltmp[0], "\t", strehltmp[1])
315  t1 = time.time()
316  print(" loop execution time:", t1 - t0, " (", n, "iterations), ", (t1 - t0) / n,
317  "(mean) ", n / (t1 - t0), "Hz")
318 
319 
320 #____________________________________________________________
321 # lib sesam
322 # sesam_class
323 # init vector fake_wfs -->
324 bool_fake_wfs = np.zeros(len(config.p_wfss), dtype=np.int32)
325 
326 bool_fake_wfs[0] = 1
327 
328 # init wfs
329 crop_factor = np.zeros(sum(bool_fake_wfs))
330 size_c = np.zeros(sum(bool_fake_wfs))
331 centers_fake_wfs = []
332 d_P = []
333 offset = np.zeros((2, 4, sum(bool_fake_wfs)))
334 size = np.zeros(sum(bool_fake_wfs))
335 bin_factor = np.zeros(sum(bool_fake_wfs))
336 d_valid_pix = []
337 #____________________________________________________________
338 
339 # pour le wfs_fake = 0
340 w = 0
341 # rebin param
342 bin_factor[w] = 3
343 #fake_wfs_param
344 bool_fake_wfs[w] = 1
345 size[w] = 800
346 #____________________________________________________________
347 
348 # import calibration and image
349 if (bool_fake_wfs[w] == 1):
350  if (config.p_wfss[w].type_wfs == 'pyrhr'):
351  centers_fake_wfs.append(pyrcalib.giveMeTheCalibs()[1]['centers'])
352 
353  nb_im = 100 # nombre d'image
354  path = '/home/sdurand/RecordPyrImages_2017_06_06_07h49/pyrImageCube.fits'
355  pyr_im_cube = import_im(nb_im, path)
356 #_____________________________________________________________
357 
358 # initialisation fake wfs
359 fake_pos = np.where(bool_fake_wfs == 1)
360 for f in range(sum(bool_fake_wfs)):
361  crop_factor[f] = size[f] - ((size[f] / bin_factor[f]) * bin_factor[f])
362  size_c[f] = size[f] - crop_factor[f]
363 
364  if (config.p_wfss[fake_pos[f]].type_wfs == 'pyrhr'):
365  offset[:, :, f] = calib_pyr(centers_fake_wfs[f], fake_pos[f],
366  bin_factor=bin_factor[f],
367  crop_factor=crop_factor[f]) # calcul offset for wfs
368  d_P.append(
369  ch.carmaWrap_obj_Float2D(ch.carmaWrap_context(), data=create_P(
370  bin_factor[f], size_c[f]))) # add wfs offset on GPU
371  else:
372  d_P.append([])
373 
374 #_____________________________________________________________
375 # valable seulmement pour wf0_fake :
376 w = 0
377 
378 if (bool_fake_wfs[w] == 1): # verif fake_wfs
379  if (config.p_wfss[w].type_wfs == 'pyrhr'): # verif fake_pyrhr
380  npup = config.p_wfss[w]._validsubsx.shape[0]
381  valid_pix = np.zeros((2, npup), dtype=np.int32)
382  d_P.append(
383  ch.carmaWrap_obj_Float2D(ch.carmaWrap_context(),
384  data=create_P(bin_factor[w], size_c[w])))
385  valid_pix[0, :] = np.int32(config.p_wfss[w]._validsubsx + offset[0, :, w].repeat(
386  config.p_wfss[w]._nvalid)) # cacul new X new validsubx
387  valid_pix[1, :] = np.int32(config.p_wfss[w]._validsubsy + offset[1, :, w].repeat(
388  config.p_wfss[w]._nvalid)) # cacul new Y new validsuby
389 
390  #d_valid_pix = ch.carmaWrap_obj_Float2D(ch.carmaWrap_context(), data=valid_pix)
391  d_valid_pix.append([
392  ch.carmaWrap_obj_Int1D(ch.carmaWrap_context(), data=valid_pix[0, :]),
393  ch.carmaWrap_obj_Int1D(ch.carmaWrap_context(), data=valid_pix[1, :])
394  ]) # add valid subpix coord in GPU
395  loop(100, d_valid_pix, d_P, offset=offset, bool_fake_wfs=bool_fake_wfs,
396  cube_im=pyr_im_cube, bin_factor=bin_factor,
397  crop_factor=crop_factor) # Run loop
398 #_______________________________________________________________
399  elif (config.p_wfss[w].type_wfs == 'sh'):
400 
401  print("TODO SH")
402  else:
403  print("Error")
404 else:
405  loop(100)
closed_loop_fake_wfs.import_im
def import_im(nb_im, path)
Definition: closed_loop_fake_wfs.py:115
closed_loop_fake_wfs.crop_im
def crop_im(im, taille_sortie)
Definition: closed_loop_fake_wfs.py:195
closed_loop_fake_wfs.loop
def loop(n, d_valid_pix=[], d_P=[], offset=[], bool_fake_wfs=np.zeros(len(config.p_wfss)), bin_factor=[], crop_factor=[], cube_im=[])
Definition: closed_loop_fake_wfs.py:238
closed_loop_fake_wfs.pyr_aquisition
def pyr_aquisition(n=0)
Definition: closed_loop_fake_wfs.py:158
closed_loop_fake_wfs.create_P
def create_P(bin_factor, size)
Definition: closed_loop_fake_wfs.py:125
closed_loop_fake_wfs.calib_pyr
def calib_pyr(centers, wfs_numbers, bin_factor=1, crop_factor=0)
Definition: closed_loop_fake_wfs.py:130
closed_loop_fake_wfs.get_slope_pyrhr
def get_slope_pyrhr(npup, valid_pixel)
Definition: closed_loop_fake_wfs.py:175
closed_loop_fake_wfs.binning_im
def binning_im(im, bin_factor)
Definition: closed_loop_fake_wfs.py:208
closed_loop_fake_wfs.binning_im_2
def binning_im_2(im, bin_factor)
Definition: closed_loop_fake_wfs.py:226