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
11 from email.mime.text
import MIMEText
12 smtp = smtplib.SMTP(
'smtp.obspm.fr')
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())
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/")
32 import carmaWrap
as ch
39 import matplotlib.pyplot
as plt
40 import hdf5_util
as h5u
41 import resDataBase
as db
42 import astropy.io.fits
as pf
45 import compassConfigToFile
as cf
46 import make_pupil
as mkP
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...")
65 simulName =
"PYR_39m_RoundPupil_FromHippo6"
67 print(
"-------------------------------------")
68 print(
"DETECTED BASH SCRIPT with parameters:")
70 print(
"-------------------------------------")
72 freq = float(sys.argv[2])
73 RON = float(sys.argv[3])
74 MODU = float(sys.argv[4])
75 gain = float(sys.argv[5])
76 magnitude = float(sys.argv[6])
77 nKL_Filt = int(float(sys.argv[7]))
78 simulName = sys.argv[8]
79 NSSP = int(sys.argv[9])
80 GPU = int(sys.argv[10])
81 comment =
"SRVsGSVsNControlledModes"
83 print((
"Freq=", freq))
85 print((
"MODU=", MODU))
86 print((
"gain=", gain))
87 print((
"magnitude=", magnitude))
88 print((
"nKL_Filt=", nKL_Filt))
90 print((
"simulName=", simulName))
92 pathResults =
"/volumes/hra/micado/" + simulName
94 dBResult = pathResults +
"/" + simulName +
".h5"
107 imat0_PATH =
"/home/fvidal/dataSimus"
108 iMatName =
"imatDiffraction_ELTPYR_35Layers.fits"
109 gainModalName =
"gains4K_ELTPYR_35Layers.fits"
110 KL2VName =
"KL2VNorm_ELTPYR_35Layers.fits"
112 iMatName = "iMat_MODU_2_ELTPUPIL.fits"
113 KL2VName = "KL2VNorm_MODU_2_ELTPUPIL.fits"
114 gainModalName = "gains4K_MODU_2_ELTPUPIL.fits"
117 ModalBasisType =
"Btt"
119 PSFWithOtherPupil =
True
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"
134 GPUs = np.array([4, 5, 6, 7], dtype=np.int32)
136 GPUs = np.array([GPU], dtype=np.int32)
137 print((
"Using GPUs: ", GPUs))
139 GPUs = np.array([4, 5, 6, 7], dtype=np.int32)
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/")
156 param_file = sys.argv[1]
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])
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
167 h5u.configFromH5(param_file, config)
169 raise ValueError(
"Parameter file extension must be .py or .h5")
171 print((
"param_file is", param_file))
173 if (hasattr(config,
"simul_name")):
174 if (config.simul_name
is None):
177 simul_name = config.simul_name
180 print((
"simul name is", simul_name))
183 if (simul_name ==
""):
187 param_dict = h5u.params_dictionary(config)
188 matricesToLoad = h5u.checkMatricesDataBase(os.environ[
"SHESHA_ROOT"] +
"/data/",
190 c = ch.carmaWrap_context(devices=GPUs)
195 def __init__(self, config, wfs, tel, atm, dms, tar, rtc):
206 hdulist = pf.open(filepath)
207 header = hdulist[0].header
208 names = np.sort(list(set(df))).tolist()
211 if (name !=
"centroider.thresh"):
213 if (type(val)
is list):
216 value += (str(v) +
" ")
217 value = value.replace(
"\n",
"")
218 elif (type(val)
is np.ndarray):
221 value += (str(v) +
" ")
222 value = value.replace(
"\n",
"")
227 if ((type(value)
is str)):
228 if (len(value) > 50):
229 print((
"warning", name,
"keyword has been cut to 100 characters"))
232 header.set(name, value,
'')
234 hdulist.writeto(filepath, clobber=
True)
238 param_dict = h5u.params_dictionary(config)
239 matricesToLoad = h5u.checkMatricesDataBase(os.environ[
"SHESHA_ROOT"] +
"/data/",
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)
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)
248 dms = ao.dm_init(config.p_dms, config.p_wfss, wfs, config.p_geom, config.p_tel)
250 tar = ao.target_init(c, tel, config.p_target, config.p_atmos, config.p_geom,
251 config.p_tel, config.p_dms)
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)
258 h5u.validDataBase(os.environ[
"SHESHA_ROOT"] +
"/data/", matricesToLoad)
260 print(
"====================")
262 print(
"====================")
263 print(
"objects initialzed on GPU:")
264 print(
"--------------------------------------------------------")
270 return wfs, tel, atm, dms, tar, rtc
274 rtc.do_control_geo(1, dms, tar, 0)
276 v = rtc.get_command(1)
277 ai = P.dot(v) * 1000.
281 def loop(n, wfs, tel, atm, dms, tar, rtc, move_atmos=True, noise=True, loopData=0,
284 print(
"----------------------------------------------------")
285 print(
"iter# | S.E. SR | L.E. SR | Est. Rem. | framerate")
286 print(
"----------------------------------------------------")
288 ph = tar.get_image(0, "se")
291 phsize = pup.shape[0]
292 npup = pupBig.shape[0] # wao.wfs.get_pyrimghr(0).shape
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]))
300 ph = tar.get_image(0,
"se")
301 PSFtarget = np.zeros((config.p_target.ntargets, ph.shape[0], ph.shape[1]))
307 slopes = np.zeros((loopData, rtc.get_centroids(0).shape[0]))
308 volts = np.zeros((loopData, rtc.get_voltage(0).shape[0]))
310 slopes = volts =
None
313 sr_se = np.zeros((n / 10, config.p_target.ntargets))
314 sr_le = np.zeros((n / 10, config.p_target.ntargets))
320 for t
in range(config.p_target.ntargets):
321 tar.atmos_trace(t, atm, tel)
323 for w
in range(len(config.p_wfss)):
324 wfs.raytrace(w,
"all", tel, atm, dms)
325 wfs.sensors_compimg(w, noise=noise)
329 if (i >= (n - loopData)):
331 s = rtc.get_centroids(0)
332 v = rtc.get_voltage(0)
333 volts[ii, :] = v.copy()
334 slopes[ii, :] = s.copy()
337 rtc.doclipping(0, -1e5, 1e5)
344 tarPhaseError = np.sqrt(np.sum(ai**2))
347 RmsErrorTot.append(tarPhaseError)
348 print(
"tarPhaseError =", tarPhaseError,
"nm rms")
349 if ((i + 1) % 10 == 0):
350 print((
"Iter#:", i + 1,
"/", n))
352 SRTmp = np.zeros(config.p_target.ntargets)
353 SRTmp2 = np.zeros(config.p_target.ntargets)
355 for t
in range(config.p_target.ntargets):
356 if (PSFWithOtherPupil):
358 SR[0] = PSFSEArray[t, :, :].max()
359 SR[1] = (PSFLEArray[t, :, :] / (i + 1)).max()
361 SR = tar.get_strehl(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()
373 numiter.append(i + 1)
377 print((
" loop execution time:", t1 - t0,
" (", n,
"iterations), ", (t1 - t0) / n,
378 "(mean) ", n / (t1 - t0),
"Hz"))
380 for t
in range(config.p_target.ntargets):
381 SR = tar.get_strehl(t)
382 PSFtarget[t, :, :] = tar.get_image(t,
"le")
385 return SRList, tar.Lambda.tolist(), sr_se.astype(int), sr_le.astype(
386 int), numiter, slopes, volts, PSFtarget, RmsErrorTot
390 colnames = h5u.params_dictionary(config)
392 "PSFFilenames":
None,
394 "rmsErrorList":
None,
401 "lambdaTarget":
None,
416 resAll = db.readDataBase(fullpath=dBResult)
417 if (
not (type(resAll) == pd.core.frame.DataFrame)):
418 print(
"Creating compass database")
419 resAll = db.createDf(list(colnames.keys()) + list(
425 config.p_wfs0.set_nxsub(NSSP)
427 decalage = int((240 - 4 - (NSSP * 2)) / 2. + NSSP / 2)
428 config.p_wfs0.set_pyr_pup_sep(decalage)
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)
433 config.p_loop.set_ittime(1 / freq)
434 config.p_wfs0.set_noise(RON)
435 config.p_loop.set_niter(niter)
437 config.p_wfs0.set_gsmag(magnitude)
439 res = pd.DataFrame(columns=list(colnames.keys()) + list(simunames.keys()))
442 if (PSFWithOtherPupil):
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)
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(
454 for target
in range(config.p_target.ntargets):
455 tar.set_pupil(target, pupELTSpiders.astype(np.float32))
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"
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)
478 print(
"Reloading imat KL2V and gains4K from files...")
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()
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)
496 cmat0, cmatKL0 = cal.computeCmatModal(imat, KL2VNorm, nfilt, gains)
497 com.set_command_matrix(cmat0)
499 print(
"Closing Loop with Imat Diffraction Limited")
502 SR, lambdaTargetList, sr_se, numiter, _, _, _ =
loop(200, wfs, tel, atm, dms, tar,
503 rtc, move_atmos=
True,
505 print(
"SR After 200 iterations of closed loop:")
508 cmat0, cmatModal0 = cal.computeCmatModal(imat, modal_basis, nfilt, gains)
509 com.set_command_matrix(cmat0)
511 print(
"Closing Loop with Imat Diffraction Limited")
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:")
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)
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"
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()
533 gainopt = np.linspace(1., 1., aoAd.Nactu - 2 - nfilt)
535 cmatT, cmatKLT = cal.computeCmatModal(imat, modal_basis, nfilt, gainopt * gain)
537 com.set_command_matrix(cmatT)
547 SR, lambdaTargetList, sr_le, sr_se, numiter, _, _ = loop(200,wfs,tel,atm,dms,tar,rtc, noise=True)
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)
568 print(
"Starting Real Loop")
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)
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())
596 dfparams = h5u.params_dictionary(config)
597 for i
in range(len(dfparams)):
598 if (type(dfparams[list(dfparams.keys())[i]])
is list):
599 if (len(dfparams[list(dfparams.keys())[i]]) == 1):
600 dfparams[list(dfparams.keys())[i]] = dfparams[list(dfparams.keys())[i]][0]
602 dfparams.update(simunames)
604 res = db.fillDf(res, dfparams)
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
622 res.srir.values[0] = SR
623 res.lambdaTarget.values[0] = lambdaTargetList
624 res.loc[0,
"gsmag"] = config.p_wfs0.gsmag
625 res.loc[0,
"gain"] = gain
628 res.loc[0,
"type_dm"] =
"pzt, tt"
631 res.loc[0,
"pixsizeInMeters"] = (config.p_tel.diam / config.p_geom.get_spupil().shape[0])
636 res.loc[0,
"simulname"] = simulName
637 res.rmsErrorList.values[0] = rmsErrorList
638 res.rmsError.values[0] = np.average(np.array(rmsErrorList))
641 print(
"Saving PSFs...")
643 for t
in range(config.p_target.ntargets):
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)
652 pf.writeto(pathResults +
"/PSFs/" + PSFName, PSFtarget[t, :, :].copy(),
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
663 hdulist = pf.open(filepath)
664 header = hdulist[0].header
666 header[
"wavelengthMic"] = tar.Lambda.tolist()[t]
667 header[
"pixsizeArcSec"] = PSFPixsize
668 hdulist.writeto(filepath, clobber=
True)
672 res.PSFFilenames.values[0] = [
"PSF NOT SAVED"]
674 res.PSFFilenames.values[0] = PSFNameList
676 resAll = db.fillDf(resAll, res)
678 resAll.to_hdf(dBResult,
"resAll", complevel=9, complib=
'blosc')
681 print(
"Simulation Done...")
683 except Exception as e:
687 sendMail("Script failed!" + str(sys.argv), str(e)+ str("\nline: ") +str(sys.exc_info()[-1].tb_lineno))
689 """Convert resAll not well formatted:
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]
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)):
707 resAll.loc[j, key] = resAll.loc[j, key][0]
709 print "Warning! Could not convert list "+key+" with 1 element to 1 element"