COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
script_PYR39m.py
1 import cProfile
2 import pstats as ps
3 import sys, os
4 sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/widgets/")
5 sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/lib/")
6 
7 #from adoptLib import computeKLModesImat, computeCmatModal
8 from shesha.util import tools
9 import numpy as np
10 import carmaWrap as ch
11 import shesha as ao
12 import time
13 import matplotlib.pyplot as plt
14 import hdf5_util as h5u
15 import resDataBase as db
16 import astropy.io.fits as pf
17 import glob
18 import pandas as pd
19 
20 print("TEST SHESHA\n closed loop: call loop(int niter)")
21 simulName = "PYR_39m"
22 pathResults = "/volumes/hra/micado/PYR39m_RoundPupil_RUN2/"
23 dBResult = "/volumes/hra/micado/PYR_39m_RoundPupil_RUN2.h5"
24 #GPUS = np.array([0, 1, 2, 3])
25 
26 if (len(sys.argv) == 1):
27  error = 'command line should be:"python -i test.py parameters_filename"\n with "parameters_filename" the path to the parameters file'
28  raise Exception(error)
29 if (len(sys.argv) == 2):
30  print("Using Internal parameters...")
31  """
32  -----------------
33  INPUTS
34  -----------------
35  """
36  freqs = [500.]
37  gainslist = [1]
38  #magnitudes=[11, 12, 13, 14, 15, 16]
39  magnitudes = [11, 15]
40  nKL_Filt = 1000
41  MODU = [5]
42  RONS = [0.5]
43 else:
44  print("DETECTED BASH SCRIPT")
45  #python $script $PARFILE $FREQ $MODU $GAIN $MAG $KLFILT
46  print(sys.argv)
47  freqs = [float(sys.argv[2])] # frequency
48  RONS = [float(sys.argv[3])] # RONS
49  MODU = [float(sys.argv[4])] # MODU
50  gainslist = [float(sys.argv[5])] # frequency
51  magnitudes = [float(sys.argv[6])] # frequency
52  nKL_Filt = float(sys.argv[7]) # frequency
53 
54 #$FREQ $NPIX $PIXSIZE $GAIN $TH $MAG $KLFILT
55 Nsimutot = len(gainslist) * len(freqs) * len(RONS) * len(MODU) * len(magnitudes)
56 
57 if (not glob.glob(pathResults)):
58  print("Results folder not found. Creating it now:")
59  tools.system("mkdir " + pathResults)
60 if (not glob.glob(pathResults + "PSFs/")):
61  print("PSFs folder not found. Creating it now:")
62 
63  tools.system("mkdir " + pathResults + "PSFs/")
64 
65 #get parameters from file
66 param_file = sys.argv[1] # par filename
67 if (param_file.split('.')[-1] == b"py"):
68  filename = param_file.split('/')[-1]
69  param_path = param_file.split(filename)[0]
70  sys.path.insert(0, param_path)
71  exec("import %s as config" % filename.split(".py")[0])
72  #sys.path.remove(param_path)
73 elif (param_file.split('.')[-1] == b"h5"):
74  sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/data/par/par4bench/")
75  import scao_sh_16x16_8pix as config
76  #sys.path.remove(os.environ["SHESHA_ROOT"]+"/data/par/par4bench/")
77  h5u.configFromH5(param_file, config)
78 else:
79  raise ValueError("Parameter file extension must be .py or .h5")
80 
81 print("param_file is", param_file)
82 
83 if (hasattr(config, "simul_name")):
84  if (config.simul_name is None):
85  simul_name = ""
86  else:
87  simul_name = config.simul_name
88 else:
89  simul_name = ""
90 print("simul name is", simul_name)
91 
92 matricesToLoad = {}
93 if (simul_name == b""):
94  clean = 1
95 else:
96  clean = 0
97  param_dict = h5u.params_dictionary(config)
98  matricesToLoad = h5u.checkMatricesDataBase(os.environ["SHESHA_ROOT"] + "/data/",
99  config, param_dict)
100 #initialisation:
101 # context
102 c = ch.carmaWrap_context(devices=np.array([0, 1, 2, 3], dtype=np.int32))
103 
104 #c.set_active_device(6)
105 
106 
107 def makeFITSHeader(filepath, df):
108  hdulist = pf.open(filepath) # read file
109  header = hdulist[0].header
110  names = np.sort(list(set(df))).tolist()
111  for name in names:
112  val = df[name][0]
113  if (type(val) is list):
114  value = ""
115  for v in val:
116  value += (str(v) + " ")
117  elif (type(val) is np.ndarray):
118  value = ""
119  for v in val:
120  value += (str(v) + " ")
121  else:
122  value = val
123  header.set(name, value, '')
124  hdulist.writeto(filepath, clobber=True) # Save changes to file
125 
126 
127 def initSimu(config, c):
128  # wfs
129  param_dict = h5u.params_dictionary(config)
130  matricesToLoad = h5u.checkMatricesDataBase(os.environ["SHESHA_ROOT"] + "/data/",
131  config, param_dict)
132  print("->wfs")
133  wfs, tel = ao.wfs_init(config.p_wfss, config.p_atmos, config.p_tel, config.p_geom,
134  config.p_target, config.p_loop, config.p_dms)
135  # atmos
136  print("->atmos")
137  atm = ao.atmos_init(c, config.p_atmos, config.p_tel, config.p_geom, config.p_loop,
138  config.p_wfss, wfs, config.p_target, rank=0, clean=clean,
139  load=matricesToLoad)
140 
141  # dm
142  print("->dm")
143  dms = ao.dm_init(config.p_dms, config.p_wfss, wfs, config.p_geom, config.p_tel)
144 
145  # target
146  print("->target")
147  tar = ao.target_init(c, tel, config.p_target, config.p_atmos, config.p_geom,
148  config.p_tel, config.p_dms)
149 
150  print("->rtc")
151  # rtc
152  rtc = ao.rtc_init(tel, wfs, config.p_wfss, dms, config.p_dms, config.p_geom,
153  config.p_rtc, config.p_atmos, atm, config.p_tel, config.p_loop,
154  do_refslp=False, clean=clean, simul_name=simul_name,
155  load=matricesToLoad, doimat=0)
156 
157  h5u.validDataBase(os.environ["SHESHA_ROOT"] + "/data/", matricesToLoad)
158 
159  print("====================")
160  print("init done")
161  print("====================")
162  print("objects initialzed on GPU:")
163  print("--------------------------------------------------------")
164  print(atm)
165  print(wfs)
166  print(dms)
167  print(tar)
168  print(rtc)
169  return wfs, tel, atm, dms, tar, rtc
170 
171 
172 def loop(n, wfs, tel, atm, dms, tar, rtc):
173  t0 = time.time()
174  print("----------------------------------------------------")
175  print("iter# | S.E. SR | L.E. SR | Est. Rem. | framerate")
176  print("----------------------------------------------------")
177  sr_se = []
178  sr_le = []
179  numiter = []
180  for i in range(n):
181  atm.move_atmos()
182 
183  if (config.p_controllers[0].type == b"geo"):
184  for t in range(config.p_target.ntargets):
185  tar.atmos_trace(t, atm, tel)
186  rtc.do_control_geo(0, dms, tar, 0)
187  rtc.apply_control(0)
188  tar.dmtrace(0, dms)
189  else:
190  for t in range(config.p_target.ntargets):
191  tar.atmos_trace(t, atm, tel)
192  tar.dmtrace(t, dms)
193  for w in range(len(config.p_wfss)):
194  wfs.raytrace(w, "all", tel, atm, dms)
195  wfs.sensors_compimg(w)
196 
197  rtc.do_centroids(0)
198  rtc.docontrol(0)
199  rtc.apply_control(0)
200 
201  if ((i + 1) % 100 == 0):
202  print("Iter#:", i + 1)
203  #for t in range(config.p_target.ntargets):
204  t = 1
205  SR = tar.get_strehl(t)
206  print("Tar %d at %3.2fMicrons:" % (t + 1, tar.Lambda[t]))
207  signal_se = "SR S.E: %1.2f " % SR[0]
208  signal_le = "SR L.E: %1.2f " % SR[1]
209 
210  print(signal_se + signal_le)
211  #print(i+1,"\t",,SR[0],"\t",SR[1])
212  sr_le.append(SR[1])
213  sr_se.append(SR[0])
214  numiter.append(i + 1)
215 
216 
217 #
218 # plt.pause(0.01)
219 # plt.scatter(numiter, sr_le, color="green", label="Long Exposure")
220 # plt.plot(numiter, sr_le, color="green")
221 # plt.scatter(numiter, sr_se, color="red", label="Short Exposure")
222 # plt.plot(numiter, sr_se, color="red")
223 
224  t1 = time.time()
225  print(" loop execution time:", t1 - t0, " (", n, "iterations), ", (t1 - t0) / n,
226  "(mean) ", n / (t1 - t0), "Hz")
227  SRList = []
228  for t in range(config.p_target.ntargets):
229  SR = tar.get_strehl(t)
230  SRList.append(SR[1]) # Saving Long Exp SR
231  return SRList, tar.Lambda.tolist(), sr_le, sr_se, numiter
232 
233 mimg = 0. # initializing average image
234 
235 SR = []
236 """
237 dictProcess, dictplot = getDataFrameColumns()
238 colnames = h5u.params_dictionary(config)
239 resAll = pd.DataFrame( columns=colnames.keys()) # res is the local dataframe for THIS data set
240 resAll = resAll.append(colnames, ignore_index=True) #Fill dataframe
241 resAll.srir = None
242 """
243 
244 colnames = h5u.params_dictionary(config) # config values internal to compass
245 #simunames = {"PSFFilenames":None, "srir":None, "lambdaTarget":None, "threshold":None, "sr_le":None, "sr_se":None, "numiter":None, "NklFilt":None, "NklTot":None, "Nkl":None, "eigenvals":None, "Nphotons":None, "Nactu":None, "Nslopes":None}# Added values computed by the simu..
246 simunames = {
247  "PSFFilenames": None,
248  "srir": None,
249  "lambdaTarget": None,
250  "nbBrightest": None,
251  "sr_le": None,
252  "sr_se": None,
253  "numiter": None,
254  "NklFilt": None,
255  "NklTot": None,
256  "Nkl": None,
257  "eigenvals": None,
258  "Nphotons": None,
259  "Nactu": None,
260  "RON": None,
261  "Nslopes": None
262 } # Added values computed by the simu..
263 
264 resAll = db.readDataBase(fullpath=dBResult) # Reads all the database if exists
265 if (not (type(resAll) == pd.core.frame.DataFrame)):
266  print("Creating compass database")
267  resAll = db.createDf(list(colnames.keys()) + list(
268  simunames.keys())) # Creates the global compass Db
269 
270 #res = db.addcolumn(res,simunames)
271 
272 #freqs = [100.,300., 500., 1000.]
273 #npixs = [4,6,8]
274 #pixsizes = [0.5,1,1.5] # in lambda/dssp
275 #gainslist = [0.1, 0.3, 0.5]
276 #ths=[0,1,2,3]
277 #magnitudes=[11.5,12.5,13.5,14.5]
278 
279 
280 #res500 = pf.get_data("/home/fvidal/res_500.fits")
281 #fig = plt.figure(num = 1)
282 #fig.show()
283 NCurrSim = 0
284 
285 for freq in freqs:
286  config.p_loop.set_ittime(1 / freq)
287  for RON in RONS:
288  config.p_wfs0.set_noise(RON)
289  for modulation in MODU:
290  rMod = modulation
291  config.p_wfs0.set_pyr_npts(
292  int(np.ceil(int(rMod * 2 * 3.141592653589793) / 4.) * 4))
293  config.p_wfs0.set_pyr_ampl(rMod)
294  for gain in gainslist:
295  config.p_controller0.set_gain(gain) # Change Gain
296  for magnitude in magnitudes:
297  NCurrSim += 1
298  config.p_wfs0.set_gsmag(magnitude)
299  res = pd.DataFrame(
300  columns=list(colnames.keys()) + list(
301  simunames.keys())) # Create Db for last result
302  print("Freq = %3.2f Hz" % (1. / config.p_loop.ittime))
303  print("Magnitude = %3.2f" % config.p_wfs0.gsmag)
304  print("Gain = %3.2f" % config.p_controller0.gain)
305 
306  wfs, tel, atm, dms, tar, rtc = initSimu(config, c) # Init Simu
307  nfilt = nKL_Filt
308  #cmat = ao.compute_cmatWithKL(rtc, config.p_controllers[0], dms, config.p_dms, config.p_geom, config.p_atmos, config.p_tel, nfilt)
309  print("Reading cMat")
310  print(1 / 0.)
311 
312  #cmat = pf.get_data(os.environ["SHESHA_ROOT"]+"/test/scripts/cmatKLGood.fits").byteswap().newbyteorder()
313  print("Setting cMat")
314  #rtc.set_cmat(0, cmat.copy().astype(np.float32))
315  rtc.set_cmat(0, cmat.copy())
316  print("Starting Loop")
317  SR, lambdaTargetList, sr_le, sr_se, numiter = loop(
318  config.p_loop.niter, wfs, tel, atm, dms, tar, rtc)
319  dfparams = h5u.params_dictionary(
320  config) # get the current compass config
321  dfparams.update(simunames) # Add the simunames params
322 
323  res = db.fillDf(res, dfparams) # Saving dictionnary config
324  res.loc[0, "NklFilt"] = nKL_Filt
325  res.loc[0, "Nkl"] = cmat.shape[0] - 2 - nKL_Filt
326  res.loc[0, "NklTot"] = cmat.shape[0] - 2
327  res.loc[0, "Nactu"] = cmat.shape[0]
328  res.loc[0, "Nslopes"] = cmat.shape[1]
329  res.loc[0, "Nphotons"] = config.p_wfs0._nphotons
330  res.loc[0, "RON"] = RON
331  res.loc[0, "Nphotons"] = config.p_wfs0._nphotons
332  #res.eigenvals.values[0] = rtc.getEigenvals(0)
333  res.srir.values[0] = SR # Saving computed values
334  res.lambdaTarget.values[0] = lambdaTargetList
335  res.loc[0, "gsmag"] = config.p_wfs0.gsmag
336  res.loc[0, "gain"] = config.p_controller0.gain
337  #res.sr_le.values[0] = sr_le
338  #res.sr_se.values[0] = sr_se
339  #res.numiter.values[0] = numiter
340  res.loc[0, "simulname"] = simulName
341  print("Saving PSFs...")
342  PSFNameList = []
343  for t in range(config.p_target.ntargets):
344  PSFtarget = tar.get_image(t, "le")
345  date = time.strftime("_%d-%m-%Y_%H:%M:%S_")
346  lam = "%3.2f" % tar.Lambda.tolist()[t]
347  lam = lam.replace(".", "_")
348  PSFName = "PYR_" + lam + "_" + date + ".fits"
349  PSFNameList.append(PSFName)
350  #PSFNameList.append("NOT SAVED")
351  pf.writeto(pathResults + "PSFs/" + PSFName, PSFtarget.copy(),
352  clobber=True)
353  lam2 = "%3.2f" % tar.Lambda.tolist()[t]
354  res.loc[0, "SR_%s" % lam2] = SR[t]
355  filepath = pathResults + "PSFs/" + PSFName
356 
357  #"Add the SR and wavelegth value at the top of the PSF header file"
358  hdulist = pf.open(filepath) # read file
359  header = hdulist[0].header
360  header["SR"] = SR[t]
361  header["wavelength"] = tar.Lambda.tolist()[t]
362  hdulist.writeto(filepath, clobber=True) # Save changes to file
363  # Adding all the parameters to the header
364  makeFITSHeader(filepath, res)
365  print("Done")
366  res.loc[0, "type_ap"] = str(res.loc[0, "type_ap"][0])
367  res.loc[0, "type"] = str(res.loc[0, "type"][0])
368  res.loc[0, "type"] = "pzt, tt"
369  res.loc[0, "npix"] = res.loc[0, "npix"][0]
370  #res.loc[0, "nbBrightest"] = res.loc[0, "nbBrightest"][0]
371  res.loc[0, "pixsize"] = res.loc[0, "pixsize"][0]
372  res.PSFFilenames.values[0] = PSFNameList
373  resAll = db.fillDf(resAll, res) # Saving in global DB
374  #resAll.to_hdf("/home/fvidal/compass/trunk/shesha/test/scripts/resultatsScripts/SH39m.h5", "resAll", complevel=9,complib='blosc')
375  resAll.to_hdf(dBResult, "resAll", complevel=9, complib='blosc')
376  #db.saveDataBase(resAll)
377 
378 print("Simulation Done...")
379 """
380 Sauver PSF dans le bon nom + directory
381  ranger... + params dans le header
382 
383  SR = np.zeros((3, len(set(resAll.gsmag))))
384  i=0
385  for mag in list(set(resAll.gsmag)):
386  SR[0,i] = resAll[resAll.gsmag == mag]["SR_1.20"].max()
387  SR[1,i] = resAll[resAll.gsmag == mag]["SR_1.65"].max()
388  SR[2,i] = resAll[resAll.gsmag == mag]["SR_2.20"].max()
389  i+=1
390 
391 
392 """
script_PYR39m.makeFITSHeader
def makeFITSHeader(filepath, df)
Definition: script_PYR39m.py:107
shesha.util
Utilities functions.
Definition: shesha/shesha/util/__init__.py:1
script_PYR39m.loop
def loop(n, wfs, tel, atm, dms, tar, rtc)
Definition: script_PYR39m.py:172
script_PYR39m.initSimu
def initSimu(config, c)
Definition: script_PYR39m.py:127