COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
script_SH39m.py
1 import cProfile
2 import pstats as ps
3 import sys, os
4 sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/widgets/")
5 
6 from shesha.util import tools
7 import numpy as np
8 import carmaWrap as ch
9 import shesha as ao
10 import time
11 import matplotlib.pyplot as plt
12 import hdf5_util as h5u
13 import resDataBase as db
14 import astropy.io.fits as pf
15 import glob
16 import pandas as pd
17 
18 print("TEST SHESHA\n closed loop: call loop(int niter)")
19 simulName = "SH_39m"
20 #pathResults="/home/fvidal/compass/shesha/test/scripts/resultatsScripts/RunSH39m/"
21 #pathResults="/opt/public/fvidal/data/RunSH39m/"
22 pathResults = "/volumes/hra/micado/RunSH39m_RoundPupil/"
23 dBResult = "SH39m_RoundPupil.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  npixs = [8]
38  pixsizes = [1] # in lambda/dssp
39  gainslist = [0.3]
40  bps = [10]
41  magnitudes = [11, 12, 13, 14, 15, 16]
42  RONS = [2, 10] # noises
43  nKL_Filt = 450
44 else:
45  print("DETECTED BASH SCRIPT")
46  print(sys.argv)
47  freqs = [float(sys.argv[2])] # frequency
48  npixs = [float(sys.argv[3])] # npixs
49  pixsizes = [float(sys.argv[4])] # pixsizes
50  gainslist = [float(sys.argv[5])] # Gains
51  bps = [float(sys.argv[6])] # nb Brightests pixels
52  magnitudes = [float(sys.argv[7])] # magnitudes
53  RONS = [float(sys.argv[8])] # noises
54  nKL_Filt = float(sys.argv[9]) # nb of KL
55 #$FREQ $NPIX $PIXSIZE $GAIN $BP $MAG $KLFILT
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)
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 = {
246  "PSFFilenames": None,
247  "srir": None,
248  "lambdaTarget": None,
249  "nbBrightest": None,
250  "sr_le": None,
251  "sr_se": None,
252  "numiter": None,
253  "NklFilt": None,
254  "NklTot": None,
255  "Nkl": None,
256  "eigenvals": None,
257  "Nphotons": None,
258  "Nactu": None,
259  "RON": None,
260  "Nslopes": None
261 } # Added values computed by the simu..
262 
263 resAll = db.readDataBase(
264  fullpath=pathResults + 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 #magnitudes=[11.5,12.5,13.5,14.5]
277 
278 
279 #res500 = pf.get_data("/home/fvidal/res_500.fits")
280 #fig = plt.figure(num = 1)
281 #fig.show()
282 Nsimutot = len(gainslist) * len(magnitudes) * len(bps) * len(RONS) * len(pixsizes) * len(
283  npixs) * len(freqs)
284 NCurrSim = 0
285 for freq in freqs:
286  config.p_loop.set_ittime(1 / freq)
287  for npix in npixs:
288  config.p_wfs0.set_npix(npix)
289  for pixsize in pixsizes:
290  pxsize = pixsize * config.p_wfs0.Lambda / (
291  config.p_tel.diam / config.p_wfs0.nxsub) * 0.206265
292  config.p_wfs0.set_pixsize(pxsize)
293  for gain in gainslist:
294  config.p_controller0.set_gain(gain) # Change Gain
295  for bp in bps:
296  config.p_centroider0.set_nmax(bp)
297  for RON in RONS:
298  config.p_wfs0.set_noise(RON)
299  for magnitude in magnitudes:
300  NCurrSim += 1
301  config.p_wfs0.set_gsmag(magnitude)
302  res = pd.DataFrame(
303  columns=list(colnames.keys()) +
304  list(simunames.keys())) # Create Db for last result
305  print("Simu #%d/%d" % (NCurrSim, Nsimutot))
306  print("Freq = %3.2f Hz" % (1. / config.p_loop.ittime))
307  print("npix = %d pixels" % config.p_wfs0.npix)
308  print("nb of Brightest pixels= %d " % bp)
309  print("%3.2f'' pixel size " % config.p_wfs0.pixsize)
310  print("Magnitude = %3.2f" % config.p_wfs0.gsmag)
311  print("RON = %3.1f e-" % RON)
312  print("Gain = %3.2f" % config.p_controller0.gain)
313 
314  wfs, tel, atm, dms, tar, rtc = initSimu(config,
315  c) # Init Simu
316  nfilt = nKL_Filt
317  cmat = ao.compute_cmatWithKL(rtc, config.p_controllers[0],
318  dms, config.p_dms,
319  config.p_geom, config.p_atmos,
320  config.p_tel, nfilt)
321 
322  rtc.set_cmat(0, cmat.copy().astype(np.float32))
323 
324  SR, lambdaTargetList, sr_le, sr_se, numiter = loop(
325  config.p_loop.niter, wfs, tel, atm, dms, tar, rtc)
326  dfparams = h5u.params_dictionary(
327  config) # get the current compass config
328  dfparams.update(simunames) # Add the simunames params
329 
330  res = db.fillDf(res, dfparams) # Saving dictionnary config
331  res.loc[0, "NklFilt"] = nKL_Filt
332  res.loc[0, "Nkl"] = cmat.shape[0] - 2 - nKL_Filt
333  res.loc[0, "NklTot"] = cmat.shape[0] - 2
334  res.loc[0, "Nactu"] = cmat.shape[0]
335  res.loc[0, "Nslopes"] = cmat.shape[1]
336  res.loc[0, "Nphotons"] = config.p_wfs0._nphotons
337  res.loc[0, "RON"] = RON
338  #res.eigenvals.values[0] = rtc.getEigenvals(0)
339  res.srir.values[0] = SR # Saving computed values
340  res.lambdaTarget.values[0] = lambdaTargetList
341  res.loc[0, "gsmag"] = config.p_wfs0.gsmag
342  res.loc[0, "gain"] = config.p_controller0.gain
343  res.loc[0, "pixsize"] = config.p_wfs0.pixsize
344  res.loc[0, "npix"] = config.p_wfs0.npix
345  res.loc[0, "nbBrightest"] = bp
346  #res.sr_le.values[0] = sr_le
347  #res.sr_se.values[0] = sr_se
348  #res.numiter.values[0] = numiter
349  res.loc[0, "simulname"] = simulName
350  print("Saving PSFs...")
351  PSFNameList = []
352  for t in range(config.p_target.ntargets):
353  PSFtarget = tar.get_image(t, "le")
354  date = time.strftime("_%d-%m-%Y_%H:%M:%S_")
355  lam = "%3.2f" % tar.Lambda.tolist()[t]
356  lam = lam.replace(".", "_")
357  PSFName = "SH_" + lam + "_" + date + ".fits"
358  PSFNameList.append(PSFName)
359  #PSFNameList.append("NOT SAVED")
360  pf.writeto(pathResults + "PSFs/" + PSFName,
361  PSFtarget.copy(), clobber=True)
362  lam2 = "%3.2f" % tar.Lambda.tolist()[t]
363  res.loc[0, "SR_%s" % lam2] = SR[t]
364  filepath = pathResults + "PSFs/" + PSFName
365 
366  #"Add the SR and wavelegth value at the top of the PSF header file"
367  hdulist = pf.open(filepath) # read file
368  header = hdulist[0].header
369  header["SR"] = SR[t]
370  header["wavelength"] = tar.Lambda.tolist()[t]
371  hdulist.writeto(filepath,
372  clobber=True) # Save changes to file
373  # Adding all the parameters to the header
374  makeFITSHeader(filepath, res)
375  print("Done")
376  res.loc[0, "type_ap"] = str(res.loc[0, "type_ap"][0])
377  res.loc[0, "type"] = str(res.loc[0, "type"][0])
378  res.loc[0, "type"] = "pzt, tt"
379  res.PSFFilenames.values[0] = PSFNameList
380  resAll = db.fillDf(resAll,
381  res) # Saving res in global resAll DB
382  #resAll.to_hdf("/home/fvidal/compass/trunk/shesha/test/scripts/resultatsScripts/SH39m.h5", "resAll", complevel=9,complib='blosc')
383  resAll.to_hdf(pathResults + dBResult, "resAll", complevel=9,
384  complib='blosc')
385  #db.saveDataBase(resAll)
386 
387 print("Simulation Done...")
388 """
389 Sauver PSF dans le bon nom + directory
390  ranger... + params dans le header
391 
392 """
shesha.util
Utilities functions.
Definition: shesha/shesha/util/__init__.py:1
script_SH39m.initSimu
def initSimu(config, c)
Definition: script_SH39m.py:127
script_SH39m.loop
def loop(n, wfs, tel, atm, dms, tar, rtc)
Definition: script_SH39m.py:172
script_SH39m.makeFITSHeader
def makeFITSHeader(filepath, df)
Definition: script_SH39m.py:107