2 Created on Wed Apr 27 09:28:23 2016
12 import carmaWrap
as ch
15 import matplotlib.pyplot
as pl
17 import hdf5_util
as h5u
19 from scipy.sparse
import csr_matrix
21 if (len(sys.argv) < 2):
22 error =
'command line should be at least:"python -i test.py parameters_filename"\n with "parameters_filename" the path to the parameters file'
23 raise Exception(error)
26 param_file = sys.argv[1]
27 if (param_file.split(
'.')[-1] == b
"py"):
28 filename = param_file.split(
'/')[-1]
29 param_path = param_file.split(filename)[0]
30 sys.path.insert(0, param_path)
31 exec(
"import %s as config" % filename.split(
".py")[0])
32 sys.path.remove(param_path)
33 elif (param_file.split(
'.')[-1] == b
"h5"):
34 sys.path.insert(0, os.environ[
"SHESHA_ROOT"] +
"/data/par/par4bench/")
35 import scao_sh_16x16_8pix
as config
36 sys.path.remove(os.environ[
"SHESHA_ROOT"] +
"/data/par/par4bench/")
37 h5u.configFromH5(param_file, config)
39 raise ValueError(
"Parameter file extension must be .py or .h5")
41 print(
"param_file is", param_file)
43 if (len(sys.argv) > 2):
44 savename = sys.argv[2]
46 savename =
"roket_default.h5"
56 if (hasattr(config,
"simul_name")):
57 if (config.simul_name
is None):
60 simul_name = config.simul_name
63 print(
"simul name is", simul_name)
66 if (simul_name == b
""):
70 param_dict = h5u.params_dictionary(config)
71 matricesToLoad = h5u.checkMatricesDataBase(os.environ[
"SHESHA_ROOT"] +
"/data/",
76 c = ch.carmaWrap_context(devices=np.array([6], dtype=np.int32))
81 wfs, tel = ao.wfs_init(config.p_wfss, config.p_atmos, config.p_tel, config.p_geom,
82 config.p_target, config.p_loop, config.p_dms)
86 atm = ao.atmos_init(c, config.p_atmos, config.p_tel, config.p_geom, config.p_loop,
87 config.p_wfss, wfs, config.p_target, rank=0, clean=clean,
92 dms = ao.dm_init(config.p_dms, config.p_wfss, wfs, config.p_geom, config.p_tel)
96 tar = ao.target_init(c, tel, config.p_target, config.p_atmos, config.p_geom,
97 config.p_tel, config.p_dms)
101 rtc = ao.rtc_init(tel, wfs, config.p_wfss, dms, config.p_dms, config.p_geom,
102 config.p_rtc, config.p_atmos, atm, config.p_tel, config.p_loop,
103 clean=clean, simul_name=simul_name, load=matricesToLoad)
106 h5u.validDataBase(os.environ[
"SHESHA_ROOT"] +
"/data/", matricesToLoad)
108 print(
"====================")
110 print(
"====================")
111 print(
"objects initialzed on GPU:")
112 print(
"--------------------------------------------------------")
119 print(
"----------------------------------------------------")
120 print(
"iter# | SE SR image | LE SR image | Fitting | LE SR phase var")
121 print(
"----------------------------------------------------")
123 error_flag =
True in [w.roket
for w
in config.p_wfss]
135 Performs the main AO loop for n interations. First, initialize buffers
136 for error breakdown computations. Then, at the end of each iteration, just
137 before applying the new DM shape, calls the error_breakdown function.
139 :param n: (int) : number of iterations
142 com : (np.array((n,nactus))) : full command buffer
144 noise_com : (np.array((n,nactus))) : noise contribution for error breakdown
146 alias_wfs_com : (np.array((n,nactus))) : aliasing estimation in the WFS direction
148 tomo_com : (np.array((n,nactus))) : tomography error estimation
150 H_com : (np.array((n,nactus))) : Filtered modes contribution for error breakdown
152 trunc_com : (np.array((n,nactus))) : Truncature and sampling error of WFS
154 bp_com : (np.array((n,nactus))) : Bandwidth error estimation on target
156 mod_com : (np.array((n,nactus))) : Commanded modes expressed on the actuators
158 fit : (float) : fitting (mean variance of the residual target phase after projection)
160 SR : (float) : final strehl ratio returned by the simulation
164 nactu = rtc.get_command(0).size
165 com = np.zeros((n, nactu), dtype=np.float32)
166 noise_com = np.zeros((n, nactu), dtype=np.float32)
167 alias_wfs_com = np.copy(noise_com)
168 wf_com = np.copy(noise_com)
169 tomo_com = np.copy(noise_com)
170 trunc_com = np.copy(noise_com)
171 H_com = np.copy(noise_com)
172 mod_com = np.copy(noise_com)
173 bp_com = np.copy(noise_com)
175 psf_ortho = tar.get_image(0,
'se') * 0.
176 Ee = np.copy(noise_com)
179 gRD = np.identity(RD.shape[0]) - config.p_controllers[0].gain * gamma * RD
184 if (config.p_controllers[0].type == b
"geo"):
185 for t
in range(config.p_target.ntargets):
186 tar.atmos_trace(t, atm, tel)
187 rtc.docontrol_geo(0, dms, tar, 0)
188 rtc.applycontrol(0, dms)
191 for t
in range(config.p_target.ntargets):
192 tar.atmos_trace(t, atm, tel)
194 for w
in range(len(config.p_wfss)):
195 wfs.sensors_trace(w,
"all", tel, atm, dms)
196 wfs.sensors_compimg(w)
202 if (error_flag
and i > -1):
205 trunc_com, bp_com, wf_com, mod_com, fit, psf_ortho, i,
208 rtc.applycontrol(0, dms)
210 if ((i + 1) % 100 == 0
and i > -1):
211 strehltmp = tar.get_strehl(0)
212 print(i + 1,
"\t", strehltmp[0],
"\t", strehltmp[1],
"\t",
213 np.exp(-strehltmp[2]),
"\t", np.exp(-strehltmp[3]))
215 print(
" loop execution time:", t1 - t0,
" (", n,
"iterations), ", (t1 - t0) / n,
216 "(mean) ", n / (t1 - t0),
"Hz")
219 SR2 = np.exp(-tar.get_strehl(0, comp_strehl=
False)[3])
220 SR = tar.get_strehl(0, comp_strehl=
False)[1]
223 return com, noise_com, alias_wfs_com, tomo_com, H_com, trunc_com, bp_com, mod_com, np.mean(
224 fit[N_preloop:]), SR, SR2, psf_ortho, Ee, Ff
229 Performs the main AO loop for n interations. First, initialize buffers
230 for error breakdown computations. Then, at the end of each iteration, just
231 before applying the new DM shape, calls the error_breakdown function.
233 :param n: (int) : number of iterations
236 com : (np.array((n,nactus))) : full command buffer
238 noise_com : (np.array((n,nactus))) : noise contribution for error breakdown
240 alias_wfs_com : (np.array((n,nactus))) : aliasing estimation in the WFS direction
242 tomo_com : (np.array((n,nactus))) : tomography error estimation
244 H_com : (np.array((n,nactus))) : Filtered modes contribution for error breakdown
246 trunc_com : (np.array((n,nactus))) : Truncature and sampling error of WFS
248 bp_com : (np.array((n,nactus))) : Bandwidth error estimation on target
250 mod_com : (np.array((n,nactus))) : Commanded modes expressed on the actuators
252 fit : (float) : fitting (mean variance of the residual target phase after projection)
254 SR : (float) : final strehl ratio returned by the simulation
256 for i
in range(0, n):
259 if (config.p_controllers[0].type == b
"geo"):
260 for t
in range(config.p_target.ntargets):
261 tar.atmos_trace(t, atm, tel)
262 rtc.docontrol_geo(0, dms, tar, 0)
263 rtc.applycontrol(0, dms)
265 for t
in range(config.p_target.ntargets):
266 tar.atmos_trace(t, atm, tel)
267 for w
in range(len(config.p_wfss)):
268 wfs.sensors_trace(w,
"all", tel, atm, dms)
269 wfs.sensors_compimg(w)
273 rtc.applycontrol(0, dms)
282 def error_breakdown(com, noise_com, alias_wfs_com, tomo_com, H_com, trunc_com, bp_com,
283 wf_com, mod_com, fit, psf_ortho, i, Ee, Ff, gamma, gRD):
285 Compute the error breakdown of the AO simulation. Returns the error commands of
286 each contributors. Suppose no delay (for now) and only 2 controllers : the main one, controller #0, (specified on the parameter file)
287 and the geometric one, controller #1 (automatically added if roket is asked in the parameter file)
288 Commands are computed by applying the loop filter on various kind of commands : (see schema_simulation_budget_erreur_v2)
290 - Ageom : Aliasing contribution on WFS direction
291 Obtained by computing commands from DM orthogonal phase (projection + slopes_geom)
293 - B : Projection on the target direction
294 Obtained as the commmands output of the geometric controller
297 Obtained by computing commands from DM parallel phase (RD*B)
299 - E : Wavefront + aliasing + ech/trunc + tomo
300 Obtained by performing the AO loop iteration without noise on the WFS
302 - F : Wavefront + aliasing + tomo
303 Obtained by performing the AO loop iteration without noise on the WFS and using phase deriving slopes
307 Note : rtc.get_err returns to -CMAT.slopes
310 noise_com : np.array((niter,nactu)) : Noise contribution
313 alias_wfs_com : np.array((niter,nactu)) : Aliasing on WFS direction contribution
316 tomo_com : np.array((niter,nactu)) : Tomographic error contribution
319 H_com : np.array((niter,nactu)) : Filtered modes error
322 trunc_com : np.array((niter,nactu)) : sampling/truncature error contribution
325 bp_com : np.array((niter,nactu)) : Bandwidth error
327 wf_com : np.array((niter,nactu)) : Reconstructed wavefront
329 mod_com : np.array((niter,nactu)) : commanded modes
331 fit : np.array((niter)) : fitting value
333 i : (int) : current iteration number
336 g = config.p_controllers[0].gain
337 Dcom = rtc.get_command(0)
338 Derr = rtc.get_err(0)
340 tarphase = tar.get_phase(0)
344 if (config.p_wfss[0].type == b
"sh"):
345 ideal_bincube = wfs.get_bincubeNotNoisy(0)
346 bincube = wfs.get_bincube(0)
347 if (config.p_centroiders[0].type == b
"tcog"
349 invalidpix = np.where(bincube <= config.p_centroiders[0].thresh)
350 ideal_bincube[invalidpix] = 0
351 rtc.setthresh(0, -1e16)
352 wfs.set_bincube(0, ideal_bincube)
353 elif (config.p_wfss[0].type == b
"pyrhr"):
354 ideal_pyrimg = wfs.get_binimg_notnoisy(0)
355 wfs.set_pyrimg(0, ideal_pyrimg)
358 if (config.p_centroiders[0].type == b
"tcog"):
359 rtc.setthresh(0, config.p_centroiders[0].thresh)
365 if (i + 1 < config.p_loop.niter):
366 noise_com[i + 1, :] = gRD.dot(noise_com[i, :]) + g * (Derr - E)
371 rtc.docentroids_geom(0)
376 if (i + 1 < config.p_loop.niter):
377 trunc_com[i + 1, :] = gRD.dot(trunc_com[i, :]) + g * (E - gamma * F)
382 rtc.docontrol_geo_onwfs(1, dms, wfs, 0)
383 rtc.applycontrol(1, dms)
384 for w
in range(len(config.p_wfss)):
385 wfs.sensors_trace(w,
"dm", tel, atm, dms)
387 wfs.sensors_compimg(0)
388 if(config.p_wfss[0].type == b"sh"):
389 ideal_bincube = wfs.get_bincubeNotNoisy(0)
390 bincube = wfs.get_bincube(0)
391 if(config.p_centroiders[0].type == b"tcog"): # Select the same pixels with or without noise
392 invalidpix = np.where(bincube <= config.p_centroiders[0].thresh)
393 ideal_bincube[invalidpix] = 0
394 rtc.setthresh(0,-1e16)
395 wfs.set_bincube(0,ideal_bincube)
396 elif(config.p_wfss[0].type == b"pyrhr"):
397 ideal_pyrimg = wfs.get_binimg_notnoisy(0)
398 wfs.set_pyrimg(0,ideal_pyrimg)
400 rtc.docentroids_geom(0)
402 Ageom = rtc.get_err(0)
403 if (i + 1 < config.p_loop.niter):
404 alias_wfs_com[i + 1, :] = gRD.dot(
405 alias_wfs_com[i, :]) + gamma * g * (Ageom)
410 tar.atmos_trace(0, atm, tel)
411 rtc.docontrol_geo(1, dms, tar, 0)
412 B = rtc.get_command(1)
417 rtc.applycontrol(1, dms)
418 tar.dmtrace(0, dms, do_phase_var=0)
419 fit[i] = tar.get_strehl(0, comp_strehl=
False)[2]
421 psf_ortho += tar.get_image(0,
'se') / niters
427 modes_filt = modes.copy() * 0.
428 modes_filt[-nfiltered - 2:-2] = modes[-nfiltered - 2:-2]
429 H_com[i, :] = Btt.dot(modes_filt)
430 modes[-nfiltered - 2:-2] = 0
431 mod_com[i, :] = Btt.dot(modes)
436 C = mod_com[i, :] - mod_com[i - 1, :]
438 bp_com[i, :] = gRD.dot(bp_com[i - 1, :]) - C
444 for w
in range(len(config.p_wfss)):
445 wfs.sensors_trace(w,
"atmos", tel, atm, dms)
446 rtc.docontrol_geo_onwfs(1, dms, wfs, 0)
447 G = rtc.get_command(1)
449 modes[-nfiltered - 2:-2] = 0
450 wf_com[i, :] = Btt.dot(modes)
452 G = mod_com[i, :] - wf_com[i, :]
453 if (i + 1 < config.p_loop.niter):
454 tomo_com[i + 1, :] = gRD.dot(tomo_com[i, :]) - g * RD.dot(G)
457 tar.set_phase(0, tarphase)
463 cgains = np.zeros(E.shape[1])
464 for k
in range(E.shape[1]):
465 cgains[k] = np.polyfit(E[:, k], F[:, k], 1)[0]
467 return np.mean(cgains)
477 IF = rtc.get_IFsparse(1).T
485 delta = IF.T.dot(IF).toarray() / N
488 Tp = np.ones((T.shape[0], T.shape[1] + 1))
490 deltaT = IF.T.dot(Tp) / N
492 tau = np.linalg.inv(delta).dot(deltaT)
496 tdt = tau.T.dot(delta).dot(tau)
497 subTT = tau.dot(np.linalg.inv(tdt)).dot(tau.T).dot(delta)
501 gdg = G.T.dot(delta).dot(G)
502 U, s, V = np.linalg.svd(gdg)
503 U = U[:, :U.shape[1] - 3]
505 L = np.identity(s.size) / np.sqrt(s)
510 Btt = np.zeros((n + 2, n - 1))
511 Btt[:B.shape[0], :B.shape[1]] = B
512 mini = 1. / np.sqrt(np.abs(TT))
515 Btt[n:, n - 3:] = mini
518 delta = np.zeros((n + T.shape[1], n + T.shape[1]))
520 delta[:-2, :-2] = IF.T.dot(IF).toarray() / N
521 delta[-2:, -2:] = T.T.dot(T) / N
524 return Btt.astype(np.float32), P.astype(np.float32)
531 Btt_filt = np.zeros((Btt.shape[0], Btt.shape[1] - nfilt))
532 Btt_filt[:, :Btt_filt.shape[1] - 2] = Btt[:, :Btt.shape[1] - (nfilt + 2)]
533 Btt_filt[:, Btt_filt.shape[1] - 2:] = Btt[:, Btt.shape[1] - 2:]
538 Dmp = np.linalg.inv(Dm.T.dot(Dm)).dot(Dm.T)
540 cmat = Btt_filt.dot(Dmp)
542 return Dm.astype(np.float32), cmat.astype(np.float32)
552 U, s, V = np.linalg.svd(DmtDm)
554 s[s.shape[0] - nfilt - 2:s.shape[0] - 2] = 0.
555 DmtDm1 = U.dot(np.diag(s)).dot(U.T)
556 Dmp = DmtDm1.dot(Dm.T)
560 return Dm.astype(np.float32), cmat.astype(np.float32)
572 def cov_cor(P, noise, trunc, alias, H, bp, tomo):
573 cov = np.zeros((6, 6))
574 cor = np.zeros((6, 6))
583 for i
in range(cov.shape[0]):
584 for j
in range(cov.shape[1]):
586 tmpi = P.dot(bufdict[str(i)])
587 tmpj = P.dot(bufdict[str(j)])
589 np.mean(tmpi * tmpj, axis=1) -
590 np.mean(tmpi, axis=1) * np.mean(tmpj, axis=1))
592 cov[i, j] = cov[j, i]
594 s = np.reshape(np.diag(cov), (cov.shape[0], 1))
597 cor[ok] = cov[ok] / np.sqrt(sst[ok])
611 IF = rtc.get_IFsparse(1)
614 tmp = (config.p_geom._ipupil.shape[0] -
615 (config.p_dms[0]._n2 - config.p_dms[0]._n1 + 1)) / 2
616 tmp_e0 = config.p_geom._ipupil.shape[0] - tmp
617 tmp_e1 = config.p_geom._ipupil.shape[1] - tmp
618 pup = config.p_geom._ipupil[tmp:tmp_e0, tmp:tmp_e1]
619 indx_pup = np.where(pup.flatten() > 0)[0].astype(np.int32)
620 dm_dim = config.p_dms[0]._n2 - config.p_dms[0]._n1 + 1
621 cov, cor =
cov_cor(P, noise_com, trunc_com, alias_wfs_com, H_com, bp_com, tomo_com)
622 psf = tar.get_image(0,
"le", fluxNorm=
False)
624 fname =
"/home/fferreira/Data/" + filename
626 "noise": noise_com.T,
627 "aliasing": alias_wfs_com.T,
628 "tomography": tomo_com.T,
629 "filtered modes": H_com.T,
630 "non linearity": trunc_com.T,
631 "bandwidth": bp_com.T,
636 "IF.indices": IF.indices,
637 "IF.indptr": IF.indptr,
640 "indx_pup": indx_pup,
646 "psfortho": np.fft.fftshift(psf_ortho),
649 "dm.xpos": config.p_dms[0]._xpos,
650 "dm.ypos": config.p_dms[0]._ypos,
654 h5u.save_h5(fname,
"psf", config, psf)
656 for k
in list(pdict.keys()):
657 h5u.save_hdf5(fname, k, pdict[k])
667 nfiltered = int(config.p_controllers[0].maxcond)
668 niters = config.p_loop.niter
670 config.p_loop.set_niter(niters + N_preloop)
672 rtc.load_Btt(1, Btt.dot(Btt.T))
674 rtc.set_cmat(0, cmat)
676 imat = rtc.get_imat(0)
678 Nact = ao.create_nact_geom(config.p_dms, 0)
692 com, noise_com, alias_wfs_com, tomo_com, H_com, trunc_com, bp_com, wf_com, fit, SR, SR2, psf_ortho, E, F =
loop(
694 noise_com = noise_com[N_preloop:, :]
695 trunc_com = trunc_com[N_preloop:, :]
696 alias_wfs_com = alias_wfs_com[N_preloop:, :]
697 H_com = H_com[N_preloop:, :]
698 bp_com = bp_com[N_preloop:, :]
699 tomo_com = tomo_com[N_preloop:, :]