COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
script_PYR39m_optimGain.py
1 #
2 """
3 Test Line
4 
5 ipython -i /home/fvidal/compass/shesha/test/scripts/script_PYR39m_optimGain.py /home/fvidal/compass/shesha/data/par/MICADO/micado_39m_PYR_ELTPupil.py 500 0.1 5 0.5 17 100 py3Test 92 5
6 """
7 
8 
9 def sendMail(message, title):
10  import smtplib
11  from email.mime.text import MIMEText
12  smtp = smtplib.SMTP('smtp.obspm.fr')
13  msg = MIMEText(title)
14  msg['From'] = 'micmac'
15  msg['To'] = 'Script micmac'
16  msg['Subject'] = message
17  smtp.sendmail('micmac@obspm.fr', ["fabrice.vidal@obspm.fr"], msg.as_string())
18 
19 
20 #try:
21 import cProfile
22 import pstats as ps
23 import sys, os
24 sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/widgets/")
25 sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/lib/")
26 sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/AOlib/")
27 sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/src/shesha_util/")
28 
29 #from adoptLib import computeKLModesImat, computeCmatModal
30 from shesha.util import tools
31 import numpy as np
32 import carmaWrap as ch
33 import shesha.config as ao
34 import shesha.sim
35 import shesha.constants as scons
36 from shesha.constants import CONST
37 
38 import time
39 import matplotlib.pyplot as plt
40 import hdf5_util as h5u
41 import resDataBase as db
42 import astropy.io.fits as pf
43 import glob
44 import pandas as pd
45 import compassConfigToFile as cf
46 import make_pupil as mkP
47 
48 if (len(sys.argv) == 1):
49  error = 'command line should be:"python -i test.py parameters_filename"\n with "parameters_filename" the path to the parameters file'
50  raise Exception(error)
51 elif (len(sys.argv) == 2):
52  print("Using Internal parameters...")
53  """
54  -----------------
55  INPUTS
56  -----------------
57  """
58  freq = 500
59  gain = 1
60  magnitude = 11
61  nKL_Filt = 450
62  MODU = 5
63  RON = 0.1
64  NSSP = 92
65  simulName = "PYR_39m_RoundPupil_FromHippo6"
66 else:
67  print("-------------------------------------")
68  print("DETECTED BASH SCRIPT with parameters:")
69  print((sys.argv))
70  print("-------------------------------------")
71 
72  freq = float(sys.argv[2]) # AO Loop frequency
73  RON = float(sys.argv[3]) # noise on the WFS measurement in electrons
74  MODU = float(sys.argv[4]) # modulation radius
75  gain = float(sys.argv[5]) # global loop gain
76  magnitude = float(sys.argv[6]) # gs magnitude
77  nKL_Filt = int(float(sys.argv[7])) # Nb KL filtered
78  simulName = sys.argv[8] # Nb KL filtered
79  NSSP = int(sys.argv[9]) # Number of ssp (pixels of pyramid)
80  GPU = int(sys.argv[10]) # GPU number
81  comment = "SRVsGSVsNControlledModes"
82 
83 print(("Freq=", freq))
84 print(("RON=", RON))
85 print(("MODU=", MODU))
86 print(("gain=", gain))
87 print(("magnitude=", magnitude))
88 print(("nKL_Filt=", nKL_Filt))
89 print(("GPU=", GPU))
90 print(("simulName=", simulName))
91 
92 pathResults = "/volumes/hra/micado/" + simulName
93 
94 dBResult = pathResults + "/" + simulName + ".h5"
95 savePSFs = True
96 PYR = True
97 
98 imatFromFile = False
99 #iMatName = "iMat39mPYR_MODU_"+str(int(MODU))+".fits"
100 #KL2VName = "KL2VNorm39mPYR_MODU_"+str(int(MODU))+".fits"
101 #gainModalName = "gains4K_MODU_"+str(int(MODU))+".fits"
102 
103 #iMatName = "iMat_MODU_5_ELTPUPIL.fits"
104 #KL2VName = "KL2VNorm_MODU_5_ELTPUPIL.fits"
105 #gainModalName = "gains4K_MODU_5_ELTPUPIL.fits"
106 
107 imat0_PATH = "/home/fvidal/dataSimus"
108 iMatName = "imatDiffraction_ELTPYR_35Layers.fits"
109 gainModalName = "gains4K_ELTPYR_35Layers.fits"
110 KL2VName = "KL2VNorm_ELTPYR_35Layers.fits"
111 """
112 iMatName = "iMat_MODU_2_ELTPUPIL.fits"
113 KL2VName = "KL2VNorm_MODU_2_ELTPUPIL.fits"
114 gainModalName = "gains4K_MODU_2_ELTPUPIL.fits"
115 """
116 
117 ModalBasisType = "Btt"
118 
119 PSFWithOtherPupil = True
120 
121 niter = 8096
122 saveCBData = True
123 nbLoopData = 512
124 """
125 simulName = "PYR_39m"
126 pathResults="/home/fvidal/dataSimus/PYR_39m_RoundPupil_RUN1/"
127 dBResult = "/home/fvidal/dataSimus/PYR_39m_RoundPupil_RUN1.h5"
128 imat0_PATH = "/home/fvidal/compass/shesha/test/scripts"
129 savePSFs = False
130 imatFromFile = False
131 """
132 
133 if (GPU == 1):
134  GPUs = np.array([4, 5, 6, 7], dtype=np.int32)
135 else:
136  GPUs = np.array([GPU], dtype=np.int32)
137 print(("Using GPUs: ", GPUs))
138 
139 GPUs = np.array([4, 5, 6, 7], dtype=np.int32)
140 #GPUs = np.array([0,1,2,3], dtype=np.int32)
141 
142 if (not glob.glob(pathResults)):
143  print("Results folder not found. Creating it now:")
144  tools.system("mkdir " + pathResults)
145 if (not glob.glob(pathResults + "/PSFs/")):
146  print("PSFs folder not found. Creating it now:")
147  tools.system("mkdir " + pathResults + "/PSFs/")
148 if (not glob.glob(pathResults + "/AODATA/")):
149  print("AODATA folder not found. Creating it now:")
150  tools.system("mkdir " + pathResults + "/AODATA/")
151 if (not glob.glob(pathResults + "/CircularBuffers/")):
152  print("CircularBuffers folder not found. Creating it now:")
153  tools.system("mkdir " + pathResults + "/CircularBuffers/")
154 
155 #get parameters from file
156 param_file = sys.argv[1] # par filename
157 if (param_file.split('.')[-1] == "py"):
158  filename = param_file.split('/')[-1]
159  param_path = param_file.split(filename)[0]
160  sys.path.insert(0, param_path)
161  exec("import %s as config" % filename.split(".py")[0])
162  #sys.path.remove(param_path)
163 elif (param_file.split('.')[-1] == "h5"):
164  sys.path.insert(0, os.environ["SHESHA_ROOT"] + "/data/par/par4bench/")
165  import scao_sh_16x16_8pix as config
166  #sys.path.remove(os.environ["SHESHA_ROOT"]+"/data/par/par4bench/")
167  h5u.configFromH5(param_file, config)
168 else:
169  raise ValueError("Parameter file extension must be .py or .h5")
170 
171 print(("param_file is", param_file))
172 
173 if (hasattr(config, "simul_name")):
174  if (config.simul_name is None):
175  simul_name = ""
176  else:
177  simul_name = config.simul_name
178 else:
179  simul_name = ""
180 print(("simul name is", simul_name))
181 
182 matricesToLoad = {}
183 if (simul_name == ""):
184  clean = 1
185 else:
186  clean = 0
187  param_dict = h5u.params_dictionary(config)
188  matricesToLoad = h5u.checkMatricesDataBase(os.environ["SHESHA_ROOT"] + "/data/",
189  config, param_dict)
190 c = ch.carmaWrap_context(devices=GPUs)
191 
192 
193 class wao_class():
194 
195  def __init__(self, config, wfs, tel, atm, dms, tar, rtc):
196  self.config = config
197  self.wfs = wfs
198  self.tel = tel
199  self.atm = atm
200  self.dms = dms
201  self.tar = tar
202  self.rtc = rtc
203 
204 
205 def makeFITSHeader(filepath, df):
206  hdulist = pf.open(filepath) # read file
207  header = hdulist[0].header
208  names = np.sort(list(set(df))).tolist()
209 
210  for name in names:
211  if (name != "centroider.thresh"):
212  val = df[name][0]
213  if (type(val) is list):
214  value = ""
215  for v in val:
216  value += (str(v) + " ")
217  value = value.replace("\n", "")
218  elif (type(val) is np.ndarray):
219  value = ""
220  for v in val:
221  value += (str(v) + " ")
222  value = value.replace("\n", "")
223 
224  else:
225  value = val
226 
227  if ((type(value) is str)):
228  if (len(value) > 50):
229  print(("warning", name, "keyword has been cut to 100 characters"))
230  #header.set(name, value[:50],'')
231  else:
232  header.set(name, value, '')
233 
234  hdulist.writeto(filepath, clobber=True) # Save changes to file
235 
236 
237 def initSimu(config, c):
238  param_dict = h5u.params_dictionary(config)
239  matricesToLoad = h5u.checkMatricesDataBase(os.environ["SHESHA_ROOT"] + "/data/",
240  config, param_dict)
241  print("->wfs")
242  wfs, tel = ao.wfs_init(config.p_wfss, config.p_atmos, config.p_tel, config.p_geom,
243  config.p_target, config.p_loop, config.p_dms)
244  print("->atmos")
245  atm = ao.atmos_init(c, config.p_atmos, config.p_tel, config.p_geom, config.p_loop,
246  config.p_wfss, wfs, config.p_target, rank=0)
247  print("->dm")
248  dms = ao.dm_init(config.p_dms, config.p_wfss, wfs, config.p_geom, config.p_tel)
249  print("->target")
250  tar = ao.target_init(c, tel, config.p_target, config.p_atmos, config.p_geom,
251  config.p_tel, config.p_dms)
252  print("->rtc")
253  rtc = ao.rtc_init(tel, wfs, config.p_wfss, dms, config.p_dms, config.p_geom,
254  config.p_rtc, config.p_atmos, atm, config.p_tel, config.p_loop,
255  do_refslp=False, clean=clean, simul_name=simul_name,
256  load=matricesToLoad, doimat=0)
257 
258  h5u.validDataBase(os.environ["SHESHA_ROOT"] + "/data/", matricesToLoad)
259 
260  print("====================")
261  print("init done")
262  print("====================")
263  print("objects initialzed on GPU:")
264  print("--------------------------------------------------------")
265  print(atm)
266  print(wfs)
267  print(dms)
268  print(tar)
269  print(rtc)
270  return wfs, tel, atm, dms, tar, rtc
271 
272 
273 def compute_modal_residuals(P, rtc, dms, tar):
274  rtc.do_control_geo(1, dms, tar, 0)
275  #self.rtc.do_control_geo_on(1, self.dms,self.tar, 0)
276  v = rtc.get_command(1)
277  ai = P.dot(v) * 1000. # np rms units
278  return ai
279 
280 
281 def loop(n, wfs, tel, atm, dms, tar, rtc, move_atmos=True, noise=True, loopData=0,
282  P=None):
283  t0 = time.time()
284  print("----------------------------------------------------")
285  print("iter# | S.E. SR | L.E. SR | Est. Rem. | framerate")
286  print("----------------------------------------------------")
287  """
288  ph = tar.get_image(0, "se")
289  pupBig = ph*0.
290 
291  phsize = pup.shape[0]
292  npup = pupBig.shape[0] # wao.wfs.get_pyrimghr(0).shape
293 
294  pupBig[(npup-phsize)/2:(npup+phsize)/2, (npup-phsize)/2:(npup+phsize)/2] = pup
295  PSFLEArray = np.zeros((config.p_target.ntargets, ph.shape[0],ph.shape[1]))
296  PSFSEArray = np.zeros((config.p_target.ntargets, ph.shape[0],ph.shape[1]))
297 
298  """
299  RmsErrorTot = []
300  ph = tar.get_image(0, "se")
301  PSFtarget = np.zeros((config.p_target.ntargets, ph.shape[0], ph.shape[1]))
302  sr_se = []
303  numiter = []
304  if (loopData):
305  if (loopData > n):
306  loopData = n
307  slopes = np.zeros((loopData, rtc.get_centroids(0).shape[0]))
308  volts = np.zeros((loopData, rtc.get_voltage(0).shape[0]))
309  else:
310  slopes = volts = None
311  ii = 0
312  jj = 0
313  sr_se = np.zeros((n / 10, config.p_target.ntargets))
314  sr_le = np.zeros((n / 10, config.p_target.ntargets))
315 
316  for i in range(n):
317  if (move_atmos):
318  atm.move_atmos()
319 
320  for t in range(config.p_target.ntargets):
321  tar.atmos_trace(t, atm, tel)
322  tar.dmtrace(t, dms)
323  for w in range(len(config.p_wfss)):
324  wfs.raytrace(w, "all", tel, atm, dms)
325  wfs.sensors_compimg(w, noise=noise)
326 
327  rtc.do_centroids(0)
328  if (loopData):
329  if (i >= (n - loopData)):
330  #print("Recording loop Data")
331  s = rtc.get_centroids(0)
332  v = rtc.get_voltage(0)
333  volts[ii, :] = v.copy()
334  slopes[ii, :] = s.copy()
335  ii += 1
336  rtc.docontrol(0)
337  rtc.doclipping(0, -1e5, 1e5)
338  rtc.apply_control(0)
339 
340  signal_le = ""
341  signal_se = ""
342  if (P is not None):
343  ai = compute_modal_residuals(P, rtc, dms, tar)
344  tarPhaseError = np.sqrt(np.sum(ai**2))
345  else:
346  tarPhaseError = 0.
347  RmsErrorTot.append(tarPhaseError)
348  print("tarPhaseError =", tarPhaseError, "nm rms")
349  if ((i + 1) % 10 == 0):
350  print(("Iter#:", i + 1, "/", n))
351  t = 0
352  SRTmp = np.zeros(config.p_target.ntargets)
353  SRTmp2 = np.zeros(config.p_target.ntargets)
354 
355  for t in range(config.p_target.ntargets):
356  if (PSFWithOtherPupil):
357  SR = list([0, 0])
358  SR[0] = PSFSEArray[t, :, :].max() # SE SR
359  SR[1] = (PSFLEArray[t, :, :] / (i + 1)).max() # LE SR
360  else:
361  SR = tar.get_strehl(t)
362  #print("Tar %d at %3.2fMicrons:" % (t+1, tar.Lambda[t]))
363  signal_se += "SR S.E %3.2fMicrons:: %1.2f " % (tar.Lambda[t], SR[0])
364  signal_le += "SR L.E %3.2fMicrons:: %1.2f " % (tar.Lambda[t], SR[1])
365  SRTmp[t] = SR[0] * 100
366  SRTmp2[t] = SR[1] * 100
367  print((signal_se + signal_le))
368  sr_se[jj, :] = SRTmp.copy()
369  sr_le[jj, :] = SRTmp2.copy()
370 
371  #sr_se.append()
372  #sr_se.append(SR[0])
373  numiter.append(i + 1)
374  jj += 1
375 
376  t1 = time.time()
377  print((" loop execution time:", t1 - t0, " (", n, "iterations), ", (t1 - t0) / n,
378  "(mean) ", n / (t1 - t0), "Hz"))
379  SRList = []
380  for t in range(config.p_target.ntargets):
381  SR = tar.get_strehl(t)
382  PSFtarget[t, :, :] = tar.get_image(t, "le")
383  SRList.append(SR[1]) # Saving Last Long Exp SR
384 
385  return SRList, tar.Lambda.tolist(), sr_se.astype(int), sr_le.astype(
386  int), numiter, slopes, volts, PSFtarget, RmsErrorTot
387 
388 
389 SR = []
390 colnames = h5u.params_dictionary(config) # config values internal to compass
391 simunames = {
392  "PSFFilenames": None,
393  "rmsError": None,
394  "rmsErrorList": None,
395  "comment": None,
396  "NCPA": None,
397  "NCPAList": None,
398  "ModalType": None,
399  "srir": None,
400  "gainModal": None,
401  "lambdaTarget": None,
402  "nbBrightest": None,
403  "sr_le": None,
404  "sr_se": None,
405  "numiter": None,
406  "NklFilt": None,
407  "NklTot": None,
408  "Nkl": None,
409  "eigenvals": None,
410  "Nphotons": None,
411  "Nactu": None,
412  "RON": None,
413  "Nslopes": None
414 } # Added values computed by the simu..
415 
416 resAll = db.readDataBase(fullpath=dBResult) # Reads all the database if exists
417 if (not (type(resAll) == pd.core.frame.DataFrame)):
418  print("Creating compass database")
419  resAll = db.createDf(list(colnames.keys()) + list(
420  simunames.keys())) # Creates the global compass Db
421 
422 # -----------------------------------------------------------------------------
423 # ----------- Replacing values from user defined variables-------------------
424 # -----------------------------------------------------------------------------
425 config.p_wfs0.set_nxsub(NSSP)
426 if (PYR):
427  decalage = int((240 - 4 - (NSSP * 2)) / 2. + NSSP / 2)
428  config.p_wfs0.set_pyr_pup_sep(decalage)
429  rMod = MODU
430  config.p_wfs0.set_pyr_npts(int(np.ceil(int(rMod * 2 * 3.141592653589793) / 4.) * 4))
431  config.p_wfs0.set_pyr_ampl(rMod)
432 
433 config.p_loop.set_ittime(1 / freq)
434 config.p_wfs0.set_noise(RON)
435 config.p_loop.set_niter(niter)
436 
437 config.p_wfs0.set_gsmag(magnitude)
438 
439 res = pd.DataFrame(columns=list(colnames.keys()) + list(simunames.keys())) # Create Db
440 wfs, tel, atm, dms, tar, rtc = initSimu(config, c) # Init COMPASS Simu!
441 
442 if (PSFWithOtherPupil):
443  #pp = config.p_geom.get_spupil().shape[0]
444  cent = config.p_geom.pupdiam / 2. + 0.5
445  oldsetting = int(config.p_tel.t_spiders)
446  config.p_tel.set_t_spiders(-1) # Enabling spiders for pupil computation
447  pupELTSpiders = mkP.make_pupil(config.p_geom.pupdiam, config.p_geom.pupdiam,
448  config.p_tel, cent, cent).astype(np.float32)
449  config.p_tel.set_t_spiders(
450  oldsetting) # Disabling spiders in case of it is used else where in GPU...
451  #PUPILPATH = "/home/fvidal/dataSimus/pupELTwithSpiders_368.fits"
452  #PUPILPATH = "/home/fvidal/dataSimus/pupELTwithSpiders.fits"
453  #PUPILPATH = "/home/fvidal/dataSimus/pupELTwithSpiders_1472.fits"
454  for target in range(config.p_target.ntargets): # Apply E-ELT pupil for each target
455  tar.set_pupil(target, pupELTSpiders.astype(np.float32))
456 
457 # ------------ ADOPT ----------------
458 ADOPTPATH = os.getenv("ADOPTPATH")
459 sys.path.insert(0, ADOPTPATH)
460 import adoptCompass as adoptComm
461 import adoptVariables as adoptVar
462 import aoCalib as cal
463 config_fileName = ADOPTPATH + "/config/ADOPT.conf"
464 wao = wao_class(config, wfs, tel, atm, dms, tar, rtc)
465 cf.returnConfigfromWao(wao, filepath=config_fileName)
466 com = adoptComm.command_class(wao, ao)
467 aoAd = adoptVar.ao_class(adoptVar.ao_attributes, adoptVar.wfs_attributes,
468  adoptVar.dm_attributes, config_fileName)
469 com.initComm(aoAd)
470 com.do_ref_slopes()
471 
472 #KL2V = com.getKL2V()
473 #
474 nfilt = nKL_Filt
475 
476 # Computing imat on diffraction limited source.
477 if (imatFromFile):
478  print("Reloading imat KL2V and gains4K from files...")
479  #print(imat0_PATH+"/"+iMatName)
480  #print(imat0_PATH+"/gains4K_MODU_"+str(int(MODU))+".fits")
481  imat = pf.get_data(imat0_PATH + "/" + iMatName)
482  modal_basis = pf.get_data(imat0_PATH + "/" + KL2VName)
483  gains4KRAW = pf.get_data(imat0_PATH + "/" + gainModalName)
484  gains4K = np.zeros(imat.shape[0] - nfilt)
485  gains4K[:-2] = gains4KRAW[:imat.shape[0] - nfilt - 2]
486  gains4K[-2:] = gains4KRAW[-2:]
487  gainopt = gains4K.copy()
488 else:
489  KL2V = com.getKL2V()
490  KL2VNorm = cal.normalizeKL2V(KL2V)
491  print("Computing Imat Diffraction Limited")
492  imat = cal.computeImatModal(com, KL2VNorm, aoAd.dm0.push4iMat, aoAd.dm1.push4iMat,
493  withTurbu=False, noise=False)
494  gains = np.linspace(1., 1., aoAd.Nactu - 2 - nfilt)
495  gains[-2:] = 1.0
496  cmat0, cmatKL0 = cal.computeCmatModal(imat, KL2VNorm, nfilt, gains)
497  com.set_command_matrix(cmat0)
498  com.close_loop()
499  print("Closing Loop with Imat Diffraction Limited")
500 
501  # Closing loop until we reach the fitting error for the given ao config + turbulence conditions (seeing ect...) but without noise and bandwidth (screen is frozen)
502  SR, lambdaTargetList, sr_se, numiter, _, _, _ = loop(200, wfs, tel, atm, dms, tar,
503  rtc, move_atmos=True,
504  noise=True)
505  print("SR After 200 iterations of closed loop:")
506 
507  if (PYR):
508  cmat0, cmatModal0 = cal.computeCmatModal(imat, modal_basis, nfilt, gains)
509  com.set_command_matrix(cmat0)
510  com.close_loop()
511  print("Closing Loop with Imat Diffraction Limited")
512 
513  # Closing loop until we reach the fitting error for the given ao config + turbulence conditions (seeing ect...) but without noise and bandwidth (screen is frozen)
514  SR, lambdaTargetList, sr_se, sr_le, numiter, _, _, _, _ = loop(
515  100, wfs, tel, atm, dms, tar, rtc, move_atmos=True, noise=True, P=P)
516  print("SR After 100 iterations of closed loop:")
517 
518  # Computing 2nd imat with this best conditions (no noise + limited by fitting)
519  imatTurbu = cal.computeImatModal(com, modal_basis, aoAd.dm0.push4iMat,
520  aoAd.dm1.push4iMat, withTurbu=True, noise=False)
521  gains4K = cal.computeOptimGainK(imat, imatTurbu, nfilt)
522 
523  date = time.strftime("_%d-%m-%Y_%H:%M:%S_")
524  gainModalName = "gains4K_" + date + ".fits"
525  iMatName = "imat_" + date + ".fits"
526  KL2VName = "modal_basis_" + date + ".fits"
527  # saving imat, modal basis, and gains...
528  pf.writeto(pathResults + "/AODATA/" + iMatName, imat)
529  pf.writeto(pathResults + "/AODATA/" + KL2VName, modal_basis)
530  pf.writeto(pathResults + "/AODATA/" + gainModalName, gains4K)
531  gainopt = gains4K.copy()
532  else:
533  gainopt = np.linspace(1., 1., aoAd.Nactu - 2 - nfilt)
534  gainopt[-2:] = 1.0
535 cmatT, cmatKLT = cal.computeCmatModal(imat, modal_basis, nfilt, gainopt * gain)
536 cmat = cmatT
537 com.set_command_matrix(cmatT)
538 com.close_loop()
539 com.resetSR()
540 
541 # ------------------------------------------------------------------------------
542 # --------------------- Modal Optim. ----------------------------------------
543 # ------------------------------------------------------------------------------
544 # Taking 2048 loop data for Optim Modal gain optim ("a la" Gendron & Lena)
545 # closing loop by adding noise + bandwidth and wait a bit that loop converge...
546 """
547 SR, lambdaTargetList, sr_le, sr_se, numiter, _, _ = loop(200,wfs,tel,atm,dms,tar,rtc, noise=True)
548 com.resetSR()
549 SR, lambdaTargetList, sr_le, sr_se, numiter, slopes, volts = loop(2048,wfs,tel,atm,dms,tar,rtc, noise=True, loopData=True)
550 V2KL =np.linalg.pinv(KL2VNorm)
551 sol = cal.recPseudoopen_loop(slopes, volts, imat, V2KL, gains4K, nfilt, 1/aoAd.Fe, aoAd.Fe)
552 gainoptCorr = cal.modalControlOptimizationopen_loopData(sol.T, cmatKL0, KL2VNorm, gmax = 1.0, Fs = aoAd.Fe, latency = 1/aoAd.Fe, BP = 1e12,ngain=200)
553 gainopt = gainopt*gainoptCorr
554 cmatOptim,_ = cal.computeCmatModal(imat, KL2VNorm, nfilt, gainopt);
555 com.set_command_matrix(cmatOptim)
556 
557 com.close_loop()
558 """
559 # ------------------------------------------------------------------------------
560 
561 #cmat = pf.get_data(os.environ["SHESHA_ROOT"]+"/test/scripts/cmatKLGood.fits").byteswap().newbyteorder()
562 #rtc.set_cmat(0, cmat.copy().astype(np.float32))
563 
564 # -----------------------------------------------------------------------------
565 # ----------- !!!!!! Starting real loop !!!!!!-------------------
566 # -----------------------------------------------------------------------------
567 
568 print("Starting Real Loop")
569 com.resetSR()
570 # com.resetDm()
571 SR, lambdaTargetList, sr_se, sr_le, numiter, slopesCB, voltsCB, PSFtarget, rmsErrorList = loop(
572  config.p_loop.niter, wfs, tel, atm, dms, tar, rtc, loopData=nbLoopData, P=P)
573 
574 if (saveCBData):
575  PYRImage = wfs.get_pyrimg(0)
576  date = time.strftime("_%d-%m-%Y_%H:%M:%S_")
577  slopesCBName = "slopesCB_" + date + ".fits"
578  voltsCBName = "voltsCB_" + date + ".fits"
579  PYRIMAGEName = "pyrImageCB_" + date + ".fits"
580  SRHistorySEName = "SRHistorySE_" + date + ".fits"
581  SRHistoryLEName = "SRHistoryLE_" + date + ".fits"
582  pf.writeto(pathResults + "/CircularBuffers/" + slopesCBName, slopesCB.copy())
583  pf.writeto(pathResults + "/CircularBuffers/" + voltsCBName, voltsCB.copy())
584  pf.writeto(pathResults + "/CircularBuffers/" + PYRIMAGEName, PYRImage.copy())
585  pf.writeto(pathResults + "/CircularBuffers/" + SRHistorySEName, sr_se.copy())
586  pf.writeto(pathResults + "/CircularBuffers/" + SRHistoryLEName, sr_le.copy())
587 
588 else:
589  slopesCBName = ""
590  voltsCBName = ""
591  PYRIMAGEName = ""
592  SRHistorySEName = ""
593  SRHistoryLEName = ""
594 
595 # ------------- Saving config and results in data frame -----
596 dfparams = h5u.params_dictionary(config) # get the current compass config
597 for i in range(len(dfparams)):
598  if (type(dfparams[list(dfparams.keys())[i]]) is list): # If list
599  if (len(dfparams[list(dfparams.keys())[i]]) == 1): # If list has only 1 element:
600  dfparams[list(dfparams.keys())[i]] = dfparams[list(dfparams.keys())[i]][0]
601 
602 dfparams.update(simunames) # Add the simunames params
603 
604 res = db.fillDf(res, dfparams) # Saving dictionnary config
605 res.loc[0, "iMatName"] = iMatName
606 res.loc[0, "ModalName"] = KL2VName
607 res.loc[0, "gainModalName"] = gainModalName
608 res.loc[0, "slopesCBName"] = slopesCBName
609 res.loc[0, "voltsCBName"] = voltsCBName
610 res.loc[0, "SRHistoryLEName"] = SRHistorySEName
611 res.loc[0, "SRHistorySEName"] = SRHistoryLEName
612 res.loc[0, "PYRIMAGEName"] = PYRIMAGEName
613 res.loc[0, "comment"] = comment
614 res.loc[0, "NklFilt"] = nKL_Filt
615 res.loc[0, "Nkl"] = imat.shape[0] - nfilt - 2
616 res.loc[0, "NklTot"] = cmat.shape[1] - 2
617 res.loc[0, "Nactu"] = cmat.shape[1]
618 res.loc[0, "Nslopes"] = cmat.shape[0]
619 res.loc[0, "Nphotons"] = config.p_wfs0._nphotons
620 res.loc[0, "RON"] = RON
621 #res.eigenvals.values[0] = rtc.getEigenvals(0)
622 res.srir.values[0] = SR # Saving computed values
623 res.lambdaTarget.values[0] = lambdaTargetList
624 res.loc[0, "gsmag"] = config.p_wfs0.gsmag
625 res.loc[0, "gain"] = gain
626 #res.loc[0, "type_ap"] = str(res.loc[0, "type_ap"][0])
627 #res.loc[0, "type_wfs"] = str(res.loc[0, "type_wfs"][0])
628 res.loc[0, "type_dm"] = "pzt, tt"
629 #res.loc[0, "npix"] = res.loc[0, "npix"][0]
630 #res.loc[0, "nbBrightest"] = res.loc[0, "nbBrightest"][0]
631 res.loc[0, "pixsizeInMeters"] = (config.p_tel.diam / config.p_geom.get_spupil().shape[0])
632 # PSF pixsize =
633 #res.sr_le.values[0] = sr_le
634 #res.sr_se.values[0] = sr_se
635 #res.numiter.values[0] = numiter
636 res.loc[0, "simulname"] = simulName
637 res.rmsErrorList.values[0] = rmsErrorList
638 res.rmsError.values[0] = np.average(np.array(rmsErrorList))
639 
640 # --------------- PSF Stuff ----------------------
641 print("Saving PSFs...")
642 PSFNameList = []
643 for t in range(config.p_target.ntargets):
644 
645  date = time.strftime("_%d-%m-%Y_%H:%M:%S_")
646  lam = "%3.2f" % tar.Lambda.tolist()[t]
647  lam = lam.replace(".", "_")
648  PSFName = "PYR_TARGET_" + str(t + 1) + "_Lambda_" + lam + "_" + date + ".fits"
649  PSFNameList.append(PSFName)
650  #PSFNameList.append("NOT SAVED")
651  if (savePSFs):
652  pf.writeto(pathResults + "/PSFs/" + PSFName, PSFtarget[t, :, :].copy(),
653  clobber=True)
654  lam2 = "%3.2f" % tar.Lambda.tolist()[t]
655  res.loc[0, "SR_%s" % lam2] = SR[t]
656  PSFPixsize = (tar.Lambda.tolist()[t] * 1e-6) / (
657  wao.config.p_tel.diam / wao.config.p_geom.get_spupil().shape[0] *
658  wao.config.p_geom.get_ipupil().shape[0]) * 206265.
659  res.loc[0, "pixsizeArcSec_%s" % lam2] = PSFPixsize
660  filepath = pathResults + "/PSFs/" + PSFName
661  if (savePSFs):
662  #"Add the SR and wavelegth value at the top of the PSF header file"
663  hdulist = pf.open(filepath) # read file
664  header = hdulist[0].header
665  header["SR"] = SR[t]
666  header["wavelengthMic"] = tar.Lambda.tolist()[t]
667  header["pixsizeArcSec"] = PSFPixsize
668  hdulist.writeto(filepath, clobber=True) # Save changes to file
669  # Adding all the parameters to the header
670  makeFITSHeader(filepath, res)
671  else:
672  res.PSFFilenames.values[0] = ["PSF NOT SAVED"]
673 print("Done")
674 res.PSFFilenames.values[0] = PSFNameList
675 
676 resAll = db.fillDf(resAll, res) # Saving in global DB
677 #resAll.to_hdf("/home/fvidal/compass/trunk/shesha/test/scripts/resultatsScripts/SH39m.h5", "resAll", complevel=9,complib='blosc')
678 resAll.to_hdf(dBResult, "resAll", complevel=9, complib='blosc')
679 #db.saveDataBase(resAll)
680 
681 print("Simulation Done...")
682 """
683 except Exception as e:
684  if(e):
685  print "Simu failed"
686  print e
687  sendMail("Script failed!" + str(sys.argv), str(e)+ str("\nline: ") +str(sys.exc_info()[-1].tb_lineno))
688 """
689 """Convert resAll not well formatted:
690 
691 
692 
693 for i in range(len(dfparams)):
694  if(type(dfparams[dfparams.keys()[i]]) is list): # If list
695  if(len(dfparams[dfparams.keys()[i]]) == 1): # If list has only 1 element:
696  dfparams[dfparams.keys()[i]] = dfparams[dfparams.keys()[i]][0]
697 
698 
699 
700 
701 for i in range(len(list(set(resAll)))):
702  if(type(resAll[list(set(resAll))[i]][0]) is list): # If list
703  if(len(resAll[list(set(resAll))[i]][0]) == 1): # If list has only 1 element:
704  key = list(set(resAll))[i]
705  for j in range(len(resAll)):
706  try:
707  resAll.loc[j, key] = resAll.loc[j, key][0]
708  except:
709  print "Warning! Could not convert list "+key+" with 1 element to 1 element"
710  pass
711 
712 
713 
714 
715 """
script_PYR39m_optimGain.initSimu
def initSimu(config, c)
Definition: script_PYR39m_optimGain.py:237
script_PYR39m_optimGain.wao_class.wfs
wfs
Definition: script_PYR39m_optimGain.py:197
shesha.util
Utilities functions.
Definition: shesha/shesha/util/__init__.py:1
shesha.constants
Numerical constants for shesha and config enumerations for safe-typing.
Definition: constants.py:1
script_PYR39m_optimGain.wao_class.tar
tar
Definition: script_PYR39m_optimGain.py:201
script_PYR39m_optimGain.wao_class.dms
dms
Definition: script_PYR39m_optimGain.py:200
script_PYR39m_optimGain.wao_class.tel
tel
Definition: script_PYR39m_optimGain.py:198
shesha.config
Parameter classes for COMPASS.
Definition: shesha/shesha/config/__init__.py:1
script_PYR39m_optimGain.wao_class.rtc
rtc
Definition: script_PYR39m_optimGain.py:202
script_PYR39m_optimGain.sendMail
def sendMail(message, title)
Test Line.
Definition: script_PYR39m_optimGain.py:9
script_PYR39m_optimGain.wao_class.config
config
Definition: script_PYR39m_optimGain.py:196
script_PYR39m_optimGain.loop
def loop(n, wfs, tel, atm, dms, tar, rtc, move_atmos=True, noise=True, loopData=0, P=None)
Definition: script_PYR39m_optimGain.py:281
script_PYR39m_optimGain.compute_modal_residuals
def compute_modal_residuals(P, rtc, dms, tar)
Definition: script_PYR39m_optimGain.py:273
script_PYR39m_optimGain.wao_class.atm
atm
Definition: script_PYR39m_optimGain.py:199
script_PYR39m_optimGain.wao_class.__init__
def __init__(self, config, wfs, tel, atm, dms, tar, rtc)
Definition: script_PYR39m_optimGain.py:195
script_PYR39m_optimGain.wao_class
Definition: script_PYR39m_optimGain.py:193
script_PYR39m_optimGain.makeFITSHeader
def makeFITSHeader(filepath, df)
Definition: script_PYR39m_optimGain.py:205